注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

Talk To Myself

Wubin's Bioinformatics Life

 
 
 
 
 

日志

 
 

生物序列基本格式fasta格式的解析  

2011-12-24 23:17:45|  分类: Script |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
文 / 屈武斌 (quwubin@gmail.com)

正常情况下,这是一个很简单的问题,因为Biopython中有现成的解析包。

当然,这里说的是非一般的情况。比如,如果我们所要解析的序列并非来自NCBI等标准数据库,而是由普通生物学者自己提供的数据,里面可能会存在各种很难想到的问题;或者我们有其他的需求。列举如下:

  1. 序列中有若干行空行
  2. >后面是空的(也就是说解析出来的id是None)
  3. 只有>行(即注释行),没有序列行
  4. 序列ID不唯一
  5. 我们项目中仅仅需要一个序列解析的模块,而不需要一个庞大的Biopython包
  6. 对解析速度的要求

以上的这些情况,是我们在开发MFEPrimer, MFEprimer-2.0, MPprimer, VizPrimer等项目中遇到的问题,作为在线免费供用户使用的生物信息学网站,在收集用户输入的引物序列时,不能保证用户的输入都像NCBI那样规范,所以有必要自己写模块来应付各种情况。

下面是我写的Python模块,代码很简单,因为fasta格式本身很简单。但是请注意,代码中没有处理ID不唯一的情况,并且逐个yield返回结果,而非一次性返回为结果列表,这对解析大容量数据库有利。

#!/usr/bin/python
# -*- coding: utf-8 -*- #
'''A Fasta format Parser return Iterator

USAGE:

import FastaIterator
infile = sys.argv[1]
records = FastaIterator.parse(open(infile))
for record in records:
print record.id
print record.desc
print record.seq
print record.size

Author: Wubin Qu <http://quwubin.sinaapp.com/>

ChangeLog

2011-12-24
* Use fh.next() to replace fh.readline()
* Now, can handle blank line in fasta record

2011-11-16
* Use iterator

2010-2-23
* v1.0 released.

'''

import sys

def parse(fh):
'''
A Fasta-format Parser return Iterator
'''
# Remove the comment and blank lines before the first record
eof = False
while True:

try:
line = fh.next()
# Note: can not be: fh.next().strip()
except StopIteration:
print >> sys.stderr, "No fasta record found, may be a empty file?"
exit()

line = line.strip()
if not line: continue

if line.startswith('>'):
break

while True:
if not line.startswith('>'):
print >> sys.stderr, "Records in Fasta files should start with '>' character"
exit()

id, sep, desc = line[1:].partition(' ')
if len(id) < 1:
print >> sys.stderr, "No id found"
exit()
seq_lines = []

try:
line = fh.next()
# Note: can not be: fh.next().strip()
except StopIteration:
print >> sys.stderr, "Error fasta format: no sequence line after '>' line"
exit()

line = line.strip()
if not line: continue

while True:

if line.startswith('>'):
break

seq_lines.append(line.replace(' ', '').replace("\r", ''))

try:
line = fh.next()
# Note: can not be: fh.next().strip()
except StopIteration:
eof = True
break

line = line.strip()
if not line: continue

yield Iterator(id, desc, ''.join(seq_lines))

if eof:
return

assert False, 'Should not reach this line'

class Iterator:
'''Create the class'''
def __init__(self, id, desc, seq):
self.id = id
self.desc = desc
self.seq = seq
self.size = len(seq)

def main ():
'''Main'''
import FastaIterator
infile = sys.argv[1]
records = FastaIterator.parse(open(infile))
for record in records:
print record.id
print record.desc
print record.seq
print record.size

if __name__ == '__main__':
main()


  评论这张
 
阅读(790)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017