Python如何解析fasta格式，並儲存為字典？

今天想到一個問題：

我如何把fasta文件讀到Python裡面，存為字典格式。

由於技術不過關，我折騰了半天還不知道如何做。所以就去Google，很多人提議用biopython。但是讀取序列其實是一件比較基礎的能力，而且任務也比較簡單， 用基礎庫的話比較通用，也能鍛煉自己的編程能力。

最後發現有一個人在stackoverflow提問，為何他寫的代碼結果不是預期那樣：

import sys

sequence = " "

fasta = {}

with open(sys.argv[1]) as file_one:

file_one_content = file_one.read()

for line in file_one_content.split("
"):

if not line.strip():

continue

if line.startswith(">"):

sequence_name = line.rstrip("
").replace(">", "")

else:

sequence = line.rstrip("
")

if sequence_name not in fasta:

fasta[sequence_name] = []

fasta[sequence_name].append(sequence)

print fasta

結果如下：

{"seq3": ["ATGTGTTTGTGTGCTCCTCCTC", "AACGTCGTGACGGGTGCGTGGTGTGTGTCCAA"], "seq2": ["ATGGGTGTGTGTGTG", "ATGTGTTTGTGTGCTCCTCCTC"],

"seq1": [" ", "ATGGGTGTGTGTGTG"]}

代碼存在問題，但是和我想要的結果類似，所以改改錯誤就行了，首先知道他到底錯在哪裡，先看看提問者的思路：

代碼逐行解釋

import sys

： 導入sys, 通過sys.argv讀取命令行參數

fasta = {}

：定義了一個字典，fasta用於存放數據。

sequence = " "

： 定義了一個含有一個空格的字元串，，sequene用於存放序列

with open(sys.argv[1]) as file_one:

:利用sys.args[1]讀取文件，然後用open打開，成為python里的file object，用於IO操作。

file_one_content = file_one.read():

就是把所有文件讀到了一個字元串里。

for line in file_one_content.split("
"):

根據換行符進行分割

第一個條件語句主要用於剔除空白行

if not line.strip():

continue

第二個條件語句，提問者的目的是對每一行進行分類，如果以」>」開頭當作序列名，否則就是序列。

if line.startswith(">"):

sequence_name = line.rstrip("
").replace(">", "")

else:

sequence = line.rstrip("
")

第三個條件語句， 提問者的意圖是判斷序列是不是在字典中，去重

if sequence_name not in fasta:

fasta[sequence_name] = []

fasta[sequence_name].append(sequence)

: 由於同一個基因序列是可以是多行存放的，所以需要追加在後面。

尋找真兇

那麼問題來了，道理是哪裡出錯了？

以下列序列為例，對代碼進行演示：

ATGGGTGTGTGTGTG

ATGTGTTTGTGTGCTCCTCCTC

AACGTCGTGACGGGTGCGTGGTGTGTGTCCAA

由於非空，跳過第一個條件語句。因為以」>」開始，那麼

sequence_name

就是seq1.因為seq1不在字典中，所以內容為空，並且append了一個空字元串。因此可以解釋為什麼seq1裡面有一個是空的。

ATGGGTGTGTGTGTG

過來的時候，不是」>」開頭，那麼就是序列了，填入到了seq1字典的值中。

由於非空，所以是序列名，添加到字典中，然後

fasta[sequence_name].append(sequence)

, 就會插入剛才的

ATGGGTGTGTGTGTG

。接著添加了自身的

ATGTGTTTGTGTGCTCCTCCTC

。

所以錯誤就是，字典賦值不應該直接放在循壞中，應該放在條件語句中。

正確的代碼

知道代碼，錯在哪裡之後，就可以根據這個邏輯寫出正確的代碼。由於file object已經時可迭代對象，不需要用file.read()一次性讀取所有字元然後分割，也不需要用file.readlines()讀取每一行。

優化後的代碼如下：

import sys

fasta = {}

seq_name = ""

sequence = ""

for line in os.open(sys.argv[1]):

if not line.strip():

continue

if line.startswith(">"):

seq_name = line

else:

sequence = line

if seq_name not in fasta:

fasta[seq_name] = []

continue

fasta[seq_name] = sequence

測試下效果:

test = """

ATGGGTGTGTGTGTG

ATGTGTTTGTGTGCTCCTCCTC

AACGTCGTGACGGGTGCGTGGTGTGTGTCCAA

"""

import sys

fasta = {}

seq_name = ""

sequence = ""

for line in test.split(sep="
"):

if not line.strip():

continue

if line.startswith(">"):

seq_name = line

else:

sequence = line

if seq_name not in fasta:

fasta[seq_name] = []

continue

fasta[seq_name] = sequence

結果如下：

當然原問題也有人回答了，感覺他寫的更加好看一點。

import sys

fasta = {}

with open(sys.argv[1]) as file_one:

for line in file_one:

line = line.strip()

if not line:

continue

if line.startswith(">"):

active_sequence_name = line[1:]

if active_sequence_name not in fasta:

fasta[active_sequence_name] = []

continue

sequence = line

fasta[active_sequence_name].append(sequence)

print fasta

想要系統學習python和免費學習資料的 可以加裙 四七四五三四九五一

喜歡這篇文章嗎？立刻分享出去讓更多人知道吧！