查看原文
其他

Python如何解析fasta格式,并储存为字典?

2017-07-06 徐洲更 生信媛

今天想到一个问题:

我如何把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

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存