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

根据gff和基因组文件提取最长转录本有很多工具，但发现没有一个适合我的，因为那些工具有各种输入文件限制，例如GetTransTool。

后来我发现gff文件里gene间的mRNA，最开始的mRNA就是最长转录本，当我只提出gene和mRNA时，第一条mRNA的长度等于基因的长度且长于后面出现的，如下图：

1687873658594.jpg

于是就想到了之前一个提取最好blast比对best hit的脚本，就可以很好提取对应的最长转录本信息：这里改编后可以提取gene和最长转录本对应的名称list，以及gene和最长转录本位置信息。

输入文件gff文件也要先处理一下：
inputfile为gff.gene.mRNA.txt，是一个根据gff文件提取的1,3,4,5,7,9列，且只保留第三列gene和mRNA的行

awk 'BEGIN{OFS="\t"}{if($3=="gene"||$3=="mRNA")print $1,$3,$4,$5,$7,$9}' gff.file |awk -F ";" '{print$1"\t"$2}' > input_file

#!/usr/bin/python
# -*- coding: utf-8 -*-
# conda activate python3

import sys

def usage():
    print('Usage: python gene_gff.py [input_file] [outfile]')

def main():
    inf = open(sys.argv[1], 'rt')
    #inputfile为gff.gene.mRNA.txt，是一个根据gff文件提取的1,3,4,5,7,9列，且只保留第三列gene和mRNA的行
    #脚本为: awk 'BEGIN{OFS="\t"}{if($3=="gene"||$3=="mRNA")print $1,$3,$4,$5,$7,$9}' gff.file |awk -F ";" '{print$1"\t"$2}' > input_file

    ouf = open("gff.gene.longest.txt", 'wt')     #保存只有gene和最长转录本的行
    gene_trans = open("gene_trans.txt", 'wt')
    longest_trans_l = open("longest_transcript.txt", 'wt')

    longest_trans_list = []          #提取最长转录本的list
    gene_list = {}                   #用字典存储gene和transcript
    flag_list = []                  
    for line in inf:
        line = line.strip()
        name = line.split("\t")[6]
        ID = line.split("\t")[5].split("=")[1]
        if name.startswith('Dbxref'):             #设置第一个条件筛选gene_list
            ouf.write(line + '\n')
            gene = ID.split("-")[1]
        elif name not in flag_list:              #设置第二个条件筛选最长转录本                
            ouf.write(line + '\n')
            transcript = ID
            longest_trans_list.append(transcript)
            gene_list[gene] = transcript
        else:
            continue
        flag_list.append(name)
    
    longest_trans_l.write('\n'.join([f'{line}'for line in longest_trans_list]))
    gene_trans.write('\n'.join([f'{key}\t{value}' for key, value in gene_list.items()])) 
    inf.close()
    ouf.close()
    gene_trans.close()
    longest_trans_l.close()

try:
    main()
except IndexError:
    usage()

之后根据gene_trans.txt可以提取最长转录本序列，genomic.gff.pep.fa是使用gffread提取的protein序列

seqkit grep -f gene_trans.txt genomic.gff.pep.fa > genome.longest.pep.fa

太久没写Python脚本，差点忘光了。。。再加上chatgpt。。。