上海锐翌生物科技有限公司

服务电话:021-51001612

邮箱:support@realbio.cn

技术课堂

教程 | 如何用cd-hit去除冗余序列?
发布时间:2017-09-07 13:34   点击率:
生信分析中经常要根据指定条件查找相似序列,比如构建多个样品间的非冗余基因集、分析样品间的相似程度等等,cd-hit这款软件就可以用较短的时间解决此类问题。其工作原理可概述为:将所有序列按照参数设定进行聚类,并将每一组聚类中的最长序列作为代表序列进行输出,同时给出每组聚类下的每个序列名可供相似度分析使用。下面我们来简单介绍一下它的使用方法。

1. 下载与安装
Cd-hit下载网址为https://github.com/weizhongli/cdhit/archive/V4.6.2.tar.gz,需要在linux系统下操作。解压压缩包后进入软件本体路径,直接输入命令:make,进行编译即可。

2. 输入文件
Cd-hit的输入文件仅有一个fasta格式文件 ,一般来说cd-hit是将几个样品的基因或蛋白序列进行聚类,所以需要将这些样品的序列汇总到一起作为输入文件,可在linux系统下通过cat命令实现:

cat a.fasta b.fasta c.fasta > all.fasta

其中a.fasta,b.fasta,c.fasta为fasta格式的三个样品基因或蛋白序列,all.fasta为汇总后的序列,在分析中作为cd-hit的输入序列。值得注意的是,在三个样品序列中不能有序列名相同的序列,否则会出现错误。因此,一般在分析时会在各样品序列名前添加样品名,这样即可避免重复。序列名是fasta文件中以“>”开头的行空格之前的内容,如图2-1中蓝色线圈出部分。

 

图2-1 fasta文件格式

3. 分析过程与重要参数介绍
Cd-hit运行时用很多参数可以进行调整设置,其运行命令为(参数仅为示例):

/path-to-cdhit-4.6.8/cd-hit -i all.fasta -o new.fa -c 0.8 -aS 0.8 -d 0

下面简单介绍一下重要的几个参数:
-i:输入文件,fasta格式。
-o:输出文件前缀,输出文件有两个,分别为fasta格式序列文件和以.clstr结尾的聚类信息文件。
-c:较短序列比对到长序列的bp与自身bp数的比值超过该数值则聚类为一组,默认为0.9。
-d:聚类信息文件中各个聚类组中序列名的长度,设为0则将取完整序列名。
-aL:控制代表序列比对严格程度的参数,默认为0,若设为0.8则表示比对区间要占到代表(长)序列的80%。
-AL:控制代表序列比对严格程度的参数,默认为99999999,若设为40则表示代表序列的非比对区间要短于40bp。
-aS:控制短序列比对严格程度的参数,默认为0,若设为0.8则表示比对区间要占到短序列的80%。
-AS:控制短序列比对严格程度的参数,默认为99999999,若设为40则表示短序列的非比对区间要短于40bp。
图3-1详解了-aL,-AL,-aS,-AS四个参数。

 
图3-1 参数图解
aL = Ra / R
AL = R - Ra
aS = Sa / S
AS = S - Sa

4. 输出文件介绍
Cd-hit有两个输出文件:一个是只含有所有代表序列(即去冗余后的序列)的fasta文件,其格式参看图2-1;另一个是以.clstr结尾的聚类信息文件,其格式如图4-1。

 
▲ 聚类信息文件

以“>”开头的是一个聚类组。每组下面按序号排列,如上图中Cluster 1组有5个聚类序列。每个聚类序列有一个百分比或“*”,百分比代表该序列与代表序列的相似度,“*”代表该序列即为代表序列。