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