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("\n"):
        if not line.strip():
            continue
        if line.startswith(">"):
            sequence_name = line.rstrip('\n').replace(">", "")
        else:
            sequence = line.rstrip('\n')
        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("\n"): 根据换行符进行分割

第一个条件语句主要用于剔除空白行

if not line.strip():
    continue

第二个条件语句，提问者的目的是对每一行进行分类，如果以”>”开头当作序列名，否则就是序列。

if line.startswith(">"):
    sequence_name = line.rstrip('\n').replace(">", "")
else:
    sequence = line.rstrip('\n')

第三个条件语句， 提问者的意图是判断序列是不是在字典中，去重

if sequence_name not in fasta:
    fasta[sequence_name] = []

fasta[sequence_name].append(sequence): 由于同一个基因序列是可以是多行存放的，所以需要追加在后面。

寻找真凶

那么问题来了，道理是哪里出错了？
以下列序列为例，对代码进行演示：

>seq1
ATGGGTGTGTGTGTG
>seq2
ATGTGTTTGTGTGCTCCTCCTC
>seq3
AACGTCGTGACGGGTGCGTGGTGTGTGTCCAA

> seq1 由于非空，跳过第一个条件语句。因为以”>”开始，那么sequence_name就是seq1.因为seq1不在字典中，所以内容为空，并且append了一个空字符串。因此可以解释为什么seq1里面有一个是空的。
ATGGGTGTGTGTGTG过来的时候，不是”>”开头，那么就是序列了，填入到了seq1字典的值中。
>seq2 由于非空，所以是序列名，添加到字典中，然后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 = '''
>seq1
ATGGGTGTGTGTGTG
>seq2
ATGTGTTTGTGTGCTCCTCCTC
>seq3
AACGTCGTGACGGGTGCGTGGTGTGTGTCCAA
'''
import sys
fasta = {}
seq_name = ''
sequence = ''
for line in test.split(sep='\n'):
    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