一个简易版提取最长转录本的方法

根据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