生物資訊指令碼練習(1) 找出fasta檔案中大於500的序列
阿新 • • 發佈:2019-01-30
最近做了一些生物資訊的指令碼練習。
這是第一個例子。
找出一個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[n]:
output.write(seq[element[n]+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()