生物信息脚本练习（1） 找出fasta文件中大于500的序列

最近做了一些生物信息的脚本练习。
这是第一个例子。
找出一个fasta文件中大于500的序列，并重定向到另一个新的文件中。
这个文件每条序列是如下的样子。

c100027.graph_c0|orf3 type=complete len=150nt loc=c100027.graph_c0:123-272:-
ATGAGGATCTTTACGCCAAATGAGGGCCTTGTTGTTGATCTTTTCAGTAAGAGACGTTGC
GGAAATATTGGCGAAAATTTAAGAAACCTAGTTAGTTTTAGCAGCCACATAGTTGGCAAG
AATTATACTATTCAAGCCATGTGGCACTGA

这是我一开始的解法：
思路是用正则匹配第一行中len的数值。找到他们之后，根据这数字算出这条序列有多少行，然后把这么些行的信息输出。

import reoutput = open('data.txt', 'w')seq = []element = []row = []with open ("d:/Zanthoxylum_Bungeanum_Maxim.Unigene.cds.fa","r") as f:    for l in f:        seq.append(l)for i in seq:    if re.match(">.*len=[5-9][0-9][0-9]nt.*",i) or re.match(">.*len=[0-9][0-9][0-9][0-9]nt.*",i):        element_number = seq.index(i)   #>500序列的标签所处的位置        element.append(element_number)                            row_number = re.search("len=\d*nt", i).group(0)        index = row_number[4:7]                row_number = (int(index)// 60) + 1   #每个>500的序列的行数        row.append(row_number)len_row = len(row)def f(n):    #输出一个完整的带序列标签的序列 ，共找到n条符合条件序列    x = 0      for i in seq:        if x <= row[n]:              output.write(seq[element[n]+x]) #输出>之后的碱基序列            x += 1    returnfor i in range(0,len_row):    f(i)output.close()

这是我的另外一种解法：

先把序列整合成一个字符串，正则找到之后整个输出，而不是每行输出。

import refw=open('use.fa','w')with open("data1.fa","r") as f:    line = f.readlines()for i in line:    if i[0]!= ">":        i  = i.strip("\n")    else:         i = i.replace(">","\n>")    fw.write(i) fw.close()ttt = []with open("use.fa","r") as f:    seq = f.readlines()    seq = seq[1:]    for i in seq:        if seq.index(i)%2 ==0:            a = re.search("len=\d+",i).group(0)            ttt.append(a[4:]) #这个列表只包含所有序列的长度值print(ttt)with open('temp.fa','w') as qq:    qq.write(seq[0])    for i in ttt:        if int(i)>500:            #qq.write(seq[ttt.index(i)])            qq.write(seq[ttt.index(i)+1])

其实python和perl一样，“答案都不止一个”
期待更好的解法

# 8.14更新。我有了更好的解法# 这种把fasta文件转化成字典的方法来自这个论坛，感谢。# http://www.biotrainee.com/forum-59-1.htmlimport osimport reos.chdir("c:/编程题")def readfasta(filename):    fa = open(filename, 'r')    res = {}    ID = ''    for line in fa:        if line.startswith('>'):                       ID = line#.strip('\n')            res[ID] = ''        else:                        res[ID] += line#.strip('\n')     return resoutput = {}res = readfasta('500.fa')regex = re.compile(r'=\d+')for k,v in res.items():    m = regex.findall(k)    for n in m:        if int(n[1:])> 500:            output[k] = vf = open("output.txt" , "w")for k,v in output.items():    f.write(k)    f.write(v)f.close()