測序數據的解析:Fastq與FastQC

二代測序平臺獲得的原始數據為fastq（或為壓縮文件fq.gz）格式，包含雙末端測序所得的正向和反向兩個文件（通常用「1」和「2」來區分），如下所示：

每一個read包含四行內容，其中第一行以@開頭，後面是reads的屬性信息，也即read名稱。中間用「:」隔開。例如上例中HISEQ為測序平臺名稱，266為測序運行run的編號，HHNWKBCXX為流通池（flowcell）編號。接下來四個數字為位置信息，2代表流通池中的第2個lane，1101代表第2個lane中的第1101個tile，10010:58789代表該read在該tile中的x：y坐標信息。1為讀取編號，雙末端一共有兩次讀取（不包含index的讀取）。N代表沒有被儀器過濾掉，也即通過了初步質量過濾；0為controlnumber，代表對照序列的鑑定情況。GGCTAC為文庫Index序列。

第二行為read序列信息。一般條件下read1裡面最前面為特異性Barcode和反向引物的序列，read2裡面最前面為正向引物的序列。

第三行以「+」開頭，跟隨者該read的名稱（一般與@後面的內容相同），可以省略，但「+」一定不能省。

第四行代表read每個鹼基的測序質量。每個鹼基對應的字符在ASCII碼中對應的十進位數字減去33即為該鹼基質量（也即Phred33體系），例如上圖中第一個鹼基的質量為D，對應的十進位數字為68（見下表），則鹼基質量為68-33=35。鹼基質量Q=-10*lgP，P為鹼基被測錯的概率。也即Q為30代表被測錯的概率為0.001，鹼基質量越高，則被測錯的概率越低。

對於新下機的原始數據，我們可以使用軟體FastQC來對其測序質量進行可視化，軟體地址：

http://www.bioinformatics.babraham.ac.uk/projects/fastqc/

FastQC使用方法如下：

其中參數-o設置結果文件輸出路徑，默認為當前路徑；-q為安靜模式運行；-t設置所使用的核數，根據伺服器情況而定，更多命令行選項使用命令fastqc -h來查看。fastqfile為原始測序數據，也可以是fq.gz壓縮文件：

#可以同時檢查正反向原始數據：fastqc -o fastqc -t 20 R1.fastq R2.fastq#對於大批量的數據，也可以用過管道命令和shell腳本進行批量處理：ls rawdata/*fq | while read id; do fastqc -o fastqc -q -t 20 $id; done#或者後臺運行模式（一次提交多任務，所以核數要調小）：ls rawdata/*fq | while read id; do nohup fastqc -o fastqc -q -t 2 $id & done#上面的while循環是很常用的一種文檔批處理方法，注意變量的運算使用的是通配符而非正則表達式。對引用標準輸入還可使用xargs函數：ls rawdata/*fq | xargs -n 1 -P 5 fastqc -o fastqc -q -t 10#有時候一個項目有大批量樣品甚至大批量文庫，需要合併來檢測質量並做報告，這時候可以使用以下命令合併序列文件：cat *1.fq > total.R1.fqcat *2.fq > total.R2.fq
打開生成的html結果報告文件，就可以看到可視化的質檢結果。在查看結果之前，我們要對自己的數據有一定的把握，例如是否已經去掉接頭，是擴增子測序數據還是鳥槍法測序數據等。基因組宏基因組鳥槍法測序數據reads比較隨機均勻，鹼基分布也會比較均勻，而擴增子數據由於兩端都有引物，以及插入片段均為16S，所以會出現很多重複序列，且鹼基分布非均勻。接下來詳細解析不同檢測結果的含義。
使用瀏覽器打開html文件後，即可看到基本的參數以及質量分析結果，結果分為綠色的"PASS"，黃色的"WARN"和紅色的"FAIL"，如下所示：
包括測序技術/平臺、reads數量、長度、GC含量等。
⑵Per base sequence quality
所有reads鹼基的測序質量統計結果。箱線圖中紅色線表示中位數，黃色是25%-75%區間，延伸線是10%-90%區間，藍線是平均數曲線。若任一位置鹼基的下四分位數低於10或中位數低於25，報"WARN"；若任一位置的下四分位數低於5或中位數低於20，報"FAIL"。
⑶Pertile sequence quality流通池中不同晶片（tile）的鹼基測序質量平均值對比，顯示了測序儀的系統差錯。熱圖中藍色部分是質量較好的點，紅色越明顯則是測序質量越低。縱坐標為tile編號，如果某tile的測序質量很低，可以考慮去除該tile的序列數據。
⑷Per sequence quality scores每條read的鹼基質量均值的統計結果。橫軸為測序質量quality，縱軸是read數目。從圖中可以容易得看出不同質量範圍內的read數量。其中當峰值也即最大read質量小於27（錯誤率0.2%）時報"WARN"，當峰值小於20（錯誤率1%）時報"FAIL"。
⑸Per Base Sequence Content
對所有reads的每一個位置，統計ATCG四種鹼基（正常情況）比例的分布情況。橫軸為鹼基位置，縱軸為百分比。正常情況下四種鹼基的出現頻率應該是接近的，而且沒有位置差異。因此好的樣本中四條線應該平行且接近。當部分位置鹼基的比例出現bias時，即四條線在某些位置紛亂交織，往往提示我們有overrepresentedsequence的汙染。當所有位置的鹼基比例一致的表現出bias時，即四條線平行但分開，往往代表文庫有bias(建庫過程或本身特點)，或者是測序中的系統誤差。當任一位置的A/T比例與G/C比例相差超過10%，報"WARN"；當任一位置的A/T比例與G/C比例相差超過20%，報"FAIL"。
統計reads的平均GC含量的分布。橫軸為GC比例，縱軸為reads數量。紅線是實際情況，藍線是理論分布（正態分布，均值不一定在50%，而是由平均GC含量推斷的）。曲線形狀的偏差往往是由於文庫的汙染或是部分reads構成的子集有偏差（overrepresentedreads）。形狀接近正態但偏離理論分布的情況提示我們可能有系統偏差。偏離理論分布的reads超過15%時，報"WARN"；偏離理論分布的reads超過30%時，報"FAIL"。
當測序儀器不能辨別某條read的某個位置到底是什麼鹼基時，就會產生「N」。對所有reads的每個位置，統計N的比例。
正常情況下N的比例是很小的，所以圖上常常看到一條直線，但放大Y軸之後會發現還是有N的存在，這不算問題。當Y軸在0%-100%的範圍內也能看到「凸起」時，說明測序系統出了問題。當任意位置的N的比例超過5%，報"WARN"；當任意位置的N的比例超過20%，報"FAIL"。⑻Sequence Length Distribution統計reads長度的分布。橫軸為片段長度，縱軸為read數量。當reads長度不一致時報"WARN"；當有長度為0的read時報"FAIL"。
⑼Sequence Duplication Levels統計序列的重複度（duplication level，也即一個文庫中某條序列的copy數），理論上大部分序列都只出現一次，低的重複度意味著高的基因組覆蓋率。測序深度越高，越容易產生一定程度的重複，這是正常的現象；但如果duplication的程度很高，就提示我們可能有bias的存在（如建庫過程中的PCR duplication）。可以想像，如果原始數據很大（事實往往如此），對所有序列的比較將會需要很大內存，所以FastQC只用前100,000條reads來進行統計，以反映全部數據中序列重複度情況。而且，大於75bp的reads只取前50bp進行比較，由於reads越長越不容易完全相同（由測序錯誤導致），所以這樣做使得重複度的統計更加嚴格。序列duplication level分布圖將會展示文庫中不同重複度的序列所佔比例，其中橫坐標是duplication levels，縱坐標是duplicated reads的比例。圖中藍色線展示了全部序列中不同重複度序列的百分比，紅線顯示的是有重複序列中不同重複度序列的百分比（所有序列的重複度減去1）。由於展示範圍的限制，重複數目大於等於10的reads會被按照區間合併統計，造成在duplicationlevel為10的時候曲線突然凸起，結果如下所示：
當非unique（也即duplication level大於1）的reads佔總數的比例大於20%時，報"WARN"；當非unique的reads佔總數的比例大於50%時，報"FAIL"。⑽Overrepresented Sequences如果有某個序列大量出現，就叫做over-represented。FastQC的標準是佔全部reads的0.1%以上。和上面的duplicate analysis一樣，為了計算方便，只取了fastq數據的前100,000條reads進行統計，所以有可能over-represented reads不全在裡面。而且大於75bp的reads也是只取50bp。統計結果以列表形式展示，當發現超過總reads數0.1%的reads時報"WARN"，當發現超過總reads數1%的reads時報"FAIL"。統計接頭序列的含量。一般測序儀自帶軟體會切去接頭序列，所以下機數據並沒有接頭序列。
如果某n個鹼基的短序列在reads中大量出現，其頻率高於統計期望的話，FastQC將其記為over-representedkmer（重複短序列）。默認的n=5，可以通過設置-k參數的值來調節n的大小，範圍是2-10。出現頻率總體上3倍於期望或是在某位置上5倍於期望的k-mer被認為是over-represented。FastQC除了列出所有over-representedkmers，還會繪製前6個k-mers的分布圖。如下圖所示我們的數據中只檢測出一個k-mer序列：
如下所示為k-mers分布圖，其中橫坐標為k-mer出現的鹼基位點，縱坐標為該位點k-mers數目：
當有出現頻率總體上3倍於期望或是在某位置上5倍於期望的k-mer時，報"WARN"；當有出現頻率在某位置上10倍於期望的k-mer時報"FAIL"。