最近做了一些生物資訊的指令碼練習。
這是第一個例子。
找出一個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()