閆占正,李玉雙
(燕山大學(xué)理學(xué)院,河北秦皇島066004)
高通量測(cè)序技術(shù)的快速發(fā)展為癌癥基因組學(xué)研究帶來了重大突破。在腫瘤與癌癥研究中,在臨床環(huán)境下獲得的腫瘤固體組織中混合有正常細(xì)胞(即腫瘤不純),使得篩選出來的癌癥數(shù)據(jù)存在偏差,如果后續(xù)不經(jīng)處理直接使用,往往會(huì)導(dǎo)致研究結(jié)果不理想甚至得到錯(cuò)誤的結(jié)果[1]。 因此,進(jìn)行腫瘤純度的估計(jì),即推測(cè)腫瘤組織中癌細(xì)胞的比例,已成為癌癥研究中十分重要的課題[2-3],并已取得重要進(jìn)展。如,ABSOLUTE 方法,利用單核苷酸多態(tài)性陣列數(shù)據(jù)(SNP array),借助最大似然估計(jì)計(jì)算腫瘤純度[2];InfiniumPurify 模 型[4-6],利 用DNA 甲 基 化 芯 片 數(shù)據(jù)[7],將位點(diǎn)劃分為超(次)甲基化位點(diǎn),構(gòu)建線性函數(shù)來估計(jì)腫瘤純度;PurityEst,利用二代測(cè)序數(shù)據(jù),通過對(duì)腫瘤樣本中突變等位基因的比例建模來估計(jì)腫瘤純度[8]。
與其他類型數(shù)據(jù)相比,DNA 甲基化數(shù)據(jù)能更穩(wěn)定和靈活地顯示腫瘤細(xì)胞和正常細(xì)胞之間的內(nèi)在差異。觀察正常樣本,發(fā)現(xiàn)其甲基化水平分布與高斯分布相似,而有限高斯混合模型可以逼近任意概率分布密度函數(shù)[9]。因此,本文利用DNA 甲基化數(shù)據(jù),借助高斯混合模型[10-11](Gmm, Gaussian mixture model)估計(jì)腫瘤純度,得到的9 種類型的腫瘤純度與ABSOLUTE 和InfiniumPurify 的結(jié)果高度一致。之所以與這2 種方法比較,是因?yàn)锳BSOLUTE 是最早且最具影響力之一的純度估計(jì)工具,常被用作評(píng)估的金標(biāo)準(zhǔn),而InfiniumPurify 與GmmPurify 使用的數(shù)據(jù)相同,所求結(jié)果更具有可比性。
癌癥和腫瘤基因圖譜TCGA 數(shù)據(jù)庫(kù)(https://portal.gdc.cancer.gov/)中包含32 種腫瘤類型,每種類型都含有個(gè)數(shù)不同的腫瘤樣本。包含正常樣本的腫瘤類型有22 種,其中擁有2 個(gè)以上匹配的正常樣本有21 種,多形性膠質(zhì)母細(xì)胞瘤(GBM)只有1 個(gè)。由于多數(shù)腫瘤類型無(或只有很少)匹配的正常樣本,本文選用文獻(xiàn)[5]中構(gòu)造的公共正常樣本集進(jìn)行操作,即在21 種腫瘤類型中每一個(gè)類型隨機(jī)抽取2個(gè)正常樣本,并與GBM 中的1 個(gè)正常樣本一起組成包含43 個(gè)樣本的公共正常樣本集。每個(gè)腫瘤或正常樣本均包含相同的485 577 個(gè)甲基化位點(diǎn),其中大約9 萬個(gè)位點(diǎn)的β 值全部缺失,有大約1 萬個(gè)位點(diǎn)的β 值部分缺失(位點(diǎn)的β 值表示該位點(diǎn)的甲基化水平,0 ≤β ≤1,其中0 表示完全未甲基化,1 表示完全甲基化)。為避免缺失值帶來的影響,將存在缺失值的位點(diǎn)全部舍棄。以包含466 個(gè)腫瘤樣本和43 個(gè)公共正常樣本的肺腺癌(lung adenocarcinoma,LUAD)為例,執(zhí)行刪除操作后剩余371 265 個(gè)位點(diǎn)。
高斯混合模型[9]是指具有以下形式的概率分布模型:
這里μk是均值,σk是方差,αk是系數(shù),φ(y|θk)是高斯分布密度,稱其為第k 個(gè)高斯模型。每個(gè)高斯模型都是一個(gè)聚類。
當(dāng)確定樣本所屬高斯模型時(shí),可以利用最大似然概率算法來求解模型參數(shù)。設(shè)樣本總數(shù)為N,樣本集合為X={x1, x2, ..., xN}。設(shè)K 個(gè)高斯分布的樣本數(shù)量分別為N1, N2, ..., NK, 第k (k∈{1,2,...,K})個(gè)高斯分布的樣本集合為L(zhǎng)(k)。則有
當(dāng)無法確定樣本所屬模型時(shí),可用期望最大化算法通過迭代進(jìn)行求解,算法分2 步: E 步,根據(jù)參數(shù)θ(n)求 得 后 驗(yàn) 概 率P(xt|X;θ);M 步,更 新xt的 期 望E(xt),然 后 用E(xt)代 替xt求 出 新 的 模 型 參 數(shù)θ(n+1)。如此迭代直到θ 趨于穩(wěn)定。
假設(shè)模型參數(shù)已知,隨機(jī)初始化一組參數(shù)θ(0),可求得后驗(yàn)概率P(xt|X;θ)。求第k 個(gè)高斯分布,分別取x1, x2, ..., xN的概率,即求樣本xt在第k 個(gè)高斯分布下的概率γ(t,k),有公式:
由式(4)~(7),可以推出:
隨機(jī)初始化一組合適參數(shù)(α(0)k, μ(0)k, σ(0)k),k∈{1,2,...,K},將樣本集合X 和初始化參數(shù)代入式(7)~(10),進(jìn)行迭代直到參數(shù)趨于穩(wěn)定,可求得樣本集合X 的高斯混合模型參數(shù)(αk, μk, σk)。
將高斯混合模型應(yīng)用于公共正常樣本集,選取不同的K 值進(jìn)行實(shí)驗(yàn)。實(shí)驗(yàn)需確保各類高斯模型存在的合理性以及各類模型之間明顯的差異性。為此,制定以下實(shí)驗(yàn)條件:最小的分類概率不小于0.1,且各類之間的均值之差不小于0.1。研究發(fā)現(xiàn),當(dāng)K≤3 時(shí),滿足實(shí)驗(yàn)條件的甲基化位點(diǎn)占絕大多數(shù)。以LUAD 為例,其位點(diǎn)總數(shù)約37 萬。當(dāng)K=2 時(shí),滿足條件的位點(diǎn)約占29%;當(dāng)K=3 時(shí),滿足條件的位點(diǎn)約占5%;當(dāng)K=4 時(shí),滿足條件的位點(diǎn)約占0.4%。故設(shè)定K=3, N=43。以肺腺癌(LUAD)為例,依次處理371 265 個(gè)樣本集合Xi={βi,1, βi,2, ...,βi,N}, i∈{1,2,...,371 265}, 其中βi,t表示第i 個(gè)位點(diǎn)在第t 個(gè)公共正常樣本下的甲基化水平值。具體為:初始化參數(shù)α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),迭代求得其解。迭代參數(shù)為:α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。依次計(jì)算,最終得到371 265 個(gè)位點(diǎn)在公共正常樣本集中的高斯混合模型參數(shù)矩陣:
差異甲基化位點(diǎn)(DMPs)是指在腫瘤樣本和正常樣本中甲基化水平具有顯著差異的位點(diǎn),可提供腫瘤純度估計(jì)的有效信息。非差異甲基化位點(diǎn)會(huì)成為噪音,從而削弱純度信息[5]。 本文依據(jù)前面計(jì)算的公共正常樣本集的高斯混合模型參數(shù)和位點(diǎn)在腫瘤樣本中的甲基化水平值,通過定義重要的統(tǒng)計(jì)量來篩選差異甲基化位點(diǎn)。
定義1設(shè)某類型的腫瘤有m 個(gè)腫瘤樣本,第i個(gè)位點(diǎn)在第j (j=1, 2, …, m)個(gè)腫瘤樣本中的信息貢獻(xiàn)率為
其中,
0≤Ri,j≤1,反映第i 個(gè)位點(diǎn)在第j 個(gè)腫瘤樣本和公共正常樣本中甲基化水平的差異程度。
定義2第i 個(gè)位點(diǎn)在腫瘤類型中的信息貢獻(xiàn)值為
其中,0≤Si≤m,反映第i 個(gè)位點(diǎn)在腫瘤樣本和公共正常樣本中甲基化水平的差異程度,Si越大,說明差異越高。 將所有位點(diǎn)按Si非遞增排序。設(shè)置閾值T,選取前T 個(gè)位點(diǎn)作為差異甲基化位點(diǎn)。
如圖1(a)所示,差異甲基化位點(diǎn)的β 值服從雙峰分布[4],其峰值位置將用于腫瘤純度估計(jì)。為提高差異甲基化位點(diǎn)的信噪比,采用文獻(xiàn)[4]中的方法,先將雙峰分布變?yōu)閱畏宸植迹簩⒉町惣谆稽c(diǎn)劃分為超甲基化位點(diǎn)和次甲基化位點(diǎn)(若位點(diǎn)在腫瘤樣本中的甲基化水平均值高于正常樣本中的甲基化水平均值,則稱該位點(diǎn)是超甲基化位點(diǎn);反之,稱為次甲基化位點(diǎn)),依據(jù)合理性假設(shè)——超甲基化位點(diǎn)與次甲基化位點(diǎn)的β 值之和等于1,將次甲基化位點(diǎn)的β 變?yōu)?-β, 超甲基化位點(diǎn)的β 保持不變。對(duì)所有變換后的β 值進(jìn)行核密度估計(jì),其峰值對(duì)應(yīng)位置的β 值即為腫瘤純度。
由于GmmPurify 與ABSOLUTE 采用的數(shù)據(jù)類型不同,為方便比較,選取來自9 種不同腫瘤類型的3 759 個(gè)腫瘤樣本,它們同時(shí)具有SNP array 和DNA 甲基化2 種數(shù)據(jù),并對(duì)不同的腫瘤類型和比較對(duì)象(ABSOLUTE 和InfiniumPurify)設(shè)置了不同的閾值T。結(jié)果表明,GmmPurify 計(jì)算得到的純度估值與 2 種方法的結(jié)果具有強(qiáng)相關(guān)性,和InfiniumPurify 的結(jié)果高度一致(見表1)。圖1(b)展示了GmmPurify 和ABSOLUTE 對(duì)196 個(gè)肺腺癌腫瘤樣本的純度估值的強(qiáng)相關(guān)性(皮爾遜相關(guān)系數(shù)為0.834);圖1(c)展示了GmmPurify 和InfiniumPurify對(duì)466 個(gè)肺腺癌腫瘤樣本的純度估值的強(qiáng)相關(guān)性(皮爾遜相關(guān)系數(shù)為0.929)。
圖1 肺腺癌(LUAD)選擇T=200 時(shí)的結(jié)果展示Fig.1 The results of lung adenocarcinoma(LUAD)with T=200
表1 最高皮爾遜相關(guān)系數(shù)表Table 1 Table of the highest Pearson correlations
此外,本文還測(cè)試了閾值T 對(duì)GmmPurify 純度估計(jì)的影響。如圖2 所示,T 過大或過小都會(huì)影響純度估計(jì)的準(zhǔn)確性。 如當(dāng)T=50 000 時(shí),GmmPurify 的 純 度 估 值 與 ABSOLUTE 和InfiniumPurify 的相關(guān)性大大降低??梢?,每種腫瘤類型都有各自的最佳閾值區(qū)間,在此區(qū)間內(nèi)得到的純度估值相對(duì)最優(yōu)且比較穩(wěn)定。如果要求所有的腫瘤類型都選擇相同的閾值,建議T=1 000,此時(shí)能夠得到令人滿意的純度估計(jì)結(jié)果(見表2)。
以上分析說明,GmmPurify 可以作為腫瘤純度估計(jì)的可選擇工具。
腫瘤純度反映了癌癥類型的內(nèi)在特性以及臨床收集樣品的精確度,對(duì)后續(xù)的癌癥數(shù)據(jù)分析有重要影響。本文提出了一種計(jì)算簡(jiǎn)單、可適用于多種腫瘤類型、具有較強(qiáng)生物學(xué)意義的腫瘤純度估算方法GmmPurify。估算結(jié)果與經(jīng)典的ABSOLUTE 和InfiniumPurify 高度一致,驗(yàn)證了方法的可行性。此外,所提出的統(tǒng)計(jì)量“信息貢獻(xiàn)值”可作為后續(xù)研究DNA 差異甲基化分析的重要工具。
圖2 9 種腫瘤類型在不同閾值下的純度估值分析Fig.2 The analysis of purity estimations for 9 tumors with different thresholds
表2 T=1 000 時(shí)的相關(guān)系數(shù)表Table 2 Table of correlations with T=1 000