一、cd-hit介绍
cd-hit 是一款用于将蛋白、核酸序列快速聚类的工具。由于宏基因组样品中可能包含相似物种,拼接结果中可能会存在一部分冗余序列,导致预测出来的基因包含冗余部分,可以通过聚类进行去冗余。
通常去冗余采用的聚类算法根据序列相似度对序列进行聚类,需要进行 all by all 的比较,例如 orthoMCL,不过这种方法非常耗时。而 cd-hit 使用一种贪婪的增量聚类方法,首先对输入的序列根据序列的长短进行排序,并从最长到最短的顺序处理它们。将最长的序列自动的分为第一类并作为第一类的代表序列,然后将剩下的序列与在其之前发现的代表性序列进行比较,根据序列相似性将其归为其中的一类或成为新的一个聚类的代表序列,如此遍历所有序列完成聚类过程。在默认方式中,序列仅和每一个聚类中的代表性序列(为这类中的最长序列)进行比较而不和这个类中的其他序列进行比对。在准确模式下,序列会和每个聚类中的所有序列进行比较然后决定是成为新的一类还是归为其中的一类中。
网址:http://weizhongli-lab.org/cd-hit/
二、软件安装
代码语言:javascript复制conda install -y cd-hit
三、使用案例
代码语言:javascript复制#对核酸聚类
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
bsub -q fat -n 12 -o %J.log -e %J.err sh cdhit.sh
seqkit stat mg.ffn mg.filter.ffn
file format type num_seqs sum_len min_len avg_len max_len
mg.ffn FASTA DNA 266,921 151,713,911 68 568.4 14,211
mg.filter.ffn FASTA DNA 255,020 145,715,723 68 571.4 14,211
#改变不同的阈值后,条数减少了
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
bsub -q fat -n 12 -o %J.log -e %J.err sh cdhit1.sh
seqkit stat mg.ffn mg.filter.ffn
file format type num_seqs sum_len min_len avg_len max_len
mg.ffn FASTA DNA 266,921 151,713,911 68 568.4 14,211
mg.filter.ffn FASTA DNA 242,836 139,449,227 68 574.3 14,211
grep '>' mg.filter.ffn >id.list #提取筛选后的id
#对氨基酸聚类
# echo "time cd-hit -i mg.faa -o mg.filter.faa -aS 0.9 -T 12 1>cdhit1.log 2>cdhit1.err" >cdhit1.sh
# bsub -q fat -n 12 -o %J.log -e %J.err sh cdhit1.sh
# 使用transeq翻译,transeq来自emboss软件包
transeq -sequence mg.filter.ffn -outseq mg.filter.faa
#用筛选后的id 提取氨基酸序列
awk '{print $1 }' id.list >id.txt
sed -i 's/>//' id.txt
seqkit grep -r -f id.txt mg.faa >mg.cdhit.faa
#seqkit也可以提取成氨基酸序列
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
写在最后:有时间我们会努力更新的。大家互动交流可以前去论坛,地址在下面,复制去浏览器即可访问,弥补下公众号没有留言功能的缺憾。
代码语言:javascript复制bioinfoer.com
有些板块也可以预设为大家日常趣事的分享等,欢迎大家来提建议。