最长mRNA、CDS、Protein（Python实现）

# !usr/bin/env python3
# -*- coding:utf-8 -*-
"""
@FileName: get_longest
@Time: 2022/3/25,19:12
@Motto: go go go 
"""
import argparse
from Bio import SeqIO  #   pip install  biopython


def read_file(file):
    t = {}  # 记录长度和序列名字
    result = {}  #这个字典用于储存最长转录本 、最长cds、最长protein
    for seq_record in SeqIO.parse(file, "fasta"): # 用biopython模块解析文件
        id = seq_record.id.rsplit(".", 1)[0]

        if id not in t:
            result[seq_record.id] = str(seq_record.seq)
            t[id] = [len(seq_record.seq), seq_record.id]
        else:
            if t[id][0] >= len(seq_record.seq):
                continue
            else:
                result.pop(t[id][1])
                result[seq_record.id] = str(seq_record.seq)
                t[id] = [len(seq_record.seq), seq_record.id]
    return result


def write(filename, res):
    with open(filename,'w') as f:
        for i, j in res.items():
            f.write(">" + i + "\n")
            f.write(j + "\n")


def main():
    parser = argparse.ArgumentParser(usage='********', description='得到最长结果')
    parser.add_argument("-i", "--input", help="input filename")
    parser.add_argument("-o", "--output", help="output filename")
    args = parser.parse_args()
    res_dict = read_file(args.input)
    write(args.output, res_dict)


if __name__ == '__main__':
    #res_dict = read_file(r'./cds.fa')
    #write(r'out_cds', res_dict)
    main()

现在手里有一个mRNA的文件，分析需要同一基因中最长mRNA。我的想法是先建立两个字典。result保存我要的结果，t辅助用来记录基因名字和序列长度。用.右边分割一次，得到基因ID序列aaaaa.01G000100，如果ID没有在t 中，result 中加入 序列名aaaaa.01G000100.t1：序列 ，t 中加入 ID：[序列名字aaaaa.01G000100.t1,序列长度]。 如果ID在t 中，比较长度，保留长的。