免疫组库（Immune Repertoire，IR） 是指个体在某个特定时间其循环系统中所有多样性的T细胞和B细胞的总和，即V(D)J序列多样性的集合。免疫组库分析是指运用高通量测序技术对个体免疫组库多样性分析，以及对T细胞和B细胞独特性分析。
关于更多免疫组库分析的相关概念可以参考https://www.nature.com/articles/s41435-022-00180-w
免疫组库分析全流程
1.测序数据标准化
因为建库质量不同我们得到的测序数据的有效序列的占比可能是是不同的且测序结果的数据大小也可能不同，为了保证分析的准确性，我们需要对测序数据进行标准化处理。
不过在大多研究中不进行标准化处理，我们在分析时可以根据自己课题不同的需求来决定是否需要标准化。
个人理解：此步标准化是为了消除建库中的误差，比如建库时间不同，样品批次不同等等。如果能保证建库操作没有误差切批次一致的情况下，可以不进行标准化处理，即认为误差是样本本身的差异。
也可以在恢复完免疫组库后在对结果进行标准化，这个后面再说。
对原始数据的标准化，我一般会先使用trimmomatic对数据进行过滤，然后使用seqkit软件抽样，保证每个样本的总reads数是一致的。
具体用法见：https://www.jianshu.com/p/a8935adebaae
#使用trimmomatic去除接头序列和低质量的序列
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数作为抽样标准
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.使用相关软件（Trust4 or Mixcr）进行组装
- 2.1TRUST4
TRUST4是哈佛Shirley Liu课题组研发的一款免疫组库分析工具，他可以直接从Bluk RNA-seq数据和scRNA-seq数据中重构免疫组库信息，当然免疫组测序也可以使用此软件分析。
文章内提到，与传统的mixcr相比，在RNAseq数据处理方面有不错的效果。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：可以使用已经比对到基因组上的bam文件
- -1,-2,-u：输入fastq文件，可以是压缩的。-u为单端测序
- -f：VDJC基因的参考基因组,包含基因坐标。 如果输入的fastq文件可以直接把IMGT的的database当做-f输入
- --ref：更详细的参考基因组，比如来自于IMGT的database，包含基因的亚类，比如IGHA*02。
两个参考文件的区别
# 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数据这个参数去掉即可
-o ./results/${name}
done
输出文件中主要关注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中_代表终止子，？代表核酸序列中出现了N。
含有barcode和umi的测序数据分析
示例：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 * *
后续分析参考：https://hbox.cool/immune-repertoire/
- 2.3MIXCR
3.使用VDJtools分析结果
基于Java的分析软件，command line运行模式。
使用vdjtools进行免疫组库分析
4.使用immunarch R包进行可视化
兼容mixcr和vdjtools的结果，可以直接可视化。
- 4.1包的简介与安装
- 4.2数据加载
- 4.3处理单细胞配对链数据
- 4.4Basic analysis and clonality
- 4.5Repertoire overlap
- 4.6Gene usage analysis
- 4.7Diversity estimation
- 4.8Track clonotypes across samples
- 4.9Annotate clonotypes
- 4.10Kmer and sequence motif analysis
5.使用R语言手动进行个性化分析
自己写代码实现，可以对免疫组库的结构更加清晰，同时对结果合理性的判断也会更加准确。
- 免疫组库多样性（Richness/Diversity/Clonality）的计算及可视化
- CDR3 length distribution analysis
- 氨基酸占比分析
- Gene usage analysis
- V-J junctions analysis
- CDR3 amino acid motif
详细代码
