FASTQ文件格式及测序文件phred质量格式判断

1. fastq格式说明

Fastq格式是一种基于文本的存储生物序列和对应碱基（或氨基酸）质量的文件格式。下面是一个Illumina平台测序的真实数据，其中包含了1条reads的信息。FASTQ格式储存的序列信息，每1条reads的信息，可以分成4行。

cat WT_rep1_BAF155.fastq | head -n 4
第1行主要储存序列测序时的坐标等信息,@ 开始的标记符号,SALLY:291:C149WACXX测序仪唯一的设备名称,2表示 lane的编号，表示1101 tail的坐标，2579 在tail中的X坐标，1951 在tail中的Y坐标,length=51表示这段read长度
@WT_rep1_BAF155.1 SALLY:291:C149WACXX:2:1101:2579:1951 length=51
第2行是测序得到的序列信息，一般用ATCGN来表示，其中N表示荧光信号干扰无法判断到底是哪个碱基。
CTGNCCAAGGTAATTTATAGATTCAATGCCATCCCCATCAAGCTACCAANG
第3行以“+”开始，可以储存一些附加信息，一般是空的
+WT_rep1_BAF155.1 SALLY:291:C149WACXX:2:1101:2579:1951 length=51
第4行储存的是质量信息，与第2行的碱基序列是一一对应的，其中的每一个符号对应的ASCII值成为phred值，可以简单理解为对应位置碱基的质量值，越大说明测序的质量越好。不同的版本对应的不同。
BCC#4ADDHHBFHIJJIIJJJIIIIJHIJIJIIJGGIJJJJIGJJJJJJ##

8fa1aecb80692b73a052f54235134976.jpg

2.FASTQ质量值的计算方法

在测序仪进行测序的时候，会自动根据荧光信号的强弱给出一个参考的测序错误概率（error probility，P）根据定义来说，P值肯定是越小越好。我们怎么储存他们呢？直接储存成小数点？比如1%储存成0.01？这肯定是不高效的，因为1个碱基的信息，占用了至少4个字符。
所以科学家们的做法想了一个办法：
1.将P取log10之后再乘以-10，得到的结果为Q。
比如，P=1%，那么对应的Q=-10*log10（0.01）=20
2.把这个Q加上33或者64转成一个新的数值，称为Phred，最后把Phred对应的ASCII字符对应到这个碱基。
如Q=20，Phred = 20 + 33 = 53，对应的符号是”5”
各版本不同Phred对应的ASCII值
在计算Q值和加上33/64的时候，不同测序仪，产生的数据不同，大概如下所示：
Solexa标准

image.png

Illumina标准

image.png

3. ASCII码与碱基质量Q值表

3e1c891fc2a933f2909237b792b50150.jpg

4.不同的质量编码格式

下面是不同版本质量得分和质量字符ASCII的关系：

be1f76ff9262475642a478d2eee305f9.jpg

phred33=-10log10（p-value）+33
phred64=-10log10（p-value）+64
根据上面的公式可以看到，这两种质量记录方式除了计算的时候加上的数字不一样，其他都是一样的。
数据处理时，有些软件会根据碱基质量得分的不同做不同处理，所以需要指定正确的编码方式。

5.判断标准

根据上面的图，我们可以总结判断标准：
ASCII值在 (33,73)，是 Sanger 的 Phred+33
ASCII值在 (33,74)，是 Illumina 1.8 的 Phred+33
ASCII值在 (59,104)，是 Solexa 的 Phred+64
ASCII值在 (64,104)，是 Illumina 1.3 的 Phred+64
ASCII值在 (67,104)，是 Illumina 1.5 的 Phred+64
进一步归纳规律：
ASCII值<59：只出现在Phred+33
ASCII值>74：只出现在Phred+64

6.phred质量格式判断

有了上面的背景知识，其实想要区分Phred+33和Phred+64就非常简单了。

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 WT_rep1_BAF155.fastq.gz
Phred33


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