李曉東，何小雨，陳瑋，李瑞琳，趙丹，祝海棟，張裕，代闖闖，陸忠華，遲學斌,3，牛北方,4*，郎顯宇*

1. 中國科學院計算機網絡信息中心，北京 100190

2. 中國科學院大學，北京 100049

3. 中國科學院計算科學應用中心，北京 100190

4. 貴州大學醫(yī)學院，貴州 貴陽 550025

引言

在對腫瘤長久以來的研究中，研究人員逐漸將它定義為一種“基因組疾病”[1]。腫瘤起源于某個上皮細胞，當它積累了足夠的體細胞突變后，會導致其自身的增殖失調，從而擴展成一個大的細胞群，細胞群中的子代細胞也保留有與祖先細胞一致的突變。能夠誘發(fā)細胞發(fā)生腫瘤變異以及促進腫瘤發(fā)展的突變可以分為兩種：抑制基因 (suppressor genes) 的失活性突變和癌基因 (oncogenes) 的激活性突變。兩種變異都可能導致攜帶它們的細胞獲得顯著的成長優(yōu)勢，從而迅速發(fā)展成腫瘤組織。

腫瘤基因組研究期望尋找到變異基因與臨床表現之間的聯系，從而識別出驅動腫瘤進化發(fā)展的變異，為病人提供有效的精準治療方案。研究發(fā)現[2]，腫瘤基因組變異主要有單核苷酸變異(Single Nucleotide Variant，SNV)，插入和缺失變異(Insertion and Deletion，Indel)、結構變異 (Structure Variant，SV) 三種類型。其中腫瘤體細胞 SNV 是一種最簡單的變異類型，在腫瘤的發(fā)生中扮演重要的角色[3]。SNV 可能導致對應的氨基酸改變，甚至所編碼蛋白質密碼子變成終止密碼子，從而影響基因的功能。通過高通量測序技術，研究人員能夠直接獲得腫瘤和正常組織樣本的原始基因序列數據。通過分析樣本基因序列，并與標準參考基因組比對，尋找與腫瘤相關的 SNV，是目前腫瘤基因研究中常用的方法。

在腫瘤序列分析中，體細胞 SNV 檢測是后續(xù)分析的基礎。尋找到變異位點集合后，通過比對已有數據庫 (如 dbSNP，COSMIC) 中的變異位點，從而尋找與疾病相關的變異基因。因此，最終分析結果的可靠性依賴于尋找到的初始變異位點集合的質量。

然而，樣品降解、遺傳異質性和組織污染 (雜質)等因素影響，對 SNV 的檢測存在較大挑戰(zhàn)。同時，測序過程中引入的測序誤差、序列重排引入的比對誤差以及堿基覆蓋度不足等，也給變異位點的檢測增加了困難。早期的一些檢測軟件如 NaiveSubtract[4]，采用統(tǒng)計學顯著性檢驗方法，分別對配對的腫瘤樣本和正常樣本進行分析，獲得樣本的 SNV 候選位點集，并從腫瘤樣本的候選位點集中去除在正常樣本中檢測到的候選位點。但是，這種簡單的“減法”沒有對兩個樣本中共有的家系多態(tài)性位點進行分析，同時也不適用于發(fā)現低等位基因頻率的變異。目前，針對體細胞突變檢測的算法，均是通過聯合分析腫瘤和正常樣本序列，從而最大程度減小因測序平臺及比對算法引入的隨機或系統(tǒng)誤差。

在對變異位點檢測時，研究人員希望檢測到的候選位點集具有較高的準確性，同時假陽性位點數目盡可能低。由于不同方法使用的模型和評估指標不同，不同軟件識別的變異位點集合存在較大的不一致性，方法準確性也存在較大差別。本文針對目前最受歡迎的四種 SNV 檢測算法 Varscan2[5]、SomaticSniper[6]、Strelka[7]和 MuTect2[8]進行測試，在已知基準 SNV突變信息的數據集上，比較四種方法的表現，同時對不同方法在不同測序深度上的結果進行分析。

1 體細胞 SNV 檢測分析

DNA 序列包含腺嘌呤核苷酸 (A)、鳥嘌呤核苷酸(G)、胞嘧啶核苷酸 (C)、胸腺嘧啶核苷酸 (T) 四種核苷酸。SNV 主要發(fā)生在 DNA 復制過程中，即某種核苷酸突變成為另外一種類型核苷酸，從而導致同一個位點上的基因存在兩個以上的狀態(tài)，稱作等位基因。在二倍體基因組中，對每一個位點，相較于參考基因組，設所有可能的等位基因型組合為 Ω={AA，AB，BB} A 表示該位點基因型與參考基因組一致，B 表示該位點基因型與參考基因組不一致。AA 為純合子參考型，AB 為雜合子突變型，BB 為純合子突變型(在人的基因組中，理論上存在 10 種組合，即 {AA，AC，AT，AG，CC，CT，CG，TT，TG，GG}，此處簡化，每個位點只包含兩種可能的基因型 )。對每一個位點，腫瘤樣本和正常樣本中聯合基因型分布如表 1所示。

1.1 變異檢測算法

本文分析的四種變異檢測算法，主要使用經過質控、比對、去重、重校驗后的 BAM 類型文件作為輸入數據。BAM 文件是 SAM 文件的二進制格式，文件中包含每個 read 片段的原始序列，堿基質量，比對到參考序列的坐標信息，比對質量等信息。算法通過聯合分析參考基因組、腫瘤和正常組織樣本基因組，計算各位點的相關特征信息，采用不同模型，確定每個位點在不同樣本中基因型，從而確定在腫瘤樣本中該位點是否為體細胞單核苷酸變異位點 (圖 1)。

1.1.1 Mutect2

Mutect2 算法基于隱馬爾可夫模型的貝葉斯方法，設計了兩個貝葉斯分類器。對變異位點的檢測分為四步：(1) 去除低質量序列數據；(2) 使用第一個貝葉斯分類器對腫瘤樣本進行變異檢測；(3) 過濾由于相關測序過程中產生的假陽性誤差，這類誤差很難被誤差模型捕獲；(4) 使用第二個貝葉斯分類器，對識別的變異位點進行區(qū)分，判斷為體細胞變異 (SNV) 或家系變異 (Germline variant)。

1.1.2 SomaticSniper

SomaticSniper 算法使用貝葉斯概率模型，基于參考序列，計算每種基因型的先驗概率分布。通過每個位點在腫瘤和正常樣本中的觀測數據，計算位點聯合基因型的后驗概率，采用對數運算轉換為每個位點的體細胞變異得分。轉換后的分值在 0-255 之間，表示位點的變異程度，用戶可指定相應的體細胞變異閾值，以過濾變異分數較低的候選位點。

圖1 假設的腫瘤-正常樣本比對示意圖，等位基因數(Allelic Counts)， 表示該位點與參考基因一致的 read 數，表示覆蓋該位點的總 read 數Fig. 1 The comparison diagram of hypothesized tumor-normal sample

1.1.3 Strelka

Strelka 首先搜索正常和腫瘤樣本 BAM 文件中發(fā)生插入缺失的序列，對其重排。算法使用重排后的序列數據，基于貝葉斯概率模型，計算出候選位點在正常樣本中最可能的基因型和正常基因型狀態(tài)下體細胞變異概率。模型中使用到的序列信息包括堿基質量和比對質量，鏈偏差，位點為變異位點的先驗概率以及正常樣本中雜合型的預期率。由于測序的正常樣本中一般包含有家系變異，腫瘤樣本包含真正體細胞變異的腫瘤組織和正常組織，因此 Strelka 使用等位基因頻率而不是二倍體基因型計算體細胞變異概率。對尋找到的原始 SNV 和 Indel 變異位點，Strelka 依據配置參數，以剔除假陽性位點。

1.1.4 Varscan2

Varscan2 采用啟發(fā)式和統(tǒng)計學算法，分別生成包含 SNV 變異和 Indel 變異的 VCF 文件。相較于其他算法，Varscan2 使用經 samtools 程序基于 BAM 文件生成的 mpileup 類型文件作為輸入[9]。在每一個位點，程序同時讀取腫瘤和正常樣本中的數據后，對位點的深度進行標準化，然后進行啟發(fā)式配對比較。當一個位點存在某種堿基出現頻率超過閾值，該位點基因型為純合型，否則為雜合型。對在腫瘤和正常組織樣本中基因型不匹配的位點，統(tǒng)計兩個樣本中參考支持型 (reference-supporting) 和變異支持型 (variantsupporting) 的序列數，進行單尾 Fisher 精確檢驗。如果檢驗的 P 值小于設定的顯著性閾值，算法依據正常樣本的基因型來確定該位點的變異類型，當正常樣本為純合型，分類為體細胞變異 (SNV)，當正常樣本為雜合性，分類為雜合性缺失 (LOH)。如果大于顯著性閾值，則為家系變異 (Germline variant)。在生成的VCF 文件中，標注每個位點的最大概率變異類型。

1.2 數據特征說明

上述四種變異檢測算法均是通過統(tǒng)計待分析位置的數據信息，選擇可用于分析的有效特征。不同的變異檢測算法使用的數據特征各有不同，根據數據特征的特點，可分為以下六種類別，如表 2 所示。

(1) 序列測序深度：正常樣本和腫瘤樣本中覆蓋某一位點的短 read 數；

(2) 堿基質量：在測序時，每測一個堿基會給出一個相應的質量值，用于衡量測序準確度，通常使用的 Phred 堿基質量值公式為，Q-score=-10×log10P，P為預估的堿基錯誤檢查率；

(3) 鏈偏差：采用雙端測序時，由于二代測序的讀長限制，某些位置的堿基可能只在正向鏈或反向鏈上被測到；

(4) 序列比對質量：將樣本序列比對到參考基因組時，由比對軟件生成，用于衡量比對準確度；

(5) 位點位置值：位點位于所在序列的位置信息。位于所在序列中間時值為 1，位于所在序列兩端值為0。越接近序列兩端，由測序和比對過程產生誤差的概率越大；

(6) 先驗概率與檢驗水平：在對腫瘤和樣本序列的比較中，檢測方法采用統(tǒng)計假設檢驗，在用戶指定的顯著水平下，判斷檢測位點的基因型是否匹配或位點是否為體細胞變異位點。

1.3 候選 SNV 位點過濾

控制假陽性率是體細胞 SNV 變異檢測方法面臨的主要難點。理想狀態(tài)下，變異識別算法應具有較高的敏感性，在最大程度發(fā)現真實變異的同時，保證檢測到的變異為假陽性變異的概率最低。上述四種檢測算法使用了多種信息特征來計算每個候選位點的體細胞變異概率。然而，由于 DNA 序列數據中的復雜因素影響，每個算法使用的特征因素并不完全相同，初始得到的候選變異位點集合中存在大量的假陽性變異位點。因此需要對候選位點進行過濾，以減少假陽性。

Varscan2 和 SomaticSniper 需要結合后續(xù)的流程，進行過濾。對初始得到的候選位點，提取位點的坐標信息，采用 bamreadcount[10]計算每個單核苷酸位點的詳細信息，剔除不滿足指定條件閾值的候選位點。Strelka 和 Mutect2 算法內置過濾流程。用戶在算法運行前，需確定相關閾值參數，算法分別輸出過濾前的變異位點集和過濾后的高可信度變異位點集。

2 實驗數據

對于突變檢測類方法性能的衡量，困難在于，沒有一個已知完全變異信息的標準數據集。因此，無法對變異識別方法準確率和召回率進行準確評估。為解決這個問題，研究人員嘗試使用兩種途徑構造可以作為衡量標準的數據集 (表 3)。

表2 四種變異檢測算法使用的特征明細Table 2 The detailed characteristics of the four detection methods

途徑 1：使用多測序平臺對數據進行多次測序。對同一個樣本，使用不同建庫方法和不同測序平臺進行測序。通過整合多個測序結果，得到高置信度的原始序列文件。美國國家標準和技術研究院 (NIST)，對人類基因組計劃中的 NA12787 (NIST RM 8398) 樣本，采用多種測序平臺和建庫方法，對樣本進行測序，獲得不同種類的數據[11]。通過不同數據間相互驗證、補充，從而獲得該樣本的各類型變異信息。該數據通?？勺鳛榧蚁底儺?(Germline variant)，插入缺失變異 (Indel) 和結構變異 (SV) 檢測及全基因組組裝的標準數據。

途徑 2：使用序列數據模型軟件，人工合成相應的變異數據。目前采用的對變異數據的模擬有兩種方法：(1) 基于參考基因組，模擬讀長和變異位點配置信息，合成突變基因組；(2) 基于已有的經過比對的正常序列數據，在其中一些混入預先設定的變異，設定的變異包括隨機變異和已發(fā)現的相關疾病變異，以模擬腫瘤組織數據。通過上述方法，研究人員能夠控制樣本的準確變異信息，以對不同算法的性能進行有效評估。由于第一種方法無法模擬真實數據中變異位點的非隨機分布，非獨立誤差和拷貝數變異等特性，目前分析上常采用第二種方法。

本文使用的測試數據來自由國際癌癥基因組協會 (ICGC) 和癌癥基因組圖譜 (TCGA) 癌癥基因組計劃聯合舉辦的變異識別挑戰(zhàn)項目[12]。正常組織的基因序列數據是通過 Illumina HiSeq 2000 測序儀產生的雙端 (paired-end) 測序數據，平均讀長約為 101bp，平均覆蓋度約為 60X。使用 GRCh37/hg19 基因組作為參考基因組，采用 BWA[13]軟件進行序列比對，通過bamSurgeons[14]軟件模擬產生對應的腫瘤基因序列數據。數據包含確定的 SNV 變異信息，可對每個檢測算法的結果進行有效評估。

表3 實驗數據描述Table 3 Experimental data description

表4 不同變異識別方法識別的相同候選位點比例Table 4 The proportion of the same candidate sites identi fi ed by different methods

3 實驗結果及分析

在對上述四種方法進行測試過程中，Mutect2 和strelka 使用算法默認的參數。SomaticSniper 設定堿基質量參數，對初始識別的候選位點，采用程序開發(fā)者推薦的過濾方法進行過濾，獲得高置信度變異位點集。Varscan2 使用基于算法默認參數，并依據數據集樣本純度，對相應參數進行調整，獲得最佳的高置信度候選位點集。詳細算法運行參數參見補充材料 1。

3.1 不同方法檢測的變異位點數量存在顯著差異

在三組正常-腫瘤全基因組樣本對上的檢測結果顯示，四種檢測算法識別的體細胞突變數量存在顯著差異。其中，strelka 檢測到的候選位點數目最少，Mutect2 和 SomaticSniper 檢測到的候選位點數目較多。算法檢測變異位點數量與樣本真實變異位點數量呈正相關 (表 4)。

為確定不同算法檢測結果的一致性，我們對每個樣本的四個候選集交叉比較，計算同時被兩種算法檢測的位點比例。如圖 2 所示，四種算法中，Mutect2、SomaticSniper 與其他算法的一致性較低，Strelka 的一致性最高。Strelka 檢測到的位點，多數也被其他方法檢測到。但同時，其他方法檢測到的位點也被 Strelka 檢測到的整體比例最低。由此說明，四種算法中 Strelka 的檢測標準最為嚴苛，檢測的候選位點數量也最少。

3.2 四種檢測算法的性能評估

通過與各樣本的真實變異位點比較，我們對各個算法的準確率、召回率進行衡量， 比較檢測算法的綜合性能。四種檢測方法中，Varscan2 和 Strelka 的綜合性能優(yōu)于另外兩種算法，SomaticSniper 的性能最低。Mutect2 損失一定的準確度，但能夠檢測到最多的真實變異。受樣本純度的影響，識別算法在 IS2 數據集的整體準確率表現相對數據集 IS1 和 IS3 更低，此時 Strelka 的準確率顯著高于其他幾種算法 (表 5)。

圖2 四種變異檢測方法檢測到的候選位對比Fig. 2 The comparison of candidate positions number detected by four different methods

表5 四種檢測算法在各個樣本集上的檢測性能評估Table 5 Four detection algorithms are evaluated on each sample set

圖3 四種檢測算法的交集位點數目和準確率Fig. 3 The intersection number and accuracy of the four detection algorithms

圖4 四種檢測算法的并集位點數目和召回率Fig. 4 The recall rate and convergence number of the four algorithm

圖5 三組樣本中真實變異位點的測序深度分布，橫軸表示位點的測序深度Fig. 5 Sequencing depth distribution of real variation sites in three groups of samples

圖6 在三組樣本上四種檢測方法的結果對比，虛線表示當前樣本低測序深度部分真實變異數，藍色折線表示檢測算法檢測到的低測序深度個數，黃色柱圖表示各算法檢測結果與真實變異一致的位點數Fig. 6 The detection effect of four detection methods at low-sequencing depth site

對四種檢測方法得到的候選位點進行統(tǒng)計發(fā)現，被所有方法均檢測到的位點數分別為 3563，4870，5768，交集中真實變異位點的準確率顯著高于單個方法 (圖 3)。四種方法的并集數分別為 8289，11511，11921，并集的召回率顯著高于任意單一方法 (圖 4)。

表6 四種檢測算法對低測序深度變異位點的檢測性能比較Table 6 The performance of different algorithms for low sequencing depth

3.3 測序深度對檢測算法準確率的影響

在各個算法模型中，序列測序深度是所有檢測算法均衡量的一個因素。提取三組樣本中真實位點分別在對應腫瘤和正常樣本中的測序深度數，在顯著水平0.05 條件下，采用 Fisher 檢驗，真實變異位點在腫瘤和正常組織中的測序深度分布無顯著差異。

在三組基準變異數據集中，我們選取在腫瘤和正常樣本中測序深度均小于平均測序深度 50% 的真實變異位點，作為低測序深度位點，其個數分別為 119，120，400，占總變異位點的比例分別為3.33%，2.77%，5.06%。分析四種檢測算法在低測序深度部分的檢測結果，與真實變異位點比對發(fā)現，SomaticSniper 檢測到的總位點數目最多，但也引入了最多的假陽性位點，Strelka 檢測到的數目最少。Mutect2 檢測到的真實低測序深度變異位點數目最多，同時引入的假陽性位點也較少 (圖 5-6)。

對低測序深度部分檢測算法的性能評估，各算法的準確率和召回率遠小于其平均水平。分析結果表明，當測序深度較低時，變異頻率較低的等位基因很難被識別到，從而對位點基因型的確定存在較大的干擾，影響算法的準確性。四種算法中，Mutect2 對低深度部分位點的檢測具有最好的準確率和召回率。當樣本數據的測序深度較低時，使用 Mutect2 會獲得較好的檢測結果 (表 6)。

4 結束語

通過對腫瘤-正常樣本基因序列分析，發(fā)現準確的體細胞點突變位點對腫瘤的研究和癌癥的個性化治療至關重要。在本文中，我們驗證了四種變異檢測方法在三組模擬的正常-腫瘤全基因組序列上的突變識別效果，從不同維度對檢測方法的性能進行了評估。盡管四種方法均檢測到大量相同的 SNV 位點，尤其是真實的變異位點，但通過對三組基準數據的總體分析表明，VarScan2 在檢測高質量 SNV方面表現最為出色，在準確率和召回率上較其他方法更為均衡；Strelka 對變異位點的識別標準最為嚴苛，故能保證較高的準確性，但會遺漏較多真實的變異位點；Mutect2 能夠檢測到最多的真實變異，但同時也會引入較多的假陽性位點，這對后續(xù)變異位點的驗證工作造成較大的干擾；SomaticSniper 與Mutect2 特點相似，但準確性上表現更低。當樣本純度較低時，Strelka 具有更好的表現。此外，樣本數據的測序深度會顯著影響各算法的檢測效果。針對低測序深度樣本，我們更推薦使用 Mutect2 進行變異位點檢測。研究人員可根據研究目的，結合各算法的不同特點，采用多種方法的組合，能夠獲得更多或更準確的候選位點集。

[1]Strausberg R L, Simpson A J G, Old L J, et al. Oncogenomics and the development of new cancer therapies[J].Nature,2004, 429(6990): 469-474

[2]Strausberg R L, Simpson A J G.Whole-genome cancer analysis as an approach to deeper understanding of tumour biology[J]. British journal of cancer, 2010,102(2):243-248

[3]Carlson B.SNPs-A shortcut to personalized medicine[J].Genetic Engineering & Biotechnology News, 2008,28(12): 12-12.

[4]DePristo M A, Banks E, Poplin R, et al.A framework for variation discovery and genotyping using next-generation DNA sequencing data[J].Nature genetics, 2011, 43(5):491-498.

[5]Koboldt D C, Zhang Q, Larson D E, et al.VarScan 2:somatic mutation and copy number alteration discovery in cancer by exome sequencing[J].Genome research, 2012,22(3): 568-576.

[6]Saunders C T, Wong W S W, Swamy S, et al.Strelka:accurate somatic small-variant calling from sequenced tumor–normal sample pairs[J].Bioinformatics, 2012,28(14): 1811-1817.

[7]Cibulskis K, Lawrence M S, Carter S L, et al. Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples[J].Nature biotechnology,2013, 31(3): 213-219.

[8]Larson D E, Harris C C, Chen K, et al.SomaticSniper:identification of somatic point mutations in whole genome sequencing data[J].Bioinformatics, 2011, 28(3):311-317.

[9]Li H.A statistical framework for SNP calling,mutation discovery,association mapping and population genetical parameter estimation from sequencing data[J].Bioinformatics,2011, 27(21): 2987-2993.

[10]Li H.Toward better understanding of artifacts in variant calling from high-coverage samples[J].Bioinformatics,2014, 30(20): 2843-2851.

[11]Zook J M, Catoe D, McDaniel J, et al. Extensive sequencing of seven human genomes to characterize benchmark reference materials[J].Scientific Data, 2016,3(25).

[12]Ewing A D, Houlahan K E, Hu Y, et al.Combining tumor genome simulation with crowdsourcing to benchmark somatic single-nucleotide-variant detection[J].Nature methods, 2015, 12(7): 623-630.

[13]Li H, Durbin R.Fast and accurate long-read alignment with Burrows–Wheeler transform[J].Bioinformatics, 2010,26(5): 589-595.

[14]Xu H, DiCarlo J, Satya R V, et al.Comparison of somatic mutation calling methods in amplicon and whole exome sequence data[J].BMC genomics, 2014,15(1): 244-253.

[15]Abeel T, Helleputte T, Van de Peer Y, et al.Robust biomarker identification for cancer diagnosis with ensemble feature selection methods[J].Bioinformatics,2009, 26(3): 392-398.

[16]Alioto T S, Buchhalter I, Derdak S, et al.A comprehensive assessment of somatic mutation detection in cancer using whole-genome sequencing[J].Nature communications,2015, 6.

[17]Guo Y, Li J, Li C I, et al. The effect of strand bias in Illumina short-read sequencing data[J].BMC genomics,2012, 13(1): 666-676.