基于高通量測序數據的插入/刪除新突變檢測方法

邢文昊， 劉永壯， 王亞東

(哈爾濱工業大學 計算機科學與技術學院， 哈爾濱 150001)

0 引 言

新突變是指并非遺傳自親代，而是在子代中首次出現的基因組變異，這種變異既可能來自于生殖細胞，也可能來源于胚胎發育早期的受精卵。大量研究表明新突變與人類疾病的發生密切相關，尤其是自閉癥、癲癇等精神類疾病[1]。隨著高通量測序技術的迅猛發展，基于高通量測序數據進行新突變檢測已經成為常規手段。目前的基于高通量測序數據的新突變檢測方法大概可以分為2類。一類是在親代和子代樣本中獨立地進行基因組變異檢測，然后通過比較子代和親代的基因型發現新突變；另一類是對親代和子代基因組數據建立聯合檢測模型直接檢測新突變[2]。例如PolyMutt[3]和Triodenovo[4]均為基于似然度的新突變檢測方法，該方法分別計算無孟德爾遺傳約束和有孟德爾遺傳約束下的最大基因型似然度，然后根據似然比確定候選位點是否是新突變，似然比越大，候選位點是新突變的可能性越高；DeNovoGear[5]是一種基于貝葉斯的新突變檢測方法，該方法利用貝葉斯公式計算候選位點是新突變的后驗概率，并根據后驗概率的大小來判斷候選位點是否是新突變；DNMFilter[2]是一種基于Gradient Boosting[6]的新突變檢測方法，該方法利用Gradient Boosting作為分類器，從子代和親代序列比對上下文選擇序列特征，該方法相比其它方法能夠顯著降低新突變的錯誤發現率。在上述4種方法中，PolyMutt、Tridenovo均能夠檢測單核苷酸新突變和插入/刪除新突變，但插入/刪除新突變的檢測準確率明顯低于單核苷酸新突變的檢測準確率；而DNMFilter只能夠檢測單核苷酸新突變而不能檢測插入/刪除新突變。

本文提出一種基于Adaboost的插入/刪除新突變檢測方法。該方法利用Adaboost作為分類模型；該方法對read比對之后的家系基因組數據進行局部重新拼接，并從局部拼接后的read比對上下文中提取序列特征；該方法利用已經驗證的插入/刪除新突變和隨機選取的錯誤插入/刪除新突變構建訓練集。最后本文將訓練的Adaboost模型應用于千人基因組計劃CEU家系基因組測序數據，并同常用的新突變檢測方法進行比較。

1 Adaboost算法

Adaboost算法[7]是一種經典的boosting算法，其主要思想是將一系列較弱的分類器結合起來，形成具有較強分類能力的組合分類器(強分類器)。

假設N個被標記的樣本數據集合{(X1,Y1), (X2,Y2), (X3,Y3)，…，(XN,YN)}為訓練集，Xi為特征向量，Yi為類別標記，定義Wk(i)為第K次循環樣本i的權重(K的取值范圍為1,2,3，…，Kmax，其中Kmax為最大循環次數)。Adaboost算法流程可表述如下。

首先，初始化W1(i)=1/N(其中i=1,2,3，…，N)。

然后，在每一輪循環中根據樣本權重Wk(i)，在訓練集中進行樣本重采樣生成子集Dk。

最后，在子集Dk上訓練弱分類器Ck，計算Ck在訓練集上的訓練誤差Ek，并且根據本輪的分類結果計算下一輪訓練的權重Wk+1(i)，研究推得其數學運算公式如下：

(1)

其中，Zk是歸一化系數；hk(xi)是分類器Ck對特征向量Xi分類的標記；αk的數學求值則可寫作如下形式：

(2)

使用更新后的權重Wk+1(i)重復上述過程，當循環次數達到最大循環次數Kmax或者是誤差率小于某一閾值時停止，樣本被劃分類別可使用以下公式做出判斷：

(3)

Adaboost算法可以有效降低分類中的偏差和方差[8]，而與其它boosting算法相比，Adaboost對過擬合有著更強的抗性[9]。選擇Adaboost算法可以有效降低假陽性錯誤，進而減少過濾結果中的假陽性位點的數量。

2 基于Adaboost的插入/刪除新突變過濾方法

基于Adaboost的插入/刪除新突變檢測方法總體流程如圖1所示，使用千人基因組計劃的家系基因組數據作為實驗數據。首先使用已有的新突變檢測工具檢測出候選的插入/刪除新突變位點，并對候選位點進行局部從頭拼接(local de novo assembly)；然后結合已經存在的真實突變位點，創建帶標簽的候選位點集合，根據候選位點從重組的家系基因組數據中提取特征構建訓練集；最后使用訓練集對模型進行調整優化，并且用測試集來驗證模型的效果。

圖1 插入/刪除新突變檢測流程圖

2.1 數據集

本文使用千人基因組計劃CEU家系全基因組測序數據進行模型構建與實驗結果分析，該數據集的BAM文件可以從此處：ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/working/20120117_ceu_trio_b37_decoy/下載得到。此外，該數據集包含了經過生物實驗驗證的53個插入/刪除新突變(6個插入突變，47個刪除突變)[5]，這些突變位點的位置信息可以在http://www.nature.com/nmeth/journal/v10/n10/extref/nmeth.2611-S2.xls獲得。

2.2 訓練集構建

本文構造訓練集應該滿足的條件是真實的新突變盡可能集中分布，假的新突變在空間分布上應該與新突變盡可能分離，并且在空間中盡可能分布廣泛。

本文構造正例集的具體方法是：首先將2.1節中經生物實驗驗證的53個真實新突變加入到正例集；由于正例集規模較小，無法滿足模型構建的需求，本文模擬部分插入/刪除新突變并將其加入到正例集中，此后接續的模擬過程就是選擇父本(或母本)和子代是野生型且母本(或父本)是雜合型的位點，將母本(或父本)與子代進行交換，即可模擬出滿足親代是野生型、子代是雜合型的新突變位點。

本文構造反例集的具體方法是：從除去經過生物實驗驗證的53個真實新突變的候選新突變位點集合中，隨機選取一定量的位點加入到反例集；從胚系基因組變異(Germline Variation)位點集合，隨機選取一定量的位點加入到反例集中。

2.3 序列特征選擇

研究中選定了55個特征來訓練模型，這些特征分為4個部分。就整體而言，第一部分是父本read的特征，第二部分是母本read的特征，第三部分是子代read特征，第四部分是公共特征。其中，前三部分特征都是相同的，所選用的特征見表1。

表1 插入/刪除新突變檢測選擇的序列特征描述信息

3 實驗結果與分析

3.1 Adaboost與Gradient boosting性能比較

Gradient boosting和Adaboost一樣，也是一種常用的boosting方法，兩者的主要區別是損失函數的不同，直觀上看前者是由梯度下降來決定，后者是由高權重的點來決定。Adaboost可以看作是Gradient boosting的特殊情況[10]。

為了測試比較兩者的運行效果，在同一組數據上分別使用Gradient boosting和Adaboost進行了實驗，如圖2所示，縱軸表示百分比，橫軸表示cutoff(cutoff是分類時使用的值，小于此值被分為假的插入/刪除突變，大于此值被劃分為真的插入/刪除突變，實驗中需要對其進行設置)，FP表示cutoff-假陽性率曲線，TP表示cutoff-真陽性率曲線。

在圖2所示的曲線中，GBFP在0～0.5下降比較快，ADATP在0.5以后下降比較快，但是所設定的cutoff都大于0.5(只有2個類別: 1個是真的插入/刪除突變，1個是假的插入/刪除突變)，在0.5以后的區域Adaboost的表現比較好，可以篩選掉大部分的假陽性的插入/刪除新突變。在cutoff-TP曲線上,盡管GBTP表現比較好，在0.97左右才出現了突然下降，即大量的真實插入/刪除新突變開始被錯分為假的插入刪除突變，而ADAFP在0.84左右就已經開始下降，但是考慮到設置cutoff時，為了防止原本數量就很少的真實的插入/刪除新突變被錯判為假的插入刪除突變，所以這個值不會過高，盡管Adaboost在cutoff-TP曲線上下降較早，但是和Gradient boosting一樣滿足精度要求，可以選擇設置在0.6～0.8之間。

3.2 基于Adaboost的插入/刪除新突變檢測模型評估

為了考察這些特征對于分類的重要性，本文選擇部分正例和反例，使用R語言中的adabag包對其進行訓練，然后提取各個特征對于分類的相對重要性，設計結果如圖3所示，在本組數據中，父本等位基因平衡性、父本和子代PValue值、母本等位基因平衡性、Homopolymer、父本平均mapping質量等因素對于分類的影響比較大，而子代鏈方向等因素影響較小。

圖2 基于Adaboost與基于Gradient boosting構建模型對比

Fig. 2 Comparison of model based on Adaboost and Gradient boosting

圖3 用于插入/刪除新突變檢測的各特征相對重要性

為了評價構造的訓練集的分類性能，對所選擇的訓練集進行主成分分析，最終繪制后則如圖4所示，前三個主成分之和達到了75%以上。其中，紅色的點代表了正例集的數據，藍色代表反例集的數據。可以看出正例樣本與反例樣本能夠明顯區分開，表明了構建的訓練集和選擇的分類特征能夠較好區分真實的插入/刪除新突變和假的插入/刪除新突變。

3.3 千人基因組計劃家系數據的實驗結果與分析

將本模型與其它常用的模型進行了對比，研究得到的對比后結果詳見表2。其中，Triodenovo是基于似然模型對新突變檢測的工具，相較于Polymutt的運行，性能要更勝一籌。DeNovoGear 是基于似然模型對新突變進行檢測的工具，GATK PhaseByTransmission用于矯正家系中的模糊的位點。使用4種方法對同一個家系基因組數據進行處理，各個方法使用默認的參數，檢出總數目中本文設計實現的DNINDELFilter具有最小的數目，而且DNINDELFilter檢測出了數據中所有的53個真實的插入/刪除新突變。從而達到了保留所有的真實插入/刪除新突變，保證真陽性率，同時盡可能刪除較多的假的插入/刪除新突變，降低假陽性率的預期目標。

圖4 訓練集數據主成分分析

表2 新突變檢測方法比較

4 結束語

本文將Adaboost算法應用在插入/刪除新突變篩選中，并且應用Adaboost算法，實現了DNINDELFilter。在此研究過程中，選擇了與插入/刪除新突變相關的55個特征作為訓練特征，對高通量測序數據進行了局部重新組裝，通過構建訓練集、繼而訓練、測試該模型、及對模型參數輔以一定調整，最后使得模型能夠在保留大部分的真實插入/刪除新突變同時，篩選掉絕大多數的假的突變，并且本文提出的模型具有一定的抗過擬合能力，篩選結果表現出較高的可信度。在本文的研究基礎之上，未來的工作可以在如下方面展開：應深入探討、并進一步優化模型中的特征，使得模型對插入/刪除新突變具有更好的篩選能力；目前該模型在插入/刪除新突變檢測中取得了良好的效果，后續可以嘗試將本文的模型做出研究改進，并應用到結構新突變篩選中。