使用TransDecoder预测CDS

TransDecoder按照其官网的说明，主要用于识别转录本序列中的潜在的编码区域，也就是预测CDS。转录本可以由RNA-Seq数据通过Trinity组装来的，也可以由RNA-Seq比对到参考基因组上构建的转录本。
最新版本的TransDecoder可在此处找到。

TransDecoder识别可能的编码区域是基于以下几个标准：

  1. a minimum length open reading frame (ORF) is found in a transcript sequence
  2. a log-likelihood score similar to what is computed by the GeneID software is > 0
  3. the above coding score is greatest when the ORF is scored in the 1st reading frame as compared to scores in the other 5 reading frames
  4. if a candidate ORF is found fully encapsulated by the coordinates of another candidate ORF, the longer one is reported. However, a single transcript can report multiple ORFs (allowing for operons, chimeras, etc)
  5. a PSSM is built/trained/used to refine the start codon prediction.
  6. optional the putative peptide has a match to a Pfam domain above the noise cutoff score

Step 1: extract the long open reading frames

TransDecoder.LongOrfs -t target_transcripts.fasta

默认情况下，TransDecoder.LongOrfs将识别至少100个氨基酸长的ORF。您可以通过’-m’参数降低此值，但可以知道，使用更短的最小长度标准，误报ORF预测的比率会急剧增加。

Step 2: (optional)

通过blast或pfam搜索鉴定与已知蛋白具有同源性的ORF。
参考下面的将同源性搜索包括为ORF保留标准部分。

Step 3: predict the likely coding regions

TransDecoder.Predict -t target_transcripts.fasta [ homology options ]
最终的文件可以在当前目录找到，也就是后缀为.pep, .cds, .gff3和.bed的文件

一般来说，我们会使用TransDecoder对无参转录组的拼接结果序列预测其CDS，所以我们可以先将拼接序列用BLAST比对nr以及swissprot蛋白数据库，然后提取其比对上的同源序列的位置来识别CDS，最后再通过TransDecoder的第一步和第三步来预测那些未比对上的序列的CDS。

以有参考基因组的转录结果GTF文件预测编码区域

此处的过程与上述过程相同，不同之处在于，我们必须首先生成一个与转录本序列相对应的fasta文件，最后，我们重新计算GFF3格式的基因组注释文件

  1. 使用基因组和transcripts.gtf文件构建transcript fasta文件，如下所示：
    util/gtf_genome_to_cdna_fasta.pl transcripts.gtf test.genome.fasta > transcripts.fasta

  2. 接下来，将成绩单结构GTF文件转换为alignment-GFF3格式的文件（之所以这样做，是因为我们的进程对gff3而不是对gtf文件进行操作-没什么大的影响）。像这样将gtf转换为alignment-gff3：
    util/gtf_to_alignment_gff3.pl transcripts.gtf > transcripts.gff3

  3. 运行上述过程以生成最佳候选ORF预测：
    TransDecoder.LongOrfs -t transcripts.fasta (optionally, identify peptides with homology to known proteins) TransDecoder.Predict -t transcripts.fasta [ homology options ]

  4. 最后，生成基于基因组的编码区域注释文件：
    util/cdna_alignment_orf_to_genome_orf.pl \ transcripts.fasta.transdecoder.gff3 \ transcripts.gff3 \ transcripts.fasta > transcripts.fasta.transdecoder.genome.gff3

输出文件说明

创建了一个工作目录（例如transcripts.transdecoder_dir /）来运行和存储管道的中间部分，其中包含：

1
2
3
4
5
6
7
8
9
10
longest_orfs.pep   : all ORFs meeting the minimum length criteria, regardless of coding potential.
longest_orfs.gff3  : positions of all ORFs as found in the target transcripts
longest_orfs.cds   : the nucleotide coding sequence for all detected ORFs

longest_orfs.cds.top_500_longest   : the top 500 longest ORFs, used for training a Markov model for coding sequences.

hexamer.scores                     : log likelihood score for each k-mer  (coding/random)

longest_orfs.cds.scores            : the log likelihood sum scores for each ORF across each of the 6 reading frames
longest_orfs.cds.scores.selected   : the accessions of the ORFs that were selected based on the scoring criteria (described at top)
1
2
3
4
transcripts.fasta.transdecoder.pep : peptide sequences for the final candidate ORFs; all shorter candidates within longer ORFs were removed.
transcripts.fasta.transdecoder.cds  : nucleotide sequences for coding regions of the final candidate ORFs
transcripts.fasta.transdecoder.gff3 : positions within the target transcripts of the final selected ORFs
transcripts.fasta.transdecoder.bed  : bed-formatted file describing ORF positions, best for viewing using GenomeView or IGV.

最终输出文件如下:

  • *.pep （是最终的候选ORF编码的蛋白序列）
  • *.cds （最终预测的CDS序列）
  • *.gff3 （是表示ORF的目标转录本的位置）
  • *.bed （用于后期的IGV可视化，以BED格式存放ORF位置信息）

一键化脚本

也有作者发布了脚本，只需要准备参考基因组的数据库位置信息以及相关软件即可自动进行分析，得出结果，非常方便快捷。
具体有需要脚本的后续可以留言～