fastq ：怎麼判斷fastq是Phred33格式還是Phred64 格式

現在的一般都是Phred33的吧。。

對于每個堿基的品質編碼标示，不同的軟體采用不同的方案，目前有5種方案： Sanger，Phred quality score，值的範圍從0到92，對應的ASCII碼從33到126，但是對于測序資料（raw read data）品質得分通常小于60，序列拼接或者mapping可能用到更大的分數。 Solexa/Illumina 1.0, Solexa/Illumina quality score，值的範圍從-5到63，對應的ASCII碼從59到126，對于測序資料，得分一般在-5到40之間； Illumina 1.3+，Phred quality score，值的範圍從0到62對應的ASCII碼從64到126，低于測序資料，得分在0到40之間； Illumina 1.5+，Phred quality score，但是0到2作為另外的标示，詳見http://solexaqa.sourceforge.net/questions.htm#illumina Illumina 1.8+

官方學習連結：

http://maq.sourceforge.net/qual.shtml   <quality value>

http://maq.sourceforge.net/fastq.shtml   <fastq format>

以下屬于轉載,僅僅作為了解：

最近在學習質控知識時, 對于品質值體系及轉換産生了一些疑問, 作了一些嘗試, 趁叢集故障, 在此總結一下品質值體系相比之前教育訓練時所學的質控内容, (我拿到的) 流程中還多了一步 phred33to64, 也就是把 .fastq 格式的資料從 Phred33 品質值體系轉換為 Phred64 品質體系, 于是先補充學習了下品質值體系:首先要從品質值說起, 測序儀器下機資料的 fastq 檔案中, 每條序列都對應了相同長度的品質值, 反映出每個堿基的準确性和可靠性, (現在主流用的) 計算公式為:Q = -10log10p而這個 p 值就是 Phred 計算出來的, 表示一個堿基被識别錯誤的可能性, Phred 一開始是一個軟體 (或者說計算方法), 對測序儀器識别到的熒光強度 (三代的不了解) 進行評估, 針對不同儀器有不同的标準表, 然後根據表中熒光強度的範圍和分辨率分析得出堿基的 p 值, 由于比較可靠, 逐漸被各大公司采納 (以上腦補翻譯自 Phred quality score 和 Phred base calling)反正, Q 值為 10 就表示這個堿基有 90% 的機率是正确的, 20 就是 99%, 40 就是 4 個 9, 很好記, 相信大家也都很熟悉但是這時候問題就來了, 因為一個堿基對應一個品質值, 可 Q 值可以是一位數也可以是兩位數, 連在一起的話就分不出哪個對應哪個堿基 (比如某兩個堿基的品質值序列為 123 , 則可能為 1|23 或 12 | 3, 當然這是個極端的例子), 此外也浪費存儲空間, 因為一個數字隻有 10 種可能性, 卻要占去一個字元位, 而一個位元組有 8 位, 理論上已經可以代表 256 種狀态, 這還沒換算一個字元要占多少位元組, 是以就會把堿基品質值轉換為相應的 ASCII 碼, 這是計算機中最基本的一套字元體系, 用來把常用符号 (共 127 個) 轉為二進制以便于機器使用, ASCII 碼表見 ASCII code chart, 這樣一個字元就可以搞定了, 很友善很省事但是這時候問題又來了, 最理想的情況當然是直接把品質值作為序号找出對應的那個 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 就行了, 這個代碼有點老, 還沒加 [35, 75] 的判斷, 我又懶得重新改判斷邏輯, 是以直接在最後加了個列印最小值最大值, 再跟上面的示例一比, 就基本确定是什麼體系了關于這兩種體系對我們實際流程的影響和分析, 請聽下回分解轉自：http://not.farbox.com/post/phred_p1