cdhit去冗余

生信喵

一、cd-hit介绍
       cd-hit 是一款用于将蛋白、核酸序列快速聚类的工具。由于宏基因组样品中可能包含相似物种，拼接结果中可能会存在一部分冗余序列，导致预测出来的基因包含冗余部分，可以通过聚类进行去冗余。
       通常去冗余采用的聚类算法根据序列相似度对序列进行聚类，需要进行 all by all 的比较，例如 orthoMCL，不过这种方法非常耗时。而 cd-hit 使用一种贪婪的增量聚类方法，首先对输入的序列根据序列的长短进行排序，并从最长到最短的顺序处理它们。将最长的序列自动的分为第一类并作为第一类的代表序列，然后将剩下的序列与在其之前发现的代表性序列进行比较，根据序列相似性将其归为其中的一类或成为新的一个聚类的代表序列，如此遍历所有序列完成聚类过程。在默认方式中，序列仅和每一个聚类中的代表性序列（为这类中的最长序列）进行比较而不和这个类中的其他序列进行比对。在准确模式下，序列会和每个聚类中的所有序列进行比较然后决定是成为新的一类还是归为其中的一类中。
       网址：http://weizhongli-lab.org/cd-hit/
二、软件安装

1. conda install -y cd-hit

复制代码

三、使用案例

1. #对核酸聚类
   
2. echo "time cd-hit-est -i mg.ffn -o mg.filter.ffn -aS 0.9 -T 12 1>cdhit.log 2>cdhit.err" >cdhit.sh
   
3. bsub -q fat -n 12 -o %J.log -e %J.err sh cdhit.sh
   
4. seqkit stat mg.ffn mg.filter.ffn
   
5.         file           format  type  num_seqs      sum_len  min_len  avg_len  max_len
   
6.         mg.ffn         FASTA   DNA    266,921  151,713,911       68    568.4   14,211
   
7.         mg.filter.ffn  FASTA   DNA    255,020  145,715,723       68    571.4   14,211
   
8. #改变不同的阈值后，条数减少了
   
9. echo "time cd-hit-est -i mg.ffn -o mg.filter.ffn -aS 0.85 -T 12 -c 0.85 1>cdhit1.log 2>cdhit1.err" >cdhit1.sh
   
10. bsub -q fat -n 12 -o %J.log -e %J.err sh cdhit1.sh
    
11. seqkit stat mg.ffn mg.filter.ffn
    
12.         file           format  type  num_seqs      sum_len  min_len  avg_len  max_len
    
13.         mg.ffn         FASTA   DNA    266,921  151,713,911       68    568.4   14,211
    
14.         mg.filter.ffn  FASTA   DNA    242,836  139,449,227       68    574.3   14,211
    
15. grep '>' mg.filter.ffn >id.list #提取筛选后的id
    
16. 
    
17. #对氨基酸聚类
    
18. # echo "time cd-hit -i mg.faa -o mg.filter.faa -aS 0.9 -T 12 1>cdhit1.log 2>cdhit1.err" >cdhit1.sh
    
19. # bsub -q fat -n 12 -o %J.log -e %J.err sh cdhit1.sh
    
20. # 使用transeq翻译，transeq来自emboss软件包
    
21. transeq -sequence mg.filter.ffn -outseq mg.filter.faa
    
22. #用筛选后的id 提取氨基酸序列
    
23. awk '{print $1 }' id.list >id.txt
    
24. sed -i 's/>//' id.txt
    
25. seqkit grep -r -f id.txt mg.faa >mg.cdhit.faa
    
26. #seqkit也可以提取成氨基酸序列
    
27. time seqkit translate mg.filter.ffn -T 11 >test.faa

复制代码

选项参数：
       -i 输入文件，fasta 格式的序列
       -o 输出文件路径和名字
       -c 相似性（clustering threshold），0.9 表示相似性大于等于 90%的为一类
       -n 两两序列进行序列比对时选择的 word size
       -d 0 表示使用 fasta 标题中第一个空格前的字段作为序列名字
       -M 16000，16GB RAM
       -T 使用的线程数
Choose of word size:
       -n 5 for thresholds 0.7 ~ 1.0
       -n 4 for thresholds 0.6 ~ 0.7
       -n 3 for thresholds 0.5 ~ 0.6
       -n 2 for thresholds 0.4 ~ 0.5