在使用python做生物信息分析,尤其是DNA序列,蛋白质序列的分析过程中,经常会需要读入外部fasta文件,传统的读入方法是以文件读入加上字符串操作实现的,例如:
seqdict={}
with open("filename.fasta",'rt') as infile:
for line in infile:
if line[0]=='>':
key=line[1:].strip()
seqdict[key]=''
else:
seqdict[key]+=line.strip()
其实上面这段代码有个小bug:如果fasta文件中第一行是空行,就会直接触发else的条件,此时由于key尚未声明,则会出现报错。修改方法则是多增加一个判断是否为空的条件:
seqdict={}
with open("filename.fasta",'rt') as infile:
for line in infile:
if line[0]=='>':
key=line[1:].strip()
seqdict[key]=''
elif len(line.strip())>0:
seqdict[key]+=line.strip()
虽然以上方法已经足够简单了,但总归要写好几行,不够快速,最近发现一个python包protfasta,可以快速读入DNA和蛋白的fasta文件
安装protfasta:
pip install protfasta
读入fasta文件
import protfasta
seq=protfasta.read_fasta("filename.fasta")
在这里,返回结果的seq
是一个字典格式的变量,其中key是序列名,value是序列
不过这个方法也有一个不好的地方,对于python中普通的文件读写,当找不到目标文件时,python会抛出FileNotFoundError
,所以对于普通的文件操作,进行异常捕获时只需要:
try:
f=open("filename.txt",'rt')
except FileNotFoundError:
pass
但protfasta还定义了另一个Error:protfasta.protfasta_exceptions.ProtfastaException
因此在这里,用户需要捕获的异常就要换个名字了:
try:
seq=protfasta.read_fasta("filename.fasta")
except protfasta.ProtfastaException:
pass
由于protfasta在遇到其他异常,比如读入fasta文件过程中发现序列有非法字符的时候,抛出的也是这个error,这样做就使得异常处理的操作对所有异常都是统一的了,而不仅仅是不能找到文件
不过一般来说如果写个脚本是自己用来分析生信数据的话也不用考虑那么多乱七八糟的东西,直接protfasta.read_fasta()读入就好了,还是很方便的
关于protfasta的更多使用方法可以参照其说明文档:
https://protfasta.readthedocs.io/en/latest/
非常简洁,推荐尝试!