56.《Bioinformatics Data Skills》之

Fastq格式作为fasta格式的延伸，提供每个碱基质量打分，常用于存储高通量测序数据。

fastq格式

fastq文件长这样：

egfr_flank.fasta

其中：

1. 以@开头，代表描述行，包括序列标志名与其它（可选）信息
2. 碱基序列，包括一行或多行
3. +号代表碱基序列结束，老式的fastq可能会重复描述行（冗余的信息可能导致文件过大）
4. 碱基质量得分，以ASCII码存储，数目与碱基一致（详细介绍见下一节）

潜在问题

1. 与fasta同样，保持你的描述行保持规范统一，常用的约定是以第一个空格将描述分为序列标志名与注释。

2. 碱基质量得分ASCII码包括了@符号，作为碱基质量得分的@符号有可能刚好出现在开头。编写脚本依赖匹配@开头作为header是易错的做法，最好使用完善的转换工具（例如readfq，后续会介绍），它的原理是使用碱基质量得分数目和碱基数目一致作为匹配header依据。

统计fasta与fastq序列数目

本节数据地址

fasta文件“>”符号只会出现在header，只需统计以“>”开头的行的数目就可以了：

$ grep -c "^>" egfr_flank.fasta
5

注意：使用引号是安全的做法，代码出错使>被解释为重定向可能损坏源文件。

当我们尝试以类似的方式统计fastq行数时，得到如下结果：

$ grep -c "^@" untreated1_chr4.fq
208779

查看文件可知每个序列占用4行，所以实际上序列的数目为行数/4（817420/4=204355）：

wc -l untreated1_chr4.fq
817420 untreated1_chr4.fq

两者结果并不一致，这是前者@开头的碱基质量得分被误认为header导致的。

当然，依赖行数的做法并不鲁棒，不同的序列可能拥有不同的行数，更好的做法是使用我们之前提到的bioawk工具(详细介绍)：

$ bioawk -c fastx "END{print NR}" untreated1_chr4.fq
204355