上次我们用prodigal对contig进行了ORF(开放阅读框)预测,可见11.windows ubuntu 子系统 contig ORF(开放阅读框)预测。-CSDN博客,接下来我们便开始构建非冗余基因集。
构建非冗余基因集有几个重要的原因:
-
简化分析: 在大规模基因组学研究中,存在大量的基因序列。去除冗余可以简化数据集,减少分析的复杂性和计算的资源消耗。
-
提高准确性: 去除冗余基因集可以减少重复的信息,从而提高分析的准确性。每个基因在非冗余基因集中只出现一次,避免了重复计算和统计分析中的偏差。
-
便于比较和注释: 非冗余基因集更容易进行基因功能注释和进化比较分析。每个基因都是唯一的,使得对其功能、结构和进化关系的研究更加清晰和准确。
上次,我们使用prodigal -i contigs.fasta -o prodigal/a1.gff -d prodigal/a1.fasta -a prodigal/a1pro.fasta -p meta ,得到
其中a1pro.fasta指定输出的蛋白质序列文件名,用于存储预测出的蛋白质序列,我们用来构建非冗余基因集。我们用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-腾讯云开发者社区-腾讯云 (tencent.com)
#安装CD-HIT
conda install -c bioconda cd-hit -y
#CD-HIT基础应用
cd-hit -i a1pro.fasta -o qurongyuprotein.fasta -c 0.95 -M 640000 -T 24 -n 5 -d 0 -aS 0.9 -g 1 -sc 1 -sf 1
-i a1pro.fasta
:输入文件,即包含蛋白质序列的 FASTA 文件。-o qurongyuprotein.fasta
:输出文件,即聚类后的蛋白质序列的 FASTA 文件。-c 0.95
:设置聚类阈值,表示两个序列之间的相似度阈值,超过该阈值的序列将被聚类在一起。在这里是 0.95,即相似度达到 95% 以上的序列会被聚类在一起。-M 640000
:设置内存限制,单位为 MB。这里设置为 640000 MB,即 640 GB。-T 24
:设置线程数,即并行处理的线程数量,这里设置为 24。-n 5
:设置聚类后的序列名字的长度,即在输出文件中每个序列的名字长度限制为 5 个字符。-d 0
:用于设置聚类时是否考虑序列的长度差异,0 表示不考虑。-aS 0.9
:设置序列比对的阈值,即序列相似性的阈值,超过该阈值的序列将被认为是相似的。这里设置为 0.9,表示相似性达到 90% 以上的序列会被认为是相似的。-g 1
:设置是否使用贪婪法进行聚类,默认为 0,表示不使用。-sc 1
:设置是否使用长度和序列相似性的联合评分来确定序列相似性,默认为 0,表示不使用。-sf 1
:设置是否按比对的序列的第一个序列的序列名字进行输出,默认为 0,表示不按序列名字进行输出。
然后得到两个输出文件qurongyuprotein.fasta qurongyuprotein.fasta.clstr
其中qurongyuprotein.fasta.clstr的内容为
>Cluster 0
0 958aa, >NODE_416_length_5114_cov_66.925948_1... *
1 55aa, >NODE_33207_length_299_cov_3.144144_2... at 96.36%
2 90aa, >NODE_47254_length_274_cov_1.076142_1... at 95.56%
3 32aa, >NODE_67526_length_259_cov_0.961538_1... at 96.88%
4 81aa, >NODE_88620_length_245_cov_1.952381_1... at 95.06%
5 79aa, >NODE_101987_length_238_cov_1.279503_1... at 96.20%
6 47aa, >NODE_112950_length_232_cov_1.483871_1... at 95.74%
其中标*的为代表序列,可以将代表序列挑取出来进行后续分析,构建非冗余基因集。