FASTA/FASTQ文件的解析看起来是一件比较容易的事情,然而如果你尝试亲自编写代码的话就会发现很难全面地考虑问题。更好的方式是使用开源的库,这些库经过了广大用户的测试,更加鲁棒,同时速度更快。
有很多工具可以实现这一目标,这里介绍一个由Li Heng开发的独立的解析FASTA/FASTQ文件的python脚本,叫做readfq.py
(下载地址)。函数的输入输出分别为:
- FASTX的文件对象(文件对象可以由
open([file])
或者标准输入流创建) - 返回值为一个元组,内容为序列的描述,序列,序列质量(FASTA文件序列质量为
None
)
这里通过一个例子来展示此脚本:统计FASTA文件碱基数目
数据下载地址:https://github.com/vsbuffalo/bds-files/tree/master/chapter-10-sequence
#!/usr/bin/env python
# -*- coding:utf-8 -* #1
# count_nucleotides.py -- 清点碱基数目
import sys
from collections import Counter #2
from readfq import readfq #3
IUPAC_BASES = "ACGTRYSWKMBDHVN-." #4
counts = Counter()
for name, seq, qual in readfq(sys.stdin): #5
counts.update(seq.upper())
for base in IUPAC_BASES:
print base + "\t" + str(counts[base])
解释:
#1:由于注释有中文,需要在头部加这行
#2:这里引入Counte
对象,此对象类似于dict
对象,增加的新特性方便统计数字
#3:通过readfq.py
文件引入readfq
函数
#4:所有可能的碱基字符,由于Counter对象没有顺序,这里列出方便按顺序展示结果
#5:这里调用readfq
函数,采用元组拆包的方式定义序列seq对象,使用update
函数统计字符的出现数目,这里将所有字母转为大写,所以soft-masked字符也会被统计(这里也可以不使用Counter而采用for循环来统计)
python脚本可以使用python ./count_nucleotides.py
的方式运行,也可以在脚本开头加#!/usr/bin/env python
配置,然后使用chmod +x count_nucleotides.py
赋予运行权限。这里以后者运行程序:
$ cat contam.fastq | ./count_nucleotides.py
A 103
C 110
G 94
T 109
R 0
Y 0
S 0
W 0
K 0
M 0
B 0
D 0
H 0
V 0
N 0
- 0
. 0
以上就是这个脚本的全部内容,此脚本有很大的改进空间,例如增加参数读取文件名而不是使用标准输入流,或者增加特性可以统计每条序列的碱基,soft-masked的碱基,CpG岛的数目,验证输入是否有效等等。这个脚本里面最复杂的部分实际上是readfq
函数,通过这个例子也说明了代码重用的重要性,我们没有必要重新实现复杂的函数,这也是专家的编程思路。