由于可變剪切一個(gè)基因得到好多個(gè)序列長(zhǎng)度不同的轉(zhuǎn)錄本，應(yīng)該挑選出序列最長(zhǎng)的最長(zhǎng)轉(zhuǎn)錄本作為數(shù)據(jù)的分析
看一下有規(guī)律的蛋白文件命名，一般是以下形式
寫(xiě)個(gè)腳本提取
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
#生成結(jié)果
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軟件
最長(zhǎng)轉(zhuǎn)錄本（Longest Transcript）指的是一個(gè)基因所產(chǎn)生的RNA分子中最長(zhǎng)的那個(gè)。這個(gè)RNA分子可能是mRNA（信使RNA），也可能是其他類(lèi)型的RNA，比如pre-mRNA（前體mRNA）、snRNA（小核RNA）或miRNA（微小RNA）等。最長(zhǎng)轉(zhuǎn)錄本的長(zhǎng)度取決于基因的轉(zhuǎn)錄調(diào)控和剪接等機(jī)制。
最長(zhǎng)ORF（Open Reading Frame）指的是一個(gè)RNA分子中具有最長(zhǎng)的連續(xù)未被中斷的編碼區(qū)域，一般指的是可以被翻譯成蛋白質(zhì)的區(qū)域。ORF通常從起始密碼子（start codon）開(kāi)始，一直延伸到終止密碼子（stop codon）結(jié)束。在一個(gè)RNA分子中可能存在多個(gè)ORF，而其中最長(zhǎng)的ORF通常被認(rèn)為是最可能被翻譯成蛋白質(zhì)的部分。
因此，最長(zhǎng)轉(zhuǎn)錄本和最長(zhǎng)ORF的區(qū)別在于，前者指的是RNA分子中整體長(zhǎng)度最長(zhǎng)的那個(gè)，而后者指的是RNA分子中具有最長(zhǎng)編碼潛力的連續(xù)區(qū)域。最長(zhǎng)轉(zhuǎn)錄本可以包含一個(gè)或多個(gè)ORF，但最長(zhǎng)ORF不一定完全重疊于最長(zhǎng)轉(zhuǎn)錄本中。研究這兩個(gè)概念可以幫助科學(xué)家更好地理解基因的轉(zhuǎn)錄和翻譯過(guò)程，以及基因產(chǎn)物（蛋白質(zhì)）的功能。
需要準(zhǔn)備文件·pep.fasta
TransDecoder通過(guò)運(yùn)行一個(gè)包含目的轉(zhuǎn)錄本序列的fasta文件來(lái)實(shí)現(xiàn)功能。簡(jiǎn)單的用法如下：
提取最長(zhǎng)的開(kāi)放閱讀框
TransDecoder.LongOrfs -t trans_wild12.fa
預(yù)測(cè)可能的編碼區(qū)
TransDecoder.Predict -t trans_wild12.fa
候選編碼區(qū)的最終集合可以在文件.transdecoder中找到。擴(kuò)展名包括.pep，.cds，.gff3和.bed。
從有參考基因組的轉(zhuǎn)錄結(jié)果GTF文件預(yù)測(cè)編碼區(qū)域：
1、需要有參轉(zhuǎn)錄組比對(duì)后拼接的轉(zhuǎn)錄本的GTF文件以及參考基因組序列：
2、將GTF文件轉(zhuǎn)化為GFF3文件
3、提取、預(yù)測(cè)
4、最后生成一個(gè)基于有參基因組的編碼區(qū)域注釋文件:
/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
輸出結(jié)果：
longest_orfs.pep : 所有達(dá)到最小長(zhǎng)度標(biāo)準(zhǔn)的ORF, 不管是否編碼
longest_orfs.gff3 : 在目的轉(zhuǎn)錄本中發(fā)現(xiàn)的所有ORF的位置
longest_orfs.cds : 所有檢測(cè)到的ORF的核酸編碼序列
longest_orfs.cds.top_500_longest : 前500個(gè)最長(zhǎng)的ORF，用于訓(xùn)練一個(gè)編碼序列的馬爾科夫模型
hexamer.scores : 每個(gè)k-mer的對(duì)數(shù)似然得分 (coding/random)
longest_orfs.cds.scores : 每個(gè)ORF同6個(gè)閱讀框間對(duì)數(shù)似然得分的總和
longest_orfs.cds.scores.selected : 根據(jù)得分標(biāo)準(zhǔn)所選出的ORF
longest_orfs.cds.best_candidates.gff3 : 轉(zhuǎn)錄本中選出的ORF的位置
最后的輸出文件在你當(dāng)前的工作目錄中。
transcripts.fasta.transdecoder.pep : 最終候選ORF的蛋白質(zhì)序列；所有較長(zhǎng)ORF中的較短的候選序列已被移除。
transcripts.fasta.transdecoder.cds : 最終候選ORF的編碼區(qū)的核酸序列。
transcripts.fasta.transdecoder.gff3 : 最終被選中的ORF在目的轉(zhuǎn)錄本中的位置
transcripts.fasta.transdecoder.bed : 用來(lái)描述ORF位置的bed格式文件
