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