Fastq格式文件儲存了生物序列的信息及其質(zhì)量信息。以電腦中的一個文件的為例
$ 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標識符，后面可跟其他描述性內(nèi)容，但序列ID與描述部分空格分開，這點很重要，是很多分析的基礎(chǔ)
第二行：序列的順序，基本不可能出現(xiàn)重復(fù)。這里是核酸。
第三行：必須+開頭，后面可以跟描述信息，此處是跟了第一行的信息。但為了節(jié)約空間，很多時候這一行只有一個+,而不跟任何內(nèi)容。
第四行：堿基或氨基酸（此處堿基）的質(zhì)量字符，對應(yīng)著第二行的堿基，反應(yīng)的是該堿基的錯誤率，所以這一行的字符數(shù)和第二行要一一對應(yīng)，否則就亂了。這里就引入了ASCII code。
所以引入下面，堿基質(zhì)量值是什么，如何獲得？怎么表示？如何轉(zhuǎn)換？
2 堿基質(zhì)量值Q值和ASCII碼之間的關(guān)系
因為第四行的編碼，開始由Phred程序開發(fā)者定義，所以一般稱為Phred quality。那堿基質(zhì)量得分怎么來的？
每合成一個堿基，即可發(fā)出一個熒光信號，該信號可以被捕捉到，并生成是是軌跡數(shù)據(jù)。不同的堿基用不同顏色標記，檢測相應(yīng)峰值即可判斷堿基。
而Phred通過計算相應(yīng)波峰參數(shù)，去查詢通過已知序列測序分析得到的一個表，即可把錯誤率轉(zhuǎn)換為質(zhì)量得分。也就是把波峰參數(shù)和質(zhì)量得分對應(yīng)起來。
堿基錯誤率與質(zhì)量得分的關(guān)系如下
也就是說，質(zhì)量值Q是測序錯誤率的對數(shù)*-10。假如錯誤率是0.01，則Q值為20。可見，錯誤率越低，其Q值越高。即Q值越高越可靠。
Q10——0.1
Q20——0.01
Q30——0.001
Q40——0.0001
那就直接以上面這種Q值表示堿基質(zhì)量不就行了嗎？為什么要用ASCII碼表示。
如果用數(shù)字表示，數(shù)字和數(shù)字之間需要有間隔符號以區(qū)分，再者會浪費存儲空間，所以可以把質(zhì)量值轉(zhuǎn)變?yōu)橄鄳?yīng)的ASCII碼，這樣就完成了把質(zhì)量數(shù)向ASCII碼的轉(zhuǎn)換，那現(xiàn)在看下ASCII碼
如果直接把Q值直接對應(yīng)ASCII碼，應(yīng)該挺方便的，但是，Q值有時會有負值，再者，看ASCII碼的0-31位都是控制字符，沒法打印和保存，能打印的從要從32位的
Space開始，所以就可以給實際的Q值加上一個固定值，這樣就可以打印出來而保存在fastq文件中了。所以下面就是固定值加多少的問題。phred33，就是加了33，也就是10變成43，查表是
+。例子中很多F,Q值是70-33=37。而當時的solexa加的是64，這就是Phred64.他們與ASCII的對應(yīng)關(guān)系如下
數(shù)據(jù)處理時，有些軟件會根據(jù)堿基質(zhì)量得分的不同做不同處理，所以需要指定正確的編碼方式，也就是需要對Phred33還是64正確判斷。不過現(xiàn)在基本都33了，但如果下載以前的數(shù)據(jù)不一定。
下面是不同版本質(zhì)量得分和質(zhì)量字符ASCII的關(guān)系
從上面可以看出，Phred33的字符使用33-73，而+64使用包括59（包括）-104之間的ASCII碼。所以，只要ASCII小于59的僅僅在Phred+33中使用，而+64的都大于59，而從表可以看到大于等于74的旨在Phred+74中使用，所以，如果軟件沒有自動判斷，根據(jù)以上特點，就可以自行判斷是Phred33還是Phred64。
3 如何判斷是Phred33還是Phred64
默認讀取1000條序列，在這1000條序列中：
- 如果有2個以上的質(zhì)量字符ASCII值小于等于58(即有兩個堿基的得分小于等于25），同時沒有任何質(zhì)量字符的ASCII值大于等于75，即判斷是Phred+33。
- 如果有2個以上的質(zhì)量字符ASCII值大于等于75(即有兩個堿基的得分大于等于10），同時沒有任何質(zhì)量字符的ASCII值小于等于58，即判斷是Phred+64。
- 如果所有質(zhì)量字符的ASCII值介于59到74之間，即判斷可能是Phred+33，但建議使用更多的序列做進一步測試（出現(xiàn)這種結(jié)果可能有兩種情況：1, Phred+33編碼，所有堿基質(zhì)量得分介于26到42之間；2，Phred+64編碼，所有堿基質(zhì)量得分介于-5到10；是前者的可能性大）。
- 如果出現(xiàn)上述3種以外的情況，建議打印出質(zhì)量字符的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
