基于rMKL-LPP方法的乳頭狀腎細(xì)胞癌多組學(xué)數(shù)據(jù)整合分型分析*

李靈梅 魏億芳 李 治 房瑞玲 崔躍華 曹紅艷，4△

【提 要】 目的 探討局部保留投影的正則化多核學(xué)習(xí)(regularized multiple kernel learning with locality preserving projections，rMKL-LPP)在乳頭狀腎細(xì)胞癌(papillary renal cell carcinoma，PRCC)多組學(xué)數(shù)據(jù)分子分型中的應(yīng)用，進(jìn)一步研究PRCC分子分型在信號(hào)通路活性和基因表達(dá)調(diào)控方面的異質(zhì)性。方法 采用rMKL-LPP方法對(duì)PRCC的mRNA、miRNA和DNA甲基化數(shù)據(jù)進(jìn)行整合，進(jìn)一步采用k-means方法聚類分型，并通過Cox回歸分析研究不同分型的預(yù)后風(fēng)險(xiǎn)。針對(duì)不同分型，進(jìn)行通路活性分析，使用差異表達(dá)分析篩選DEmRNAs(differentially expressed mRNAs)，DEmiRNAs(differentially expressed miRNAs)和DMGs(differentially methylated genes)，并對(duì)三者的重合基因進(jìn)行GO(gene ontology)富集分析，最后使用相關(guān)及生存分析篩選可能受DNA甲基化或miRNA調(diào)控且影響患者生存的基因。結(jié)果 PRCC患者分為三型，不同亞型在通路活性和基因表達(dá)方面均有差異。篩選出10條活性存在差異的通路；1185個(gè)DEmRNAs，13個(gè)DEmiRNAs及416個(gè)DMGs，其中36個(gè)重合基因富集于有統(tǒng)計(jì)學(xué)差異的8個(gè)GO生物項(xiàng)。相關(guān)分析發(fā)現(xiàn)，ABL2可能受hsa-miR-107調(diào)控，13個(gè)基因可能受DNA甲基化調(diào)控。生存分析表明，ZNF135和RBPMS2可能與患者生存結(jié)局相關(guān)。結(jié)論 rMKL-LPP能夠有效識(shí)別PRCC亞型，篩選出的通路及潛在生物標(biāo)志物，可為PRCC針對(duì)性治療提供依據(jù)。

乳頭狀腎細(xì)胞癌(papillary renal cell carcinoma，PRCC)是腎細(xì)胞癌(renal cell carcinoma，RCC)的第二常見亞型，占RCC病例的10%～15%[1]。PRCC具有很強(qiáng)的異質(zhì)性，組織學(xué)上通常分為兩型：Type Ⅰ和Type Ⅱ。其中，Type Ⅰ型屬于低級(jí)別腫瘤，Type Ⅱ型屬于高級(jí)別腫瘤[2]。Type Ⅱ較Type Ⅰ異質(zhì)性更強(qiáng)，預(yù)后更差[3]，可分化為高度惡性的RCC肉瘤樣型[4]。該組織學(xué)分型常用于傳統(tǒng)臨床對(duì)PRCC患者進(jìn)行預(yù)后評(píng)估，然而同一組織類型且臨床分期相近的患者即使采用相同的治療手段，其治療效果和預(yù)后結(jié)局亦相差較大，評(píng)估效果并不理想。隨著測(cè)序技術(shù)的發(fā)展，積累了大量組學(xué)數(shù)據(jù)，從組學(xué)層面研究癌癥分子分型隨之興起。整合多組學(xué)數(shù)據(jù)不僅能夠捕獲PRCC不同組學(xué)的異質(zhì)性，同時(shí)還可獲得組學(xué)間的關(guān)聯(lián)信息[5]，從多層面揭示疾病的復(fù)雜調(diào)控機(jī)制。如何利用組學(xué)數(shù)據(jù)對(duì)患者精準(zhǔn)分型，為治療方案的選擇及預(yù)后評(píng)估提供幫助，進(jìn)而實(shí)現(xiàn)精準(zhǔn)治療，是PRCC臨床治療的重要發(fā)展方向。

TCGA研究組[3](2016)采用COCA方法[6](cluster-of-clusters analysis)綜合PRCC患者的miRNA/mRNA、拷貝數(shù)變異、蛋白表達(dá)數(shù)據(jù)及DNA甲基化數(shù)據(jù)，首次對(duì)PRCC患者進(jìn)行了分子分型。COCA是一種兩步聚類法，首先基于不同數(shù)據(jù)類型的聚類結(jié)果構(gòu)建一個(gè)二進(jìn)制矩陣，然后輸入該矩陣進(jìn)行一致性聚類，得到一個(gè)綜合不同數(shù)據(jù)集的全局聚類結(jié)構(gòu)。然而，COCA屬于后期整合方法，在對(duì)每個(gè)組學(xué)數(shù)據(jù)進(jìn)行單獨(dú)聚類時(shí)，易損失較弱的數(shù)據(jù)信號(hào)[7]，而且在組合不同數(shù)據(jù)的聚類結(jié)構(gòu)時(shí)未能考慮不同組學(xué)對(duì)分型的貢獻(xiàn)[8]。而基于多核學(xué)習(xí)[9]的方法，將不同n×pi的組學(xué)數(shù)據(jù)分別轉(zhuǎn)換為n×n的樣本相似矩陣，通過學(xué)習(xí)優(yōu)化，得到最優(yōu)樣本相似矩陣的線性組合，能夠反映不同類型數(shù)據(jù)的權(quán)重，在多組學(xué)數(shù)據(jù)整合分型中獨(dú)具優(yōu)勢(shì)。Speicher等[10]將多核學(xué)習(xí)與局部保留投影降維方法[11](locality preserving projections，LPP)結(jié)合，提出了rMKL-LPP。rMKL-LPP具有以下特點(diǎn)：(1)可基于數(shù)據(jù)類型靈活選擇核函數(shù)；(2)樣本相似矩陣的權(quán)重即為不同組學(xué)數(shù)據(jù)的貢獻(xiàn)度，反映了不同組學(xué)對(duì)分型的貢獻(xiàn)；(3)每個(gè)數(shù)據(jù)類型可設(shè)置多個(gè)核函數(shù)，避免了核參數(shù)設(shè)定的局限性。此外，Rappoport和Shamir[7]研究不同整合方法在10種TCGA癌癥分型中的應(yīng)用時(shí)，指出rMKL-LPP相比其他方法，更能有效識(shí)別出與臨床特征及生存率顯著相關(guān)的分子亞型。

因此，本文采用rMKL-LPP算法，整合PRCC患者mRNA、miRNA及DNA甲基化數(shù)據(jù)進(jìn)行分型，并尋找不同分型的重要通路及差異表達(dá)基因，為實(shí)現(xiàn)乳頭狀腎細(xì)胞癌不同分型的針對(duì)性治療提供參考。

數(shù)據(jù)和方法

1.數(shù)據(jù)來源

使用R包TCGAbiolinks[12]下載PRCC的mRNA、miRNA、DNA甲基化及臨床數(shù)據(jù)，進(jìn)行ID匹配后，得到表達(dá)矩陣：56493×219的mRNA矩陣、1881×219的miRNA矩陣及485577×219的DNA甲基化矩陣，其中行表示每個(gè)組學(xué)數(shù)據(jù)的特征，列表示樣本。數(shù)據(jù)預(yù)處理方法如下：(1)對(duì)啟動(dòng)子區(qū)域CpG甲基化位點(diǎn)進(jìn)行注釋，啟動(dòng)子區(qū)域?yàn)檗D(zhuǎn)錄起始位點(diǎn)2kbp內(nèi)的區(qū)域[13]，進(jìn)一步去除性染色體上的啟動(dòng)子區(qū)CpG位點(diǎn)。(2)刪除缺失比例大于30%的特征，用KNN(k-nearest neighbors)算法填補(bǔ)剩余缺失值，并對(duì)mRNA和miRNA數(shù)據(jù)進(jìn)行l(wèi)og2轉(zhuǎn)換。最終得到16534個(gè)mRNA，437個(gè)miRNA和49022個(gè)DNA甲基化位點(diǎn)。

2.分析方法

多核學(xué)習(xí)降維(multiple kernel learning for dimensionality reduction，MKL-DR)方法[14]通過使用核函數(shù)，將不同數(shù)據(jù)集映射到高維空間并進(jìn)行集成，然后通過降維算法將集成結(jié)果映射到低維空間，進(jìn)行后續(xù)分析[15]。rMKL-LPP在MKL-DR的基礎(chǔ)上，采用LPP進(jìn)行降維，同時(shí)為了避免優(yōu)化問題中的過擬合，加入了正則約束項(xiàng)。方法原理如下：

(1)多核學(xué)習(xí)

多核學(xué)習(xí)將M個(gè)給定的基本核函數(shù){k1，…，kM}線性組合，通過優(yōu)化權(quán)重系數(shù)得到一個(gè)融合核K，如公式(1)所示。

(1)

其中Km表示基本核函數(shù)，βm是核函數(shù)Km的權(quán)重系數(shù)。

(2)局部保留投影降維LPP

LPP是一種基于圖嵌入框架的無監(jiān)督方法，旨在尋找最優(yōu)投影向量v，使得經(jīng)v映射后，樣本在優(yōu)化空間中仍然能夠保持高維空間中的近鄰關(guān)系。v根據(jù)圖保留準(zhǔn)則(graph-preserving criterion)進(jìn)行優(yōu)化：

(2)

(3)

(4)

(5)

其中xi和xj表示第i和j個(gè)樣本，元素wij構(gòu)成相似矩陣W，元素dij組成約束矩陣D，Nk(i)和Nk(j)為數(shù)據(jù)點(diǎn)i和j的最近鄰數(shù)。

(3)引入正則約束的優(yōu)化

(6)

其中α是一維情況下的投影向量，Ki為集成空間，β是核函數(shù)的權(quán)重向量。對(duì)于多維數(shù)據(jù)，將針對(duì)投影矩陣A=[α1…αp]進(jìn)行優(yōu)化，并采用坐標(biāo)下降法交替對(duì)A與β進(jìn)行迭代優(yōu)化，直到達(dá)到收斂或最大迭代次數(shù)。若從優(yōu)化A開始，令所有核矩陣權(quán)重β數(shù)值相同，且總和為1；若從優(yōu)化β開始，AAT應(yīng)初始化為I。

(4)k-means聚類

rMKL-LPP通過LPP將集成結(jié)果投影到低維空間，進(jìn)一步采用k-means方法進(jìn)行聚類分型，并根據(jù)輪廓系數(shù)選擇最優(yōu)分型數(shù)。

(5)rMKL-LPP參數(shù)選擇

3.基于PRCC分型結(jié)果的下游分析

(1)Cox回歸分析

控制初始診斷年齡、性別和病理診斷分期等情況下，采用Cox回歸分析對(duì)PRCC患者的分型結(jié)果進(jìn)行預(yù)后評(píng)估。

(2)通路活性分析

利用progeny軟件包[16]對(duì)不同亞型進(jìn)行通路活性分析，并采用非參數(shù)檢驗(yàn)篩選活性存在差異的通路，篩選標(biāo)準(zhǔn)為Padj<0.01。

(3)差異基因篩選

采用Kruskal-Wallis秩和檢驗(yàn)篩選DEmRNAs、DEmiRNAs及DMGs，閾值設(shè)為Padj<0.01；進(jìn)一步采用超幾何分布檢驗(yàn)[17]篩選在每個(gè)分型上富集的特征，篩選標(biāo)準(zhǔn)為Padj<0.01。為選擇最具代表性的特征，要求特征在該分型中至少有2/3的樣本發(fā)生改變，同時(shí)至少在一個(gè)其他分型中少于1/3樣本發(fā)生改變，按此標(biāo)準(zhǔn)選出的特征即為最終的差異基因。

(4)GO富集分析

利用miRWalk[18]在線工具預(yù)測(cè)DEmiRNAs的靶基因，進(jìn)一步采用clusterprofile R包[19]對(duì)DEmRNAs、DEmiRNAs靶基因及DMGs的重合基因進(jìn)行富集分析。

(5)相關(guān)分析

分別對(duì)DEmRNAs與DMGs的重合基因，DEmRNAs與DEmiRNAs靶基因的重合基因進(jìn)行Pearson相關(guān)分析，根據(jù)相關(guān)系數(shù)r和P值篩選出可能受DNA甲基化負(fù)調(diào)控的基因，篩選標(biāo)準(zhǔn)為-1.0

(6)基因生存分析

對(duì)可能受DNA甲基化或miRNA調(diào)控的基因進(jìn)行生存分析。根據(jù)基因表達(dá)水平將患者分為高表達(dá)組和低表達(dá)組，應(yīng)用Kaplan-Meier生存分析篩選與患者生存相關(guān)的基因，篩選標(biāo)準(zhǔn)為P<0.05。

結(jié) 果

1.PRCC患者分型結(jié)果評(píng)價(jià)

采用rMKL-LPP對(duì)219名PRCC患者的組學(xué)數(shù)據(jù)進(jìn)行整合分析，最優(yōu)分型數(shù)為4(圖1)，生存曲線見圖2，不同分型患者的生存率存在差異(χ2=89.566，P<0.0001)。經(jīng)Cox回歸分析發(fā)現(xiàn)，Cluster2和Cluster3生存率差異無統(tǒng)計(jì)學(xué)意義(χ2=0.050，P=0.823)。因此，將Cluster2和Cluster3合并成為一個(gè)新的Cluster2，三組基本資料見表1，生存曲線見圖3。結(jié)果顯示，Cluster3患者相比于其他兩型，預(yù)后差，說明基于rMKL-LPP的分型與PRCC患者生存相關(guān)。

圖1 PRCC分型結(jié)果的三維圖

圖2 PRCC分為四型的生存曲線圖

圖3 PRCC分為三型的生存曲線圖

表1 PRCC患者分型的基本資料

在校正協(xié)變量的情況下，研究不同分型對(duì)預(yù)后的影響，即分型作為自變量，生存時(shí)間和生存狀態(tài)作為因變量，擬合Cox回歸模型，結(jié)果如表2，預(yù)后最差的Cluster3患者的死亡風(fēng)險(xiǎn)是Cluster1的47.731倍，Cluster2患者的死亡風(fēng)險(xiǎn)是Cluster1的6.143倍；病理分期差異有統(tǒng)計(jì)學(xué)意義，Ⅳ期患者死亡風(fēng)險(xiǎn)是Ⅰ期患者的20.351倍。

表2 219例PRCC患者的Cox回歸分析結(jié)果

2.通路活性分析

對(duì)PRCC亞型進(jìn)行通路活性分析，存在差異的10條通路如圖4所示，其中TGF-β、EGFR、NF-Kβ、MAPK、Hypoxia、TNF-α和PI3K通路在Cluster3中活性最高；通路Wnt和VEGF在Cluster2中活性最高，Estrogen通路在Cluster1中活性最高。不同分型通路活性的差異也在一定程度上反映了PRCC不同亞型的異質(zhì)性。

圖4 PRCC不同亞型的差異通路

3.差異基因篩選及分析

(1)差異基因的篩選結(jié)果

篩選出1185個(gè)DEmRNAs，其中上調(diào)626個(gè)，下調(diào)559個(gè)；459個(gè)差異甲基化位點(diǎn)映射到416個(gè)DMGs，包括111個(gè)高甲基化基因和305個(gè)低甲基化基因；篩選出13個(gè)DEmiRNAs，其中2個(gè)上調(diào)，11個(gè)下調(diào)。圖5依次為DEmRNAs，DMGs及DEmiRNAs表達(dá)熱圖，從圖中可以明顯看出這些特征在不同亞型中的表達(dá)差異。

圖5 不同亞型中差異基因表達(dá)熱圖

對(duì)13個(gè)DEmiRNAs進(jìn)行預(yù)測(cè)得到36個(gè)靶基因，通過對(duì)DEmRNAs，DMGs以及36個(gè)DEmiRNAs靶基因進(jìn)行聯(lián)合分析，發(fā)現(xiàn)DEmiRNAs靶基因與DEmRNAs有1個(gè)重合基因，DMGs與DEmRNAs有35個(gè)重合基因(圖6)。

圖6 差異基因的韋恩圖

(2)GO富集分析

為進(jìn)一步驗(yàn)證基于rMKL-LPP分型的生物學(xué)意義，對(duì)聯(lián)合分析得到的重合基因進(jìn)行富集分析。36個(gè)重合基因富集于8個(gè)GO生物項(xiàng)，見圖7，基因與GO生物項(xiàng)的關(guān)系如圖8所示。GO富集分析可從生物過程(biological process，BP)、分子功能(molecular function，MF)和細(xì)胞組成(cellular component，CC)等三部分對(duì)基因及基因產(chǎn)物進(jìn)行注釋。8個(gè)GO生物項(xiàng)主要體現(xiàn)在生物過程和細(xì)胞組成兩個(gè)方面。圖7中8個(gè)GO生物項(xiàng)縱軸自上而下依次為中胚層發(fā)育、色氨酸分解過程、含吲哚化合物分解代謝過程、吲哚烷基胺分解過程、色氨酸代謝過程、吲哚烷基胺代謝過程、頂端質(zhì)膜與細(xì)胞的頂端部分。圖中實(shí)心圓的大小表示富集于該通路基因的數(shù)量。

圖7 GO通路分析圖

圖8 基因與8個(gè)GO生物項(xiàng)的網(wǎng)絡(luò)關(guān)系圖

(3)相關(guān)分析

對(duì)DEmRNAs與DMGs重合基因進(jìn)行相關(guān)分析，最終得到13個(gè)存在相關(guān)關(guān)系的基因，即可能受DNA甲基化調(diào)控的基因，如圖9所示?；蛟诓煌瑏喰捅憩F(xiàn)出不同的相關(guān)關(guān)系，如ZNF135僅在Cluster1中存在相關(guān)關(guān)系；而RBPMS2僅在Cluster3中存在相關(guān)關(guān)系。對(duì)DEmRNAs與DEmiRNAs靶基因的重合基因進(jìn)行相關(guān)性分析，發(fā)現(xiàn)ABL2僅在Cluster2中與hsa-miR-107存在負(fù)相關(guān)關(guān)系(圖9)?；蛑g的相關(guān)關(guān)系表明基因間可能存在生理學(xué)調(diào)控作用。

圖9 重合基因相關(guān)關(guān)系熱圖

(4)基因生存分析

對(duì)可能受DNA甲基化或miRNA調(diào)控的基因進(jìn)行生存分析，最終得到2個(gè)可能影響PRCC患者預(yù)后的基因，如圖10所示。在Cluster1中，ZNF135低表達(dá)組患者的總生存期低于高表達(dá)組；在Cluster3中，RBPMS2高表達(dá)的患者較低表達(dá)的患者預(yù)后更差。

圖10 基因ZNF135與RBPMS2的生存曲線圖

討 論

本文采用rMKL-LPP方法對(duì)PRCC多組學(xué)數(shù)據(jù)進(jìn)行整合分型，將PRCC患者分為三型，不同分型在通路的活性、基因表達(dá)調(diào)控方面均有差異。基于分型得到的潛在生物標(biāo)記物(基因或信號(hào)通路)，將為PRCC針對(duì)性的干預(yù)治療提供重要的參考依據(jù)。

PRCC患者分為三型，Cluster1與Cluster2型PRCC發(fā)病年齡在60～65歲，且男性居多，與大多數(shù)病例的高發(fā)年齡、性別構(gòu)成基本吻合[20]。預(yù)后最差的Cluster3型患者初始診斷年齡偏小，且在女性中更為常見。結(jié)合三個(gè)亞型來看，發(fā)病年齡越早的患者預(yù)后越差，而且不同性別的患者高發(fā)年齡可能不一致。對(duì)此，臨床上應(yīng)予以重視，多關(guān)注小于50歲的患者，同時(shí)加強(qiáng)對(duì)女性患者的篩查和評(píng)估，及時(shí)進(jìn)行干預(yù)。本研究女性樣本含量較小，有關(guān)PRCC年齡及性別的差異仍需進(jìn)一步研究。

不同亞型信號(hào)通路活性的差異可幫助理解PRCC異質(zhì)性的分子基礎(chǔ)。通路TGF-β、EGFR、NF-Kβ、MAPK、Hypoxia、PI3K和TNF-α在Cluster3中活性最高。其中TGF-β可通過誘導(dǎo)患者上皮間質(zhì)轉(zhuǎn)換來促進(jìn)RCC發(fā)展[21]。EGFR通路在腎臟發(fā)育過程中起關(guān)鍵作用，可能是PRCC一個(gè)潛在的治療方向[22]。通路NF-Kβ和MAPK可調(diào)節(jié)RCC細(xì)胞增殖、侵襲和遷移[23-24]。Hypoxia通路作為腎癌的主要驅(qū)動(dòng)因素被廣泛研究，與透明細(xì)胞腎細(xì)胞癌關(guān)系密切[25-26]。PI3K在RCC中可加速細(xì)胞周期，促進(jìn)腫瘤細(xì)胞遷移[27]，據(jù)此，可通過抑制PI3K/AKT信號(hào)通路來阻止RCC的生長(zhǎng)和轉(zhuǎn)移[28]。TNF-α通路可為癌細(xì)胞的激活、分化、侵襲和增殖提供信號(hào)，促進(jìn)癌癥發(fā)展[29-30]。可見，這些通路的異常激活可能與Cluster3不良預(yù)后相關(guān)。此外，Wnt通路的異常激活可促進(jìn)RCC的轉(zhuǎn)移和惡化[31]。而VEGF和Estrogen通路可影響腎癌的發(fā)生發(fā)展[32-33]，其對(duì)于PRCC靶向治療的意義有待進(jìn)一步挖掘。

本研究基于PRCC分子分型得到三個(gè)可能受DNA甲基化或miRNA調(diào)控且影響患者生存的基因ABL2、ZNF135與RBPMS2，這三個(gè)差異分子靶標(biāo)與PRCC的關(guān)系尚不明確，但有研究發(fā)現(xiàn)它們與其他癌癥有關(guān)。ABL2是一種原癌基因，可參與調(diào)控細(xì)胞生長(zhǎng)、侵襲和遷移等過程[34]，與腎細(xì)胞癌的發(fā)生發(fā)展相關(guān)[35]。ZNF135編碼一種轉(zhuǎn)錄抑制蛋白，在透明細(xì)胞腎細(xì)胞癌、宮頸癌與乳腺癌等多種癌癥中高度甲基化[36]，本研究發(fā)現(xiàn)其在RPCC中也高度甲基化，具體作用機(jī)制有待進(jìn)一步闡明。RBPMS2的高表達(dá)與胃腸道間質(zhì)瘤有密切聯(lián)系[37]，還可促進(jìn)乳腺癌的發(fā)生發(fā)展[38]。這些基因?qū)RCC的預(yù)后有一定的預(yù)測(cè)價(jià)值，有望作為PRCC未來藥物治療的潛在靶點(diǎn)，而且基因間的調(diào)控作用也值得進(jìn)一步研究。

綜上所述，本文基于rMKL-LPP方法對(duì)PRCC多組學(xué)數(shù)據(jù)的整合分析，能夠有效地識(shí)別亞型，為PRCC的分型研究提供了新的思路。識(shí)別出的PRCC亞型在信號(hào)通路活性、基因表達(dá)及調(diào)控方面均存在差異，這有助于進(jìn)一步理解不同分型發(fā)生發(fā)展的潛在分子機(jī)制。此外，篩選出的潛在生物標(biāo)志物將為PRCC治療和預(yù)后評(píng)估提供一定的理論依據(jù)和臨床指導(dǎo)。