一个免疫组库的实战

前面我带领大家通过IMGT数据库认知免疫组库，而且也一起[从IMGT数据库下载免疫组库相关fasta序列](https://mp.weixin.qq.com/s/zPN2F3aKxeqDtwuWvClJqg)，免疫组库重要的研究对象就是分成BCR的IGH,IGK,IGL这3类，以及TCR的TRA,TRB,TRD,TRG，它们各自都有V，D（可选），J，C基因。

已经预告了有一个免疫组库的实战，现在终于有时间来带领大家搞定它。

来自于文章；https://www.tandfonline.com/doi/full/10.1080/2162402X.2019.1644110
数据：https://www.ncbi.nlm.nih.gov/bioproject/PRJEB33490

熟读文献，理解免疫组库背景

前面我们提到了，BCR有IGH,IGK,IGL这3类，而TCR有TRA,TRB,TRD,TRG，它们各自都有V，D（可选），J，C基因。

而本研究采取多重PCR来扩增TCR里面的TRB，测序采用的是 Illumina MiSeq sequencer (600–cycle, single indexed, paired-end run). 数据分析使用的是MiXCR这个工具。

主要的关注点也是 random V(D) J recombination 产生的CDR3多样性，其中下游分析使用 tcR 包，当然了，其它优秀R包，比如 VDJtools 也可以做同样的分析。

多重PCR来扩增TCR里面的TRB

下载文献的免疫组库测序数据

同样的，需要熟悉GEO和SRA数据库：

获得文献里面的数据集里面的样本的数据库里面的ID列表：

$head id.list
ERR3445007
ERR3445008
ERR3445009
ERR3445010
ERR3445011
ERR3445012
ERR3445013
ERR3445014
ERR3445015
ERR3445016

使用下面脚本批量下载：

mkdir -p ~/data/public/TRB
cd ~/data/public/TRB
cat id.list |while read id;do ( ~/biosoft/sratoolkit/sratoolkit.2.8.2-1-centos_linux64/bin/prefetch $id -X 100G );done

注意：这个sratoolkit版本好像经常有问题，根据 id.list里面存储的id批量下载后的sra文件夹，存储到sra文件夹，文件大小如下：

$ls -lh sra/*|cut -d" " -f 5-|head
 23M Apr 25 17:57 sra/ERR3445007.sra
 27M Apr 25 17:57 sra/ERR3445008.sra
 20M Apr 25 17:58 sra/ERR3445010.sra
 40M Apr 25 17:58 sra/ERR3445011.sra
 20M Apr 25 17:59 sra/ERR3445012.sra
 19M Apr 25 18:00 sra/ERR3445014.sra
 44M Apr 25 18:04 sra/ERR3445016.sra
 32M Apr 25 18:04 sra/ERR3445017.sra
 20M Apr 25 18:05 sra/ERR3445018.sra
 31M Apr 25 18:06 sra/ERR3445019.sra

可以看到免疫组库测序数据的文件都很小哦。如果你的sratoolkit有问题，或者prefetch命令下载sra文件速度太慢，可以参考：使用ebi数据库直接下载fastq测序数据 , 需要自行配置好，然后去EBI里面搜索到的 fq.txt 路径文件：

项目地址是： https://www.ebi.ac.uk/ena/browser/view/PRJEB33490

脚本如下：

# conda activate download 
# 自己搭建好 download 这个 conda 的小环境哦。
cat fq.txt |while read id
do
ascp -QT -l 300m -P33001 \
-i ~/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh \
era-fasp@$id .
done
# nohup bash step1-aspera.sh 1>step1-aspera.log 2>&1 &

这个脚本会根据你在EBI里面搜索到的 fq.txt 路径文件，来批量下载fastq测序数据文件。

两个下载方式都可以，取决于你的网络情况，不同国家地区不同的选择，人啊要学会灵活应变！有一些学员在b站发弹幕吐槽我前面的ngs教学课程的sratoolkit代码失败，其实是它们自己的网络问题，我的课程录制较早，那个时候也没有想到过大家居然连数据都无法下载，这一点只能说请大家见谅。

我们的最新教程，都是走：使用ebi数据库直接下载fastq测序数据 , 这个教程来直接下载fastq文件啦。

然后把sra文件，批量作为fastq文件，这个时候，可以加入样本信息，所以需要制作sra文件和样本对应表格，如下：

$head id2sample.txt
healthy1 ERR3445007
healthy2 ERR3445008
healthy3 ERR3445009
healthy4 ERR3445010
healthy5 ERR3445011
healthy6 ERR3445012
healthy7 ERR3445013
healthy8 ERR3445014
healthy9 ERR3445015
healthy10 ERR3445016

log日志如下：

healthy1 ERR3445007
Read 72065 spots for sra/ERR3445007.sra
Written 72065 spots for sra/ERR3445007.sra

healthy2 ERR3445008
Read 87649 spots for sra/ERR3445008.sra
Written 87649 spots for sra/ERR3445008.sra

healthy4 ERR3445010
Read 76007 spots for sra/ERR3445010.sra
Written 76007 spots for sra/ERR3445010.sra

healthy5 ERR3445011
Read 203464 spots for sra/ERR3445011.sra
Written 203464 spots for sra/ERR3445011.sra

得到的测序文件如下：

$ls -lh raw_fq/*|cut -d" " -f 5-|head
 18M Apr 27 19:05 raw_fq/healthy10_1.fastq.gz
 24M Apr 27 19:05 raw_fq/healthy10_2.fastq.gz
 15M Apr 27 19:06 raw_fq/healthy11_1.fastq.gz
 18M Apr 27 19:06 raw_fq/healthy11_2.fastq.gz
7.7M Apr 27 19:06 raw_fq/healthy12_1.fastq.gz
 11M Apr 27 19:06 raw_fq/healthy12_2.fastq.gz
 12M Apr 27 19:06 raw_fq/healthy13_1.fastq.gz
 16M Apr 27 19:06 raw_fq/healthy13_2.fastq.gz
 11M Apr 27 19:06 raw_fq/healthy14_1.fastq.gz
 14M Apr 27 19:06 raw_fq/healthy14_2.fastq.gz

可以看到，每个样本还是有两个测序数据的，意味着是双端测序。前面其实也提到了的，测序采用的是 Illumina MiSeq sequencer (600–cycle, single indexed, paired-end run)。而且他们的样本的数据量很稳定。

如果你参考使用ebi数据库直接下载fastq测序数据，那么就无需经过sra文件转fastq啦，本来就是下载fastq.gz 文件。如下所示：

 7.8M May 23 09:43 ERR3445007_1.fastq.gz
 11M May 23 09:43 ERR3445007_2.fastq.gz
 9.8M May 23 09:43 ERR3445008_1.fastq.gz
 14M May 23 09:43 ERR3445008_2.fastq.gz
 7.0M May 23 09:44 ERR3445009_1.fastq.gz
 10M May 23 09:44 ERR3445009_2.fastq.gz
 7.2M May 23 09:44 ERR3445010_1.fastq.gz
 9.3M May 23 09:44 ERR3445010_2.fastq.gz

不管是采用何种公共数据的下载方式，最后都是拿到fastq文件，走免疫组库分析流程哦。这个文章：We used next-generation immunosequencing to study T cell receptor (TCR) repertoire metrics on 346 blood samples from healthy individuals and cancer patients producing a dataset with around 8.8 million TCR reads.

而且免疫组库的每个样本数据量很小，几百个样本合起来的数据量都比不上一个转录组测序。

构建免疫组库文件夹及软件环境

这个时候我们有3个策略：

MiXCR is a bioinformatic tool that processes B- or T-cell immune repertoire data from raw sequences to quantified clonotypes
IgBLAST can analyze both BCR and TCR. This is done by annotating V, D, J regions, CDR3 identification and mutational analysis of the variable regions using the BLAST search algorithm.
IMGT/HighV-QUEST identifies the V, D and J genes and alleles by alignment with the germline receptor gene and allele sequences of the IMGT germline database

我们都写了教程，如下：

使用igblast进行免疫组库分析
使用MiXCR进行免疫组库分析
因为IMGT/HighV-QUEST是网页工具，我们略过

这里我们选择MiXCR进行讲解，因为这个示例文章也是这个MiXCR流程。而且MiXCR是MiLaboratory开发的，他们实验室出品了3款：MiXCR, MIGEC, VDJtools ，都是值得推荐的软件哦。

免疫组库测序数据的过滤

我们首先认识了免疫组库测序数据，知道了免疫组库测序数据的一些特性，并不能简单的看fastqc报告，就判断我们的免疫组库测序数据不合格。

不过，任何测序数据的过滤都是有一个共性，去除低质量序列，以及去除接头，我们这里仍然是选择走 trim_galore 流程。参考：使用igblast进行免疫组库分析

source activate qc 
# conda install -c bioconda flash
mkdir -p raw clean merge blast MiXCR
cat config |while read id;do (trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -o clean raw/${id}*gz);done
# 这里的config文件里面是本次免疫组库文献里面的数据样本ID

过滤之后的fq文件如下：

 7.1M May 28 10:26 clean/ERR3445007_1_val_1.fq.gz
 8.4M May 28 10:26 clean/ERR3445007_2_val_2.fq.gz
 9.0M May 28 10:26 clean/ERR3445008_1_val_1.fq.gz
 11M May 28 10:26 clean/ERR3445008_2_val_2.fq.gz
 6.4M May 28 10:27 clean/ERR3445009_1_val_1.fq.gz
 7.6M May 28 10:27 clean/ERR3445009_2_val_2.fq.gz
 6.7M May 28 10:27 clean/ERR3445010_1_val_1.fq.gz
 8.2M May 28 10:27 clean/ERR3445010_2_val_2.fq.gz
 14M May 28 10:27 clean/ERR3445011_1_val_1.fq.gz
 16M May 28 10:27 clean/ERR3445011_2_val_2.fq.gz

然后把PE数据使用FLASH合并，仍然是参考：使用igblast进行免疫组库分析，代码如下：

# flash raw/ERR3445170_1.fastq.gz raw/ERR3445170_2.fastq.gz -M 250 -o clean
# conda install -c bioconda flash
cat config |while read id;do (flash clean/${id}*gz -M 250 -o merge/${id} );done

得到的文件如下：

 39M May 28 11:37 ERR3445007.extendedFrags.fastq
 47M May 28 11:37 ERR3445008.extendedFrags.fastq
 34M May 28 11:37 ERR3445009.extendedFrags.fastq
 41M May 28 11:37 ERR3445010.extendedFrags.fastq
 86M May 28 11:37 ERR3445011.extendedFrags.fastq
 34M May 28 11:37 ERR3445012.extendedFrags.fastq
 30M May 28 11:37 ERR3445013.extendedFrags.fastq
 40M May 28 11:37 ERR3445014.extendedFrags.fastq
 54M May 28 11:37 ERR3445015.extendedFrags.fastq
 91M May 28 11:37 ERR3445016.extendedFrags.fastq

因为igblast比对接受fa文件输入，所以需要首先把fq进行格式转换，这里我们采用seqkit软件的fq2fa子命令哈：

cd merge
# conda install -c bioconda seqkit 
ls *.extendedFrags.fastq|while read id;do ( seqkit fq2fa $id > ${id%%.*}.fa );done

得到文件如下：

 22M May 28 11:44 ERR3445007.fa
 26M May 28 11:44 ERR3445008.fa
 19M May 28 11:44 ERR3445009.fa
 23M May 28 11:44 ERR3445010.fa
 49M May 28 11:44 ERR3445011.fa
 19M May 28 11:44 ERR3445012.fa
 17M May 28 11:44 ERR3445013.fa
 22M May 28 11:44 ERR3445014.fa
 30M May 28 11:44 ERR3445015.fa
 52M May 28 11:44 ERR3445016.fa

有意思的是，中间有几个fq文件转fa文件失败，不过我们的样本数量太多了，而且这个免疫组库系列教程反正看的人不多，我就懒得去检查失败的几个样本是什么问题了。（有可能是批量下载失败，因为写这个教程的时候我人在中国大陆地区）

走igblastn流程进行比对

完全参考：使用igblast进行免疫组库分析，需要自己配置好igblast软件，以及下载好数据库文件并且构建好索引。跑代码就非常简单了：

cd blast/
cat ../config |while read id;
do 
echo $id 
igblastn \
-germline_db_V ~/biosoft/igblast/database/human_TRBV.fa \
-germline_db_D ~/biosoft/igblast/database/human_TRBD.fa \
-germline_db_J ~/biosoft/igblast/database/human_TRBJ.fa \
-domain_system imgt \
-organism human \
-ig_seqtype TCR \
-auxiliary_data ~/biosoft/igblast/optional_file/human_gl.aux \
-show_translation \
-outfmt 7 \
-num_threads 4 \
-query ${id}.fa \
-out ${id}.blast.output.7
done

这样就批量把免疫组库测序数据过滤整理好的fa文件比对到了IMGT数据库的TRB的V,D,J基因啦。

也可以走MiXCR流程

完全参考：使用MiXCR进行免疫组库分析，首先下载解压MiXCR软件压缩包：

mkdir -p ~/biosoft/MiXCR/
cd ~/biosoft/MiXCR/
wget https://github.com/milaboratory/mixcr/releases/download/v3.0.13/mixcr-3.0.13.zip
unzip mixcr-3.0.13.zip

然后就可以使用MiXCR软件，走免疫组库分析流程，只需要是clean的reads组成的fq文件即可，MiXCR软件会自动检测属于BCR的IGH,IGK,IGL这3类，还是TCR的TRA,TRB,TRD,TRG，的V，D（可选），J，C基因。

批量运行的代码如下：

cd MiXCR/ 
mixcr=${HOME}/biosoft/MiXCR/mixcr-3.0.13/mixcr
echo $mixcr 
# conda activate qc
cat ../config |while read id;
do 
echo $id 
ls ../clean/${id}*gz 

$mixcr analyze amplicon \
 -s hsa \
 --starting-material dna \
 --5-end v-primers --3-end c-primers \
 --adapters adapters-present \
 ../clean/${id}*gz $id 
done

免疫组库下游分析

解析igblastn比对结果

非常复杂，需要单独写一个推文。

解析MiXCR流程的结果文件

也非常复杂，需要单独写一个推文。

高级可视化技巧

更加复杂，需要单独写多个推文。

明码标价

如果是一个免疫组库项目，走igblastn或者MiXCR流程只需要800即可，但是如果要提供结果解读服务或者个性化的统计可视化，需要具体谈。

生信菜鸟团

欢迎去论坛biotrainee.com留言参与讨论，或者关注同名微信公众号biotrainee