fastq格式文件及phred33的判断

Fastq格式文件储存了生物序列的信息及其质量信息。以电脑中的一个文件的为例

$ zcat SRR8518347.fq.gz |head

@SRR8518347.1 1 length=150
GTATTTTCTAAGACGTAATGCTACTTTACACTTAACAGACTACAGTATAGTGTAAACGTAACGGTTATGTCCACTAGGAAACCAAAATATTCGTGTTACTTGCTTTACTGTACTTGTTTTATTTCAGTAGTCGGGAACTAAACTCACAGT
+SRR8518347.1 1 length=150
AAA,,F<F,FKKAFFF<FFFKF<FA,KKFFKKKKKKKKKF,AKFK,A7<FKK,AFFKKFKKKA,AFF,FKAKKKA<FFKAFKKKKKA,F<KKKKFFKKKKK,FKKAKAKKKKF7FFFAKFAFFKK7,,FF<F,,<,AA<KAFKKKKKAF,
@SRR8518347.2 2 length=150
CCAGCTAATTTTTGTATTACTAGTAGAGACGGGGTTTCACCAGGTTGGTCAGGATGGTCCTGAACTCCTGAGCTCAGGCGATCCACCCACCTCAGCCTCCCAAAGTGGTAGGATTACAGGCGTGAGATCGGAAGAGCACACGTCTGAACT
+SRR8518347.2 2 length=150
AAAFAAFFFKKFFKKKKKKKKKKKKFKKA7FKKK(FFFKKKFFFKK,A,7F,FFKK,FK7AKKFKKK,FKA,AFF7<AA<AAFFFKFK<KKK,FKFFKKKKKKKKA,,7FKKKFKKK7FAFKAA7,A7FF,FFAK7<FAAFKKKK7FFF7

1 格式说明

第一行：必须@开头，紧跟唯一的序列的ID标识符，后面可跟其他描述性内容，但序列ID与描述部分空格分开，这点很重要，是很多分析的基础 第二行：序列的顺序，基本不可能出现重复。这里是核酸。 第三行：必须+开头，后面可以跟描述信息，此处是跟了第一行的信息。但为了节约空间，很多时候这一行只有一个+,而不跟任何内容。 第四行：碱基或氨基酸（此处碱基）的质量字符，对应着第二行的碱基，反应的是该碱基的错误率，所以这一行的字符数和第二行要一一对应，否则就乱了。这里就引入了ASCII code。

fastq format 1

fastq format 2

所以引入下面，碱基质量值是什么，如何获得？怎么表示？如何转换？

2 碱基质量值Q值和ASCII码之间的关系

因为第四行的编码，开始由Phred程序开发者定义，所以一般称为Phred quality。那碱基质量得分怎么来的？

Chromas软件展示的一个DNA 序列质量结果

每合成一个碱基，即可发出一个荧光信号，该信号可以被捕捉到，并生成是是轨迹数据。不同的碱基用不同颜色标记，检测相应峰值即可判断碱基。 而Phred通过计算相应波峰参数，去查询通过已知序列测序分析得到的一个表，即可把错误率转换为质量得分。也就是把波峰参数和质量得分对应起来。 碱基错误率与质量得分的关系如下

Phred quality score

也就是说，质量值Q是测序错误率的对数*-10。假如错误率是0.01，则Q值为20。可见，错误率越低，其Q值越高。即Q值越高越可靠。 Q10——0.1 Q20——0.01 Q30——0.001 Q40——0.0001 那就直接以上面这种Q值表示碱基质量不就行了吗？为什么要用ASCII码表示。 如果用数字表示，数字和数字之间需要有间隔符号以区分，再者会浪费存储空间，所以可以把质量值转变为相应的ASCII码，这样就完成了把质量数向ASCII码的转换，那现在看下ASCII码

ASCII码

如果直接把Q值直接对应ASCII码，应该挺方便的，但是，Q值有时会有负值，再者，看ASCII码的0-31位都是控制字符，没法打印和保存，能打印的从要从32位的Space开始，所以就可以给实际的Q值加上一个固定值，这样就可以打印出来而保存在fastq文件中了。 所以下面就是固定值加多少的问题。phred33，就是加了33，也就是10变成43，查表是+。例子中很多F,Q值是70-33=37。而当时的solexa加的是64，这就是Phred64.他们与ASCII的对应关系如下

Q score vs Phred

数据处理时，有些软件会根据碱基质量得分的不同做不同处理，所以需要指定正确的编码方式，也就是需要对Phred33还是64正确判断。不过现在基本都33了，但如果下载以前的数据不一定。 下面是不同版本质量得分和质量字符ASCII的关系

不同测序标记中的Phred的使用

从上面可以看出，Phred33的字符使用33-73，而+64使用包括59（包括）-104之间的ASCII码。所以，只要ASCII小于59的仅仅在Phred+33中使用，而+64的都大于59，而从表可以看到大于等于74的旨在Phred+74中使用，所以，如果软件没有自动判断，根据以上特点，就可以自行判断是Phred33还是Phred64。

3 如何判断是Phred33还是Phred64

默认读取1000条序列，在这1000条序列中：

• 1. 如果有2个以上的质量字符ASCII值小于等于58(即有两个碱基的得分小于等于25），同时没有任何质量字符的ASCII值大于等于75，即判断是Phred+33。
• 1. 如果有2个以上的质量字符ASCII值大于等于75(即有两个碱基的得分大于等于10），同时没有任何质量字符的ASCII值小于等于58，即判断是Phred+64。
• 1. 如果所有质量字符的ASCII值介于59到74之间，即判断可能是Phred+33，但建议使用更多的序列做进一步测试（出现这种结果可能有两种情况：1, Phred+33编码，所有碱基质量得分介于26到42之间；2，Phred+64编码，所有碱基质量得分介于-5到10；是前者的可能性大）。
• 1. 如果出现上述3种以外的情况，建议打印出质量字符的ASCII值人工判断。 脚本如下

$ cat fq_qual_type.sh 

less $1 | head -n 1000 | awk '{if(NR%4==0) printf("%s",$0);}' \
| od -A n -t u1 -v \
| awk 'BEGIN{min=100;max=0;} \
{for(i=1;i<=NF;i++) {if($i>max) max=$i; if($i<min) min=$i;}}END \
{if(max<=126 && min<59) print "Phred33"; \
else if(max>73 && min>=64) print "Phred64"; \
else if(min>=59 && min<64 && max>73) print "Solexa64"; \
else print "Unknown score encoding"; \

运行

sh fq_qual_type.sh example.fq