K-means聚類算法在腫瘤基因變異識別中的應用

葉 驍

(復旦大學計算機科學技術學院智能信息處理重點實驗室 上海 200433)

0 引 言

二代測序技術(NGS)，又名高通量測序技術，是一項通過隨機打斷人體DNA樣本，再進行短片段測序，從而推測出樣本基因的實驗技術手段[1]。由二代測序數據產生的讀段再經過測序數據分析流程，包括讀段匹配、重復片段檢測、讀段重新匹配以及堿基質量修正等操作后進入下游的一系列基因突變分析，如單核苷酸多態性(SNP)、結構變異(SV)、拷貝數變異(CNV)和其他變異分析[2]。

之前的研究主要集中在單樣本正常細胞的測序數據分析，近年來，隨著二代測序技術的不斷發展和完善，檢測成本和錯誤率不斷下降，人們把更多的目光投向了致癌基因研究這一領域，如前身為癌癥基因圖譜(The Cancer Genome Atlas)的國際癌癥基因協會(International Cancer Genome Consortium)[3]。對于癌細胞基因突變的檢測有利于人們更多地了解癌癥基因型與表型之間的關聯，從而提高臨床診斷的效率，加速新型藥物的研發。無論是系統性錯誤如二代測序技術本身測序片段長度短、檢測堿基不精確、PCR擴增時引入重復序列，還是隨機錯誤如下游數據分析中存在多個匹配位置以及INEDL附近難以匹配精確等問題都使得突變檢測困難重重。加之癌癥樣本不同于單樣本，其中存在著癌細胞樣本純度、癌細胞倍性及多種癌細胞雜合等問題[4]。癌細胞樣本純度表現為癌細胞群中混雜有正常細胞，這使得基因型不再是傳統的AA/BB/AB三類，頻率也不再是0/0.5/1這三個離散的數值。研究表明在極低純度的癌細胞中，變異等位基因的頻率可能只有5%，在深度較低的測序數據中將無法被觀測到。癌細胞倍性是指癌細胞中有可能出現單倍體或者多倍，從而引發雜合性丟失。多種癌細胞雜合指出現多種癌細胞的亞種，亞種之間基因型存在差別等問題。上述三類問題使得各種用于正常細胞突變檢測的工具在面對癌細胞測序數據時表現得不是很好。

目前，在癌細胞突變檢測方面已經有不少用于檢測SNVs和INDEL的工具。如比較主流的應用等位基因頻率建立貝葉斯模型的Somatic-Sniper[5]和計算各種基因型后驗概率來得到某一位點最可能基因型進而判斷是否發生變異的Mutect2[6]。但這些模型普遍包含大量參數需要使用者調整，因而對于各類情景適應性不好。

本文旨在設計一種利用無監督聚類方法檢測配對樣本基因突變的算法，摒棄以往算法中各類參數在不同場景中調節困難的缺點，同時仍能達到較高的精確率和召回率。

1 K-means聚類算法簡介

K-means算法的輸入是聚類個數k，大小為n的樣本集，輸出為滿足方差最小的k個樣本集中點聚成的類，滿足同一聚類中的樣本相似度較高，而不同聚類中的樣本相似度較小[7]。算法步驟為：

(1) 從維度空間任意選擇k個點作為初始聚類中心。

(2) 根據每個聚類樣本的特征值，計算每個樣本與這些中心點的距離，并根據最小距離重新劃分樣本點所屬類。

(3) 重新計算每個類的均值作為新的類中心。

(4) 計算測度函數，當滿足一定條件，如函數收斂時，則算法終止；如果條件不滿足則回到步驟(2)。

2 算法設計

2.1 工作流程

保存著讀段信息的fastq文件(fastq是一種存儲了生物序列以及相應的質量評價的文本格式)經二代測序標準流程處理后得到bam格式的文件。我們利用samtools[8]提供的讀段堆砌功能，分別將tumor.bam、normal.bam轉化為pileup格式的文件(pileup格式文件相當于把每條染色體都豎起來，將每條reads也豎起來平行的匹配到基因組上)tumor.pu、normal.pu。pileup文件格式提供了單個位點所處染色體、染色體上位置、參考基因上的堿基型、堆砌到此位點上堿基的集合、堆砌堿基質量的集合、堆砌堿基所處讀段質量的集合等信息。我們先根據同時出現在癌細胞與正常細胞中這一規則篩選出候選位點，再利用pileup文件信息抽取7個特征，并進行標準化處理。由于正負樣本數量相差巨大，在數據過濾步驟中我們選擇兩個最為重要的特征FISHER_PVAL_T_VS_N和VAF_T_VS_N根據分位數過濾掉大部分數據，再將剩下的位點用K-means模型聚類。關于FISHER_PVAL_T_VS_N和VAF_T_VS_N我們會在2.3節中詳細說明。通過對比聚類結果中兩個類中心的特征值選擇正負樣本。變異識別算法的工作流程如圖1所示。

2.2 使用堆砌算法生成候選變異集

讀段堆砌(pileup)的意思是，對于基因組上某一個位點，將所有覆蓋了該位點的讀段在該點處進行對齊，有利于我們的算法計算模型所需要的相關信息，更進一步地判斷該點是否發生了變異，發生的是何種變異。

言不在重，貴在實?！捌秸Z近人”，《摘編》用極其凝練、高度概括的話語提綱挈領、大開大合?！墩帯罚ǘ┲刑岬健爸赂徊恢赂?，關鍵看干部”，一句話強調了人的重要性，指出了脫貧攻堅要注重干部隊伍和人才隊伍的建設?！啊鏊乓妰赏饶唷?。扶貧工作必須務實，脫貧過程必須扎實，扶真貧、真扶貧，脫貧結果必須真實……”《摘編》（六）這段話既拉近了與群眾的距離，又開門見山直奔主題，強調了實事求是的思想路線和務實扎實的工作作風。

堆砌過程中將逐個遍歷全基因組或者基因組上指定區域的所有位點，在對于單核苷酸突變的研究中，每個位點上可能出現4種變異。舉例來說，以‘A’腺嘌呤脫氧核糖核苷酸作為參考堿基，則其可能突變為‘C’胞嘧啶脫氧核糖核苷酸，‘G’鳥嘌呤脫氧核糖核苷酸或‘T’胸腺嘧啶脫氧核糖核苷酸3中。理論上還有可能發生INDEL——短基因片段的刪除或插入，CNV基因片段拷貝數變異等變異，由于該模型僅針對單堿基突變，故不做研究。對于每一個同時出現在癌細胞和正常細胞中的基因位點，我們都將利用讀段堆砌的信息進行特征計算。堆疊的可視化展示如圖2所示。

圖2 堆疊的可視化展示

2.3 變異特征

在查閱相關文獻和實驗驗證后，我們總結出7個在檢測腫瘤突變上有意義的特征，對于所有位點，分別提取以下7個特征：

(1) VAF_T癌細胞中等位基因頻率。癌細胞等位基因頻率是衡量該點堿基是否發生突變的一項重要指標，同正常細胞中發生變異的基因頻率不同(通常只有0、0.5、1三個數值)，癌細胞中的基因頻率可能為0～1間的任意值。但是其分布仍然呈現一定的規律，因此僅考慮癌細胞中的等位基因頻率也是有意義的。

(2) VAF_T_VS_N癌細胞和正常細胞等位基因頻率之比。前文中提到，配對樣本突變檢測目的就是要找出只出現在癌細胞而不出現在正常細胞中的突變，將AGCT 4種堿基在癌和正常細胞中的頻率相除可以很明顯地看出頻率上的差別，我們將它列為一項重要的特征。

(3) MEAN_MQVQ_T癌細胞位點上讀段匹配質量均值。讀段匹配質量均值代表該堿基所在讀段(READS)的匹配好壞。如果匹配質量較差，則很可能是一個誤匹配位點，那么后續的計算過程也是無效的。

(4) MEAN_BQVQ_T癌細胞位點上堿基質量均值。堿基質量代表從測序儀上讀到該堿基的正確率，和匹配質量相同，我們需要把這一項也納入聚類模型的特征中。

(5) FR_RATIO_T癌細胞中變異基因正反鏈數目之比。之前的研究表明，如果變異基因只位于一條鏈(正鏈或反鏈)上，那么很有可能是實驗帶來的誤差。反之，如果該基因在正反鏈上分布較均等，那么就很有可能是變異基因。

(6) FISHER_PVAL_T_VS_N癌細胞和正常細胞中的變異。堿基數的Fisher精確檢驗P值,Fisher精確檢驗是一種假設檢驗，通常用來檢驗一次隨機實驗的結果是否支持對于某個隨機實驗的假設。

將某一位點上正常細胞和癌細胞里出現的參考堿基數和變異堿基數列成2×2的列聯表，如表1所示，進行Fisher精確檢驗，來判定兩類細胞中的變異堿基數是否存在顯著差異。這是許多傳統算法采用的一個重要標準，我們也將它列為一項特征。

表1 Fisher精確檢驗的列聯表

(7) BQP癌細胞中變異等位基因與參考基因堿基質量的秩和檢驗P值。P值越低代表了alt基因和ref基因的質量差異越大，alt也就越有可能是測序誤讀。

上述特征在以往傳統的貝葉斯模型算法中也常常作為重要的參數，因此我們有理由認為這些特征是有意義的。

3 實驗分析

3.1 實驗數據集

模擬實驗數據中，腫瘤樣本通過隨機采樣人類樣本基因NA12981的序列替換NA12878的序列片段生成，被替換的堿基也就是發生了體細胞突變的堿基，而原來NA12878中的突變則為生殖細胞突變。由此我們得到了一批可以用于評估算法性能的模擬數據。

真實實驗數據來源于2017年全國腫瘤體細胞基因突變高通量測序檢測生物信息學分析室間質量評價活動所提供Illumina HiSeq2500高通量測序平臺的原始靶向測序數據。原始數據為fastq文件格式，分為腫瘤細胞和正常細胞兩組實驗數據。經標準二代測序分析流程處理后得到bam格式的文件。讀段匹配過程使用bwa進行匹配，隨后對產生的bam文件應用GATK[8]工具包做了標準的預處理，包括標記重復讀段、局部重匹配和堿基質量再校正等。將本文算法與目前主流的3款腫瘤突變檢測工具進行了對比，包括Mutect2、Somatic-Sniper和JointSNVMix[10]。三者使用的都是官方默認的參數配置或是官方推薦的“最佳實踐”方案。

3.2 評價標準

分別選擇Mutect2的Tlod，Somatic-Sniper的SSC，JointSNVMix的Prob和K-means中距離類中心的距離作為算法的輸出分數，據此生成PR曲線，并在Tlod>30,SSC>50,Prob>0.9時計算出各算法的Precision、Recall和F-Score。

(1)

(2)

(3)

PR曲線以Precison作為縱軸，以Recall作為橫軸，設定不同的閾值將算法輸出分數映射為正負樣本，整張圖上曲線越向右上方凸往往算法性能越優。

3.3 實驗結果

4種算法在模擬數據上的結果如圖3所示。可以看出，在模擬數據上本文算法在準確率和召回率上均占優勢。

圖3 模擬數據PR曲線

在有標簽的IlluminaB17數據上的變異識別結果如圖4、圖5和表2所示?？梢钥闯?，采用無監督聚類的算法較以往傳統算法可以檢測出更多的變異位點，同時保證較高的查準率。

圖5 IlluminaB17數據PR曲線

指標工具Mutect2Somatic-SniperJointSNVMixK-meansPrecision0.469 40.732 10.666 70.890 4Recall0.302 60.539 50.184 20.855 3F-Score0.368 00.621 20.288 60.872 4

在沒有標簽的肺癌患者肺部細胞數據上,如圖6和表3所示?？梢园l現，用本文方法找出的位點與其他3種方法的二者交集2-cross、三者交集3-cross均表現較好。K-means與其他3種方法的交集數量都比較大，說明該算法可以找出以往算法所認同的大部分變異。

圖6 肺癌數據上算法交集數量對比

表3 肺癌數據算法2-cross數量對比

4 結 語

本文提出了一種利用無監督聚類來識別單核苷酸變異的腫瘤突變識別算法。不同于以往算法一貫采用貝葉斯模型，我們通過結合高效產生候選變異集合的讀段堆砌算法與無監督的聚類算法，自動而準確地識別變異位點，得以同時實現召回率、精度和速度三大指標上的良好表現。將其及另外三個變異識別器應用于來深度測序樣本數據，其結果證明了本文算法具有良好的效果。

本文算法中使用的特征均為標準化后的特征，在實際使用中由于數據集的差異，不同特征的重要性存在差異，為每個特征選擇適當的權重或可明顯提升算法的性能。權重的選擇既可以基于經驗估計，也可以迭代取值進行實驗求得最優的參數。