宏基因组分析初识-基本分析流程

士后阶段从湿实验转到干实验后，做事方式改变了很多，主要表现在做事目标开始变得十分明确，不再佛系。“学书未成先习剑，用剑无功再读书”，大概最能描述我现在的做事节奏。年龄大了，觉得不能再慢悠悠的一点一点打基础咯。

所以，开始第一次分析宏基因组数据前，我可以说没有完整的看过一篇文献。而是直接找了各个测序公司的项目结题报告，快速了解整个过程是怎样后，再根据自己的判断确定按照哪些流程来做分析，从实战中再去学习。

初步尝试的分析流程及注意事项总结如下：

1.原始数据质量检测：

文件完整性判断：拿到测序文件后，需要看文件是否完整，因为文件在拷贝或下载过程中可能会出现问题。当然这个不是用肉眼来看的。测序公司返回的数据除了压缩的reads数据文件文件外，还会给一个含有文件校验和的文本文件（我拿到的文件都是用的md5校验和）。需要使用md5sum -c *.txt来查看所得文件是否完整。

Clean reads文件（fastq格式）质量检测：这里一般使用fastqc对得到的fastq文件进行检测，根据Q20、Q30、adaptor、barcode等的结果判断是否还需要自己进行质控与序列修剪。我拿到过测序公司返回的很多fastq文件，它们的质控都做的非常好，基本不要自己再做处理。不过，以防万一，fastqc看看clean reads的质量仍然是必须的。

2.宏基因组组装：

能用于宏基因组组装的软件有很多，如CLCwork bench，metavelvet，megahit，ibda-ud，spades等。通过文献调研，我最终选用了spades，其实想要做好的话，最好是比较两到三个软件的拼接效果，择其善者而用之。spades的用法非常简单，可用于单菌拼接和宏基因组拼接。不过有几点需要注意：对于宏基因组拼接，它不支持--careful参数和--cov-cutoff参数。所以对于宏基因组拼接，它的命令和参数其实是十分固定的；此外，官方manual建议，对于PE150的测序策略，建库长度最好是在350-500bp，组装kmer长度建议为21,33,55,77，对于PE250的测序策略，建库长度最好是在550-700bp，组装kmer长度建议为21,33,55,77，99,127；另外，spades拼接完成后没有拼接效果的总结文件，这时候需要借助软件quast对拼接结果进行统计与评价，包括统计N50，总长度，GC含量等。Spades没有对contigs长度做过滤的参数，所以它拼接得到的fasta结果文件一般会比其他软件大，这样的好处是后续做物种分析，损失的信息会很少。如果自己对contigs最低长度有要求，可以写个程序脚本提取fasta结果文件中大于该长度的序列。

3.orf预测与基因组分箱：

拼接完成后，可以同时做两件事:orf预测与基因组分箱。

orf预测的软件也不少。我查阅的各个测序公司的结题报告中，用MetaGene的最多，但是这个软件输出结果单一，只有一个文件，包含orf在contigs上的起始位置，正负链，是基于数据库预测还是密码子直接预测等信息。后续的orf序列提取和翻译需要自己写脚本完成。我尝试着做了一遍，很是麻烦，而且存在orf起始密码子或者终止密码子不全的问题。如果编程零基础，不建议使用这个软件。最终，我选用了prodigal来预测orf。Prodigal预测结果文件十分丰富，包括orf的核酸序列文件、氨基酸序列文件、gff文件等，需要注意的是用于宏基因组orf预测的时候得加-pmeta参数。

基因组分箱：也就是常说的genomebinning，这一步很多公司的结题报告中都不会去做。但是其实对于做功能研究的大佬们来说，这一步非常有必要，可能能从乱七八糟的一堆序列中发现有意思的darkmatter。如果有很强分子生物学和生化背景，后期的验证也不会很难。这一步我用了maxbin，输入文件除了contigs外，还可以加上reads文件，这样输出结果除了有每个bin的完整度、序列GC含量以及maker基因外，还会包含每个bin在样品中的相对丰度信息。接下来可以对完整度较高的bin做功能基因组分析、代谢途径重构等。

4.非冗余基因集的构建：

将预测出来的orf序列，使用CD-HI T软件，进行聚类（聚类相似度≥0.95，覆盖度≥0.9），每个类取最长的orf序列作为代表序列，构建非冗余基因序列集合（感觉类似于扩增子分析中得到的reference序列文件）。这个软件比较简单，根据按照软件手册输入命令就可以。

5.基因丰度计算：

得到非冗余基因集后，可以将cleanreads回帖比对到非冗余基因集上，计算每个基因在各样品中的丰度。可用的软件有bedtools、SOAPAligner等，不过使用较多的是SOAPAligner。SOAPAligner是BGI开发一款快速序列比对软件，很适合reads的快速回帖比对。但是这个软件好像还有不少bug，甚至参数先后顺序不对都会导致程序运行失败，另外，它对bash的原生命令支持也不好。试了很多次，下面这样的顺序可以运行成功：soap -D数据库索引文件-a reads1.fq -b reads1.fq-o输出文件1（双端reads比对上的orf及其reads书）-2输出文件2（仅单端reads比对上orf）-u未比对上的orf序列id -p线程数。对我们计算基因丰度有用的文件是输出文件1和输出文件二。但是具体怎么计算，网上的各种博客什么都避而不谈。我找了几十篇文献才在原始文献的supplementary material中追溯到具体方法（doi:10.1038/nature11450）。但是，进行基因丰度计算前，我们先要写个程序脚本从结果文件中统计每个orf比对上的reads数以及每个orf的长度。如果编程零基础，可以尝试用其他软件来进行基因丰度统计。

6.基因注释：

得到基因丰度表后，需要知道这些基因的功能，才能分析生物样品表型与基因之间的关系，这就对基因进行功能注释。简言之，基因功能注释是将基因序列比对回数据库，根据序列相似度预测基因功能。目前最常用的比对软件是DIAMOND，其对短序列（reads）的比对速度比blast快两万倍，比对结果和blast的重复度有80~90%。对于长序列，需要采取sensitive模式，比对速度比blast快两千多倍，重复度达94%以上。需要注意的是diamond只能用氨基酸序列建库。常用的基因注释库包括NR，CAZy，KEGG，GO，COG等。当然，如果要分析特定基因，也可以用其他库。注释完成后，我们可以知道功能基因在样品中的丰度，从而推测功能基因在对各样品表型的贡献，甚至发现基因的新功能。此外，我觉得功能基因的丰度信息可以为接下来做宏转录组进行初步的筛选与指导。

7.物种丰度分析：

除了根据基因丰度分析基因功能外，还可以根据物种丰度分析特定物种的功能。

物种丰度分析有一个神器---metaphlan，必须推荐。它可以直接以reads为输入文件，得到物种的丰度信息并画出样品间物种丰度聚类热图。其姐妹软件包括graphlan（样品内）等。Metaphlan依赖数据库为软件自带数据库。

另外，非冗余基因集比对后的daa结果文件可以输入到megan软件，并通过NR库对应的分类学信息数据库获得物种注释，然后使用物种对应的基因丰度总和计算该物种的丰度（这一步总觉得有问题，未找到原始文献）。

可以将两种方法做出的物种丰度结果作比较，评估结果。

8. 其他：

以上七步应该是基本的分析步骤，也就是说，拿到原始数据了以后，先将上面的步骤过一遍，对自己的数据有个预判。然后再结合自己的样品特征、生物学表型数据，做进一步分析，比如PCA因子关联、网络分析等，将数据信息与生物学功能联系起来，这才是最终目的。