免疫組庫（Immune Repertoire，IR）是指?jìng)€(gè)體在某個(gè)特定時(shí)間其循環(huán)系統(tǒng)中所有多樣性的T細(xì)胞和B細(xì)胞的總和，即V(D)J序列多樣性的集合。免疫組庫分析是指運(yùn)用高通量測(cè)序技術(shù)對(duì)個(gè)體免疫組庫多樣性分析，以及對(duì)T細(xì)胞和B細(xì)胞獨(dú)特性分析。
mixcr允許多種數(shù)據(jù)輸入，特異性免疫組庫測(cè)序或RNA-seq測(cè)序都可以作為輸入。
mixcr工作流程
實(shí)例展示
下載8個(gè)新冠樣本8個(gè)健康樣本RNA-seq的單端測(cè)序作為輸入。
批量下載數(shù)據(jù)：可參考RNA-seq轉(zhuǎn)錄組差異分析及可視化
#下載數(shù)據(jù)
prefetch -O /path to out/ --option-file SRR_Acc_List.txt
#將NCBI下載的數(shù)據(jù)轉(zhuǎn)為fastq格式
cat SRR_Acc_List.txt|while read id;do fastq-dump --gzip -O /path to out/ /path to input/${id}.sra;done
#path to out:輸出路徑
#path to input:輸入文件路徑
#mixcr non-enriched
#single-read
cat SRR_Acc_List.txt|while read name;do mkdir ${name};done
cat SRR_Acc_List.txt|while read name;do mixcr analyze shotgun -s hsa --starting-material rna --only-productive --impute-germline-on-export --report ./report/${name}.report ./rawdata/${name}.fastq.gz ./results/${name}/${name};done
#運(yùn)行時(shí)間比較長(zhǎng)，16個(gè)樣本大約用時(shí)35h。
#推薦使用多線程分開提交任務(wù)運(yùn)行。
mixcr輸出結(jié)果
得到mixcr的輸出結(jié)果后vdjtools與immunarch（R包）都可以用來進(jìn)行可視化輸出。
1.準(zhǔn)備工作
mixcr與vdjtools是基于java平臺(tái)開發(fā)的處理從原始序列到定量克隆型的大量免疫組數(shù)據(jù)的免疫分析軟件，在使用前要確保java環(huán)境是ok的。
最低配置
- Any Java-enabled platform (Windows, Linux, Mac OS X)
- Java version 7 or higher (download from Oracle web site)
- 1–16 Gb RAM (depending on number of clones in the sample)
檢查環(huán)境并下載軟件
官網(wǎng)下載 Java Runtime Environment，jre是java的運(yùn)行環(huán)境。
java -version #檢查java環(huán)境是否ok
官網(wǎng)下載mixcr的壓縮包并解壓，并添加至環(huán)境變量。
若未添加至環(huán)境變量，每次運(yùn)行時(shí)需使用此代碼
java -jar path_to_mixcr\mixcr.jar ... #path_to_mixcr = 解壓mixcr的地址
2.MIXCR
MIXCR分析流程主要包含三個(gè)程序步驟，分別是align，assemble，export。
- align：將測(cè)序結(jié)果比對(duì)到TCR或BCR的V,D,J和C的參考基因上。
- assemble：將上一步得到的結(jié)果進(jìn)行組裝，可以得到特定的基因區(qū)域序列，如CDR3序列。
- export：將比對(duì)結(jié)果或組裝結(jié)果輸出，得到可閱讀的txt文件。
MIXCR支持輸入單端或雙端的測(cè)序文件，可使用原始測(cè)序文件也可使用質(zhì)控后的文件，質(zhì)控方法和RNA-seq差異分析中的方法一致，使用fastp或者trimmomatic進(jìn)行質(zhì)控，參考RNA-seq轉(zhuǎn)錄組差異分析及可視化。
3.一站式流程使用
MIXCR軟件提供了兩種analysis模式，包含了上述三種參數(shù)，一站式進(jìn)行分析。
-
analyze amplicon： for analysis of
targetedTCR/IG library amplification (5’RACE, Amplicon, Multiplex, etc). -
analyze shotgun： for analysis of
random fragments(RNA-Seq, Exome-Seq, etc).
#environment：Linux
# for enriched targeted TCR/IG libraries
mixcr analyze amplicon
-s hsa \ # or mmu
--starting-material dna \#or rna
--5-end v-primers \#or no-v-primers
--3-end j-primers \#or j-c-intron-primers or c-primers
--adapters adapters-present \#or no-adapters
R1.fastq.gz R2.fastq.gz ./path/mixcr_result/sample1
| Option?????????????????????????? | Description |
|---|---|
-s,--species |
hsa (or HomoSapiens), mmu (or MusMusculus), rat (currently only TRB, TRA and TRD are supported), or any species from IMGT ? library, if it is used (see here import segments). |
--starting-material |
rna or dna
|
--5-end |
no-v-primers（e.g. 5’RACE with template switch oligo ）or v-primers
|
--3-end |
j-primers or j-c-intron-primers or c-primers
|
--adapters |
adapters-present or no-adapters
|
The following parameters are optional:
| Option???????????????????????????????????? | Default | Description |
|---|---|---|
--report |
analysis_name.report |
Report file. |
--receptor-type |
xcr |
By default, all T- and B-cell receptor chains are analyzed.Possible values are:xcr (all chains), tcr, bcr, tra, trb, trg, trd, igh, igk, igl. |
--contig-assembly |
false |
Whether to assemble full receptor sequences (assembleContigs). This option may slow down the computation. |
--impute-germline-on-export |
false |
Use germline segments (printed with lowercase letters) for uncovered gene features. |
--region-of-interest |
CDR3 |
MiXCR will use only reads covering the whole target region; reads which partially cover selected region will be dropped during clonotype assembly. All non-CDR3 options require long high-quality paired-end data. See Gene features and anchor points for details. |
--only-productive |
false |
Filter out-of-frame and stop-codons in export |
--align |
Additional parameters for align step specified with double quotes (e.g –align “–limit 1000” –align “-OminSumScore=100”) | |
--assemble |
Additional parameters for assemble step specified with double quotes (e.g –assemble “-ObadQualityThreshold=0”). | |
--assembleContigs |
Additional parameters for assembleContigs step specified with double quotes. | |
--export |
Additional parameters for exportClones step specified with double quotes. |
# for non-enriched or random fragments
mixcr analyze shotgun
-s hsa # or mmu
--starting-material dna #or rna
[OPTIONS] input_file1 [input_file2] analysis_name
| Option???????????????????????????????????? | Default | Description |
|---|---|---|
--report |
analysis_name.report |
Report file. |
--receptor-type |
xcr |
Dedicated receptor type for analysis. By default, all T- and B-cell receptor chains are analyzed. MiXCR has special aligner kAligner2, which is used when B-cell receptor type is selected. Possible values for --receptor-type are: xcr (all chains), tcr, bcr, tra, trb, trg, trd, igh, igk, igl. |
--contig-assembly |
false |
Whether to assemble full receptor sequences (assembleContigs). This option may slow down the computation. |
--impute-germline-on-export |
false |
Use germline segments (printed with lowercase letters) for uncovered gene features. |
--only-productive |
false |
Filter out-of-frame and stop-codons in export |
--assemble-partial-rounds |
2 |
Number of consequent assemblePartial executions. |
--do-not-extend-alignments |
Do not perform extension of good TCR alignments. | |
--align |
Additional parameters for align step specified with double quotes (e.g –align “–limit 1000” –align “-OminSumScore=100”) | |
--assemblePartial |
Additional parameters for assemblePartial step specified with double quotes. | |
--extend |
Additional parameters for extend step specified with double quotes. | |
--assemble |
Additional parameters for assemble step specified with double quotes (e.g –assemble “-ObadQualityThreshold=0”). | |
--assembleContigs |
Additional parameters for assembleContigs step specified with double quotes. | |
--export |
Additional parameters for exportClones step specified with double quotes. |
上述兩種分析模式可滿足大多要求，更多參數(shù)細(xì)節(jié)參考mixcr用戶手冊(cè)
4.逐級(jí)分析
沒有特殊的分析需求這部分可略過。
4.1 Alignment
mixcr align [options] input_file1 [input_file2] output_file.vdjca
| options????????????????? | default value | description |
|---|---|---|
-r,--report
|
輸出文件名稱 | |
-l,--loci
|
ALL | 比對(duì)的靶點(diǎn)，可由,連接多個(gè)。IGH,IGL,IGK, TRA, TRB, TRG, TRD, IG (for all immunoglobulin loci),TCR (for all T-cell receptor loci), ALL (for all loci) . |
-s,--specoes
|
HomoSapiens | 同上個(gè)表格 |
-p,--parameters
|
default |
defaultorrna-seq,rna-seq是專門為分析rna-seq數(shù)據(jù)而優(yōu)化的。 |
-i,--diff-loci
|
接受與V和J基因不同loci的比對(duì)(在默認(rèn)情況下，這種比對(duì)被丟棄)。 | |
-t |
number of available CPU cores | 線程數(shù) |
-n,--limit
|
只處理輸入文件的前-n個(gè)序列。 |
|
-Oparameter = value |
比對(duì)區(qū)域，具體見下面表格 |
Alignment parameters
| Parameters???? | Default value | Description |
|---|---|---|
vParameters.geneFeatureToAlign |
VRegion |
region in V gene which will be used as target in align |
dParameters.geneFeatureToAlign |
DRegion |
region in D gene which will be used as target in align |
jParameters.geneFeatureToAlign |
JRegion |
region in J gene which will be used as target in align |
cParameters.geneFeatureToAlign |
CRegion |
region in C gene which will be used as target in align |
-OvParameters.geneFeatureToAlign是比對(duì)過程中的一個(gè)重要參數(shù)，有三個(gè)常用value。
-
– VRegion（默認(rèn)值）：如果你的免疫組庫是使用多重PCR構(gòu)建的5端文庫則使用這個(gè)選項(xiàng)。 -
– VTranscript：以RNA為起始材料進(jìn)行的非模板特異性擴(kuò)增，例如5'RACE。使用這個(gè)選項(xiàng)有助于提高測(cè)序信息5'端的利用率，有助于提高V基因的準(zhǔn)確性。 -
– VGene：DNA為起始材料，5'端的信息（包括V區(qū)內(nèi)含子，leader sequence，5'UTR）應(yīng)該存在于測(cè)序信息中。
如果想獲得TCR或BCR的全部克隆型（包括FR和CDR區(qū)域），使用- VTranscript或者- VGene。
4.1.1 分析RNA-seq數(shù)據(jù)
Analysis of RNA-Seq data performed with -p rna-seq option is almost equivalent to the following set of aligners parameters:
-
(most important) turned off floating bounds of V and J alignments:
–-OvParameters.parameters.floatingLeftBound=false
–-OjParameters.parameters.floatingRightBound=false -
higher thresholds:
–-OvParameters.parameters.absoluteMinScore=80(was 40)
–-OjParameters.parameters.absoluteMinScore=70(was 40)
–-OminSumScore=200(was 120; see below) -
more strict scoring for all alignments (V, J, C):
–-OxParameters.parameters.scoring.gapPenalty=-21
–-OxParameters.parameters.scoring.subsMatrix=’simple(match=5,mismatch=-12)’
4.2 Assemble clones
這一步使用到上一步得到的.vdjca文件，提取特定的基因序列，最后輸出.clns文件。
mixcr assemble [options] alignments.vdjca output.clns
| options???????????????????????? | default value | description |
|---|---|---|
-r --report
|
Report file name. | |
-t --threads
|
number of available CPU cores | Number of processing threads. |
-a, --write-alignments
|
Save initial alignments and alignments <> clones mapping in the resulting .clna file. |
|
-Oparameter=value |
Overrides default value of assembler parameter (see next subsection). |
-OassemblingFeatures參數(shù)用來設(shè)置需要組裝的區(qū)域。默認(rèn)值是CDR3。如果要組裝全序列，使用VDJRegion選項(xiàng)。
mixcr assemble -OassemblingFeatures="[V5UTR+L1+L2+FR1,FR3+CDR3]" alignments.vdjca output.clns
(note:assemblingFeatures must cover CDR3).
Other global parameters are:
| Parameter???? | Default value | Description |
|---|---|---|
minimalClonalSequenceLength ?? |
12 |
Minimal length of clonal sequence |
badQualityThreshold |
20 |
Minimal value of sequencing quality score: nucleotides with lower quality are considered as “bad”. If sequencing read contains at least one “bad” nucleotide within the target gene region, it will be deferred at initial assembling stage, for further processing by mapper. |
maxBadPointsPercent |
0.7 |
Maximal allowed fraction of “bad” points in sequence: if sequence contains more than maxBadPointsPercent “bad” nucleotides, it will be completely dropped and will not be used for further processing by mapper. Sequences with the allowed percent of “bad” points will be mapped to the assembled core clonotypes. Set -OmaxBadPointsPercent=0 in order to completely drop all sequences that contain at least one “bad” nucleotide. |
qualityAggregationType |
Max |
Algorithm used for aggregation of total clonal sequence quality during assembling of sequencing reads. Possible values: Max (maximal quality across all reads for each position), Min (minimal quality across all reads for each position), Average (average quality across all reads for each position), MiniMax (all letters has the same quality which is the maximum of minimal quality of clonal sequence in each read). |
minimalQuality |
0 |
Minimal allowed quality of each nucleotide of assembled clone. If at least one nucleotide in the assembled clone has quality lower than minimalQuality, this clone will be dropped (remember that qualities of reads are aggregated according to selected aggregation strategy during core clonotypes assembly; see qualityAggregationType). |
addReadsCountOnClustering |
false |
Aggregate cluster counts when assembling final clones: if addReadsCountOnClustering is true, then all children clone counts will be added to the head clone; thus head clone count will be a total of its initial count and counts of all its children. Refers to further clustering strategy (see below). Does not refer to mapping of low quality sequencing reads described above. |
example：
mixcr assemble -ObadQualityThreshold=10 alignments.vdjca output.clns
mixcr assemble -OmaxBadPointsPercent=0 alignments.vdjca output.clns #In order to prevent mapping of low quality reads (filter them off) one can set maxBadPointsPercent to zero
4.3Export
將比對(duì)結(jié)果或組裝結(jié)果輸出，得到可閱讀的txt文件。
#
mixcr exportAlignments [options] alignments.vdjca alignments.txt
mixcr exportClones [options] clones.clns clones.txt
5.結(jié)果解讀
參考mixcr官方手冊(cè)export部分。
