使用MiXCR进行免疫组库分析
其实我不是很想写这个免疫组库专题了，阅读量太低，估计认真跟下去也不会很多。
前面我带领大家通过IMGT数据库认知免疫组库，而且也一起从IMGT数据库下载免疫组库相关fasta序列，免疫组库重要的研究对象就是分成BCR的IGH,IGK,IGL这3类，以及TCR的TRA,TRB,TRD,TRG，它们各自都有V，D（可选），J，C基因。
接下来又认识了免疫组库测序数据，知道了免疫组库测序数据的一些特性，并且使用igblast进行免疫组库分析了，但是那个是初步的比对，虽然找到每一个测序片段的V，D（可选），J，C基因，并且拿到 CDR3序列，中间步骤有点多，常规的测序数据过滤就算了，还需要把PE数据合并，fastq格式转为fasta格式，而且本身igblast软件就很难使用，数据库文件构建也繁琐。其实也有包装好的一站式流程。
免疫组库测序数据可以使用MiXCR进行分析，比如文章 T cell receptor next-generation sequencing reveals cancer-associated repertoire metrics and reconstitution after chemotherapy in patients with hematological and solid tumors ，就提到了：
NGS and de-multiplexing was performed on an Illumina MiSeq sequencer (600–cycle, single indexed, paired-end run). Analysis of the TRB locus was computed using the MiXCR analysis tool
MiXCR是MiLaboratory开发的，他们实验室出品了3款：MiXCR, MIGEC, VDJtools ，在免疫组库领域都很出名。其中MiXCR工作流程示意图如下：
安装MiXCR软件
因为MiXCR有二进制版本，所以推荐直接下载它并且解压即可使用啦！
https://github.com/milaboratory/mixcr/releases/download/v3.0.13/mixcr-3.0.13.zip （2020年05月05日查看的最新版） 文档：https://mixcr.readthedocs.io/en/latest/index.html
不过MiXCR是基于java平台开发，所以启动MiXCR时候，需要保证java环境是ok的，而且要求java8哦。
如果是Linux或者macOS，那么非常简单，直接使用里面的mixcr命令即可调用该软件。
如果是Windows平台呢，需要使用 java -Xmx4g -Xms3g -jar 来调用下载到的MiXCR有二进制版本里面的jar文件哦。
我的是macOS操作系统，软件下载解压后如下所示：
(base) jmzengdeMBP:mixcr-3.0.13 jmzeng$ pwd/Users/jmzeng/biosoft/MiXCR/mixcr-3.0.13(base) jmzengdeMBP:mixcr-3.0.13 jmzeng$ tree -h.├── [1.4K] LICENSE├── [ 96] libraries├── [4.0K] mixcr└── [9.3M] mixcr.jar
1 directory, 3 files一站式流程使用
MiXCR软件提供一个 analyze模式，可以一站式完成(align, assemblePartial, extend, assemble, assembleContigs and export) 这些分析，需要认真读 https://mixcr.readthedocs.io/en/latest/analyze.html 说明书，其提供的例子如下：
mixcr analyze amplicon -s <species> \ --starting-material <startingMaterial> \ --5-end <5End> --3-end <3End> \ --adapters <adapters> \ [OPTIONS] input_file1 [input_file2] analysis_name在前面我们认识的免疫组库测序数据，是人类的，MiSeq测序仪，PE300测序策略，TRB，DNA测序。所以我们构建的命令是：
mixcr=/Users/jmzeng/biosoft/MiXCR/mixcr-3.0.13/mixcr$mixcr analyze amplicon \ -s hsa \ --starting-material dna \ --5-end v-primers --3-end c-primers \ --adapters adapters-present \ raw/ERR3445170_1.fastq.gz raw/ERR3445170_2.fastq.gz ERR3445170因为我设置了 --adapters adapters-present，所以这个软件也会帮忙去除我们测序数据里面的接头序列，所以我这里直接使用了raw测序数据fq文件进行分析，无需经过前面igblastn的trim_galore过程。
很轻松就出来了结果，看起来调用了我Mac全部计算资源
运行日志如下；
NOTE: report file is not specified, using ERR3445170.report to write report.Alignment: 2.1%Alignment: 15.4% ETA: 00:00:38Alignment: 27.9% ETA: 00:00:17Alignment: 39.5% ETA: 00:00:15Alignment: 51.5% ETA: 00:00:12Alignment: 63% ETA: 00:00:09Alignment: 74.2% ETA: 00:00:06Alignment: 85.4% ETA: 00:00:03Alignment: 96.9% ETA: 00:00:00============= Report ==============Analysis time: 29.52sTotal sequencing reads: 84833
Successfully aligned reads: 83421 (98.34%)Paired-end alignment conflicts eliminated: 3027 (3.57%)Alignment failed, no hits (not TCR/IG?): 850 (1%)Alignment failed because of absence of V hits: 4 (0%)Alignment failed because of absence of J hits: 507 (0.6%)No target with both V and J alignments: 49 (0.06%)Alignment failed because of low total score: 2 (0%)Overlapped: 81097 (95.6%)Overlapped and aligned: 80179 (94.51%)Alignment-aided overlaps: 1313 (1.64%)Overlapped and not aligned: 918 (1.08%)V gene chimeras: 29 (0.03%)J gene chimeras: 4 (0%)TRB chains: 83418 (100%)IGH chains: 3 (0%)Realigned with forced non-floating bound: 10098 (11.9%)Realigned with forced non-floating right bound in left read: 3 (0%)Realigned with forced non-floating left bound in right read: 3 (0%)Initialization: progress unknownAssembling initial clonotypes: 100%Building clones: 25.2%Writing clones: 0%============= Report ==============Analysis time: 2.33sFinal clonotype count: 2116
Average number of reads per clonotype: 35.18Reads used in clonotypes, percent of total: 74448 (87.76%)Reads used in clonotypes before clustering, percent of total: 76478 (90.15%)Number of reads used as a core, percent of used: 75378 (98.56%)Mapped low quality reads, percent of used: 1100 (1.44%)Reads clustered in PCR error correction, percent of used: 2030 (2.65%)Reads pre-clustered due to the similar VJC-lists, percent of used: 33 (0.04%)Reads dropped due to the lack of a clone sequence, percent of total: 4236 (4.99%)Reads dropped due to low quality, percent of total: 116 (0.14%)Reads dropped due to failed mapping, percent of total: 2590 (3.05%)Reads dropped with low quality clones, percent of total: 1 (0%)Clonotypes eliminated by PCR error correction: 832Clonotypes dropped as low quality: 1Clonotypes pre-clustered due to the similar VJC-lists: 27TRB chains: 2116 (100%)软件内置了免疫组库数据库，所以会根据提供的物种名字，自己去比对BCR的IGH,IGK,IGL这3类，以及TCR的TRA,TRB,TRD,TRG，它们各自都有V，D（可选），J，C基因。因为我们这个测试数据是TRB的，所以只有TRB的结果出来，如下：
2117 ERR3445170.clonotypes.ALL.txt 1 ERR3445170.clonotypes.IGH.txt 1 ERR3445170.clonotypes.IGK.txt 1 ERR3445170.clonotypes.IGL.txt 1 ERR3445170.clonotypes.TRA.txt 2117 ERR3445170.clonotypes.TRB.txt 1 ERR3445170.clonotypes.TRD.txt 1 ERR3445170.clonotypes.TRG.txt理解MiXCR软件的输出结果
如下：
输出的列比较多：
1cloneId 2cloneCount 3cloneFraction 4targetSequences 5targetQualities 6allVHitsWithScore 7allDHitsWithScore 8allJHitsWithScore 9allCHitsWithScore 10allVAlignments 11allDAlignments 12allJAlignments 13allCAlignments 14nSeqFR1 15minQualFR1 16nSeqCDR1 17minQualCDR1 18nSeqFR2 19minQualFR2 20nSeqCDR2 21minQualCDR2 22nSeqFR3 23minQualFR3 24nSeqCDR3 25minQualCDR3 26nSeqFR4 27minQualFR4 28aaSeqFR1 29aaSeqCDR1 30aaSeqFR2 31aaSeqCDR2 32aaSeqFR3 33aaSeqCDR3 34aaSeqFR4 35refPoints其实比较容易理解，分成V，D（可选），J，基因的比对情况，分成FR1-4以及CDR1-3的序列情况详情。实际上它可以把这拆分成为多个文件，这样方便浏览
第一部分是免疫组库核心序列
因为V,D,J,C基因具有一定程度的相似性，所以并不需要全部的几百个碱基序列才能确定是V,D,J,C基因的哪些组成。
cloneIdcloneCountcloneFractiontargetSequences020112.00.27014829142488717TGTGCCAGCAGCTTGGCCCGAGAAGGGATTGAAAACACCATATATTTT19264.00.12443584784010316TGTGCCAGCAGCCGCCTGGACCTAGGATTTTACGAGCAGTACTTC25037.00.06765796260477111TGTGCCAGCAATGTACAGCTCCTATAACACTGAAGCTTTCTTT32710.00.036401246507629485TGCAGCGGGCAAGGGGGGGAGAGGGGAGCAGATACGCAGTATTTT41851.00.024862991618310768TGTGCCAGCAGTTTATCTCAAGGGACGCTATACAATGAGCAGTTCTTC51084.00.014560498603051794TGCAGCGTTGTACCCGATAGCACAGATACGCAGTATTTT61082.00.014533634214485279TGCGCCAGCAGCTTGACAAACACCATATATTTT7890.00.01195465291209972TGCAGTGTGTGGTGGACCGGGGGGTATGGCTACACCTTC8538.00.007226520524392865TGTGCCTGGAGTCCCGGGTTTGGGGATTCACCCCTCCACTTT命令是：cut -f 1-4 ERR3445170.clonotypes.TRB.txt|head
第二部分是V,D,J,C基因详情
主要是方便统计该样本的V,D,J,C基因占比
allVHitsWithScoreallJHitsWithScoreallVAlignmentsallJAlignmentsTRBV5-5*00(886.1)TRBJ1-3*00(197)549|567|585|0|18||90.026|42|70|32|48||80.0TRBV6-5*00(872.3)TRBJ2-7*00(175.7)519|530|556|0|11||55.024|39|67|30|45||75.0TRBV6-5*00(863.7)TRBJ1-1*00(217.9)519|529|556|0|10||50.022|40|68|25|43||90.0TRBV29-1*00(942.1)TRBJ2-3*00(204.2)700|707|734|0|7||35.024|41|69|28|45||85.0TRBV27*00(923.6)TRBJ2-1*00(200.4)556|573|593|0|17||85.024|42|70|30|48||90.0TRBV29-1*00(955.5)TRBJ2-3*00(230.1)700|710|734|0|10||50.019|41|69|17|39||110.0TRBV5-1*00(693.5)TRBJ1-3*00(195.9)554|569|590|0|15||75.026|42|70|17|33||80.0TRBV20-1*00(935.1)TRBJ1-2*00(208.9)759|766|793|0|7||35.025|40|68|24|39||75.0TRBV30*00(904.5)TRBJ1-6*00(160.6)787|799|821|0|12||60.028|45|73|25|42||85.0命令是：cut -f 6-13 ERR3445170.clonotypes.TRB.txt|head|cut -f 1,3,5,7
第3，4部分大家其实仅仅是关系CDR3区域
其中第3 部分是CDR3区域的核苷酸序列，第4部分是CDR3区域的氨基酸序列。
nSeqCDR3aaSeqCDR3TGTGCCAGCAGCTTGGCCCGAGAAGGGATTGAAAACACCATATATTTTCASSLAREGIENTIYFTGTGCCAGCAGCCGCCTGGACCTAGGATTTTACGAGCAGTACTTCCASSRLDLGFYEQYFTGTGCCAGCAATGTACAGCTCCTATAACACTGAAGCTTTCTTTCASNVQL_YNTEAFFTGCAGCGGGCAAGGGGGGGAGAGGGGAGCAGATACGCAGTATTTTCSGQGGERGADTQYFTGTGCCAGCAGTTTATCTCAAGGGACGCTATACAATGAGCAGTTCTTCCASSLSQGTLYNEQFFTGCAGCGTTGTACCCGATAGCACAGATACGCAGTATTTTCSVVPDSTDTQYFTGCGCCAGCAGCTTGACAAACACCATATATTTTCASSLTNTIYFTGCAGTGTGTGGTGGACCGGGGGGTATGGCTACACCTTCCSVWWTGGYGYTFTGTGCCTGGAGTCCCGGGTTTGGGGATTCACCCCTCCACTTTCAWSPGFGDSPLHF命令是：cut -f 24,33 ERR3445170.clonotypes.TRB.txt|head
也可以自己构建数据库
因为这个软件是内置了免疫组库数据库，所以会根据提供的物种名字，自己去比对BCR的IGH,IGK,IGL这3类，以及TCR的TRA,TRB,TRD,TRG，它们各自都有V，D（可选），J，C基因。
可能会出现版本更新不及时，或者你研究的物种，并不在软件考虑范围内的情况，那么你可能是需要自行构建哦。
参考https://mixcr.readthedocs.io/en/latest/importSegments.html
文末友情宣传
生信爆款入门-全球听（买一得五）（第5期）（可能是最后一期）你的生物信息学入门课 (必看！)数据挖掘第3期（两天变三周，实力加量），医学生/临床医师首选技能提高课 生信技能树的2019年终总结 ，你的生物信息学成长宝藏 2020学习主旋律，B站74小时免费教学视频为你领路，还等什么，看啊！！！