3.2.2 cd-hit去除冗余序列
不建议用cd-hit来进行序列聚类。虽然速度很快，但不具备序列去重功能，我两个一样的序列居然被分到了不同的cluster。亲测。
一、CD-HIT简介
去冗余，也可以叫做相似序列的聚类。
CD-HIT stands for Cluster Database at High Identity with Tolerance. The program (cd-hit) takes a fasta format sequence database as input and produces a set of ‘non-redundant’ (nr) representative sequences as output. In addition cd-hit outputs a cluster file, documenting the sequence ‘groupies’ for each nr sequence representative.
cd-hit现在已经形成了一个软件包，可以对单个数据集进行去冗余，包括DNA/RNA序列和蛋白序列，也可以对两个数据集进行比较。
网址：
二、安装与使用
2.1 安装
- 下载：http://www.bioinformatics.org/cd-hit/
- 运行环境：Linux
它的使用也很简单，下载软件包后解压缩，进入该文件夹，“make"就可以了。
git clone https://github.com/weizhongli/cdhit.git
cd cdhit
make
#修改环境变量
vim /etc/profile
export PATH=/data/software/cdhit:$PATH
source /etc/profile
2.2 使用
输入的文件是fasta格式的序列文件，通过序列比对聚类（Cluster）的方法去除冗除、相似的序列，最后输出一个非冗除（non-redundant，nr）的序列文件。 另外，还有一个序列比对的结果。
其使用也简单，如下：
./cd-hit -i inputfile -o outputfile -c threshold -n wordLength
./cd-hit -i inputfile -o outputfile -c 0.97 -n 5
输入文件inputfile
>1
AGGGTGGGAACAGAGTGACCGAGGGGGCAGCCTTGGG
>2
CAGGGGGAAGACCGATGGGCCCTTGGTGGAAGCTGAA
>3
CCAGGGGGAAGACCGATGGGCCCTTGGTGGAAGCTGA
>4
CCAGGGGGAAGACCGATGGGCCCTTGGTGGAGGCTGA
>5
AGGGCGGGAACAGAGTGAACGATGGGGCATCCTTGGG
>5-2
AGGGCGGGAACAGAGTGAACGATGGGGCATCCTTGGG
>5-3
AGGGCGGGAACAGAGTGAACGATGGGGCATCCTTGGG
会输出两个文件：
一个是outputfile.clstr文件，会记录那些序列被聚类成一团
>Cluster 0
0 514aa, >2... *
>Cluster 1
0 513aa, >4... *
>Cluster 2
0 511aa, >3... *
>Cluster 3
0 502aa, >1... *
>Cluster 4
0 477aa, >5... *
1 477aa, >5-2... at 100.00%
2 477aa, >5-3... at 100.00%
一个是outputfile的fasta文件，每一个中选择一个序列
>1
AGGGTGGGAACAGAGTGACCGAGGGGGCAGCCTTGGG
>2
CAGGGGGAAGACCGATGGGCCCTTGGTGGAAGCTGAA
>3
CCAGGGGGAAGACCGATGGGCCCTTGGTGGAAGCTGA
>4
CCAGGGGGAAGACCGATGGGCCCTTGGTGGAGGCTGA
>5
AGGGCGGGAACAGAGTGAACGATGGGGCATCCTTGGG
CD-HIT在线服务 CD-HIT也有在线运行的服务。
网址：http://weizhong-lab.ucsd.edu/cdhit_suite/cgi-bin/
三、算法说明
基本思路：
- 首先对所有序列按照其长度进行排序，然后从最长的序列开始，形成第一个序列类，
- 然后依次对序列进行处理，如果新的序列与已有的序列类的代表序列的相似性在cutoff以上则把该序列加到该序列类中，否则形成新的序列类。
之所以快主要是两个方面的原因：
- 一个是使用了word过滤方法，即如果两条序列之间的相似性在80%（假设序列长度为100），那么它们至少有60个相同的长度为2的word，至少有40个相同的长度为3的word，至少有20个相同的长度为4的word。基于这个原则，在处理新的序列的时候，如果新的序列与已有序列的相同word的长度不能满足这些要求则不需要进行比对了，这极大的降低了时间消耗；
- 另外一个速度快的原因是使用了index table，可以很快的计算序列之间相同word的数目。
#当序列相似性在80%时，有20个位点是有差异的，极端的情况就是这20个位点对应的长度为2的字符串都不一样，因此是40个不一样，当有更多的不一样时，两条序列的相似性不可能在80%；同理，如果这20个位点对应的长度为4的字符串都不一样，则有80个不一样。
四、讨论
尽管很好用也很快，但是也需要注意其缺点：
- 它不能保证同一个序列类中的序列的相似性都在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的匹配。
4.1 CD-HIT sequence identity
CD-HIT identity = (number of alignment columns containing matching letters) / (length of shorter sequence)
该定义在实践中存在许多问题。 最大的问题是，gaps并不算作差异。 在极端情况下，100％ identit有可能比对结果中有gap。
CD-HIT定义alignment参数的选择比BLAST更敏感。 这是因为减少的空位罚分倾向于产生具有更多同一性的“gappier”比对，这通过添加更多的gap而成为可能，尤其是对于核苷酸序列。 这总是会增加CD-HIT定义的identity，因为gap不算作差异。
考虑这个简单的例子：
根据BLAST和CD-HIT，该比对具有6/10 = 60％的同一性。 现在假设间隙惩罚减少以允许以下对齐，其具有两个更多identity：
现在BLAST identity是8/12 = 67％，并且CD-HIT identity是8/10 = 80％。 请注意，CD-HIT identity 变化的更多。
CD-HIT定义倾向于在给定的identity阈值下产生更少的簇（更大的平均大小），但这些更多是低质量的。
ps
大致明白了这个聚类的思路，但具体是怎么一个设定阈值达到什么聚类效果，以及具体是怎么聚的还需要后续的关注和理解
五、我的案例
./cd-hit -i inputfile -o outputfile -c 0.97 -n 5
-d ： length of description in .clstr file, default 20。if set to 0, it takes the fasta defline and stops at first space
参考资料
