RNA-Seq能替代WES完成外顯子的變異檢測嗎?二代測序的四類Read重複如何產生?|《解螺旋技術交流圈》精華第2期

來自專欄解螺旋的礦工22 人贊了文章

1. RNA-Seq是否可以替代WES完成對外顯子的變異檢測？這不但省去用探針做外顯子捕獲這個步驟，代價更小；而且，它還在融合基因檢測上還更有優勢？

這是一個非常好的問題。我的回答是：RNA-Seq不能代替WES完成外顯子的變異檢測，原因如下：

(1). 轉錄本不是全部的外顯子。由於基因通過可變剪切出不同的轉錄本，實現多能性。那麼，沒被該轉錄本包括的外顯子就丟失了；

(2).轉錄本數據在基因上的覆蓋深度是極度不均勻的。不同基因的表達量不同，有些很高，有些甚至沒有。進行變異檢測的時候，這種不均勻性會極大影響變異結果的有效檢出。導致很多發現的變異可能都是那些高表達但是卻很可能不具備什麼關注點的基因上。如果這時你還是想獲得更多的變異，那麼到頭來還是得花更多的錢加大測序深度；

(3).目前對轉錄本數據進行變異檢測，還是一個偏於補充性質的分析。RNA-Seq的目的主要還是集中在基因表達方面，以及尋找差異表達基因和融合基因上。對於變異檢測，這類數據中也肯定可以發現，但假陰一定是很高的，比如低表達的基因，甚至是在這個組織（或者樣本）中不表達的基因，你就無法有效檢出它基因組上的變異了。另外，由於目前的二代測序系統並不能對RNA中的U鹼基進行識別，因此，RNA測序的時候需要先反轉錄為cDNA，這個過程會為RNA的變異檢測帶來一定程度的假陽結果；

(4).變異檢測範圍有限。使用RNA-Seq數據很難發現除單鹼基變異（SNV）之外的其他突變（比如Indel）。

2.為什麼二代測序的原始數據中會出現Read重複現象？

要搞清楚這個read重複（duplicate）的問題，我想我們需要從NGS數據的產出過程說起，具體來說如下：

1. 基因組DNA提取；

2. DNA隨機打斷，最常用的是超聲打斷；

3. 對被打斷的DNA片段進行末端修復（通常是3加A），然後在兩端加接頭，選擇特定長度的片段文庫進行PCR擴增（通過PCR的擴增會選！擇！性！地提高加上了接頭的文庫分子數量）；

4.文庫上機與測序晶元（Flowcell）上的引物結合，經過橋式PCR擴增，在晶元上形成測序所需的cluster；

5. 進行SBS測序，光學信號捕獲，生成序列。

我們一般認為第1步DNA提取出來的是完整的基因組，打斷則是完全隨機的——通常來說也確實如此。

在第3步，PCR擴增時，同一個DNA片段會產生多個相同的拷貝，第4步測序的時候，這些來源於同！一！個！拷貝的DNA片段會結合到Fellowcell的不同位置上，生成完全相同的測序cluster，然後被測序出來，這些相同的序列就是duplicate。這是duplicate的第一個來源，也是主要來源，稱為PCR duplicates（PCR重複）。

同樣，在第4步，生成測序cluster的時候，某一個cluster中的DNA序列可能搭到旁邊的另一個cluster的生成位點上，又再重新長成一個相同的cluster，這也是序列duplicate的另一個來源，這個現象在Illumina HiSeq4000之後的Flowcell中會有這類Cluster duplicates，這是第二類duplicate（如下圖）。

在第5步中，某些cluster在測序的時候，捕獲的熒光亮點由於光波的衍射，導致形狀出現重影（如同近視散光一樣），導致它可能會被當成兩個熒光點來處理。這也會被讀出為兩條完全相同的reads，這是第三類duplicate，稱之為Optical duplicates（光學重複）；

以上三種比較常見，還有第四種，稱為Sister duplicates，這是比較特殊的一個情況。它是文庫分子的兩條互補鏈同時都與Flowcell上的引物結合分別形成了各自的cluster被測序，最後產生的這對reads是完全反向互補的。比對到參考基因組時，也分別在正負鏈的相同位置上，在有些分析中也會被認為是一種duplicates。

另外，據說 NextSeq 平台上還出現過由於熒光信號捕獲相機移動位置不夠，導致 tile 邊緣被重複拍攝，每次採樣區域的邊緣由於重複採樣而出現了duplicates，下圖中藍色點代表 duplicates，可以看到在tile的左右兩側明顯富集。

以上，除了NextSeq的情況之外，所有這些不同類型的duplicates都各有特點。比如PCR duplicate的特點是隨機分布於Flowcell表面；而cluster duplicates和optical duplicates 的特點是它們都來自Flowcell上位置相鄰的cluster。Cluster的位置一般都會被記錄在原始測序fastq文件@Sequence-id那一行中。

這些Read重複都會一定程度上導致一些鹼基信號被錯誤地拉高或者減低，會對後續分析帶來干擾，特別是在WGS和WES分析時都需要去除。如果測序過程沒什麼特殊問題或者原因，那麼，測序數據的duplicate比例一般都在10%以下。

另外，PCR duplicates可以通過PCR-free來避免。並且PCR本身還會帶來一些其他的問題，比如擴增過程自帶了一定的偏向性，這會損失一定的測序隨機性，使得某些序列信息被擴大或者減小。所以，只要DNA起始量足夠，那麼我們就應該盡量採用PCR Free的方式來建庫。

3.GATK4不能進行多線程？

當我們本地跑GATK4的時候你會發現，它竟然沒有多線程的功能，這和GATK3很不同，在3中我們可以用-nt或者-nct設定多線程，但是4卻沒有類似的參數，這是為啥呢？

其實，關於這個問題已經有研究者問過GATK的團隊。簡單來說回答是，沒有！如果要使用多核來跑流程，那麼節點配置好Spark，用GATK4中的Spark功能模塊（如，HaplotypeCallerSpark）就可以了，讓Spark來幫你完成多線程。

我想可能很少有人知道以前GATK3中的多線程功能的效果其實並不好，而且還容易出問題。可能也是由於這方面的原因，GATK團隊這一次在4中就乾脆放棄了自己實現多線程的想法，直接使用現成的Spark來完成這個調度。另外，值得一提的是在GATK4中跑並行任務的最好做法是採用WDL和Cromwell相結合的方式。

話雖如此，但GATK團隊實際上還是留下了唯一的一個例外！那就是HaplotypeCaller中最消耗計算資源的模塊——pariHMM，這個是可以本地單獨多線程的！通過「--native-pair-hmm-threads」這個參數來設置，它默認是4，功能有些隱蔽！

4.Ti/Tv比率能說明變異的什麼特徵嗎？

Ti/Tv（轉換和顛換的比例）的值，實際上是物種演化的過程中在基因組上留下來的序列選擇標記，是對自然選擇的一種反應，這兒值在物種中具有一定的穩定性。因此，Ti/Tv的值常常會被我們作為一個評判變異的質控結果好壞的一個較為宏觀的指標。

而且一般來說，在基因組上C->T的轉換比較多，也就是胞嘧啶（C）容易自發轉為胸腺嘧啶（T），這是因為基因組上的C在甲基化的修飾下容易發生C->T的轉變。

另外，需要注意的是Ti/Tv畢竟是一個比值，它的大小和我們計算基因組上的哪一個區域有一定關係。比如對人類基因組而言，如果是全基因組區域，Ti/Tv比值一般在2.0x-2.1左右；而對於外顯子區域，Ti/Tv比值是3.0左右；對於那些novel variants來說，由於還沒受到明顯的自然選擇，它們的Ti/Tv的值會相對較低，可能會低至1.5左右。對其他區域而已，Ti/Tv範圍也應該在1.5-3.0之間。

有關Ti/Tv的更多內容建議看我之前討論變異數據應該如何質控的一篇文章。

5.最新人類參考序列的兩個版本NCBI（GRCh38）和UCSC（hg38）有什麼差異？

?可能大部分人都認為這兩份參考序列是一致的，但實際上它們存在著細微的差別。主要有以下四點：

(1). UCSC的hg38相比於NCBI的GRCh38缺少chrEBV（Epstein-Barr virus）序列。EBV本身不是人所有的，但由於很多細胞系在培養的過程中都需要藉助EBV，因此對於許多通過細胞系測序而來的數據中（比如海拉細胞系），就會混有這個序列。在NCBI的新版本中已經加入EBV序列，但是UCSC並沒有相應的更新；

(2). UCSC版本缺少decoy序列，這是參考序列中沒有，但卻是其他人群特有的人類序列，作為參考序列不應該缺少；另外還缺少HLA分型序列；

(3). 簡併鹼基的表達方式不同。NCBI的GRCh38中存在極少量的簡併鹼基，一共有94個，對於巨大的人類基因組序列來說，這個數字可以忽略，這些簡併鹼基在UCSC的版本中則用N代表；

(4). Y染色體上存在較大差別。NCBI上的GRCh38，在Y染色體上的兩個PAR區域（pseudoautosomal region，偽染色體區域）用N來代替了，而UCSC並沒有做類似的處理。這樣會導致它們在Y染色體序列中存在明顯的差異，進行序列比對時結果也會有所不同。

那麼這兩份我們應該用哪一份？答案是NCBI的版本。它也是GATK bundle所使用的版本，並且NCBI對其更新和維護的節奏明顯好於UCSC版本，細節也做得更好，一般我們會直接在GATK的bundle中下載。不久前UCSC也把NCBI GRCh38直接接入到它的基因組瀏覽器了（上圖）。

技術交流圈往期精華

• RNA-seq原始數據質控後，是否要合併PE和SE的比對結果
• 我是解螺旋的礦工，我熱愛生命科學
• 該如何自學入門生物信息學

推薦閱讀：