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

Option1 SOAPaligenr

图片.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等标准化的思路较为类似，先标准化基因的长度，然后再除以标准化的值的总和。还是挺简单的嘛。

图片.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列就是我们想要的标准化后的基因丰度了

图片.png

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

下面是刚创建个人公众号，会定时更新R、linux、python，组学方面的学习内容，请多多支持呦~~

image