零代码计算Contigs与Genes丰度！

在宏基因组分析中，一个最常见的任务就是计算contigs和genes的丰度，这是我们后续定量分析的基础。原理大家都懂，无非就是将reads序列map到contigs或者genes序列上，根据map到的reads数量或碱基数目计算丰度。然而实际操作起来可能是比较麻烦的，也需要自己写一些脚本。今天，我为大家分享一个不需要写代码的contigs和genes丰度计算方法。

宏基因组分析Pipeline

测序数据的解析：Fastq与FastQC

测序数据的质控：Trimmomatic！

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

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

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

免组装宏基因组群落分析

GraPhlAn物种谱可视化

宏基因组编码基因预测

宏基因组bins质量评估

Contigs与genes丰度计算

更新中……

Contigs与genes丰度计算是类似的，只不过目标序列不同，这里以contigs为例。组装所得到的contig的丰度可以通过将全部质控后的reads序列map到拼接结果中，统计落到每个contig中的全部reads数目作为contig的丰度，用每条Contig上映射的reads的总碱基数除以该Contig长度来计算每个Contig的覆盖度。同时全部序列的比对率也是评估拼接结果的一个重要手段，较好的拼接结果通常会有较高的比对率，大于50%。

reads的mapping可通过Bowtie2进行，首先根据拼接的contigs构建新的Index，如下所示：

bowtie2-build --threads 20 contigs.fasta contig

接下来将宏基因组测序的全部reads映射到拼接得到的Contigs上，每个reads至多只能分配到一条Contigs上：

bowtie2 -p 20 -x contig -1 sample_clean_reads_1.fq -2 sample_clean_reads_2.fq -S contig.sam --fast

使用samtools工具将sam文件转化为bam文件：

samtools view -bS --threads 20 contig.sam > contig.bam

对bam文件按照比对的位置坐标对reads进行排序：

samtools sort contig.bam -o contig.reads.sorted.bam --threads 20

此bam文件中便储存了reads的mapping结果，接下来一般是自己写脚本来解析。这里我们也可以借助CheckM来实现。要计算coverage首先需要准备bam的index文件，如下所示：

samtools index contig.reads.sorted.bam

运行结束后生成会生成伴随的index文件contig.reads.sorted.bam.bai，其与对应的sorted bam放在同一路径。CheckM是一个宏基因组bins评估工具，这时候我们可以把所有的contigs序列文件当作一个“bin”放到bins文件夹中。接下来使用CheckM计算coverage：

mkdir bins
cp contigs.fasta bins
checkm coverage -x fasta -m 20 -t 20 bins contigs_coverage.out contig.reads.sorted.bam

结果包含contig序列ID、所在的bin的ID、coverage等信息，如下所示：

当然这里的Bin Id是没有意义的，我们只需要第一列以及Coverage信息即可。接下来可以将reads数转化为相对丰度，也即除以总reads数和序列长度，类似于RNA-seq中的RPKM标准化方法。假如是多样品混合拼接，只需要将每一个样品的reads数据重复上面操作，最后根据contig id进行整合。

软件资源：

Bowtie2:

http://bowtie-bio.sourceforge.net/bowtie2/

Samtools:https://github.com/samtools/

CheckM:参见往期文章基因组质量评估