根據(jù)gff和基因組文件提取最長轉(zhuǎn)錄本有很多工具，但發(fā)現(xiàn)沒有一個適合我的，因為那些工具有各種輸入文件限制，例如GetTransTool。

后來我發(fā)現(xiàn)gff文件里gene間的mRNA，最開始的mRNA就是最長轉(zhuǎn)錄本，當我只提出gene和mRNA時，第一條mRNA的長度等于基因的長度且長于后面出現(xiàn)的

gff.file

因此，策略就是：先處理gff文件只提取gene和RNA的行（看自己需求，只要mRNA或者包含非編碼RNA的也可以），然后用Python寫腳本得到最長轉(zhuǎn)錄本的list，使用gffread提取的或者下載的所有氨基酸序列pep.fa文件，利用seqkit grep提取該文件中的最長轉(zhuǎn)錄本序列。

輸入文件gff文件也要先處理一下：
inputfile為gff.gene.mRNA.txt，是一個根據(jù)gff文件提取只要gene和RNA的行且只要1,3,4,5,7,9列，并將最后一列（第9列）的轉(zhuǎn)錄本ID和gene的ID根據(jù)分號分隔（只要分隔的前兩列后面的東西都不要）。

##不只是mRNA，還有非編碼RNA也包含在內(nèi)，根據(jù)你的注釋文件而定#
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如下圖所示，一個基因?qū)芏嘈修D(zhuǎn)錄本的文件：

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和最長轉(zhuǎn)錄本的行 

    longest_trans_list = []  # 提取最長轉(zhuǎn)錄本的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] = []      # 將每個基因單獨創(chuàng)建一個轉(zhuǎn)錄本列表
                longest_trans_list[gene] = transcript     #判斷gene沒有重復的，即出現(xiàn)的第一條RNA是最長轉(zhuǎn)錄本
                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():
            # 將基因和轉(zhuǎn)錄本列表以制表符分隔寫入文件
            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

基因及其對應的最長轉(zhuǎn)錄本longest_gene_trans.txt:

longest_gene_trans.txt

基因及其所有轉(zhuǎn)錄本gene_trans.txt：

gene_trans.txt

之后根據(jù)longest_gene_trans.txt可以提取最長轉(zhuǎn)錄本序列，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