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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

得到clean 序列后，利用比對(duì)軟件，將clean 序列比對(duì)到參考基因組上去。比對(duì)簡(jiǎn)單的說就是對(duì)每一個(gè)序列， 找到它在參考基因組上對(duì)應(yīng)的的位置。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

答：是的，利用WES檢測(cè)純合和雜合的CNV都比較困難。這是因?yàn)槿祟惖耐怙@子都比較短，一般都小200bp。另外，WES測(cè)序的測(cè)序深度不均一， 用reads count檢測(cè)CNV的話也會(huì)比較困難。

5.老師，因?yàn)橥怙@子短于200，所以不好看CNV這個(gè)不明白。是因?yàn)槿钡奶郙APPING不上了嗎。是缺的比如100bp  剩下80bp是不能比對(duì)到參考基因嗎？

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

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

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