根据gff和基因组文件提取最长转录本有很多工具，但发现没有一个适合我的，因为那些工具有各种输入文件限制，例如GetTransTool。
后来我发现gff文件里gene间的mRNA，最开始的mRNA就是最长转录本，当我只提出gene和mRNA时，第一条mRNA的长度等于基因的长度且长于后面出现的
gff.file
因此，策略就是：先处理gff文件只提取gene和RNA的行（看自己需求，只要mRNA或者包含非编码RNA的也可以），然后用Python写脚本得到最长转录本的list，使用gffread提取的或者下载的所有氨基酸序列pep.fa文件，利用seqkit grep提取该文件中的最长转录本序列。
输入文件gff文件也要先处理一下：
inputfile为gff.gene.mRNA.txt，是一个根据gff文件提取只要gene和RNA的行且只要1,3,4,5,7,9列，并将最后一列（第9列）的转录本ID和gene的ID根据分号分隔（只要分隔的前两列后面的东西都不要）。
##不只是mRNA，还有非编码RNA也包含在内，根据你的注释文件而定#
awk 'BEGIN{FS=OFS="\t"}{if($3=="gene"||$3=="mRNA"||$3=="lnc_RNA")print $1,$3,$4,$5,$7,$9}' gff.file |awk '{ n=split($NF, arr, ";"); $NF=""; for (i=1; i<NF; i++) printf "%s\t", $i; print arr[1] "\t" arr[2] }' FS='\t' OFS='\t' > input_file
处理完的input file如下图所示，一个基因对应很多行转录本的文件：
input file
下面是运行脚本gff_trans_finder.py
#!/usr/bin/python
# -*- coding: utf-8 -*-
# conda activate python3
import sys
def usage():
print('Usage: python gff_trans_finder.py [input_file]')
def main():
inf = open(sys.argv[1], 'rt')
ouf = open("gff.gene.longest.txt", 'wt') #保存只有gene和最长转录本的行
longest_trans_list = [] # 提取最长转录本的list
gene_list = {} # 用字典存储gene和transcript
flag_list = []
longest_trans_list = {}
for line in inf:
line = line.strip()
name = line.split("\t")[6]
ID = line.split("\t")[5].split("=")[1]
type = line.split("\t")[1]
if name.startswith('Parent'): #提取基因名
gene = name.split("-")[1]
transcript = ID
if gene not in gene_list:
ouf.write(line + '\n') # 将符合条件的行写入 ouf 文件
gene_list[gene] = [] # 将每个基因单独创建一个转录本列表
longest_trans_list[gene] = transcript #判断gene没有重复的，即出现的第一条RNA是最长转录本
gene_list[gene].append(type)
gene_list[gene].append(transcript) #判断gene是重复的，放入同一个gene列表里
else:
ouf.write(line + '\n') # 将符合条件的行写入 ouf 文件
with open('longest_gene_trans.txt', 'wt') as ouf:
for gene, transcripts in longest_trans_list.items():
output_line = gene + "\t" + transcripts + "\n"
ouf.write(output_line)
with open('gene_trans.txt', 'wt') as ouf:
for gene, transcripts in gene_list.items():
# 将基因和转录本列表以制表符分隔写入文件
output_line = gene + "\t" + "\t".join(transcripts) + "\n"
ouf.write(output_line)
inf.close()
ouf.close()
try:
main()
except IndexError:
usage()
利用上面脚本处理完后，得到如下图所示的一系列文件gff.gene.longest.txt
gff.gene.longest.txt
基因及其对应的最长转录本longest_gene_trans.txt:
longest_gene_trans.txt
基因及其所有转录本gene_trans.txt：
gene_trans.txt
之后根据longest_gene_trans.txt可以提取最长转录本序列，genomic.gff.pep.fa是使用gffread提取的protein序列
cut -f 2 longest_gene_trans.txt > longest_transcript.txt
seqkit grep -f longest_transcript.txt genomic.gff.pep.fa > genome.longest.pep.fa
