今天想到一个问题：
我如何把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
