很多情况下由于内存限制等原因无法将多个样本混合在一起拼接,这样每个样品单独拼接、预测获得的基因集在合并分析时可能会有很多冗余。要构建多个样品、多个项目的非冗余基因集,需要根据一定的相似度阈值对不同样品的基因序列进行聚类。常用的软件有CD-HIT、MMseqs、Linclust等。
宏基因组分析Pipeline
宏基因组基因集去冗余:CD-HIT
更新中……
CD-HIT(http://weizhongli-lab.org/cd-hit/)是一个广泛使用的蛋白或核酸序列比较聚类工具,其将所有序列按照参数设定进行聚类,并将每一组聚类中的最长序列作为代表序列进行输出,同时给出每组聚类下的每个序列名可供相似度分析使用。CD-HIT的基本思路是首先对所有序列按照其长度进行排序,然后从最长的序列开始,形成第一个序列类,然后依次对序列进行处理,如果新的序列与已有的序列类的代表序列的相似性在cutoff以上则把该序列加到该序列类中,否则形成新的序列类。
CD-HIT速度快主要是两个方面的原因:一个是使用了word过滤方法,即如果两条序列之间的相似性在80%(假设序列长度为100),那么它们至少有60个相同的长度为2的word,至少有40个相同的长度为3的word,至少有20个相同的长度为4的word。基于这个原则,在处理新的序列的时候,如果新的序列与已有序列的相同word的长度能满足这些要求则不需要进行比对了,这极大的降低了时间消耗;另外一个速度快的原因是使用了index table,可以很快的计算序列之间相同word的数目。
尽管CD-HIT速度快使用方便,但是也需要注意其缺点:
①它不能保证同一个序列类中的序列的相似性都在threshold之上,因为每次比对都是用新序列与序列类的代表序列进行,这就有可能使得序列类中除了代表序列外其他序列之间的相似性在threshold之下。比如A是代表序列,B与A的相似性大于0.95,C与A的相似性也大于0.95,但是这并不能保证B与C的相似性也大于0.95。
②它不能保证一个序列类的病毒与另外一个序列类中的病毒的相似性也在threshold之上,原因还是在于用代表序列代表了整个序列类。
③基于word filter的方法使得使用每个长度的word能够处理的冗余性水平有限,如使用长度为2的word只能够得到相似性在50%以上的序列,长度为3的word只能够得到相似性在66.7%以上的序列类,类似的,长度为5的word只能够得到相似性在80%以上的序列。在实际应用的时候需要注意选择的word长度与threshold的匹配。
CD-HIT可以在GitHub下载,安装方法如下所示:
CD-HIT有两个主程序:
cd-hit:(cd-hit-est)将相似的蛋白聚类成聚类簇。
cd-hit-2d:(cd-hit-est-2d)比较两个数据库,并识别数据库2中与数据库1相似的序列。
cd-hit的命令参数如下所示:
下面以6个宏基因组为例进行分析(一般设置95%的identity、90%的短序列coverage)[1]:
运行结束后,会生成非冗余的fasta序列文件,以及后缀为clstr的聚类簇列表,如下所示:
其中标*的为代表序列,可以将代表序列挑取出来进行后续分析。
参考文献:
[1] NielsenH B, Almeida M, Juncker A S, et al. Identification and assembly of genomes andgenetic elements in complex metagenomic samples without using referencegenomes[J]. Nature biotechnology, 2014, 32(8): 822-828.
END