各位老師好，很高興今天有機會和大家一起來探討全外顯子數(shù)據(jù)分析中的相關(guān)知識和技術(shù)。

顧名思義，全外顯子測序是對基因組中所有外顯子進行測序分析的方法。該方法首先利用DNA或RNA探針，將全基因組中外顯子區(qū)域DNA進行捕獲，然后對捕獲的DNA進行PCR擴增，最后對擴增后的產(chǎn)物進行高通量測序。雖然外顯子只占基因組的2%，卻含有85%的已知致病變異， 但測序和分析價格卻遠比全基因組測序要低得多。所以目前很多科研和臨床選用全外顯子測序，而不是全基因組測序。

和WGS相比，WES有以下優(yōu)勢：

1.測序和數(shù)據(jù)存儲的花費低。盡管WES一般比WGS測序深度高（100X vs 30X), WES的測序序列主要在占基因組2%的外顯子區(qū)域，使得總的測序文件比WGS的測序文件小得多。平均測序深度為100X的全外顯子測序文件的大小是6GB左右，而30X的全基因組測序文件的大小一般是90GB左右;

2.由于數(shù)據(jù)文件小，數(shù)據(jù)分析的速度就快，數(shù)據(jù)分析的花費也低； 

3.由于測序深度高，WES更容易檢測出罕見變異；

4.跟非編碼區(qū)相比，蛋白編碼區(qū)的研究更透測，注釋度更高，所以WES數(shù)據(jù)分析中的變異能夠更好地被解讀。 

盡管和WGS相比，WES存在很多優(yōu)勢，但也有不足之處：

1.WES的測序深度不均一，一些區(qū)域測序深度過高造成浪費，也有一些區(qū)域測序深度過低而無法檢測到存在的變異。比如一些SNP密集的區(qū)域，富集過程中RNA或DNA探針無法和該區(qū)域雜交導致無法有效捕獲，該區(qū)域的測序深度就會過低；

2.由于WES只測外顯子區(qū)域（人類的外顯子區(qū)域長度一般小于200bp），因而無法有效檢測CNV，結(jié)構(gòu)重組和其他大的結(jié)構(gòu)變異；

3.由于沒有捕獲探針，新的編碼基因無法捕獲。

在高通量測序和分析中有三個常用的質(zhì)量分數(shù)，了解這些質(zhì)量分數(shù)將有助于對外顯子數(shù)據(jù)分析的理解。 第一個是堿基質(zhì)量分數(shù)，它是測序過程中堿基被測序儀錯誤識別的概率。 另一個是比對質(zhì)量分數(shù)， 它是一個序列在參考基因組上比對正確與否的可信度。還有一個是變異質(zhì)量分數(shù)，它是我們檢測到的變異是否是生物變異的可信度。變異質(zhì)量分數(shù)是變異位點所有序列比對質(zhì)量分數(shù)的均方根。

由美國哈佛-麻省理工博德研究所 (Broad institute of MIT and Harvard) 開發(fā)的基因組分析工具包GATK（ genome analysis toolkit)是目前常用的NGS數(shù)據(jù)分析軟件。今天我們以GATK的best practices 為列來介紹全外顯子數(shù)據(jù)分析的流程。WES 數(shù)據(jù)分析大致分為三步：對WES數(shù)據(jù)進行前處理；利用前處理后的數(shù)據(jù)進行變異檢測，最后對檢測到的變異進行篩選和注釋。

第一部分：WES數(shù)據(jù)的前處理流程。

WES數(shù)據(jù)分析的前處理包括對原始fastq形式的序列進行質(zhì)量控制，將清理后的fastq形式的序列比對到參考基因組，對比對到參考基因組的序列進行多步質(zhì)量控制，最終得到可以用來檢測變異的序列。

一般情況下，我們從測序公司拿到的WES數(shù)據(jù)，是原始數(shù)據(jù)狀態(tài)（raw data)。 拿到數(shù)據(jù)后的第一步一定要查看fastq文件是否完整，是否存在錯誤，以及數(shù)據(jù)質(zhì)量的好壞。常用的軟件是FASTQC. 根據(jù)FASTQC的結(jié)果報告，如果fastq文件有錯誤，就無法進行下面的分析。如過fastq文件沒有錯誤，就對原始數(shù)據(jù)進行清理，包括去除低質(zhì)量或者過短序列，剪切末端低質(zhì)量的堿基，去除接頭污染以及可能的外源污染序列。 常用的軟件包括Trimmamatic, cutadapt 和fastq_mcf等。清理之后得到clean 序列。對于NGS數(shù)據(jù)質(zhì)量控制的知識和技術(shù)，請參考賽福基因微課程的第四講。

得到clean 序列后，利用比對軟件，將clean 序列比對到參考基因組上去。比對簡單的說就是對每一個序列， 找到它在參考基因組上對應的的位置。

理論上，將序列比對到參考基因組上非常簡單，就是找到每個序列在參考基因組上的位置并記錄下來。但事實上，很多情況下，測序數(shù)據(jù)和參考基因組間存在著各種差異，包括SNV，INDEL，CNV， SV等，這使得比對變得非常復雜。比對軟件需要運用復雜的算法，將序列比對到可信度最高的一個或者少數(shù)幾個位置上去。這個可信度就是我們前面提到的比對質(zhì)量分數(shù)（mapping quality ）。另外，如這里的示意圖，由于參考基因組中存在重復序列等原因，有些序列可以同等地比對到不同的位置，這也增加比對的難度。

將序列比對到參考基因組常用的軟件包括BWA(Burrows-Wheeler Aligner), Bowtie2, Blat, SOAP 等。 比對之后的結(jié)果文件以SAM 的格式存儲。 如果要對比對結(jié)果進行其他一些操作， 一般需要將SAM文件轉(zhuǎn)化成壓縮的二進制格式的BAM文件。

如果我們想直觀地查看比對結(jié)果，可以使用可視化軟件。常見的可視化軟件有IGV，tablet等。用IGV或tablet可視化比對結(jié)果時，就必須要用加過索引的BAM文件。這里顯示的是一個WES數(shù)據(jù)在ARX基因區(qū)的比對結(jié)果。我們可以看到，外顯子區(qū)域以外也有比對的序列。 一般情況下，這些序列不要去除，更大范圍內(nèi)更多的比對序列有利于下游更精確的變異檢測。在IGV里，我們可以對感興趣的區(qū)域進行放大。

這個是我們對上圖中中間的那個外顯子區(qū)進行放大的結(jié)果。在IGV里default的設置下，比對到參考基因組上的序列有不同的顏色?；疑硎驹撔蛄性诨蚪M上的比對結(jié)果是單一的，而且和序列對中另一個序列間的距離在正常范圍內(nèi)。白色說明該序列比對到2個或者2個以上的位置。如果reads是其他顏色，說明與該序列配對的序列比對到了其它的染色體， 或者兩個序列間的距離大于或小于正常的插入序列長度。

另外，可以用samtools或者bamtools查看比對結(jié)果的統(tǒng)計值。 samtools的flagstat軟件，可以查看總的比對率， 成功比對的序列里多少個是序列對，多少是單個的序列等。Samtools 的idxstat軟件可以查看每個染色體的長度，該染色體上比對的序列數(shù)以及沒有比對的序列數(shù)。

原始比對結(jié)果里含有一些比對質(zhì)量較低或者錯誤比對的序列，為了提高變異識別的敏感性和精確性，要對原始比對結(jié)果進行質(zhì)量控制。 這些質(zhì)量控制包括：去除重復序列，對可疑的比對區(qū)域進行重新比對（Local Realignment）和重新校正堿基質(zhì)量分數(shù)。 對于這些質(zhì)量控制的相關(guān)知識和技術(shù)，請參考賽?；蛭⒄n堂第四講。比對結(jié)果經(jīng)過質(zhì)控后，就可以進入WES數(shù)據(jù)分析的第二步了，即檢測樣品中的變異。

第二部分：全外顯子數(shù)據(jù)變異檢測流程。

比對結(jié)果通過質(zhì)量控制后，我們可以開始變異識別的流程。 也就是說找出比對結(jié)果里，樣品數(shù)據(jù)和參考基因組不同的位點，并計算這些位點的基因型。　

在變異識別軟件盡可能多地識別真正的生物變異時，也會識別一些非生物變異，即由比對或者測序錯誤帶來的數(shù)據(jù)和參考基因組之間的差異，這些被識別的差異稱為假陽性變異。在得到原始變異文件后，需要盡可能地去除所有非生物學變異。把真正的生物學變異從來自于系統(tǒng)誤差的非生物學變異分離出來是一個比較棘手的問題。GATK 的VQSR （Variant Quality Score Recalibration）是目前最好的去除非生物學變異的軟件。

VQSR的理論基礎是： 變異軟件識別的每個原始變異都有一套對應的注釋參數(shù)。根據(jù)這些參數(shù)做聚類分析，真正的變異趨向于聚集在一起。而這些聚類形成的組群遵循高葉斯分布。根據(jù)已知的變異建立高葉斯混合模型，然后根據(jù)建立的模型對可能的新變異進行評估，進而去除假陽性的變異。

這個是2011年發(fā)表在nature genetics上的一篇文章利用高葉斯模型評估新變異的例子。左圖是利用HapMap SNP作為訓練數(shù)據(jù)建立的模型，右圖是根據(jù)左圖的模型對新SNP評估的結(jié)果。紫色的部分和左圖中橘色的部分類似，很有可能是error， 是假陽性SNP， 所以要過濾掉。 綠色部分和左圖中紅色部分類似，很有可能是真正的新發(fā)現(xiàn)SNP。

在做VQSR時，SNP和INDEL要分別進行質(zhì)量校正。但可以將VQSR連續(xù)run兩次，而不必將SNP和INDEL分離、分別run VQSR之后再合并。以原始的包括SNP和INDEL的文件作為輸入，在run第一步時，只對SNP進行質(zhì)量分時校正，而INDEL不受影響。第二步對INDEL做質(zhì)量分數(shù)校正，而SNP不受影響。

VQSR是目前最好的識別并過濾假陽性變異的方法，但只有原始變異集足夠大的情況下，VQSR才能正常運行并過濾假陽性變異。在原始變異集不夠大，不能夠運行VQSR的情況下，可以通過對注釋參數(shù)設定閾值來手工過濾非生物學變異。這種情況下，就必須先分離SNP和INDEL，分別手工過濾后再重新合并。

從前面IGV可視化比對結(jié)果的圖里可以看出，比對到參考基因組的序列不只是在外顯子區(qū)，很多在外顯子以外的區(qū)域。那么變異軟件識別的變異里有些也會分步在外顯子區(qū)域以外。利用bedtools和捕獲區(qū)域文件，可以去除外顯子區(qū)以外的變異，只保留外顯子區(qū)內(nèi)的變異。

在得到外顯子區(qū)域內(nèi)高可信度的變異后，我們可以選取可靠的遺傳致病數(shù)據(jù)庫OMIM、HGMD、Clinvar、dbSNP等作為主要檢索數(shù)據(jù)庫對變異進行注釋，然后結(jié)合樣本來源的遺傳信息或科研需要，篩選變異。這些內(nèi)容我們會在下一節(jié)課里進行講述。

感謝各位老師在百忙之中抽出時間收看及收聽我們這次課程。歡迎各位老師加入我們的基因大課堂共同探討、共同交流。

1.什么是gVCF？gVCF和普通的VCF有什么區(qū)別？

答：VCF是Variant Call Format的縮寫。是標準化的存儲SNP，INDEL，和結(jié)構(gòu)變異的text文件。VCF里只含有變異的位點信息。  gVCF是genomic VCF。是VCF的一種。gVCF 和VCF最大的區(qū)別是gVCF里含有基因組（或者感興趣的區(qū)間）所有位點的堿基信息，不論該位點是否存在變異。

2.在數(shù)據(jù)量多大的情況下不可以用VQSR，而只能用手工過濾？

答：一般情況下，全外顯子測序樣品少于30個，或者全基因組測序樣品少于2個，VQSR都不能很好的運行。即使運行通過不出錯，結(jié)果也不可信。

3.怎樣知道bam文件是否已經(jīng)排序了？

答：可以用samtools查看bam文件的header。 如果在header里有SO:coordiate的字樣，就是已經(jīng)sort過的。

4.我先想問下WES為啥不能有效檢測CNV，這個不是特別理解。這個人是純合CNV或雜合CNV都不能有效檢測嗎？

答：是的，利用WES檢測純合和雜合的CNV都比較困難。這是因為人類的外顯子都比較短，一般都小200bp。另外，WES測序的測序深度不均一， 用reads count檢測CNV的話也會比較困難。

5.老師，因為外顯子短于200，所以不好看CNV這個不明白。是因為缺的太多MAPPING不上了嗎。是缺的比如100bp  剩下80bp是不能比對到參考基因嗎？

答：不是mapping不上。而是說一般的CNV或者其他的SV都比較大，一般大于外顯子。

6.有reads  您覺得counts數(shù)不準。 老師，這個測序深度不均一性是測序建庫技術(shù)操作的問題 ，還是每個人的個體化差異導致不均一的現(xiàn)象。比如特定位置A基因是每個人的樣本都會那不均一嗎，還是有的人均一有的人不均一？

答：不均一的原因個體DNA和測序應該都有，看不同的情況。比如我們曾經(jīng)分析過一個WES, 在一個本應該檢測到變異的基因，我們無論如何檢測不到?？梢暬l(fā)現(xiàn)該區(qū)域內(nèi)沒有reads。而查看這段基因序列，發(fā)現(xiàn)是100% 的G。而對測序儀來說，如果G含量超過80%，就很難成功測序。而如果個體某段DNA含有比較多的SNP，捕獲探針就無法很好和DNA雜交，而不能有效的捕獲。