宏基因组之基因丰度计算

简介: 目前有两种方式可计算宏基因组基因的丰度,一种是基于比对,比如bwa,bowtie,soapaligner等主流的比对软件,另一种是不比对快速估计基因丰度,可以用近俩年来流行的salmon软件,这个软件在转录组的数据比对中也经常用到,可以直接计算出原始的Counts值和标准化的TPM值,此外由于是基于非比对,计算的速度得到很大的提升,同时也节省了部分的内存(少了庞大的sam/bam文件),可以说是多快好省,但是目前的高分文章中也还是不少用基于比对的方法去计算宏基因组的基因丰度的,下面我就分别简单介绍一下基于比对的soapaligner和不比对快速估计的samlon俩个软件的操作流程!!

目前有两种方式可计算宏基因组基因的丰度,一种是基于比对,比如bwa,bowtie,soapaligner等主流的比对软件,另一种是不比对快速估计基因丰度,可以用近俩年来流行的salmon软件,这个软件在转录组的数据比对中也经常用到,可以直接计算出原始的Counts值和标准化的TPM值,此外由于是基于非比对,计算的速度得到很大的提升,同时也节省了部分的内存(少了庞大的sam/bam文件),可以说是多快好省,但是目前的高分文章中也还是不少用基于比对的方法去计算宏基因组的基因丰度的,尤其是Soapaligner比对软件,是由华大公司开发的Soap系列的工具,首先在速度上具有较大的优势,同时也很好的解决了在小内存情况下将短序列比对到参考基因组上的问题,作用也不容小觑。当然也有不少研究人员在前期质控环节进行去宿主也是不错的选择~下面我就分别简单介绍一下基于比对的soapaligner和不比对快速估计的samlon俩个软件的操作流程!!

Option1 SOAPaligenr

5d6e43500c3fce1a838b59a8b355902.png

安装与比对


1. 安装  
conda install soap
2. 构建索引  
2bwt-builder contig.fa 
会产生13个索引文件
3.双端序列进行比对  
soap -p 6 -r 2 -m 200 -x 400 -M 4 -a ./cleandata/B1.clean.fq1.gz -b ./sample1/B1.clean.fq2.gz -D contig.index -o B1_PE.soap -2 B1_SE.soap -u B1_UN.fasta

这个软件的参数有点多,-a,-b分别输入双端序列, -D 索引文件的地址,-o,输出能比对上组装后的参考序列且reads成对存在的比对文件,-2,能比对上组装后的参考序列但reads不成对存在的比对文件;-u,不能比对上组装后的参考序列。至于-m,-x参数比较主观,看过很多文献都有不同的设定,这里贴出一篇文献供参考(DOI:10.3389/fmicb.2017.01546),文献中的参数为−m 4 −r 2 −m 100 −x 1000!!

标准化

看过十几篇文献在基因丰度标准化这一块都引用到了一篇2012年发表在Nature上的一篇文献(qin2012,doi:10.1038/nature11450),看过之后在该文章的补充材料里发现了计算基因丰度的步骤,如下图俩步骤,看过之后感觉其实和转录组中的RPKM/TPM等标准化的思路较为类似,先标准化基因的长度,然后再除以标准化的值的总和。还是挺简单的嘛。

b4cd4c1460e54ad498796e67fd17a4c.png

我们也按照这个思路根据soapaligner比对的结果来进行标准化,在结果文件中我们能得到三个文件,分别为:B1_PE.soap:能比对上组装后的参考序列且reads成对存在的比对文件;B1_SE.soap: 能比对上组装后的参考序列但reads不成对存在的比对文件;sample1_UN.fasta:不能比对上组装后的参考序列,通常为fasta格式;我们只需要比对上的文件B1_PE.soap和B1_SE.soap即可。这两个文件的第8列中比对上的基因,我们将其取出来sort|uniq -c 就可以计算出基因的原始copy数目,也就是比对上特定基因的数目。之后再计算出比对基因的长度,也就是你组装去冗余之后的非冗余基因集的各基因长度,这步可以写个小脚本或者用现成的轮子seqkit/seqtk都可以很快的解决,最后利用上面的公式进行标准化就行了~  还需要注意的一点,没有比对上的基因我们就用0代替他的原始丰度就好了。

Option2 Samlon

这种方法就更加简单了,也是大致两步走,建立索引,序列比对,只不过时间将大大缩短,这种方法网上的教程有很多,我这里贴出我的代码大家可以比较一下使用~ 建议大家使用时候都安装最新版本的软件,github地址为:https://github.com/COMBINE-lab/salmon/releases/


#1.建立索引
salmon index -t B1_NR100nl.fasta -p 9 -k 31 -i ./index
#2.比对
salmon quant --validateMappings -i ./index -l A -p 3  --meta -1 ../clean1_${b}.fastq -2 ../clean2_${b}.fastq -o ${b}.quant #基因定量

注意的就是宏基因组数据加上-p meta这个参数,其他参数默认就好~,得到的结果有个quant.sf文件,丰度信息就在里面。为5列的信息,第4列就是我们想要的标准化后的基因丰度了

1b1e254bb3f6b93a2e34bae74830afb.png

OK,拿到基因丰度文件之后,我们就可以用来进行后续分析,比如计算eggnog数据库比对后的功能基因丰度~这个留到下次再说吧!!

相关文章
|
5月前
|
安全
三维基因组|基因组结构 (2)
三维基因组|基因组结构 (2)
55 0
|
5月前
R语言偏最小二乘回归PLS回归分析制药产品化学制造过程数据、缺失值填充、变量重要性
R语言偏最小二乘回归PLS回归分析制药产品化学制造过程数据、缺失值填充、变量重要性
|
5月前
|
芯片
基因测序的原理是什么
基因测序的原理是什么
|
10月前
微分方程——药物在体内分布的房室模型
微分方程——药物在体内分布的房室模型
124 0
|
算法 芯片
DNA测序原理:illumina和Pacbio对比介绍
DNA测序原理:illumina和Pacbio对比介绍
|
存储 Python
候选基因如何分析?
候选基因如何分析?
|
机器学习/深度学习 安全 数据挖掘
用于 DNA 测序的机器学习模型,理论上可以解码任何测序读数中所编码的数据值
用于 DNA 测序的机器学习模型,理论上可以解码任何测序读数中所编码的数据值
163 0
用于 DNA 测序的机器学习模型,理论上可以解码任何测序读数中所编码的数据值
|
数据采集
宏基因组数据 数据预处理
宏基因组下载数据过程
278 2
|
Shell Python Perl
宏基因组之基因预测
12年有篇BMC的文献对几款预测的软件做了评估,其实参考大多数的文献中最常见的俩个软件也就是Prodigal和Metagenemark这俩个软件,分析过程中我这俩个软件都感受一下,现在将过程记录一下~~有兴趣的话可以看看这篇文献哦。
241 0
|
算法 索引 Python
宏基因组之基因组装
宏基因组组装,即把短的reads拼装成连续的序列contig,再根据PE等关系将contig拼装成scaffold。与单个基因组组装不同,宏基因组组装最终得到的是环境样品中全部微生物的混合scaffold。理想情况下一条scaffold对应一个物种的全基因组。但由于序列太短或者覆盖度不够,很难拼出一条完整的基因组。针对高通量测序数据,出现了多种拼接算法和软件。
499 0