1. 原始測(cè)序數(shù)據(jù)的質(zhì)控
2. read比對(duì),，排序和去除重復(fù)序列
3. Indel區(qū)域重（“重新”的“重”）比對(duì)
4. 堿基質(zhì)量值重校正
5. 變異檢測(cè)
6. 變異結(jié)果質(zhì)控和過濾

       在這個(gè)圖中,，我把WGS數(shù)據(jù)分析流程的各個(gè)步驟和關(guān)系都畫下來了,。這個(gè)流程雖然只針對(duì)于人,，但對(duì)于其它二倍體生物來說,，同樣具有借鑒價(jià)值。這6個(gè)步驟,，接下來我也會(huì)進(jìn)行詳細(xì)介紹,，在本篇文章中我們首先介紹原始測(cè)序數(shù)據(jù)的質(zhì)控。

認(rèn)識(shí)測(cè)序數(shù)據(jù)——數(shù)據(jù)質(zhì)控的意義

       在第1節(jié)測(cè)序技術(shù)中,，我們已經(jīng)知道現(xiàn)在的NGS測(cè)序,，以illumina為首基本都是運(yùn)用邊合成邊測(cè)序的技術(shù)。堿基的合成依靠的是化學(xué)反應(yīng),，這使得堿基鏈可以不斷地從5’端一直往3’端合成并延伸下去,。但在這個(gè)合成的過程中隨著合成鏈的增長，DNA聚合酶的效率會(huì)不斷下降,，特異性也開始變差,，這就會(huì)帶來一個(gè)問題——越到后面堿基合成的錯(cuò)誤率就會(huì)越高^【注】，這也是為何當(dāng)前NGS測(cè)序讀長普遍偏短的一個(gè)原因,。

【注】：有時(shí)候測(cè)序儀在剛開始進(jìn)行合成反應(yīng)的時(shí)候也會(huì)由于反應(yīng)還不夠穩(wěn)定,，同樣會(huì)帶來質(zhì)量值的波動(dòng)，不過這個(gè)波動(dòng)一般都在高質(zhì)量值區(qū)域（如下圖）,。

       測(cè)序數(shù)據(jù)的質(zhì)量好壞會(huì)影響我們的下游分析,。但不同的測(cè)序平臺(tái)其測(cè)序錯(cuò)誤率的圖譜都是有差別的。因此,，非常建議在我們分析測(cè)序數(shù)據(jù)之前先搞清楚如下兩個(gè)地方：

• 原始數(shù)據(jù)是通過哪種測(cè)序平臺(tái)產(chǎn)生的,，它們的錯(cuò)誤率分布是怎么樣的，是否有一定的偏向性和局限性,，是否會(huì)顯著受GC含量的影響等,；
• 評(píng)估它們有可能影響哪些方面的分析；

       第一點(diǎn)是我們認(rèn)識(shí)數(shù)據(jù)質(zhì)量的第一步,，也是我們一定要去知道的地方,。除了看官方的資料之外,，最好的做法是自己分析。

       雖然隨著NGS測(cè)序數(shù)據(jù)變得越來越普遍,，整體的測(cè)序質(zhì)量和錯(cuò)誤率分布情況大家也都了解一些,。在實(shí)際的工作中，我也常常發(fā)現(xiàn)很多人其實(shí)并不十分關(guān)心這個(gè)數(shù)據(jù)到底長啥樣,，拿到之后,，就直接跑過濾流程，也不管這些參數(shù)或者工具是否真的是合適的,，更加不看看過濾后的數(shù)據(jù)和過濾前到底有什么不同,。盡管，大多數(shù)情況下問題不大,，但是我想跟大家說的是： 認(rèn)識(shí)你的數(shù)據(jù),，不要相信你的工具！ 這樣也能夠更好地避開很多不必要的坑,。

       因此,，在本文中我將談?wù)勗撊绾胃玫厝フJ(rèn)識(shí)一個(gè)測(cè)序數(shù)據(jù)，而不只是簡單地告訴大家一個(gè)質(zhì)控流程,，當(dāng)然,，這也是本篇文章將要進(jìn)行介紹的內(nèi)容。至于第二點(diǎn)其實(shí)需要視情況而定,，例如你的測(cè)序深度是多少,，檢測(cè)變異的時(shí)候，變異的位點(diǎn)是否過于集中在read的末尾,，比對(duì)的時(shí)候是否會(huì)出現(xiàn)了一定的正反鏈偏向性等諸如此類的問題,；又或者我們?cè)谶M(jìn)行基因組序列組裝的時(shí)候，由于對(duì)read的中出錯(cuò)的堿基更加敏感,，因此往往需要進(jìn)行更嚴(yán)格的切除,，不然會(huì)由于這些錯(cuò)誤的堿基消耗更大的計(jì)算資源和時(shí)間。

       那么說到底該如何認(rèn)識(shí)一個(gè)原始的測(cè)序數(shù)據(jù)（fastq data）呢,？一般我們可以從如下幾個(gè)方面來分析：

• read各個(gè)位置的堿基質(zhì)量值分布
• 堿基的總體質(zhì)量值分布
• read各個(gè)位置上堿基分布比例,，目的是為了分析堿基的分離程度
• GC含量分布
• read各位置的N含量
• read是否還包含測(cè)序的接頭序列
• read重復(fù)率，這個(gè)是實(shí)驗(yàn)的擴(kuò)增過程所引入的

       以上,，這幾個(gè)地方都弄明白了,，那么這個(gè)數(shù)據(jù)的基本情況也就差不多都清楚了。

       我們首先來說說read各位置的堿基質(zhì)量分布,。在第2節(jié)里面,，我們已經(jīng)知道該如何通過簡單的Python代碼計(jì)算出read的質(zhì)量值和堿基的測(cè)序錯(cuò)誤率了。但對(duì)于成千上萬的read來說，這樣做并不合適,，我們需要更直觀的表達(dá)方式——畫出來,，正所謂 一圖勝千言！目前也有很多現(xiàn)成的工具可以高效地來完成這樣的事情,，比如用得最廣的FastQC,，它是一個(gè)java程序，能夠用于給出測(cè)序數(shù)據(jù)的QC報(bào)告,，報(bào)告中會(huì)同時(shí)給出上述幾個(gè)方面的數(shù)據(jù)圖,，并提示原來的數(shù)據(jù)可能還存在著哪些問題。它可以很好地幫助我們理解測(cè)序數(shù)據(jù)的質(zhì)量情況,，但缺點(diǎn)就是 圖,！太！丑,！

      在做read質(zhì)量值分析的時(shí)候，F(xiàn)astQC并不單獨(dú)查看具體某一條read中堿基的質(zhì)量值,，而是將Fastq文件中所有的read數(shù)據(jù)都綜合起來一起分析,。下圖是一個(gè)測(cè)序質(zhì)量非常好的read各位置堿基質(zhì)量分布圖（如下圖）。

      這個(gè)圖的橫軸是read上堿基的位置,，縱軸是堿基質(zhì)量值,。在這個(gè)例子中，read的長度是126bp（來自HiSeq X10的測(cè)序結(jié)果）,，這應(yīng)該算是比較長的二代測(cè)序序列了,。我們可以看到read上的每一個(gè)位置都有一個(gè)黃色的箱型圖表示在該位置上所有堿基的質(zhì)量分布情況。除了最后一個(gè)堿基之外,，其他的堿基質(zhì)量值都基本都在大于30,，而且波動(dòng)很小，說明質(zhì)量很穩(wěn)定,，這其實(shí)是一個(gè)非常高質(zhì)量的結(jié)果,。而且我們可以看到圖中質(zhì)量值的分布都在綠色背景（代表高質(zhì)量）的區(qū)域。

       那如果是質(zhì)量很差的結(jié)果看起來會(huì)是怎么樣的呢,？我手邊一時(shí)找不到這樣的數(shù)據(jù),，就在網(wǎng)上找到了一個(gè)代替品，樣子如下：

      在這個(gè)圖中我們可以明顯看到,，read各個(gè)位置上的堿基質(zhì)量分布波動(dòng)都比較大,，特別從第18個(gè)堿基往后全部出現(xiàn)了大幅度的波動(dòng)，而且很多read的堿基質(zhì)量值都掉到非常低（紅色）的區(qū)域中了,，說明這個(gè)數(shù)據(jù)的測(cè)序結(jié)果真的非常差,，有著大量不及格的read。最好的情況是重新測(cè)序,，但如果不得不使用這個(gè)數(shù)據(jù),，就要把這些低質(zhì)量的數(shù)據(jù)全都去除掉才行,，同時(shí)還需留意是否還存在其他的問題，但不管如何都一定會(huì)丟掉很大一部分的數(shù)據(jù),。

       除了上面read各位置的堿基質(zhì)量值分布之外,，F(xiàn)astQC還會(huì)為我們計(jì)算其他幾個(gè)非常有價(jià)值的統(tǒng)計(jì)結(jié)果，包括：

1）堿基總體質(zhì)量值分布,，只要大部分都高于20,，那么就比較正常。

       在第2節(jié)里面我也提到了關(guān)于Q20和Q30的比例是我們衡量測(cè)序質(zhì)量的一個(gè)重要指標(biāo),。這其實(shí)也是從這里來進(jìn)行體現(xiàn)的,，一般來說，對(duì)于二代測(cè)序,，最好是達(dá)到Q20的堿基要在95%以上（最差不低于90%）,，Q30要求大于85%（最差也不要低于80%）。

2）read各個(gè)位置上堿基比例分布

這個(gè)是為了分析堿基的分離程度,。何為堿基分離,？我們知道AT配對(duì)，CG配對(duì),，假如測(cè)序過程是比較隨機(jī)的話（隨機(jī)意味著好）,，那么在每個(gè)位置上A和T比例應(yīng)該差不多，C和G的比例也應(yīng)該差不多,，如上圖所示,，兩者之間即使有偏差也不應(yīng)該太大，最好平均在1%以內(nèi),，如果過高,，除非有合理的原因，比如某些特定的捕獲測(cè)序所致,，否則都需要注意是不是測(cè)序過程有什么偏差,。

3）GC含量分布圖

       GC含量指的是G和C這兩種堿基占總堿基的比例。二代測(cè)序平臺(tái)或多或少都存在一定的測(cè)序偏向性,，我們可以通過查看這個(gè)值來協(xié)助判斷測(cè)序過程是否足夠隨機(jī),。對(duì)于人類來說，我們基因組的GC含量一般在40%左右,。因此,，如果發(fā)現(xiàn)GC含量的圖譜明顯偏離這個(gè)值那么說明測(cè)序過程存在較高的序列偏向性，結(jié)果就是基因組中某些特定區(qū)域被反復(fù)測(cè)序的幾率高于平均水平,，除了覆蓋度會(huì)有偏離之后,，將會(huì)影響下游的變異檢測(cè)和CNV分析。

4）N含量分布圖

       N在測(cè)序數(shù)據(jù)中一般是不應(yīng)該出現(xiàn)的，如果出現(xiàn)則意味著,，測(cè)序的光學(xué)信號(hào)無法被清晰分辨,，如果這種情況多的話，往往意味著測(cè)序系統(tǒng)或者測(cè)序試劑的錯(cuò)誤,。

5）接頭序列

       在第1節(jié) 測(cè)序技術(shù)里面我們提到了在測(cè)序之前需要構(gòu)建測(cè)序文庫,，測(cè)序接頭就是在這個(gè)時(shí)候加上的，其目的一方面是為了能夠結(jié)合到flowcell上,，另一方面是當(dāng)有多個(gè)樣本同時(shí)測(cè)序的時(shí)候能夠利用接頭信息進(jìn)行區(qū)分,。當(dāng)測(cè)序read的長度大于被測(cè)序的DNA片段^【注】時(shí)，就會(huì)在 read的末尾測(cè)到這些接頭序列（如下圖）,。一般的WGS測(cè)序是不會(huì)測(cè)到這些接頭序列的,，因?yàn)闃?gòu)建WGS測(cè)序的文庫序列（插入片段）都比較長，約幾百bp,，而read的測(cè)序長度都在100bp-150bp這個(gè)范圍,。不過在進(jìn)行一些RNA測(cè)序的時(shí)候，由于它們的序列本來就比較短,，很多只有幾十bp長（特別是miRNA）,，那么就很容易會(huì)出現(xiàn)read測(cè)通的現(xiàn)象，這個(gè)時(shí)候就會(huì)在read的末尾測(cè)到這些接頭序列,。

【注】這些DNA片段也常被我們稱之為“插入片段”

       最后，這些被測(cè)到的接頭序列和低質(zhì)量堿基一樣都是需要在正式分析之前進(jìn)行切除的read片段,。

       當(dāng)我們看完了上面的這些結(jié)果之后就可以比較清楚地了解一個(gè)測(cè)序數(shù)據(jù)的概況了,。

      那么，說了這么多,，上述提到的FastQC該怎么用呢,？

      FastQC的安裝非常簡單，我們可以通過網(wǎng)頁搜索或者直接到它的主頁上下載最新的版本,。

       FastQC的運(yùn)行非常簡單,，直接在終端通過命令行是最有效直接的,，下面我給出一個(gè)例子：

       命令比較簡單，這里 唯一值得注意的地方就是 -o 參數(shù)用于指定FastQC報(bào)告的輸出目錄,，這個(gè)目錄需要事先創(chuàng)建好,，如果不指定特定的目錄，那么FastQC的結(jié)果會(huì)默認(rèn)輸出到文件untreated.fq的同一個(gè)目錄下,。它輸出結(jié)果只有兩個(gè),，一個(gè)html和一個(gè).zip壓縮包。

       我們可以直接通過瀏覽器打開html，就可以看到FastQC給出的所有結(jié)果,，zip壓縮包解壓后,，從中我們也可以在對(duì)應(yīng)的目錄下找到所有的QC圖表和Summary數(shù)據(jù)。

       除了上述用法之外,，F(xiàn)astQC支持同時(shí)輸入多個(gè)fq文件（或者以通配符的形式輸入fq）,，當(dāng)我們的fq文件比較多時(shí)，這種用法會(huì)比較方便,，如：

這個(gè)鏈接是FastQC官網(wǎng)給出的一個(gè)Online報(bào)告的模板,，方便參考。

切除測(cè)序接頭序列和read的低質(zhì)量序列

       前面關(guān)于如何認(rèn)識(shí)fq數(shù)據(jù)的事情已經(jīng)說完了,，接下來是我們本篇文章中最后的一個(gè)重點(diǎn)——去除測(cè)序接頭和低質(zhì)量序列,！

       當(dāng)我們理解了fq數(shù)據(jù)之后，做這些過濾就不會(huì)很難,，你也完全可以自己編寫工具來進(jìn)行個(gè)性化的過濾,。目前也已有很多工具用來切除接頭序列和低質(zhì)量堿基，比如SOAPnuke,、cutadapt,、untrimmed等不下十個(gè)，但這其中比較方便好用的是Trimmomatic（也是一個(gè)java程序）,、sickle和seqtk,。Trimmomatic的好處在于，它不但可以用來切除illumina測(cè)序平臺(tái)的接頭序列,，還可以去除由我們自己指定的特定接頭序列,，而且同時(shí)也能夠過濾read末尾的低質(zhì)量序列，sickle和seqtk只能去除低質(zhì)量堿基,。具體的原理就是通過滑動(dòng)一定長度的窗口,，計(jì)算窗口內(nèi)的堿基平均質(zhì)量，如果過低,，就直接 往后全部切除,，注！意,！不是挖掉read中的這部分低質(zhì)量序列,，而是像切菜一樣，直接從低質(zhì)量區(qū)域開始把這條read后面的所有其它堿基全,！部,！剁！掉,！否則就是在人為改變實(shí)際的基因組序列情況,。

       如果下機(jī)的fq數(shù)據(jù)中不含有這些測(cè)序接頭,，那么我們除了trimmomatic之外，也可以直接使用sickle（同時(shí)支持PE和SE數(shù)據(jù)）或者seqtk（僅支持SE）,，這兩個(gè)處理起來會(huì)更快,，消耗的計(jì)算資源也更少。

現(xiàn)在我們說回如何用Trimmomatic構(gòu)造序列過濾流程,。

      首先是安裝Trimmomatic,。我們可以到它的官網(wǎng)上獲取最新的版本，下載打包好的binary即可,，如果打算看它具體的代碼,，可以在github上找到。

下載后,，直接解壓,，目錄下的trimmomatic-*.jar（我下載的是0.36版本）就是執(zhí)行程序，可以直接使用java來運(yùn)行,。

同個(gè)目錄下還有一個(gè)名為adapters的文件夾,，這個(gè)文件夾中的內(nèi)容對(duì)于我們?nèi)コ宇^序列來說非常重要。其中默認(rèn)存放的是illumina測(cè)序平臺(tái)的接頭序列（fasta格式）,，在實(shí)際的使用過程中,，如果需要去除接頭，我們需要明確指定對(duì)應(yīng)的序列作為輸入?yún)?shù),。

那么這些接頭序列具體該如何選擇呢,？一般來說，目前的HiSeq系列和MiSeq系列用的都是TruSeq3,，TruSeq2是以前GA2系列的測(cè)序儀所用的,，已經(jīng)很少見了。這些信息都可以在illumina的官網(wǎng)上找到,，至于具體該用PE（Pair End）還是SE（Single End）就按照具體的測(cè)序類型進(jìn)行選擇就ok了。如果用的不是illumina測(cè)序平臺(tái),，那么我們也可以按照adapters文件夾下的這些文件的格式做一個(gè)新的接頭序列,，然后再作為參數(shù)傳入。不過在自定義接頭序列的時(shí)候,，命名時(shí)有一些小的細(xì)節(jié)需要注意,，可以參考Trimmomatic的主頁文檔（The Adapter Fasta），這里就不展開了,。

       Trimmomatic有兩種運(yùn)行模式：PE和SE,。顧名思義，PE就是對(duì)應(yīng)Pair End測(cè)序的,，SE則是對(duì)應(yīng)Single End測(cè)序的,。

下面我分別給出例子來進(jìn)行說明：

PE模式,，HiSeq PE測(cè)序：

SE模式，HiSeq SE測(cè)序：

我們可以看到PE和SE,，顧名思義,，分別代表了 'PE模式’和'SE’模式。

同時(shí)需要明確指明質(zhì)量值體系是Phred33還是Phred64,，默認(rèn)是Phred64,，這需要特別注意，因?yàn)槲覀儸F(xiàn)在的測(cè)序數(shù)據(jù)基本都是Phred33的了,，所以一定要指定這個(gè)參數(shù),。

剩下的就是輸入的fq和輸出的fq，可以用-basein和-baseout指定,，也可以不用（如上例子）,，以及被過濾掉的fq要輸出到文件。細(xì)心的讀者可能已經(jīng)發(fā)現(xiàn)這里PE和SE有一個(gè)區(qū)別：

在SE模式中,，是不需要指定文件來存放被過濾掉的read信息的,，后面直接就接Trimmer信息！這是需要注意到的一個(gè)地方,。

關(guān)于后面的Trimmer信息,，規(guī)定了很多切除接頭序列和低質(zhì)量序列的細(xì)節(jié)，我挑重點(diǎn)的說,，具體如下：

• ILLUMINACLIP,，接頭序列切除參數(shù)。LLUMINACLIP:TruSeq3-PE.fa:2:30:10（省掉了路徑）意思分別是：TruSeq3-PE.fa是接頭序列,，2是比對(duì)時(shí)接頭序列時(shí)所允許的最大錯(cuò)配數(shù),；30指的是要求PE的兩條read同時(shí)和PE的adapter序列比對(duì)，匹配度加起來超30%,，那么就認(rèn)為這對(duì)PE的read含有adapter,，并在對(duì)應(yīng)的位置需要進(jìn)行切除^【注】 。10和前面的30不同,，它指的是,，我就什么也不管，反正只要這條read的某部分和adpater序列有超過10%的匹配率,，那么就代表含有adapter了,，需要進(jìn)行去除；

【注】測(cè)序的時(shí)候一般只會(huì)測(cè)到一部分的adapter,，因此read和adaper對(duì)比的時(shí)候肯定是不需要要求百分百匹配率的,，上述30%和10%其實(shí)是比較推薦的值。
SLIDINGWINDOW,，滑動(dòng)窗口長度的參數(shù),，SLIDINGWINDOW:5:20代表窗口長度為5,，窗口中的平均質(zhì)量值至少為20，否則會(huì)開始切除,；

• LEADING,，規(guī)定read開頭的堿基是否要被切除的質(zhì)量閾值；
• TRAILING,，規(guī)定read末尾的堿基是否要被切除的質(zhì)量閾值,；
• MINLEN，規(guī)定read被切除后至少需要保留的長度,，如果低于該長度,，會(huì)被丟掉。

此外,，另一個(gè)值得注意的地方是,，Trimmomatic的報(bào)錯(cuò)給出的提示信息都比較難以定位錯(cuò)誤問題（如下圖），但這往往都只是參數(shù)用沒設(shè)置正確所致,。

小結(jié)

數(shù)據(jù)質(zhì)控的內(nèi)容終于都講完了,，接下來第4節(jié)就是主流程的構(gòu)建。