教程 | 如何用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 java_apple登录 后端java实现最终版-程序员宅基地
文章浏览阅读298次。import com.alibaba.fastjson.JSONArray;import com.alibaba.fastjson.JSONObject;import com.auth0.jwk.Jwk;import com.helijia.appuser.modules.user.vo.AppleCredential;import com.helijia.common.api.model.Api..._com.auth0.jwk.jwk
NLP学习记录（六）最大熵模型MaxEnt_顺序潜在最大熵强化学习(maxent rl)-程序员宅基地
文章浏览阅读4.7k次。原理在叧掌握关于未知分布的部分信息的情况下，符合已知知识的概率分布可能有夗个，但使熵值最大的概率分布最真实地反映了事件的的分布情况，因为熵定义了随机变量的不确定性，弼熵值最大时，随机变量最不确定，最难预测其行为。最大熵模型介绍我们通过一个简单的例子来介绍最大熵概念。假设我们模拟一个翻译专家的决策过程，关于英文单词in到法语单词的翻译。我们的翻译决策模型p给每一个单词或短语分配一..._顺序潜在最大熵强化学习(maxent rl)
计算机毕业设计ssm科研成果管理系统p57gs系统+程序+源码+lw+远程部署-程序员宅基地
文章浏览阅读107次。计算机毕业设计ssm科研成果管理系统p57gs系统+程序+源码+lw+远程部署。springboot基于springboot的影视资讯管理系统。ssm基于SSM高校教师个人主页网站的设计与实现。ssm基于JAVA的求职招聘网站的设计与实现。springboot校园头条新闻管理系统。ssm基于SSM框架的毕业生离校管理系统。ssm预装箱式净水站可视化信息管理系统。ssm基于SSM的网络饮品销售管理系统。
Caused by: org.xml.sax.SAXParseException; lineNumber: 38; columnNumber: 9; cvc-complex-type.2.3: 元素_saxparseexception; linenumber: 35; columnnumber: 9-程序员宅基地
文章浏览阅读1.6w次。不知道大家有没有遇到过与我类似的报错情况，今天发生了此错误后就黏贴复制了报错信息“Caused by: org.xml.sax.SAXParseException; lineNumber: 38; columnNumber: 9; cvc-complex-type.2.3: 元素 'beans' 必须不含字符 [子级], 因为该类型的内容类型为“仅元素”。”然后就是一顿的百度啊， 可一直都没有找到..._saxparseexception; linenumber: 35; columnnumber: 9; cvc-complex-type.2.3:
计算机科学与技术创新创业意见,计算机科学与技术学院大学生创新创业工作会议成功举行...-程序员宅基地
文章浏览阅读156次。(通讯员 粟坤萍 2018-04-19)4月19日，湖北师范大学计算机科学与技术学院于教育大楼学院会议室1110成功召开大学生创新创业工作会议。参与本次会议的人员有党总支副书记黄海军老师，创新创业学院吴杉老师，计算机科学与技术学院创新创业活动指导老师，15、16、17级各班班主任及学生代表。首先吴杉老师介绍了“互联网+”全国大学生创新创业大赛的相关工作进度，动员各级班主任充分做好“大学生创新创业大..._湖北师范 吴杉
【Android逆向】爬虫进阶实战应用必知必会-程序员宅基地
文章浏览阅读1.1w次，点赞69次，收藏76次。安卓逆向技术是一门深奥且充满挑战的领域。通过本文的介绍，我们了解了安卓逆向的基本概念、常用工具、进阶技术以及实战案例分析。然而，逆向工程的世界仍然在不断发展和变化，新的技术和方法不断涌现。展望未来，随着安卓系统的不断更新和加固，逆向工程将面临更大的挑战。同时，随着人工智能和机器学习技术的发展，我们也许能够看到更智能、更高效的逆向工具和方法的出现。由于篇幅限制，本文仅对安卓逆向技术进行了介绍和案例分析。
随便推点
Python数据可视化之环形饼图_数据可视化绘制饼图或圆环图-程序员宅基地
文章浏览阅读1.1k次。制作饼图还需要下载pyecharts库，Echarts 是一个由百度开源的数据可视化，凭借着良好的交互性，精巧的图表设计，得到了众多开发者的认可。随着学习python的热潮不断增加，Python数据可视化也不停的被使用，那我今天就介绍一下Python数据可视化中的饼图。在我们的生活和学习中，编程是一项非常有用的技能，能够丰富我们的视野，为各行各业的领域提供了新的角度。环形饼图的制作并不难，主要是在于数据的打包和分组这里会有点问题，属性的标签可以去 这个网站进行修改。图中的zip压缩函数，并分组打包。_数据可视化绘制饼图或圆环图
SpringMVC开发技术~5~基于注解的控制器_jsp/servlet到controller到基于注解的控制器-程序员宅基地
文章浏览阅读325次。1 Spring MVC注解类型Controller和RequestMapping注释类型是SpringMVC API最重要的两个注释类型。基于注解的控制器的几个优点：一个控制器类可以控制几个动作，而一个实现了Controller接口的控制器只能处理一个动作。这就允许将相关操作写在一个控制器类内，从而减少应用类的数量基于注解的控制器的请求映射不需要存储在配置文件中，而是使用RequestM..._jsp/servlet到controller到基于注解的控制器
利用波特图来满足动态控制行为的要求-程序员宅基地
文章浏览阅读260次，点赞3次，收藏4次。相位裕量可以从增益图中的交越频率处读取（参见图2）。使用的开关频率、选择的外部元件（例如电感和输出电容），以及各自的工作条件（例如输入电压、输出电压和负载电流）都会产生巨大影响。图2所示为波特图中控制环路的增益曲线，其中提供了两条重要信息。对于图2所示的控制环路，这个所谓的交越频率出现在约80 kHz处。通过使用波特图，您可以查看控制环路的速度，特别是其调节稳定性。图2. 显示控制环路增益的波特图（约80 kHz时，达到0 dB交越点）。图3. 控制环路的相位曲线，相位裕量为60°。
Glibc Error: `_obstack@GLIBC_2.2.5‘ can‘t be versioned to common symbol ‘_obstack_compat‘_`_obstack@glibc_2.2.5' can't be versioned to commo-程序员宅基地
文章浏览阅读1.8k次。Error: `_obstack@GLIBC_2.2.5’ can’t be versioned to common symbol '_obstack_compat’原因：https://www.lordaro.co.uk/posts/2018-08-26-compiling-glibc.htmlThis was another issue relating to the newer binutils install. Turns out that all was needed was to initi_`_obstack@glibc_2.2.5' can't be versioned to common symbol '_obstack_compat
基于javaweb+mysql的电影院售票购票电影票管理系统(前台、后台)_电影售票系统javaweb-程序员宅基地
文章浏览阅读3k次。基于javaweb+mysql的电影院售票购票电影票管理系统(前台、后台)运行环境Java≥8、MySQL≥5.7开发工具eclipse/idea/myeclipse/sts等均可配置运行适用课程设计，大作业，毕业设计，项目练习，学习演示等功能说明前台用户：查看电影列表、查看排版、选座购票、查看个人信息后台管理员：管理电影排版，活动，会员，退票，影院，统计等前台：后台：技术框架_电影售票系统javaweb
分分钟拯救监控知识体系-程序员宅基地
文章浏览阅读95次。分分钟拯救监控知识体系本文出自：http://liangweilinux.blog.51cto.com0 监控目标我们先来了解什么是监控，监控的重要性以及监控的目标，当然每个人所在的行业不同、公司不同、业务不同、岗位不同、对监控的理解也不同，但是我们需要注意，监控是需要站在公司的业务角度去考虑，而不是针对某个监控技术的使用。监控目标1.对系统不间断实时监控:实际上是对系统不间..._不属于监控目标范畴的是 实时反馈系统当前状态