利用Python解决生物问题-根据指定基因ID提取序列

• 今天花了挺久时间写的一个序列提取的小程序，运行成功了，但可能在效率和实现方面存在不足，以后再改进，并希望大佬们提供宝贵的指导意见以及思路

准备文件

1.存放基因id号的txt文件
2.某物种的全部蛋白序列

生成文件

生成所需基因的序列文件

代码实现一

实现思路：

1.将所需要的基因ID存放于列表中，gene_list
2.将全部序列的fasta文件按行存放于列表中,all_seq_list
3.获取对应基因的序列在all_seq_list列表中的索引值(start和end)，并存放于列表中，seq_need_index
4.根据索引值，将需要基因的序列提取出来存放于列表中，seq_need_list
5.将seq_need_list的内容写入新文件中

###需要两个初始文件：指定id号的txt文件；序列fasta文件
###需要创建一个存储需要id序列的新文件
def extract_seq(geneid,allseq,seq_need):
    geneid_list = []
    #将基因id存放在列表中
    with open(geneid,'r') as f1:
        for each in f1:
            each = '>' + each
            geneid_list.append(each)
    #print(geneid_list)
    
    ###将全部序列存放于列表中
    all_seq_list = []
    with open(allseq,'r') as f2:
        for eachline in f2:
            all_seq_list.append(eachline)
    
    ###获取对应id序列在all_seq_list的索引值
    seq_need_index = []
    for each in geneid_list:
        start_index = all_seq_list.index(each)
        seq_need_index.append(start_index)
        for i in all_seq_list[start_index + 1:]:
            if '>' in i:
                end_index = all_seq_list.index(i)
                seq_need_index.append(end_index)
                break
    #print(seq_need_index)

    ###根据索引信息存放于列表中
    seq_need_list = []
    for n in range(0,len(seq_need_index),2):
        sindex = seq_need_index[n]
        eindex = seq_need_index[n+1]
        for m in range(sindex,eindex):
            seq_need_list.append(all_seq_list[m])
    #print(seq_need_list)


    ###将所需要的序列写入新文件中
    with open(seq_need,'w') as f3:
        for a in seq_need_list:
            f3.writelines(a)
        
              
geneid = input('请输入所需ID的文件名：')
allseq = input('请输入序列信息的文件名：')
seq_need = input('请输入生成文件的文件名：')
extract_seq(geneid,allseq,seq_need)
                

代码实现二-来自梁院士优化版改进

实现思路-无需获取索引值

1.将所需要的基因ID存放于列表中，gene_list
2.将全部序列的fasta文件按行存放于字典中seq_dict，若每行开头为‘>’,则为键，其余每行自动添加到对应键的值中

代码实现

###需要两个初始文件：指定id号的txt文件；序列fasta文件
###需要创建一个存储需要id序列的新文件
def extract_seq(geneid,allseq,seq_need):
    #将基因id存放在列表中
    with open(geneid,'rt') as f1:
        geneid_list=f1.read().strip("\n").split("\n")#俩种方法的前后顺序不能反，因为split()后返回列表
    #print(geneid_list)

    ###将全部序列和id存放于字典中
    seqdict={}
    allidlist = []
    with open(allseq,'rt') as infile:
        for line in infile:
            if line[0]=='>':
                key=line.strip()
                seqdict[key]=''
                allidlist.append(line[1:].strip())
            elif len(line.strip())>0:
                seqdict[key]+=line#若为line.strip()，最后序列没有换行符为1行


    ###获取需要的序列bing写入新文件中,注释掉的部分没有按照基因顺序提取
    #with open(seq_need,'wt') as f3:
        #for seqname in seqdict:
            #if seqname in geneid_list:
                #f3.write(f"{seqname}\n{seqdict[seqname]}\n")

    with open(seq_need,'wt') as outfile:
        for eachid in geneid_list:
            if eachid in allidlist:
                outfile.write(f'>{eachid}\n{seqdict[">"+eachid]}')
        
              
geneid = input('请输入所需ID的文件名：')
allseq = input('请输入序列信息的文件名：')
seq_need = input('请输入生成文件的文件名：')
extract_seq(geneid,allseq,seq_need)