fastaq文件中reads数的统计

在汇总测序的基本信息时经常需要从fastaq文件中计算某个样本的原始的read数量。

fastaq文件的一般格式为:

@FCD056DACXX:3:1101:2163:1959#TCGCCGTG/1
TCCGATAACGCTCAACCAGAGGGCTGCCAGCTCCGATCGGCAGTTGCAACCCATTGGCCGTCTGAGCCAGCAACCCCGGA
+
gggiiiiiiiiiiiiiiiiiiiiiiiiiigggggeeecccccc^bcbcccccccbccccc]aaccbbccc^R^^acccc_

每四行一条记录, 每条记录以@符号开始,后面跟着的是该reads的标识符,第二行为对应的核酸序列,第三行为‘+’,第四行为序列的质量信息,以字符对应的ASCII码来表示对应碱基的测序质量。因此,想要知道有多少个reads只要统计有多少个只包含+的行就知道有多少个reads了。如下:

cat your.fastq.gz | grep -c '^+$'

有些人给出的方法是:cat your.fastq.gz | grep -c ‘+’  或 cat your.fastq.gz | grep -c ‘^+’ 这两种方法都是不严谨的,因为第四行的质量信息中也是包含了‘+’的,如果只是统计包含‘+’或以‘+’开头的行会多统计出一部分数据来,尤其是第一种情况。因此,正确的方式应该是计数只包含‘+’的行。此外要注意的是,如果是pair-end的数据,一个样本是包含了两个文件的,要放在一起计数。这和后面统计质控后的Clean reads以及计算后面比对reads的比例有关系。