Closed ixxmu closed 2 years ago
网罗了一大波生物信息学数据分析方面的工程师,是时候官宣咱们的ngs数据处理业务列表:
扫描下面二维码即可添加微信咨询!
如果是一个免疫组库项目,走igblastn或者MiXCR流程只需要800元人民币即可,但是如果要提供结果解读服务或者个性化的统计可视化,需要具体谈。都会提供全套代码。不管是公共数据集还是你自己的实验测序数据,一样的费用!我们会代替你跑如下所示的流程:
前面我带领大家通过IMGT数据库认知免疫组库,而且也一起从IMGT数据库下载免疫组库相关fasta序列,免疫组库重要的研究对象就是分成BCR的IGH,IGK,IGL这3类,以及TCR的TRA,TRB,TRD,TRG,它们各自都有V,D(可选),J,C基因。
已经预告了有一个免疫组库的实战,现在终于有时间来带领大家搞定它。
前面我们提到了,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 也可以做同样的分析。
同样的,需要熟悉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 路径文件:
脚本如下:
# 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进行讲解,因为这个示例文章也是这个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文件失败,不过我们的样本数量太多了,而且这个免疫组库系列教程反正看的人不多,我就懒得去检查失败的几个样本是什么问题了。(有可能是批量下载失败, 因为写这个教程的时候我人在中国大陆地区)
完全参考:使用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软件压缩包:
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
非常复杂,需要单独写一个推文。
也非常复杂,需要单独写一个推文。
更加复杂,需要单独写多个推文。
网罗了一大波生物信息学数据分析方面的工程师,是时候官宣咱们的ngs数据处理业务列表:
扫描下面二维码即可添加微信咨询!
如果是一个免疫组库项目,走igblastn或者MiXCR流程只需要800元人民币即可,但是如果要提供结果解读服务或者个性化的统计可视化,需要具体谈。都会提供全套代码。不管是公共数据集还是你自己的实验测序数据,一样的费用!我们会代替你跑如下所示的流程:
前面我带领大家通过IMGT数据库认知免疫组库,而且也一起从IMGT数据库下载免疫组库相关fasta序列,免疫组库重要的研究对象就是分成BCR的IGH,IGK,IGL这3类,以及TCR的TRA,TRB,TRD,TRG,它们各自都有V,D(可选),J,C基因。
已经预告了有一个免疫组库的实战,现在终于有时间来带领大家搞定它。
前面我们提到了,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 也可以做同样的分析。
同样的,需要熟悉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 路径文件:
脚本如下:
# 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进行讲解,因为这个示例文章也是这个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文件失败,不过我们的样本数量太多了,而且这个免疫组库系列教程反正看的人不多,我就懒得去检查失败的几个样本是什么问题了。(有可能是批量下载失败, 因为写这个教程的时候我人在中国大陆地区)
完全参考:使用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软件压缩包:
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
非常复杂,需要单独写一个推文。
也非常复杂,需要单独写一个推文。
更加复杂,需要单独写多个推文。
https://mp.weixin.qq.com/s/-VK20oj473Kz7TUHVuaUnA