一篇文章說清楚基因組結構性變異檢測的方法

weixin_34116110發表於2018-07-23

這源自於我2013年的一篇舊文，另外，那篇文章也已被轉載得有些走樣了，我決定要重新組織。雖然主體沒有大的變化，但是內容在這裡作了較多的翻新。對於初學者來說，如果你之前看過了我的那篇舊文，還是建議你再看一次，或許可以有新的發現，如果你從未看過，那麼請不要輕易錯過。

人類基因組中的變異和人類的演化、疾病風險等方面都有著密切的聯絡。當前二代短讀長高通量測序技術（NGS），雖然能夠讓測序成本大大降低，但這種短讀長的測序方法也給基因組的變異檢測（特別是結構性變異檢測）帶來了不小的挑戰。SNP和Indel大家應該都見得比較多了，因此在這篇文章裡我將主要討論常見結構性變異的檢測方法和有關軟體以及它們的一些優缺點。

變異的分類

在開始之前，有必要先梳理一下人類基因組上的變異種類，按照目前業界的看法可以分為如下三個大類：

單鹼基變異，即單核苷酸多型性(SNP)，最常見也最簡單的一種基因組變異形式；
很短的Insertion 和 Deletion，也常被我們合併起來稱為Indel。主要指在基因組某個位置上發生較短長度的線性片段插入或者刪除的現象。強調線性的原因是，這裡的插入和刪除是有前後順序的與下述的結構性變異不同。Indel長度通常在50bp以下，更多時候甚至是不超過10bp，這個長度範圍內的序列變化可以通過Smith-Waterman 的區域性比對演算法來準確獲得，並且也能夠在目前短讀長的測序資料中較好地檢測出來；
基因組結構性變異（Structure Variantions，簡稱SVs），這篇文章的重點，通常就是指基因組上大長度的序列變化和位置關係變化。型別很多，包括長度在50bp以上的長片段序列插入或者刪除（Big Indel）、串聯重複（Tandem repeate）、染色體倒位（Inversion）、染色體內部或染色體之間的序列易位（Translocation）、拷貝數變異（CNV）以及形式更為複雜的嵌合性變異。

圖1. 結構性變異的不同種類

值得一提的是，研究人員對基因組的結構性變異發生興趣，主要還是由於在研究中發現：

SVs對基因組的影響比起SNP更大，一旦發生往往會給生命體帶來重大影響，比如導致出生缺陷、癌症等；
有研究發現，基因組上的SVs比起SNP而言，更能代表人類群體的多樣性特徵；
稀有且相同的一些結構性變異往往和疾病（包括癌症）的發生相互關聯甚至還是其直接的致病誘因。比如，《我不是藥神》電影中提到的慢粒白血病，它就和基因組的結構性變異直接相關。它是由於細胞中的9號染色體長臂與22號染色體長臂相互易位，導致ABL基因和BCR基因融合，形成了一個會導致ABL異常表達的小型染色體（稱為費城染色體）發生的。這是一個典型的結構性變異致癌的例子。

圖2. 9號染色體和22號染色體長臂發生易位，形成"費城染色體"，常導致慢性粒細胞白血病。沒錯！也就是《我不是藥神》中可以用”格列衛“治療的一種白血病

類似的例子還有很多。不過應該注意的地方是，大多數的結構性變異並不直接與疾病的發生相關聯，很多也不一定會致病，它還與周圍環境的響應或者其他的一些表型多型性有關。

一個東西一旦引起多方關注，很快就會紅起來，SVs也是這樣的。近年來，在NGS技術快速發展的大背景下，SNP和Indel這類比較容易發現的變異的檢測演算法已經做得七七八八了，大家就開始選擇新的未解問題——研究熱點，所以在這些情況的加持下，人類基因組上的結構性變異就逐漸被全面而又集中地研究了起來，相應的變異檢測工具和演算法也越來越多，如下附圖3，在omictools上就有超過130個與SVs有關的工具。

除了測序的方法之外，一些專用的基因晶片也可以用於檢測一部分結構性變異，好處是成本低，但是侷限大，而且晶片技術實際上只對大片段的序列刪除這類變異比較靈敏，無法做到全面，同時做不到單鹼基精度的斷點檢測。

圖3. 數量眾多的SVs檢測工具

工具雖然不少，但其實歸納起來，基於NGS資料的變異檢測演算法並不多。總的來說主要有以下4種不同的策略和方法，分別是：

Read Pair，一般稱為Pair-End Mapping，簡稱RP或者PEM；
Split Read，簡稱SR；
Read Depth，簡稱RD，也有人將其稱為RC——Read Count的意思，它與Read Depth是同一回事，顧名思義都是利用read覆蓋情況來檢測變異的方法；
序列從頭組裝（de novo Assembly，簡稱AS）的方法。

圖4. 四種類不同的演算法進行結構性變異的檢測分類，行是結構性變異型別，列是演算法

接下來我將對這四種不同的方法以及它們各自的特點逐一進行展開介紹。但是，具體要用什麼軟體，引數要怎麼設定，命令要怎麼跑，這寫操作性的內容我就不在這篇文章中展開了，大家可以看其參考文件（如果它的參考文件沒有寫明白，那麼多半說明這個軟體做的不怎麼樣，慎用）。這篇文章的目的是想讓大家對這些SVs方法有一個巨集觀上的認識，讓你知道原來基因組上結構性變異檢測是這樣的。

基於Read Pair（RP），也就是Pair-end Mapping（PEM）的方法

圖5. 一般的RP流程

這是RP方法的一個主要框架，用的比較多。它的思想是這樣的：我們知道PE測序的兩條read（通常稱為read1和read2），它們原本就是來自於同一個序列片段的（如下圖，對此不清楚的同學也可以檢視這一篇文章）。

圖6. Pair-End 測序

因此，Read1和Read2之間存在著客觀的物理關聯（它們就是在一起的），而且它們之間的距離——圖中這個淡藍色序列片段（通常稱為插入片段）的長度，稱為插入片段長度（Insert size）。一般來說我們是無法直接獲得每一對read1和read2之間真實的插入片段長度的，但通過序列比對，計算它們彼此之間比對位置上的距離卻可以間接獲得這個長度——BAM檔案的第9列記錄的就是這個值（這裡是BAM的格式和關於使用Pysam處理BAM檔案的文章[https://mp.weixin.qq.com/s/c0O5qHBnybZNKCERLzYbJQ]中我還舉了提取插入片段長度的例子），如下圖：

圖7. 在BAM檔案中可以方便地得到Insert size資訊

這個插入片段長度的分佈是RP方法進行變異檢測的一個關鍵資訊。我們知道，在測序之前需要先用超聲或者酶切的辦法把原始的DNA序列進行打斷處理，然後再挑選某一個長度（比如說500bp）的DNA片段來上機測序。在這個過程中，雖然我們都希望挑選的序列是一樣長的，但這肯定是不可能的。事實上，我們得到的片段長度通常都會圍繞某個期望值（比如我們這裡的500bp）做上下波動。如果我們把這些按照比對位置計算獲得的“插入片段長度”提取出來做個分佈圖，會看到它的樣子通常如下：

圖8. Insert size 長度分佈

沒錯！如無意外你見到的一般都是類似這樣的一個鐘形線——正態分佈，理論上它也應該是一個正態分佈。而且由於基因組上存在變異的序列畢竟是少數，所以，這個分佈本身就能夠反映“插入片段長度”在絕大多數情況下的真實分佈情形，那麼對於那些小部分不能反映真實情況（偏離分佈中心的片段長度）的是什麼呢？是隱含的變異訊號！

因為，如果插入片段長度有異常，它實際上包含的意思是，組成read1和read2的這個序列片段和參考基因組相比存在著序列上的變異。舉個例子，如果我們發現它這個計算出來的插入片段長度與正態分佈的中心相比大了200bp（假設這個200bp已經大於3個標準差了），那麼就意味著參考基因組比read1和read2所在的片段要長200bp，通過類似這樣的方式，我們就可以發現read1和read2所在的序列片段相比與參考基因組而言發生了200bp的刪除（Deletion）。RP除了可以利用異常插入片段長度的資訊進行線性變異（特指Deletion和Insertion）的發現之外，通過比對read1和read2之間的序列位置關係，還能夠發現更多非線性的序列變異。比如，序列倒置（Inversion），因為，按照PE的測序原理，read1和read2與參考基因組相比對，正好是一正一負，要麼是read1比上正鏈，read2比上負鏈，要麼是反過來，而且read1和read2都應處於同一個染色體上，如果不是這種現象，那麼就很可能是序列的非線性結構性變異所致，比如前者是序列倒置（Inversion），後者是序列易位（Translocation）等。

總的來說，RP通過利用比對距離、read1和read2之間的位置關係這兩個重要資訊實現了基因組上多種結構性變異的檢測。

利用這樣的原理，RP方法理論上能夠檢測到的變異型別可以包含如下6種：

序列刪除（Deletion）
序列插入（Dnsertion）
序列倒置（Inversion）
染色體內部和染色體之間的易位（intra- and inter-chromosome translocation）
序列串聯重複/倍增，也就是常說的Tandem duplications
序列在基因組上的散在重複（Interspersed duplications）

圖9. 適合RP策略檢測的6類結構性變異，藍綠色的小塊代表的是測序Read

應該說檢測的範圍還是比較廣的，但對於RP存在缺陷的地方需要指出：

第一，對於Deletion的檢測，由於要求插入片段長度的變化要具有統計意義上的顯著性，所以它所能檢測到的片段長度就會受插入片段長度的標準差（SD）所影響。簡單來說，就是越大的序列刪除越偏離正常的長度中心，才越容易被檢測到。RP對於大長度Deletion（通常是大於1kbp）比較敏感，準確性也高，而50bp-200bp的這個範圍內的變異，由於基本還處於2倍標準差以內，在統計檢驗上它的變化就顯得不那麼顯著，所以這通常就成了它的一個檢測 暗區。
第二，它所能檢測的Insertion序列，長度無法超過插入片段的長度。因為，如果這個Insertion序列很長——舉個極端的例子——整個插入片段都是Insertion序列，那麼你會發現read1和read2根本就不會比對上基因組，它在基因組上一點訊號都沒有，你甚至都不知道有這個序列的存在。另外，Insertion的檢測精度也同時受限於插入片段長度的標準差。

Delly和Breakdancer都是應用RP方法的SVs檢測軟體，應該也是當前用得比較廣的軟體，可能相比於Breakdancer，Delly用得還要更多一些。此外，Delly不僅在進行結構性變異檢測的時候應用了RP，它還同時使用了SR的方法。

Split Read（分裂read，簡稱SR）

SR這個方法，起源於Sanger長片段測序資料變異檢測的場景，而且如果單條序列足夠長（>400bp），那麼它還可以有效地檢測出MEI（Mobile Element Insertion）。NGS迅速發展起來之後，測序序列變成了短序列，不過SR同樣被應用過來了。

插播一句，媒體常常把NGS被翻譯為“下一代”測序技術，實際上這是一種不明所以的翻譯，說是二代測序技術也就算了，這個“下一代”是啥？是還沒來嗎？現在都到三代了。還不如說是測序的2G和3G時代。所以，我還是更喜歡用NGS這三個字母。

在短序列資料中 SR演算法的核心也是對非正常PE比對資料的利用。這裡有必要和上文RP（Read-Pair）中提到的非正常PE比對做個區分。RP中的非正常比對，通常是read1和read2在距離或者位置關係上存在著不正常的情形，而它的一對PE read都是能夠“無傷”地進行比對的；但SR一般是指這兩條PE的read，有一條能夠正常比對上參考基因組，但是另一條卻不行的情形。這個時候，比對軟體（比如BWA）會嘗試把這條沒能夠正常比上基因組的read在插入片段長度的波動範圍內，使用更加寬鬆的Smith-Waterman區域性比對方法，嘗試搜尋這條read最終可能比對得上的位置。如果這條read有一部分能夠比上，那麼BWA會對其進行軟切除——soft-clip（CIGAR序列中包含S的那些read，圖9中也有這類比對情況），標記能夠比上的子序列（但不能比上的序列還是留在原read中，這也是軟切除的含義）。但這個過程有時候可能不會太順利，甚至會發生多次切除再比對的情況，所以，你會看到一條read有時候竟然有很多個Soft-clip的比對結果（這個不是指多比對的情況，而是多次切除的情況），而發生這些情況的read就是SR方法的用武之地。並且soft-clip保留原序列的方式對於後續應用SR很重要，因為，它們往往不會只是依賴原有BWA的比對結果，而是會對這條read進行重新區域性比對（如果沒有保留那麼資訊的丟失就會導致大量的假陰），比如Pindel就是這麼做的，下文詳述。

除了Pindel，Delly也是應用SR方法進行變異檢測的軟體，不過Pindel在早期千人基因組計劃中用的最多，Delly是後起之秀，當然現在有代表性的工具還有不少，比如lumpy、SVseq2等。

圖10，展示了SR方法所能檢測到的變異型別以及它是如何用這些Split reads的訊號來進行結構性變異檢測的，其中的紅色線段就是被split出來並重新比對的序列。可以看出在不同的變異型別中，這些split的訊號是不同的，而且它們和RP相比也很不同，SR的read通常都會被“撕裂”出來，而RP則是“無傷”完整的read，大家注意對比這個區別。

圖10. SR所能夠檢測的變異型別和方法特點

以Pindel為例子，再說一下其利用SR進行變異檢測的思想原理。

首先，在獲得了單端唯一比對到基因組上的PE read之後，Pindel會將不能正常比上的那條read切開成2或者3小段。然後，分別按照使用者設定的最大deletion長度重新進行比對，並獲得最終的比對位置和比對方向，而斷點位置的確定就根據soft-clipped的結果來獲得。

Pindel理論上能夠檢測所有長度範圍內的Deletion以及小片段的Insertion（<50bp），Inversion，Tandem duplication和一些large insertion也都能夠找到。

SR和RP有一個共同之處，同樣會利用比對的方向去判定相關的變異，它倆雖然方式和方法有區別，但所檢測到的變異是存在重疊和互補關係的。SR的一個優勢在於，它所檢測到的SVs斷點能精確到單個鹼基，但是也和大多數的RP方法一樣，無法解決複雜結構性變異的情形。

而且 對於SR來說，它要求測序的read要更長才能體現它的優勢，比如上文我提到的，如果有400bp以上的長度，那麼一些MEI或者Alu序列的變異都能檢測到，read太短，許多變異都會不可避免地被漏掉，它的檢測功效在基因組的重複區域也會比較差。

Read Depth

Read Depth（有時也叫Read Count）簡稱RD，是目前解決基因組拷貝數變異檢測（Copy number variantion，簡稱CNV）的主要方法。拷貝數變異通常包括，序列丟失和序列重複或倍增兩個大類，在腫瘤基因組資料分析中用的比較多。其實CNV實質上是序列Deletion或Duplication，是可以歸類於Deletion和Insertion這個大的分類的，只是由於它的發生有著其獨特的特點，而且往往還比較長，所以也就習慣了獨立區分。

RD的原理基於read覆蓋深度。全基因組測序（WGS）得到的覆蓋深度呈現出來的是一個泊松分佈——因為基因組上任意一個位點被測到的機率都是很低的——是一個小概率事件，在很大量的測序read條件下，其覆蓋就會呈現一個泊松分佈，如下圖。

圖11. 典型的WGS覆蓋深度分佈圖

目前有兩種利用Read depth資訊檢測CNV的策略。一種是，通過檢測樣本在參考基因組上read的深度分佈情況來發現CNV，這類適用於單樣本，也是用的比較多的一個方法；另一種則是通過識別並比較兩個樣本在基因組上存在丟失和重複倍增的區域，以此來獲得彼此相對的CNV，適用於case-control模型的樣本，或者腫瘤樣本Somatic CNV的發現，這有點像CGH晶片的方法。

CNVnator使用的是第一種策略，同時也廣泛地被用於檢測大的CNV，當然還有很多冷門的軟體，這裡就不再列舉了；CNV-seq使用的則是第二種策略。

基於其原理，RD的方法可以很好地檢測一些大的Deletion或者Duplication事件，但是對於小的變異事件就無能為力了（如圖9所示）。

基於 de novo assembly

其實從上面看下來，SVs檢測最大的難點實際上是read太短導致的。就因為read太短，我們不能夠在比對的時候橫跨基因組重複區域；就因為read太短，很多大的Insertion序列根本就沒能夠看到資訊；就因為read太短，比對才那麼糾結，我們才需要用各種數學模型來猜測這個變異到底應該是什麼等等。

那麼既然如此，就去想辦法加長read的長度啊！辦法有兩個：三代長read測序和序列從頭組裝（de novo assembly）。

目前不太好評價三代測序和從頭組裝哪個更好，它們各有優勢也各有各的問題。比如，三代測序錯誤率高和對Indel錯誤的引入會帶來比較大的糾錯成本，而組裝則對資料量有比較高的要求，而且裝重複區域也有困難等。

但從理論上來講，三代測序和de novo assembly 的方法應該要算是基因組結構性變異檢測上最有效的方法，它們都能夠檢測所有型別的結構性變異。以組裝為例子，就目前來說，它是大長度 Insertion和很多複雜結構性變異的最好檢測方法（這是所有基於短序列比對方法的盲區）。

我們在丹麥人國家基因組專案中用的就是這方法，其中的基於長序列的SVs檢測演算法也是我主要開發的，效果還不錯。我們發了三篇文章（第一篇Pilot1在《Nature Communication》上，方法學一篇在《GigaScience》上，以及最後的主文章在《Nature》上），應該說是做出了有史以來最完整的人群基因組結構性變異圖譜，填補了以前很多地方上的空白（圖13），這個專案從2015年到目前為止加上我們和丹麥三所高校應該已經發了10來篇文章了。

圖12. 柱子的白色部分是其它檢測演算法所發現不了的SVs暗區，但我們在利用組裝的方法在丹麥人基因組中解決了

但就像我前面說到的，序列從頭組裝要裝的好，還是比較棘手的，對高等植物和脊椎動物來說，也是如此。最主要的原因在於，這些物種基因組上所存在的重複性序列和序列的雜合會嚴重影響組裝的質量，除去資金成本，這也在很大程度上阻礙了利用組裝的方法在基因組上進行變異檢測的應用。

最後，結構性變異檢測的方法這麼多，我該如何選擇

image

上面這個表列出了目前主流工具所適合的SVs，可以根據實際情況進行選擇。

其實通過上面對四種不同SVs檢測策略的比較也可以發現，小長度範圍內的變異以及較長的deletion，問題不大，但對於大多數的Insertion和更復雜的結構性變異情況，當前的檢測軟體基本都沒法還解決。Assembly應是當前全面獲得基因組上各種變異的最好方法，但是目前的侷限卻也發生在Assembly本身，若是基因組沒能裝得好，後面的變異檢測就更是無從說起。從目前的情況看，de novo assembly的方法並不能很快進入實際的應用。因此，暫且不提assembly，其餘的三種策略都各有各的優勢，從目前的結果看，並沒有哪一款軟體能夠一次性地將基因組上的各種不同情況變異型別都獲得。因此就目前短reads高通量測序技術來說，最合適的方案應是結合多個不同的策略，將結果合併在一起，這樣可以最大限度地將FN和FP降低。Speedseq和HugeSeq應該是這方面的一個典範。

文末，再來一個各類SVs檢測方法大比拼的大表供大家參考，很有價值不要錯過哦，它來自於這篇文章，如下：

各類SV檢測演算法大比拼

參考文獻

DePristo, M. a et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nature genetics43, 491–8.
Alkan, C., Coe, B. P. & Eichler, E. E. Genome structural variation discovery and genotyping. Nature reviews. Genetics12, 363–76.
Mills, R. E. et al. Mapping copy number variation by population-scale genome sequencing. Nature470, 59–65.
Africa, W. A map of human genome variation from population-scale sequencing. Nature467, 1061–73.
Korbel, J. O. et al. PEMer: a computational framework with simulation-based error models for inferring genomic structural variants from massive paired-end sequencing data. Genome biology10, R23.
Ye, K., Schulz, M. H., Long, Q., Apweiler, R. & Ning, Z. Pindel: a pattern growth approach to detect break points of large deletions and medium sized insertions from paired-end short reads. Bioinformatics (Oxford, England)25, 2865–71.

你還可以看

我的微信公眾號：解螺旋的礦工