基于高斯混合模型的腫瘤純度估計

閆占正，李玉雙

（燕山大學理學院，河北秦皇島066004）

高通量測序技術的快速發展為癌癥基因組學研究帶來了重大突破。在腫瘤與癌癥研究中，在臨床環境下獲得的腫瘤固體組織中混合有正常細胞（即腫瘤不純），使得篩選出來的癌癥數據存在偏差，如果后續不經處理直接使用，往往會導致研究結果不理想甚至得到錯誤的結果[1]。 因此，進行腫瘤純度的估計，即推測腫瘤組織中癌細胞的比例，已成為癌癥研究中十分重要的課題[2-3]，并已取得重要進展。如，ABSOLUTE 方法，利用單核苷酸多態性陣列數據(SNP array)，借助最大似然估計計算腫瘤純度[2]；InfiniumPurify 模 型[4-6]，利 用DNA 甲 基 化 芯 片 數據[7]，將位點劃分為超（次）甲基化位點，構建線性函數來估計腫瘤純度；PurityEst，利用二代測序數據，通過對腫瘤樣本中突變等位基因的比例建模來估計腫瘤純度[8]。

與其他類型數據相比，DNA 甲基化數據能更穩定和靈活地顯示腫瘤細胞和正常細胞之間的內在差異。觀察正常樣本，發現其甲基化水平分布與高斯分布相似，而有限高斯混合模型可以逼近任意概率分布密度函數[9]。因此，本文利用DNA 甲基化數據，借助高斯混合模型[10-11](Gmm, Gaussian mixture model)估計腫瘤純度，得到的9 種類型的腫瘤純度與ABSOLUTE 和InfiniumPurify 的結果高度一致。之所以與這2 種方法比較，是因為ABSOLUTE 是最早且最具影響力之一的純度估計工具，常被用作評估的金標準，而InfiniumPurify 與GmmPurify 使用的數據相同，所求結果更具有可比性。

1 數據及其處理

癌癥和腫瘤基因圖譜TCGA 數據庫(https://portal.gdc.cancer.gov/)中包含32 種腫瘤類型，每種類型都含有個數不同的腫瘤樣本。包含正常樣本的腫瘤類型有22 種，其中擁有2 個以上匹配的正常樣本有21 種，多形性膠質母細胞瘤(GBM)只有1 個。由于多數腫瘤類型無（或只有很少）匹配的正常樣本，本文選用文獻[5]中構造的公共正常樣本集進行操作，即在21 種腫瘤類型中每一個類型隨機抽取2個正常樣本，并與GBM 中的1 個正常樣本一起組成包含43 個樣本的公共正常樣本集。每個腫瘤或正常樣本均包含相同的485 577 個甲基化位點，其中大約9 萬個位點的β 值全部缺失，有大約1 萬個位點的β 值部分缺失（位點的β 值表示該位點的甲基化水平，0 ≤β ≤1，其中0 表示完全未甲基化，1 表示完全甲基化）。為避免缺失值帶來的影響，將存在缺失值的位點全部舍棄。以包含466 個腫瘤樣本和43 個公共正常樣本的肺腺癌(lung adenocarcinoma，LUAD)為例，執行刪除操作后剩余371 265 個位點。

2 GmmPurify 方法

2.1 高斯混合模型及應用

高斯混合模型[9]是指具有以下形式的概率分布模型：

這里μk是均值，σk是方差，αk是系數，φ(y|θk)是高斯分布密度，稱其為第k 個高斯模型。每個高斯模型都是一個聚類。

當確定樣本所屬高斯模型時，可以利用最大似然概率算法來求解模型參數。設樣本總數為N，樣本集合為X={x1, x2, ..., xN}。設K 個高斯分布的樣本數量分別為N1, N2, ..., NK, 第k (k∈{1,2,...,K})個高斯分布的樣本集合為L(k)。則有

當無法確定樣本所屬模型時，可用期望最大化算法通過迭代進行求解，算法分2 步: E 步，根據參數θ(n)求 得 后 驗 概 率P(xt|X;θ)；M 步，更 新xt的 期 望E(xt)，然 后 用E(xt)代 替xt求 出 新 的 模 型 參 數θ(n+1)。如此迭代直到θ 趨于穩定。

假設模型參數已知，隨機初始化一組參數θ(0)，可求得后驗概率P(xt|X;θ)。求第k 個高斯分布，分別取x1, x2, ..., xN的概率，即求樣本xt在第k 個高斯分布下的概率γ(t,k)，有公式：

由式(4)～(7)，可以推出：

隨機初始化一組合適參數(α(0)k, μ(0)k, σ(0)k)，k∈{1,2,...,K}，將樣本集合X 和初始化參數代入式(7)～(10)，進行迭代直到參數趨于穩定，可求得樣本集合X 的高斯混合模型參數(αk, μk, σk)。

將高斯混合模型應用于公共正常樣本集，選取不同的K 值進行實驗。實驗需確保各類高斯模型存在的合理性以及各類模型之間明顯的差異性。為此，制定以下實驗條件：最小的分類概率不小于0.1，且各類之間的均值之差不小于0.1。研究發現，當K≤3 時，滿足實驗條件的甲基化位點占絕大多數。以LUAD 為例，其位點總數約37 萬。當K=2 時，滿足條件的位點約占29%；當K=3 時，滿足條件的位點約占5%；當K=4 時，滿足條件的位點約占0.4%。故設定K=3, N=43。以肺腺癌(LUAD)為例，依次處理371 265 個樣本集合Xi={βi,1, βi,2, ...,βi,N}, i∈{1,2,...,371 265}, 其中βi,t表示第i 個位點在第t 個公共正常樣本下的甲基化水平值。具體為：初始化參數αi,1=1/3，μi,1= min(Xi)，σi,1=1，αi,2=1/3，μi,2= mean(Xi)，σi,2=1，αi,3=1/3， μi,3= max(Xi)，σi,3=1。將X1代入高斯混合模型，由式(7)～(10)，迭代求得其解。迭代參數為：α1,1=0.450，μ1,1=0.033，σ1,1=0.006；α1,2=0.479，μ1,2=0.032，σ1,2=0.006；α1,3=0.071，μ1,3=0.303，σ1,3=0.149。依次計算，最終得到371 265 個位點在公共正常樣本集中的高斯混合模型參數矩陣：

2.2 差異甲基化位點篩選

差異甲基化位點(DMPs)是指在腫瘤樣本和正常樣本中甲基化水平具有顯著差異的位點，可提供腫瘤純度估計的有效信息。非差異甲基化位點會成為噪音，從而削弱純度信息[5]。 本文依據前面計算的公共正常樣本集的高斯混合模型參數和位點在腫瘤樣本中的甲基化水平值，通過定義重要的統計量來篩選差異甲基化位點。

定義1設某類型的腫瘤有m 個腫瘤樣本，第i個位點在第j (j=1, 2, …, m)個腫瘤樣本中的信息貢獻率為

其中，

0≤Ri,j≤1，反映第i 個位點在第j 個腫瘤樣本和公共正常樣本中甲基化水平的差異程度。

定義2第i 個位點在腫瘤類型中的信息貢獻值為

其中，0≤Si≤m，反映第i 個位點在腫瘤樣本和公共正常樣本中甲基化水平的差異程度，Si越大，說明差異越高。 將所有位點按Si非遞增排序。設置閾值T，選取前T 個位點作為差異甲基化位點。

2.3 腫瘤純度估計

如圖1(a)所示，差異甲基化位點的β 值服從雙峰分布[4]，其峰值位置將用于腫瘤純度估計。為提高差異甲基化位點的信噪比，采用文獻[4]中的方法，先將雙峰分布變為單峰分布：將差異甲基化位點劃分為超甲基化位點和次甲基化位點（若位點在腫瘤樣本中的甲基化水平均值高于正常樣本中的甲基化水平均值，則稱該位點是超甲基化位點；反之，稱為次甲基化位點），依據合理性假設——超甲基化位點與次甲基化位點的β 值之和等于1，將次甲基化位點的β 變為1－β, 超甲基化位點的β 保持不變。對所有變換后的β 值進行核密度估計，其峰值對應位置的β 值即為腫瘤純度。

3 結果及分析

由于GmmPurify 與ABSOLUTE 采用的數據類型不同，為方便比較，選取來自9 種不同腫瘤類型的3 759 個腫瘤樣本，它們同時具有SNP array 和DNA 甲基化2 種數據，并對不同的腫瘤類型和比較對象（ABSOLUTE 和InfiniumPurify）設置了不同的閾值T。結果表明，GmmPurify 計算得到的純度估值與 2 種方法的結果具有強相關性，和InfiniumPurify 的結果高度一致（見表1）。圖1(b)展示了GmmPurify 和ABSOLUTE 對196 個肺腺癌腫瘤樣本的純度估值的強相關性（皮爾遜相關系數為0.834）；圖1(c)展示了GmmPurify 和InfiniumPurify對466 個肺腺癌腫瘤樣本的純度估值的強相關性（皮爾遜相關系數為0.929）。

圖1 肺腺癌（LUAD）選擇T=200 時的結果展示Fig.1 The results of lung adenocarcinoma（LUAD）with T=200

表1 最高皮爾遜相關系數表Table 1 Table of the highest Pearson correlations

此外，本文還測試了閾值T 對GmmPurify 純度估計的影響。如圖2 所示，T 過大或過小都會影響純度估計的準確性。 如當T=50 000 時，GmmPurify 的 純 度 估 值 與 ABSOLUTE 和InfiniumPurify 的相關性大大降低?？梢?，每種腫瘤類型都有各自的最佳閾值區間，在此區間內得到的純度估值相對最優且比較穩定。如果要求所有的腫瘤類型都選擇相同的閾值，建議T=1 000，此時能夠得到令人滿意的純度估計結果（見表2）。

以上分析說明，GmmPurify 可以作為腫瘤純度估計的可選擇工具。

4 結 論

腫瘤純度反映了癌癥類型的內在特性以及臨床收集樣品的精確度，對后續的癌癥數據分析有重要影響。本文提出了一種計算簡單、可適用于多種腫瘤類型、具有較強生物學意義的腫瘤純度估算方法GmmPurify。估算結果與經典的ABSOLUTE 和InfiniumPurify 高度一致，驗證了方法的可行性。此外，所提出的統計量“信息貢獻值”可作為后續研究DNA 差異甲基化分析的重要工具。

圖2 9 種腫瘤類型在不同閾值下的純度估值分析Fig.2 The analysis of purity estimations for 9 tumors with different thresholds

表2 T=1 000 時的相關系數表Table 2 Table of correlations with T=1 000