美文网首页RNA_SEQ
转录组重建系统发育(四)使用cdhit去冗余和使用transde

转录组重建系统发育(四)使用cdhit去冗余和使用transde

作者: 惊鸿影 | 来源:发表于2021-09-25 19:30 被阅读0次

    1.cd-hit的安装和使用

    cd-hit 是用于蛋白质序列或核酸序列聚类的工具,根据序列的相似度对序列进行聚类以去除冗余的序列,一般用于构建非冗余的数据集用于后续的实验分析。

    通常来说,根据序列相似度对序列进行聚类,首先想到的可能是通过计算两两序列之间的相似度对序列进行聚类,这样需要进行all by all的比较,相对来说比较费时,而 cd-hit 软件可以避开all by all比较的问题,大大缩短了运行时间。cd-hit 是一种贪婪的增量聚类方法,首先对输入的序列根据序列的长短进行排序,并从最长到最短的顺序处理它们。将最长的序列自动的分为第一类并作为第一类的代表序列,然后将剩下的序列与在其之前发现的代表性序列进行比较,根据序列相似性将其归为其中的一类或成为新的一个聚类的代表序列,如此遍历所有序列完成聚类过程。 在默认方式中,序列仅和每一个聚类中的代表性序列(为这类中的最长序列)进行比较而不和这个类中的其他序列进行比对。 在准确模式下,序列会和每个聚类中的所有序列进行比较然后决定是成为新的一类还是归为其中的一类中。
    cd-hit下载和安装
    也可以使用conda进行安装

    conda install -c  bioconda cdhit
    

    使用很简单,输入目标文件名和输出文件名以及一些简单的参数就像能运行。

    cd-hit-est -i ./SRS7102777.Trinity -o ./SRS7102777_cd-hit-est.fasta -c 0.98a  -d 0 -T 8 -M 64000
    

    核苷酸序列使用cd-hit-est命令,氨基酸序列使用cd-hit命令。
    -i 输入文件,fasta格式的序列
    -o 输出文件路径和名字
    -c 相似性(clustering threshold),0.9表示相似性大于等于90%的为一类
    -d 0表示使用 fasta 标题中第一个空格前的字段作为序列名字
    -M 64000,64GB RAM
    -T 使用的线程数
    之后就能获得去冗余后的文件。
    去冗余后就能使用transdecoder进行开放阅读框预测了。

    2.transdecoder的安装和使用

    软件说明书:https://github.com/TransDecoder/TransDecoder/wiki

    TransDecoder 基于以下标准识别可能的编码序列:

    1.在转录本序列中需要能够找到一个(满足)最小(限定)长度的ORF;
    2.对数似然数得分大于0。(与GeneID软件计算得到的得分相类似);
    3.第一阅读框的对数似然数打分同其它5个阅读框比较为最大值时;
    4.如果候选的ORF完全被包含在其它候选ORF的框架内,那么报告最长的ORF。否则,一个单独的转录本会得到多个ORF的报告。(考虑到有操纵子、嵌合体等情况);
    5.作为可选项,预测出的多肽在Pfam domain库中存在比对分值高于得分阈值之上的。

    TransDecoder通过运行一个包含目的转录本序列的fasta文件来实现功能。简单的用法如下:

    # Step 1: 提取最长的开放阅读框
    TransDecoder.LongOrfs -t target_transcripts.fasta  -O .
    默认情况下,TransDecoder.LongOrfs将识别长度至少为100个氨基酸的开放阅读框。你可以通过-m参数来降低这个值,但是要知道随着最小长度的变短,ORF预测的假阳性率迅速增长。
    # Step 2: (可选)
    可选地,可以通过blast或者pfam搜索已知蛋白的同源序列来识别ORF。见下面的Including homology searches as ORF retention criteria章节。
    # Step 3: 预测可能的编码区
    TransDecoder.Predict -t target_transcripts.fasta -O .
    如果转录本按照有义链进行了调整,那么使用-S标记来只检查top链。完整的使用信息见后面。
    候选编码区的最终集合可以在文件.transdecoder中找到。扩展名包括.pep,.cds,.gff3和.bed。
    
    输出文件如下:
    longest_orfs.pep : 所有达到最小长度标准的ORF, 不管是否编码
    longest_orfs.gff3 : 在目的转录本中发现的所有ORF的位置
    longest_orfs.cds : 所有检测到的ORF的核酸编码序列
    longest_orfs.cds.top_500_longest : 前500个最长的ORF,用于训练一个编码序列的马尔科夫模型
    hexamer.scores : 每个k-mer的对数似然得分 (coding/random)
    longest_orfs.cds.scores : 每个ORF同6个阅读框间对数似然得分的总和
    longest_orfs.cds.scores.selected : 根据得分标准所选出的ORF
    longest_orfs.cds.best_candidates.gff3 : 转录本中选出的ORF的位置
    然后,最后的输出文件在你当前的工作目录中。
    transcripts.fasta.transdecoder.pep : 最终候选ORF的蛋白质序列;所有较长ORF中的较短的候选序列已被移除。
    transcripts.fasta.transdecoder.cds : 最终候选ORF的编码区的核酸序列。
    transcripts.fasta.transdecoder.gff3 : 最终被选中的ORF在目的转录本中的位置
    transcripts.fasta.transdecoder.bed : 用来描述ORF位置的bed格式文件,最好用GenomeView或IGV来查看。
    

    其中pep文件就是我们需要的预测的氨基酸序列信息,pep文件内容如下:

    image.png

    文件中:header行包含了蛋白的ID信息、原始转录本ID信息、type信息、长度、正负链、打分信息、ORF坐标信息

    相关文章

      网友评论

        本文标题:转录组重建系统发育(四)使用cdhit去冗余和使用transde

        本文链接:https://www.haomeiwen.com/subject/bfhggltx.html