教程 | 如何用cd-hit去除冗余序列？_cdhit-程序员宅基地

0.简介

生信分析中经常要根据指定条件查找相似序列，比如构建多个样品间的非冗余基因集、分析样品间的相似程度等等，cd-hit这款软件就可以用较短的时间解决此类问题，可以对单个数据集进行去冗余，包括DNA/RNA序列和蛋白序列，也可以对两个数据集进行比较。其工作原理可概述为：将所有序列按照参数设定进行聚类，并将每一组聚类中的最长序列作为代表序列进行输出，同时给出每组聚类下的每个序列名可供相似度分析使用。下面我们来简单介绍一下它的使用方法。

1.下载与安装

网址：http://cd-hit.org ；http://www.bioinformatics.org/cd-hit/ ；https://github.com/weizhongli/cdhit/archive/V4.6.2.tar.gz；

这是一个在linux系统下使用的工作，我们可以给自己的电脑装一个双系统或者在windows下使用linux的虚拟机。然后我们可以执行下面的命令进行解压（注意我们要将路径先切换到安装包所在的文件目录下，或者在执行命令时使用完整路径）。

gzip -d cdhit-4.2.tar.gz

然后进入到解压后的文件夹(我解压后的文件夹为cdhit-4.2，同样要注意我们的文件路径问题，如果上面使用的是完整路径，最好这里也使用完整路径，比如我使用完整路径是‘cd /home/zpf/cdhit-4.2’)

cd cdhit-4.2

最后编译一下就可以了，执行make

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文件中以“>”开头的行空格之前的内容，如下图中蓝色线圈出部分。

3.输出文件

Cd-hit有两个输出文件：一个是只含有所有代表序列（即去冗余后的序列）的fasta文件，其格式参看图1；另一个是以.clstr结尾的聚类信息文件，其格式如图2。

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

4.去除冗余的基本思路

首先对所有序列按照其长度进行排序，然后从最长的序列开始，形成第一个序列类，然后依次对序列进行处理，如果新的序列与已有的序列类的代表序列的相似性在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个不一样。

5.使用方法和参数介绍

Cd-hit运行时用很多参数可以进行调整设置，其运行命令为（参数仅为示例）在刚才编译的文件路径下执行：

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。

下图详解了-aL，-AL，-aS，-AS四个参数。

aL = Ra / R

AL = R - Ra

aS = Sa / S

AS = S - Sa

6.cdhit的缺点

1 它不能保证同一个序列类中的序列的相似性都在threshold之上，因为每次比对都是用新序列与序列类的代表序列进行，这就有可能使得序列类中除了代表序列外其他序列之间的相似性在threshold之下。比如A是代表序列，B与A的相似性大于0.95，C与A的相似性也大于0.95，但是这并不能保证B与C的相似性也大于0.95.

2 它不能保证一个序列类的病毒与另外一个序列类中的病毒的相似性也在threshold之上，原因还是在于用代表序列代表了整个序列类。

3 基于word filter的方法使得使用每个长度的word能够处理的冗余性水平有限，如使用长度为2的word只能够得到相似性在50%以上的序列，长度为3的word只能够得到相似性在66.7%以上的序列类，类似的，长度为5的word只能够得到相似性在80%以上的序列。在实际应用的时候需要注意选择的word长度与threshold的匹配。

本文链接：https://blog.csdn.net/weixin_41481113/article/details/83824743

原作者删帖不实内容删帖广告或垃圾文章投诉

智能推荐

docker删除重装，以及极其重要的/etc/docker/key.json文件-程序员宅基地

文章浏览阅读8.9k次，点赞2次，收藏2次。先说以下我为何要删除docker的原因吧：因为我感觉Docker Hub有点慢，就配置了阿里云的镜像加速器，可是按阿里云的官方配置完后我怎么感觉它更慢了，对比昨天配置阿里的maven镜像仓库后快到起飞的速度，我认为是此次配置没有生效。多次确认新加入的/etc/docker/demon.json文件无误后又多次systemctl了未果，随即我怀疑阿里给出的以下方案中的“修改”的/etc/dock...

空间物理——概述_空间物理概论-程序员宅基地

文章浏览阅读1.9k次，点赞3次，收藏4次。文章目录空间物理的研究对象太阳风能量向地球传输的三种方式和所需要的时间太阳内部结构、太阳活动太阳内部结构太阳活动太阳风速度从太阳表面到地球轨道附近变化参考空间物理的研究对象大气层：10KM以上，分成平流层、中层、低热层、热层、逃逸层电离层：60-90KM以上，一直到1000KM左右，部分电离气体，中性成风碰撞的影响不可忽略地球磁层：完全电离的气体，1000KM以上，可忽略碰撞，有太阳风和..._空间物理概论

BQ4050学习笔记（二）-程序员宅基地

文章浏览阅读2.9k次，点赞5次，收藏25次。BQ4050学习笔记（二）永久失效：如果发⽣严重故障，该设备可以永久禁⽤电池组。永久故障检查（IFC 和 DFW 除外）可以通过设置Settings:Enabled PFA、 Settings:Enabled PF B、 Settings:Enabled PF C 和Settings:Enabled PF D 中的相应位单独启⽤或禁⽤。所有永久在设置ManufacturingStatus()[PF]之前，故障检查（IFC 和 DFW 除外）被禁⽤。当任何PFStatus()位置位时，器件进⼊ PER_bq4050

SpringCloudAlibaba-Nacos注册中心的使用_spring-cloud-alibaba 使用nacos 2.1版本-程序员宅基地

文章浏览阅读152次。第二步：填写配置文件参数，这里定义了一个名字为application-user-dev.yaml的配置，使用的是YAML格式。DataID : 非常重要，可以看做是配置的文件的名字，在程序中拉取配置文件的时候需要指定Data ID。如果不使用默认的public命名空间，那么需要指定namespace配置为要使用的命名空间的Id值。第一步：打开Nacos监控面板，进入配置列表，新增一个user服务的配置文件。进入配置列表，切换到新建立的命名空间，创建配置文件。修改Nacos，添加命名空间。_spring-cloud-alibaba 使用nacos 2.1版本

CVE-2023-21716 Microsoft Word远程代码执行漏洞Poc_cve-2023-21716复现-程序员宅基地

文章浏览阅读293次。受害者打开python代码生成的RTF文件，RTF解析器解析恶意代码，触发堆溢出，Microsoft Word会闪退，用户其它Word中未保存的内容会丢失。_cve-2023-21716复现

c语言程序设计让a变成A,c语言程序设计a期末模拟试题.docx-程序员宅基地

文章浏览阅读451次。文件排版存档编号：[UYTR-OUPT28-KBNTL98-UYNN208]文件排版存档编号：[UYTR-OUPT28-KBNTL98-UYNN208]C语言程序设计A期末模拟试题C语言程序设计A期末模拟试题一一、单项选择题(每小题2分，共20分)由C++目标文件连接而成的可执行文件的缺省扩展名为( )。A. cpp B. exe C. obj D. li..._c语言如何将a转换成a

随便推点

利用beef和msf实现session远程命令_msf如何切换一个session-程序员宅基地

文章浏览阅读534次。笔记beef启动 beef 的方法msfbeef工具目录 /usr/share/beef-xss配置文件 config.yaml启动 beef 的方法1.beef-xss2./usr/share/beef-xss/beef(使用前需要修改默认的用户名称和密码)Web 管理界面 http://127.0.0.1:3000/ui/panelShellcode http://127.0.0.1:3000/hook.js测试页面 http://127.0.0.1:3000/demos/butch_msf如何切换一个session

python判断丑数_丑数问题及变种小结-程序员宅基地

文章浏览阅读503次。丑数问题及变种小结声明1 判断丑数因子只包含2,3,5的数称为丑数(Ugly Number)，习惯上把1当作第一个丑数面试lintcode 517 ugly numbersegmentfault剑指offer 面试题34 丑数数组解法：参考剑指offer，将待判断目标依次连续整除2，3，5，若是最后获得1，证实该数为丑数；优化/*** 依次整除2，3，5判断(2，3，5顺序判断时间最优)* htt..._编写python来证明一个数是丑数

python自动化测试实战 —— WebDriver API的使用_python webdriver api-程序员宅基地

文章浏览阅读1.9k次，点赞30次，收藏11次。Selenium 简介： WebDriver是Selenium Tool套件中最重要的组件。Selenium 2.0之后已经将Selenium和WebDriver进行合并，作为一个更简单、简洁、有利于维护的API提供给测试人员使用。它提供了一套标准的接口，可以用多种编程语言调用，并且和浏览器进行交互。 WebDriver可以对浏览器进行控制，包括输入URL，点击按钮，填写表单，滚动页面，甚至是执行JavaScript代码。同时，它也能够获取网页中的信息，如文本，标签，属_python webdriver api