提取最长转录本

由于可变剪切一个基因得到好多个序列长度不同的转录本，应该挑选出序列最长的最长转录本作为数据的分析

看一下有规律的蛋白文件命名，一般是以下形式

前面是基因，P后面是转录本

写个脚本提取

import sys
import re


filename = sys.argv[1]

seqlist = {}

i = 0
with open(filename,'r') as r:
  lines=r.readlines()
  for l in lines:
    i += 1
    linecontent = l.strip()

    if linecontent.startswith(">"):
        if i > 1:
            name = linecontent[1:]
            id = name[0: name.rfind(".",1)]
            if id not in seqlist.keys():
                seqlist[id] = {}
            seqlist[id][name] = ""

        elif i == 1:
            name = linecontent[1:]
            id = name[0: name.rfind(".",1)]
            seqlist[id] = {}
            seqlist[id][name] = ""
    else:
        seqlist[id][name] = seqlist[id][name] + linecontent

maxseq = {}

for k,v in seqlist.items():
    d=0
    for k1,v1 in v.items():
        d +=1
        
        if d ==1:
            maxseq[k1] = v1
            maxlen = len(v1)
            k_last = k1
        else:
            if len(v1) > maxlen:
                maxseq[k1] = v1
                del maxseq[k_last]
                k_last = k1
 
        
#生成结果
result1 = filename + ".max"
result2 = filename + ".max.list"
result3 = filename + ".max.len"
with open(result1,'w') as w1:
    with open(result2,'w') as w2:
        with open(result3,'w') as w3:
            for key,value in maxseq.items():
                w1.write(">" + key + "\n" + value + "\n")
                w2.write(key + "\n")
                w3.write(key + "\t" + str(len(value))  +"\n")

这个是提取的结果每个蛋白的长度信息

使用TransDecoder软件

最长转录本（Longest Transcript）指的是一个基因所产生的RNA分子中最长的那个。这个RNA分子可能是mRNA（信使RNA），也可能是其他类型的RNA，比如pre-mRNA（前体mRNA）、snRNA（小核RNA）或miRNA（微小RNA）等。最长转录本的长度取决于基因的转录调控和剪接等机制。
最长ORF（Open Reading Frame）指的是一个RNA分子中具有最长的连续未被中断的编码区域，一般指的是可以被翻译成蛋白质的区域。ORF通常从起始密码子（start codon）开始，一直延伸到终止密码子（stop codon）结束。在一个RNA分子中可能存在多个ORF，而其中最长的ORF通常被认为是最可能被翻译成蛋白质的部分。
因此，最长转录本和最长ORF的区别在于，前者指的是RNA分子中整体长度最长的那个，而后者指的是RNA分子中具有最长编码潜力的连续区域。最长转录本可以包含一个或多个ORF，但最长ORF不一定完全重叠于最长转录本中。研究这两个概念可以帮助科学家更好地理解基因的转录和翻译过程，以及基因产物（蛋白质）的功能。

需要准备文件·pep.fasta

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

提取最长的开放阅读框

TransDecoder.LongOrfs -t trans_wild12.fa 

运行界面

预测可能的编码区

TransDecoder.Predict -t trans_wild12.fa

候选编码区的最终集合可以在文件.transdecoder中找到。扩展名包括.pep，.cds，.gff3和.bed。

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

1、需要有参转录组比对后拼接的转录本的GTF文件以及参考基因组序列：
2、将GTF文件转化为GFF3文件
3、提取、预测
4、最后生成一个基于有参基因组的编码区域注释文件:

/public/home/fengting/tools/TransDecoder-v5.5.0/util/gtf_genome_to_cdna_fasta.pl sbam/$id.gtf /public/home/fengting/task/5.12anno/masked/${id}.fa.masked > tra/transcripts.${id}.fasta

/public/home/fengting/tools/TransDecoder-v5.5.0/util/gtf_to_alignment_gff3.pl sbam/$id.gtf > tra/transcripts.${id}.gff3

TransDecoder.LongOrfs -t tra/transcripts.${id}.fasta
TransDecoder.Predict -t tra/transcripts.${id}.fasta

/public/home/fengting/tools/TransDecoder-v5.5.0/util/cdna_alignment_orf_to_genome_orf.pl \
     trans_wild12.fa.transdecoder.gff3 \
     transcripts.wild12.gff3 \
     trans_wild12.fasta > out/trs.${id}.gff3

输出结果：
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格式文件