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

最近做了一些生物信息的脚本练习。

这是第一个例子。

找出一个fasta文件中大于500的序列，并重定向到另一个新的文件中。

这个文件每条序列是如下的样子。

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

ATGAGGATCTTTACGCCAAATGAGGGCCTTGTTGTTGATCTTTTCAGTAAGAGACGTTGC

GGAAATATTGGCGAAAATTTAAGAAACCTAGTTAGTTTTAGCAGCCACATAGTTGGCAAG

AATTATACTATTCAAGCCATGTGGCACTGA

这是我一开始的解法：

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

import re
output = 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
:
output.write(seq[element
+x]) #输出>之后的碱基序列
x += 1
return
for i in range(0,len_row):
f(i)
output.close()

这是我的另外一种解法：

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

import re
fw=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.html 
import os
import re
os.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 res
output = {}
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] = v

f = open("output.txt" , "w")
for k,v in output.items():
f.write(k)
f.write(v)
f.close()