全網(wǎng)最全免疫組庫分析流程+代碼復現(xiàn)（一）

免疫組庫（Immune Repertoire，IR） 是指個體在某個特定時間其循環(huán)系統(tǒng)中所有多樣性的T細胞和B細胞的總和，即V(D)J序列多樣性的集合。免疫組庫分析是指運用高通量測序技術(shù)對個體免疫組庫多樣性分析，以及對T細胞和B細胞獨特性分析。
關(guān)于更多免疫組庫分析的相關(guān)概念可以參考https://www.nature.com/articles/s41435-022-00180-w

免疫組庫分析全流程

1.測序數(shù)據(jù)標準化

因為建庫質(zhì)量不同我們得到的測序數(shù)據(jù)的有效序列的占比可能是是不同的且測序結(jié)果的數(shù)據(jù)大小也可能不同，為了保證分析的準確性，我們需要對測序數(shù)據(jù)進行標準化處理。
不過在大多研究中不進行標準化處理，我們在分析時可以根據(jù)自己課題不同的需求來決定是否需要標準化。

個人理解：此步標準化是為了消除建庫中的誤差，比如建庫時間不同，樣品批次不同等等。如果能保證建庫操作沒有誤差切批次一致的情況下，可以不進行標準化處理，即認為誤差是樣本本身的差異。

也可以在恢復完免疫組庫后在對結(jié)果進行標準化，這個后面再說。
對原始數(shù)據(jù)的標準化，我一般會先使用trimmomatic對數(shù)據(jù)進行過濾，然后使用seqkit軟件抽樣，保證每個樣本的總reads數(shù)是一致的。
具體用法見：http://www.itdecent.cn/p/a8935adebaae

#使用trimmomatic去除接頭序列和低質(zhì)量的序列
cat name.txt |while read id
do
trimmomatic PE -phred33 \
./rawdata/${id}_R1_001.fastq.gz ./rawdata/${id}_R2_001.fastq.gz \
./rawdata/clean_rawdata/${id}_paired_R1.fastq.gz ./rawdata/clean_rawdata/${id}_unpaired_R1.fastq.gz \
./rawdata/clean_rawdata/${id}_paired_R2.fastq.gz ./rawdata/clean_rawdata/${id}_unpaired_R2.fastq.gz \
ILLUMINACLIP:/home/heboxiao/project/conda3/envs/RNAseq/share/trimmomatic-0.39-2/adapters/TruSeq2-PE.fa:2:30:10:1:true \
LEADING:3 \
TRAILING:3 \
SLIDINGWINDOW:4:15 \
MINLEN:50
done
#使用seqkit提取有效序列，建庫時添加了一段固定序列
cat name.txt |while read id
do
zcat ./rawdata/clean_rawdata/${id}_paired_R1.fastq.gz | seqkit grep -s -r -i -p ^ATGCATCGGATCTTCAGCATGA -o ./rawdata/clean_${id}_R1.fastq.gz
seqkit seq ./rawdata/clean_${id}_R1.fastq.gz -n -i > tmp.txt
seqkit grep -f tmp.txt ./rawdata/clean_rawdata/${id}_paired_R2.fastq.gz -o ./rawdata/clean_${id}_R2.fastq.gz
seqkit stats ./rawdata/clean_${name}_R1.fastq.gz >> ./stats.txt
done
#查看stats.txt,以最小的樣本reads數(shù)作為抽樣標準
cat name.txt |while read id
do
zcat ./rawdata/clean_${id}_R1.fastq.gz | seqkit sample -n 1500000 -o ./rawdata/normal/${id}_R1.fastq.gz
seqkit seq ./rawdata/normal/${id}_R1.fastq.gz -n -i > tmp.txt
seqkit grep -f tmp.txt ./rawdata/clean_${id}_R2.fastq.gz -o ./rawdata/normal/${id}_R2.fastq.gz
done

2.使用相關(guān)軟件（Trust4 or Mixcr）進行組裝

TRUST4是哈佛Shirley Liu課題組研發(fā)的一款免疫組庫分析工具，他可以直接從Bluk RNA-seq數(shù)據(jù)和scRNA-seq數(shù)據(jù)中重構(gòu)免疫組庫信息，當然免疫組測序也可以使用此軟件分析。
文章內(nèi)提到，與傳統(tǒng)的mixcr相比，在RNAseq數(shù)據(jù)處理方面有不錯的效果。https://doi.org/10.1038/s41592-021-01142-2

Tcr Receptor Utilities for Solid Tissue (TRUST) is a computational tool to analyze TCR and BCR sequences using unselected RNA sequencing data, profiled from solid tissues, including tumors. TRUST4 performs de novo assembly on V, J, C genes including the hypervariable complementarity-determining region 3 (CDR3) and reports consensus of BCR/TCR sequences. TRUST4 then realigns the contigs to IMGT reference gene sequences to report the corresponding information. TRUST4 supports both single-end and paired-end sequencing data with any read length.

用法說明

Usage: ./run-trust4 [OPTIONS]
        Required:
                -b STRING: path to bam file 
                -1 STRING -2 STRING: path to paired-end read files
                -u STRING: path to single-end read file
                -f STRING: path to the fasta file coordinate and sequence of V/D/J/C genes
                
Optional:
                --ref STRING: path to detailed V/D/J/C gene reference file, such as from IMGT database. (default: not used).
  • -b：可以使用已經(jīng)比對到基因組上的bam文件
  • -1,-2,-u：輸入fastq文件，可以是壓縮的。-u為單端測序
  • -f：VDJC基因的參考基因組,包含基因坐標。 如果輸入的fastq文件可以直接把IMGT的的database當做-f輸入
  • --ref：更詳細的參考基因組，比如來自于IMGT的database，包含基因的亞類，比如IGHA*02。

兩個參考文件的構(gòu)建

兩個參考文件的區(qū)別
#  GRCm38_bcrtcr.fa
>TRBV1 6 40891229 40891885 +
GGATTCTCTTCTCTTGCCTGATGCCCTGCATGCCCCACAGAGATAGAGAGAACCTGAGGTCTCAGAGATGTGGCAGTTTTGCATTCTGTGCCTCTGTGTACTCATGGCTTCTGTGGCTACAGACCCCACAGTGACTTTGCTGGAGCAAAACCCAAGGTGGCGTCTGGTAC
>TRBV2 6 41047340 41047995 +
GCACAAGAAACTGTCTGAGGAGACCCAAGAGGTAGCTGGGTATGGTCTAGGTCTGGGAGCCTATGCCCTGACAAGCATGAAGGGCAAAGAGGAAGTGTGAGCTGAAACCATCTACACAACAGGGAGCATCCAAGCAAAGCATTTTACAGGTCAGCAAAAGGCACCCAGAC
.........

#  mouse_IMGT+C.fa
>IGHA*02
GGTCCTACTCCTCCTCCTCCTATTACTATTCCT..................TCCTGCCAG
CCCAGCCTGTCACTGCAGCGGCCAGCTCTTGAGGACCTGCTCCTG......GGTTCAGAT
GCCAGCATCACATGTACTCTGAATGGCCTGAGAAATCCT...GAGGGAGCTGTCTTCACC
TGGGAGCCCTCCACTGGGAAGGAT............GCAGTGCAGAAGAAAGCTGCGCAG
AAT
>IGHA*03
............GAGTCTGCGAGAAATCCCACCATCTACCCACTGACACTCCCACCAGCT
CTGTCA............AGTGACCCAGTGATAATCGGCTGCCTGATTCACGATTACTTC
CCTTCC...GGCACGATGAATGTGACCTGGGGAAAGAGTGGGAAGGATATA.........
...ACCACCGTAAACTTCCCACCTGCCCTGGCCTCTGGG..................GGA
CGGTACACCATGAGCAGCCAGTTGACCCTGCCAGCTGTCGAGTGC......CCAGAAGGA
.......

更多的Usage請參考：https://github.com/liulab-dfci/TRUST4

這里演示免疫組測序的分析方法

Bluk測序

cat name.txt|while read name
do
~/project/biosoft/TRUST4/run-trust4 \
-1 ./rawdata/${name}_R1_001.fastq.gz \
-2 ./rawdata/${name}_R2_001.fastq.gz \
-t 8 \
-f /home/heboxiao/project/biosoft/TRUST4/mouse/GRCm38_bcrtcr.fa \
--ref /home/heboxiao/project/biosoft/TRUST4/mouse/mouse_IMGT+C.fa \
--repseq \ #如果是RNA-seq數(shù)據(jù)這個參數(shù)去掉即可
-o ./results/${name}
done

輸出文件中主要關(guān)注xxx_report.tsv，格式如下：

#count frequency CDR3nt CDR3aa V D J C cid cid_full_length
216544 1.025334e-01 TATCTCTGTGCTTTGGGGAATTATGGGAGCAGTGGCAACAAGCTCATCTTT YLCALGNYGSSGNKLIF TRAV13-1*02 . TRAJ32*01 . assemble6 0
182189 8.626656e-02 TATCTCTGTGCTCCTCTCGGGATACAACAAACTCACTTTT out_of_frame TRAV13-1*01 . TRAJ11*01 . assemble4 0
181785 8.607514e-02 TACTACTGCGCTCTGAGTGATCCAGGCACTGGGAGTAACAGGCTCACTTTT YYCALSDPGTGSNRLTF TRAV6N-6*01 . TRAJ28*01 . assemble5 0

需要注意的是頻率是IG和TR單獨分開計算的，在CDR3中_代表終止子，？代表核酸序列中出現(xiàn)了N。

含有barcode和umi的測序數(shù)據(jù)分析

示例：10x空間的免疫組庫分析

cat name.txt|while read name
do
~/project/biosoft/TRUST4/run-trust4 \
-1 ./rawdata/${name}_R1_001.fastq.gz \
-2 ./rawdata/${name}_R2_001.fastq.gz \
-f /home/heboxiao/project/biosoft/TRUST4/mouse/GRCm38_bcrtcr.fa \
--ref /home/heboxiao/project/biosoft/TRUST4/mouse/mouse_IMGT+C.fa \
--read1Range 28 -1 \
--read2Range 0 -1 \
--barcode ./rawdata/${name}_R1_001.fastq.gz \
--barcodeRange 0 15 + \
--UMI ./rawdata/${name}_R1_001.fastq.gz \
--umiRange 16 27 + \
--repseq \
--outputReadAssignment \
-o ./results/${name}
done

會多出來一個xxxx_barcode_report.tsv文件，每個barcode會選取兩個CDR3鏈，格式如下：

#barcode cell_type chain1 chain2 secondary_chain1 secondary_chain2
AATACCGGAGGGCTGT abT TRBV2*01,*,TRBJ2-2*01,*,TGTGCCAGCAGCCAAGACCCCGTAAACACCGGGCAGCTCTACTTT,CASSQDPVNTGQLYF,1.00,AATACCGGAGGGCTGT_1089,0.00,0 * TRBV5*01,T
TGCATGTGGTAATCTA abT TRBV14*01,TRBD1*01,TRBJ1-1*01,*,TGTGCCAGCAGCCGGGGACAAACAAACACAGAAGTCTTCTTT,CASSRGQTNTEVFF,1.00,TGCATGTGGTAATCTA_7760,100.00,0 * *
GAGCTATGGGCATATC abT TRBV14*01,TRBD1*01,TRBJ1-1*01,*,TGTGCCAGCAGCCGGGGACAAACAAACACAGAAGTCTTCTTT,CASSRGQTNTEVFF,1.00,GAGCTATGGGCATATC_8159,100.00,0 * *

后續(xù)分析參考：https://hbox.cool/immune-repertoire/

3.使用VDJtools分析結(jié)果

基于Java的分析軟件，command line運行模式。
使用vdjtools進行免疫組庫分析

4.使用immunarch R包進行可視化

兼容mixcr和vdjtools的結(jié)果，可以直接可視化。

5.使用R語言手動進行個性化分析

自己寫代碼實現(xiàn)，可以對免疫組庫的結(jié)構(gòu)更加清晰，同時對結(jié)果合理性的判斷也會更加準確。

  • 免疫組庫多樣性（Richness/Diversity/Clonality）的計算及可視化
  • CDR3 length distribution analysis
  • 氨基酸占比分析
  • Gene usage analysis
  • V-J junctions analysis
  • CDR3 amino acid motif
    詳細代碼
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容