最近在學習質控知識時, 對於質量值體系及轉換產生了一些疑問, 作了一些嘗試, 趁叢集故障, 在此總結一下

質量值體系

相比之前培訓時所學的質控內容, (我拿到的) 流程中還多了一步 phred33to64, 也就是把 .fastq 格式的資料從 Phred33 質量值體系轉換為 Phred64 質量體系, 於是先補充學習了下質量值體系:

首先要從質量值說起, 測序儀器下機資料的 fastq 檔案中, 每條序列都對應了相同長度的質量值, 反映出每個鹼基的準確性和可靠性, (現在主流用的) 計算公式為:

Q = -10log₁₀p

而這個 p 值就是 Phred 計算出來的, 表示一個鹼基被識別錯誤的可能性, Phred 一開始是一個軟體 (或者說計算方法), 對測序儀器識別到的熒光強度 (三代的不瞭解) 進行評估, 針對不同儀器有不同的標準表, 然後根據表中熒光強度的範圍和解析度分析得出鹼基的 p 值, 由於比較可靠, 逐漸被各大公司採納 (以上腦補翻譯自 

反正, Q 值為 10 就表示這個鹼基有 90% 的概率是正確的, 20 就是 99%, 40 就是 4 個 9, 很好記, 相信大家也都很熟悉

但是這時候問題就來了, 因為一個鹼基對應一個質量值, 可 Q 值可以是一位數也可以是兩位數, 連在一起的話就分不出哪個對應哪個鹼基 (比如某兩個鹼基的質量值序列為 123 , 則可能為 1|23 或 12 | 3, 當然這是個極端的例子), 此外也浪費儲存空間, 因為一個數字只有 10 種可能性, 卻要佔去一個字元位, 而一個位元組有 8 位, 理論上已經可以代表 256 種狀態, 這還沒換算一個字元要佔多少位元組, 因此就會把鹼基質量值轉換為相應的 ASCII 碼, 這是計算機中最基本的一套字元體系, 用來把常用符號 (共 127 個) 轉為二進位制以便於機器使用, ASCII 碼錶見 

但是這時候問題又來了, 最理想的情況當然是直接把質量值作為序號找出對應的那個 ASCII 碼, 比如質量值為 40 就換成十進位制 40 對應的 ASCII 碼,可惜質量值根據測序儀公司的標準不同, 範圍也各不相同, 基本都包括了 0 至 40 的區間, 甚至還可能是負值, 這就沒法愉快地玩耍了, 而且人家 ASCII 碼錶也不配合, 0 到 31 對應的都是些控制字元 (比如回車, 退格), 根本不適合列印和儲存, 可列印的都得從 32 號排起 (參見 ASCII printable code chart), 所以各家測序儀器公司就把質量值再加上某個固定值作為 ASCII 碼轉換成了可列印字元從而儲存在 FASTQ 檔案中

可是問題還沒有完, 這個固定值是多少好呢, 各家公司是競爭對手, 怎麼可能你用什麼我也用什麼, 所以 Sanger 公司加了 33, 也就是質量值為 0 就轉換成 ASCII 碼 33, 查表可知為 !, 也即從可列印的字元開始 (排除了空格), 這就是現在所謂的 Phred33 體系, 當時的 Solexa (後來被 Illumina 收購) 公司就偏不用 33 (此處為個人腦補), 偏要加個 64, 這樣質量值為 0 就用 @ 表示, 後面從 1 開始的就依次對應了 ABCD, 於是就成了 Phred64 體系, 至於當時三巨頭的另一家測序儀公司 454 Life Sciences (後被 Roche 收購) 就更絕, 人家從鹼基開始就不用 ACTG 表示, 直接整了個 ColorSpace 體系出來, 根本不和你們玩, (話說 Color Space 這玩意兒曾經把我狠狠地坑了好久), 當然後來大家也不跟 454 玩了, 最後他也就沒得玩了

回到質量值體系, 這樣就由 Sanger 公司和 Illumina 公司產生了 Phred33 體系和 Phred64 體系, 兩家互相拗著, 這就辛苦了寫生物資訊分析軟體的人, 兩種質量體系都要考慮, 當然好一點的軟體都是有引數接受體系型別的, 更好一點的軟體就會自動判斷體系型別進行對應轉換

本來就這樣結束的話也算是個圓滿的故事了, 可是墨菲他老人家不高興了 (Murphy's law), 所以在 2011 年, Illumina 公司表示他們又要改成 Phred33 體系了 (Upcoming changes in CASAVA), 真是大(wo)快(le)人(ge)心(qu)啊! 這麼來回一折騰, 結果還是回到了最初的起點, 與老基友相擁而泣了, Phred33 一統江湖! (當然實際上現在 Ilumina 的 Phred33 和最初 Sanger 的 Phred33 還是有點區別的, 詳見後文)

扯了這麼多, 發現寫了半天才剛交代完故事背景, 剩下的部分就等有空再寫了, 最後上點乾貨:

先是 wikipedia 上非常直觀的示例, 可以看出各家公司各個版本的質量值體系

  SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS.....................................................
  ..........................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX......................
  ...............................IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII......................
  .................................JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ......................
  LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL....................................................
  !"#$%&'()*+,-./0123456789:;<=>[email protected][\]^_`abcdefghijklmnopqrstuvwxyz{|}~
  |                         |    |        |                              |                     |
 33                        59   64       73                            104                   126
  0........................26...31.......40                                
                           -5....0........9.............................40 
                                 0........9.............................40 
                                    3.....9.............................40 
  0.2......................26...31........41                              

 S - Sanger        Phred+33,  raw reads typically (0, 40)
 X - Solexa        Solexa+64, raw reads typically (-5, 40)
 I - Illumina 1.3+ Phred+64,  raw reads typically (0, 40)
 J - Illumina 1.5+ Phred+64,  raw reads typically (3, 40)
     with 0=unused, 1=unused, 2=Read Segment Quality Control Indicator (bold) 
     (Note: See discussion above).
 L - Illumina 1.8+ Phred+33,  raw reads typically (0, 41)

然後光看也沒什麼用啊, 就有人寫了個小函式, 用來分析 fastq 檔案 (壓縮或未壓縮均可), 程式碼如下:

fqtype () {
        less $1 | head -n 999 | 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<=74 && min<59) print "Phred+33"; \
        else if(max>73 && min>=64) print "Phred+64"; \
        else if(min>=59 && min<64 && max>73) print "Solexa+64"; \
        else print "Unknown score encoding"; \
        print "( " min ", " max, ")";}'
}

寫在 shell 配置檔案裡, source 然後 fqtype <fastq_file> 就行了, 這個程式碼有點老, 還沒加 [35, 75] 的判斷, 我又懶得重新改判斷邏輯, 所以直接在最後加了個列印最小值最大值, 再跟上面的示例一比, 就基本確定是什麼體系了

關於這兩種體系對我們實際流程的影響和分析, 請聽下回分解

感覺寫的很好理解，特此謝謝。