为了更精确的定量：宏基因组gene丰度分析工具的比较

之前常有做宏基因组的朋友问我，为什么他们计算基因丰度获得的结果中，有些基因的丰度为零。理论上所有的contig序列均由reads拼装而得，而基因作为contig序列上一个区域，不该没有reads比对上。其实，这些丰度为零的基因反映了宏基因组gene丰度计算一个很容易犯的错误。非常抱歉的是，我在之前文章零代码计算contigs与genes丰度一文中，并没有及时认识并纠正这个错误，现在亡羊补牢，希望对大家能有所帮助！

宏基因组分析Pipeline

测序数据的解析：Fastq与FastQC

测序数据的质控：Trimmomatic！

宏基因组reads筛选：去除宿主序列

测序数据的组装：常用软件工具

宏基因组多样品的混合组装

Contigs与genes丰度计算

宏基因组gene丰度不同工具比较

免组装宏基因组群落分析

GraPhlAn物种谱可视化

宏基因组编码基因预测

宏基因组bins质量评估

宏基因组binning: Metabat

更新中……

在宏基因组项目中，微生物的DNA被随机打断成短的序列片段，经过测序获得reads，然后使用组装工具将reads进行拼接复原DNA序列，也即contig序列。要想获得contigs的测序深度（即depth），我们可以使用Bowtie2将reads回帖到contigs，再根据回贴情况分析平均depth。

基因序列为contigs上预测的编码区域。因此理论上，同一条contig序列上的genes应该具有相似的深度。由于在contigs回贴时，我们很容易根据bam文件获得一条contig不同区域的depth，再结合gene的起止位置，我们可以获得一条contig上所有genes的真实depth信息：

图中红线为该contig的平均depth。可以看到，虽然gene的depth并不严格等于contig的平均depth，但总体上围绕着该值呈近似的高斯分布。如果你具备小小的编程基础，完全可以根据这些depth信息进行标准化，从而获得基因的丰度。不过，可能很多朋友直接对gene进行回贴：

直接对gene进行回贴的结果出人意料，获得的gene深度呈现了一定的长度依赖性，长度超过1000bp的gene深度比较正常，而1000bp以下的gene其depth随着长度的减小快速降低，甚至有些短gene的深度为零。实际上，很大一部分的原核生物的基因长度小于1000bp，这样获得的depth信息会给后续的分析带来很大的困扰。

究其原因，Bowtie2默认为ene-to-end的回帖方法，也即要求reads必须全部比对到gene序列，而gene仅为contig中间区域，gene两端部分匹配的reads会被丢弃，导致回帖率降低。极端情况下，假如gene的长度小于reads长度，其depth肯定为零。根据这个原理，我们可以推断gene的depth与contig的depth关系：

第二幅图中的红色曲线即为上述公式的图像，可以看到实际情况正好符合我们的推测。Bowtie2也支持local的比对方式（添加参数--local），获得的结果如下所示：

可以看到，结果大为改善，但仍不甚理想，可以明显看到短gene深度的降低。接下来我们试一下其他回贴工具，首先是BBmap：

结果和Bowtie2的默认模式一样糟糕，接下来是BWA的MEM算法，该算法也支持剪接比对，结果如下所示：

可以看到，结果非常完美，与我们自己计算获得的真实depth几乎一样，我们可以比较BWA-MEM和Bowtie2-local的结果：

可以看到，BWA-MEM的比对结果要大大好于Bowtie2-local，能很好地还原gene的实际depth结果。因此，假如你一定要采用gene回贴的方法计算gene的depth和丰度，非常推荐BWA-MEM。

最近也有人推荐使用一款基于非比对方法的快速gene丰度计算工具Salmon，该工具来自转录组和宏转录组领域，速度非常快，它在宏基因组的结果如下：

该工具可以直接给出标准化后的TPM，因此受到很多人追捧。然而其对宏基因组的分析结果非常诡异，短gene的丰度急剧升高，甚至超过正常水平一个数量级。因此，十分不建议在宏基因组的分析中使用Salmon。

当然，一条contig序列的depth不一定是均匀的。当我们提取DNA时，一些微生物的DNA可能正在复制中，此时复制起点附近已经复制完毕，而复制终点附近还未被复制。由于大多数细菌是典型的θ型复制，因此可能出现一个细菌基因组中一半depth较高、一半depth较低：

这也是一条contig上gene丰度产生变化的主要原因之一。

附BWA-MEM的安装方法：

git clone https://github.com/lh3/bwa.git
cd bwa
make

使用方法：

#首先对参考序列构建idex：
bwa index -p 83_armatimo_gene 83_armatimo_gene.fna
#使用BWA-MEM进行比对：
bwa mem -t 20 83_armatimo_gene 83_clean_1.fq.gz 83_clean_2.fq.gz
bwa mem -t 20 83_armatimo_gene 83_clean_1.fq.gz 83_clean_2.fq.gz > 83_armatimo_gene.sam
#接下来的处理与bowtie2相同。

—END—