孟令豪，章　琳，厲力華

(杭州電子科技大學(xué)生命信息與儀器工程學(xué)院，浙江 杭州 310018)

?

基于多維基因組學(xué)的卵巢癌亞型分析

孟令豪，章琳，厲力華

(杭州電子科技大學(xué)生命信息與儀器工程學(xué)院，浙江 杭州 310018)

摘要：卵巢癌預(yù)后較差且個體差異很大，有必要從多維基因組學(xué)的角度來理解卵巢癌復(fù)雜的致癌機制，以期獲得導(dǎo)致卵巢癌亞型間預(yù)后差異的分子機制.鑒于從基因組學(xué)數(shù)據(jù)中所識別的預(yù)后特征往往不具備較好的一致性，提出一種多維基因組學(xué)中應(yīng)用通路活性對卵巢癌亞型進行區(qū)分的方法.首先，通過DNA拷貝數(shù)、DNA甲基化及miRNA等基因表達調(diào)控因素與基因表達變化相結(jié)合的方法來識別卵巢癌預(yù)后相關(guān)的14個通路.其次，基于卵巢癌預(yù)后相關(guān)的通路活性在不同亞型中有不同的模式，得到4種存在著生存差別的穩(wěn)定亞型，并提取了其中預(yù)后最差的亞型區(qū)別于其它亞型的TGF-β等顯著差異通路.最后，獨立數(shù)據(jù)驗證顯示，該組顯著差異通路的活性在驗證集預(yù)后最差的亞型中也發(fā)生了一致地改變.識別預(yù)后顯著差別的亞型及其內(nèi)在一致的分子機制可以對卵巢癌的預(yù)后預(yù)測及治療提供一定的參考.

關(guān)鍵詞：多維基因組學(xué)；生存分析；通路活性；聚類分析；卵巢癌亞型

0引言

卵巢癌是較為常見的婦科惡性腫瘤之一，同時也是致死率最高的婦科腫瘤[1].卵巢癌預(yù)后較差并且不同患者的預(yù)后差別較大，目前影響卵巢癌預(yù)后的分子機制仍不明確.研究表明，卵巢癌可能在分子層面上存在著多種亞型.文獻[2]基于基因表達譜發(fā)現(xiàn)了包括間葉型(mesenchymal)在內(nèi)的4種轉(zhuǎn)錄組亞型，然而這4種亞型間并沒有顯著的生存差異.文獻[3]在4種轉(zhuǎn)錄亞型的間葉型基礎(chǔ)上，通過整合基因表達譜、拷貝數(shù)變異、DNA甲基化和miRNA表達等數(shù)據(jù)信息，進一步發(fā)現(xiàn)了具有顯著生存差異的2種亞型.因此，對于癌癥這種復(fù)雜疾病的研究應(yīng)考慮遺傳組和表觀遺傳組上的改變，而不是僅僅基于基因表達譜數(shù)據(jù)來理解癌癥亞型[4].同時，大規(guī)模的癌癥基因組計劃已經(jīng)積累了大量癌癥樣本的全基因組數(shù)據(jù)，為全面解讀癌癥發(fā)生發(fā)展的分子機制提供了前所未有的機會.

目前，已有一些研究人員利用高通量基因表達數(shù)據(jù)，得到了與生存顯著相關(guān)的基因標志并用其對卵巢癌預(yù)后亞型進行區(qū)分[5-7]，然而不同研究得到的基因標志往往不具有一致性.比如，文獻[5]利用功能類評分方法得到了88個與預(yù)后相關(guān)的基因標志，然而這些基因和文獻[7]利用SPLASH方法得到的115個與預(yù)后相關(guān)的基因標志沒有交疊.事實上，這類分析忽略了基因間的相互關(guān)系，基因通常是通過網(wǎng)絡(luò)或通路上的相互關(guān)聯(lián)來發(fā)揮作用的.為了解決這個問題，許多研究利用通路信息來獲取與癌癥預(yù)后相關(guān)的分子標志.比如，文獻[5]對卵巢癌分析得到的17個通路與文獻[6]得到的48個通路有16個是重疊的.進一步研究表明，基于通路信息的分類器通常能達到或超過基于基因的分類器的分類效果，并且能更好地解釋卵巢癌亞型相關(guān)的生物學(xué)機制[8-9].因此，為了提取卵巢癌預(yù)后的一致特征并由此鑒別出預(yù)后有顯著差別的亞型，本文提出了一種基于多維基因組學(xué)和通路活性的亞型劃分方法.通過對卵巢癌數(shù)據(jù)的獨立驗證表明，本方法可以識別出存在著預(yù)后顯著差異的卵巢癌亞型及可能導(dǎo)致預(yù)后差異的內(nèi)在通路的活性差異.同時，通路活性變化與亞型預(yù)后的對應(yīng)關(guān)系在獨立驗證中呈現(xiàn)一致性，這對有針對性地治療卵巢癌患者有積極作用.

1數(shù)據(jù)與方法

1.1數(shù)據(jù)的采集與預(yù)處理

本文采用的卵巢癌的多維基因組和臨床數(shù)據(jù)來源于TCGA數(shù)據(jù)庫.選取了同時記錄生存信息和多維基因組(基因表達、miRNA表達、DNA甲基化、拷貝數(shù))數(shù)據(jù)的卵巢癌樣本556例.其中，基因表達譜經(jīng)過lowess標準化處理，原始的拷貝數(shù)改變由CBS分割后經(jīng)GISTIC2.0軟件提取獲得，DNA甲基化數(shù)據(jù)采用的是映射到基因組后的beta值.選用的驗證集來自GEO數(shù)據(jù)庫的GSE9891[10]，共有可用卵巢癌樣本224例，采集了臨床數(shù)據(jù)和基因表達譜數(shù)據(jù).

1.2預(yù)后相關(guān)的特征選擇及打分

Cox回歸模型對于樣本的生存時間無分布要求，且可以應(yīng)用于刪失數(shù)據(jù).首先，利用Cox回歸模型，針對mRNA表達變化、miRNA表達變化、DNA甲基化異常和DNA拷貝數(shù)變異，分別篩選出其中對生存時間有顯著影響的風(fēng)險因子(wald檢驗，p≤0.05).對于上述篩選得到的風(fēng)險因子，依據(jù)以下條件進一步篩選出候選的基因列表：1)mRNA表達變化與生存顯著相關(guān)；2)影響mRNA表達的miRNA表達、DNA甲基化或DNA拷貝數(shù)變異也與生存顯著相關(guān).

其次，針對篩選出的基因列表，利用ConsensusPathDB在線工具[11]進行富集分析，得到相應(yīng)的KEGG通路.最后，采用生存相關(guān)的候選基因所富集到的通路的活性作為區(qū)分卵巢癌亞型的特征，通路的活性值打分用通路內(nèi)基因表達值的綜合性度量[8](本文采用平均值)來表示.

1.3基于通路活性的亞型分型

采用非監(jiān)督聚類分析中經(jīng)典的K-means方法來區(qū)分卵巢癌亞型.其中，K-means聚類中存在著聚類數(shù)目的估計問題，采用常用的Gap統(tǒng)計量的方法來估計[12].Gap統(tǒng)計量方法中的關(guān)鍵問題是如何選擇合適的參考分布.在多維情形下，選取參考分布的方法可以利用奇異值分解方法得到原數(shù)據(jù)的主成分，在服從均勻分布的主成分數(shù)據(jù)中生成參考分布，因為其產(chǎn)生的參考數(shù)據(jù)集往往更接近原始數(shù)據(jù)集.

G(k)≥G(k+1)-sk+1

(1)

1.4亞型穩(wěn)定性檢驗

為了評價聚類分析得到的類簇是否具有聚類穩(wěn)定性，本文采用集合間相似性的度量Jaccard系數(shù)來進行衡量.Jaccard系數(shù)定義為兩個集合交集元素的數(shù)目和該兩個集合并集元素的數(shù)目之比，對于兩個類簇A和B進行相似性度量，表示為：

(2)

通常，Jaccard系數(shù)的值是在0～1之間.本文中，用Jaccard系數(shù)的bootstrap分布來估計聚類類簇的穩(wěn)定性[13]，當(dāng)Jaccard系數(shù)的均值小于0.60時，認為該類簇是相當(dāng)不穩(wěn)定的；當(dāng)處于0.60～0.75之間時，認為該類簇可以作為從數(shù)據(jù)發(fā)掘的有效模式；當(dāng)大于0.85時，認為該類簇具有高穩(wěn)定性且是有意義的.

1.5提取亞型間顯著差異的特征通路

RankProduct方法是一種基于倍數(shù)變化的秩的簡單無參數(shù)的差異特征提取方法[14].此方法不要求數(shù)據(jù)集滿足特定的分布，采用了簡單而嚴謹?shù)慕y(tǒng)計方法來確定每個特征的顯著性水平，并允許靈活控制錯誤發(fā)現(xiàn)率.

(3)

2結(jié)果與分析

2.1卵巢癌預(yù)后特征及其相關(guān)通路的選取

為了篩選mRNA表達、miRNA表達、DNA甲基化和拷貝數(shù)變異數(shù)據(jù)中與卵巢癌預(yù)后顯著相關(guān)的特征，使用單變量Cox回歸模型來進行回歸分析.在5%的顯著性水平下，選取了1 710個mRNA表達、131個miRNA表達、1 449個DNA甲基化位點和2 901個基因拷貝數(shù)變異特征.對于上述選取的特征，進一步利用前述的規(guī)則篩選出候選的基因列表，得到570個預(yù)后相關(guān)的基因.使用ConsensusPathDB在線工具對所選的570個基因進行富集分析，得到14個預(yù)后相關(guān)的KEGG通路，如表1所示.

表1　預(yù)后相關(guān)的KEGG通路

2.2基于通路活性的卵巢癌亞型及其穩(wěn)定性驗證

利用上述得到的14個通路，以通路內(nèi)基因表達值的均值作為通路的活性打分，對卵巢癌樣本進行聚類以區(qū)分亞型.在確定樣本的最優(yōu)聚類數(shù)時，利用Gap統(tǒng)計量結(jié)合K-means聚類的方法，生成100個參考樣本集，取k=1,2,…,10，分別計算G(k)，聚類數(shù)目k所對應(yīng)的Gap統(tǒng)計量如圖1所示.圖1(a)中，當(dāng)訓(xùn)練集中聚類數(shù)k=4時，Gap統(tǒng)計量滿足最優(yōu)聚類數(shù)目條件，可以確定訓(xùn)練集樣本聚為4種亞型.圖1(b)中，當(dāng)驗證集中聚類數(shù)k=4時，Gap統(tǒng)計量也滿足最優(yōu)聚類數(shù)目條件，可以確定驗證集樣本也聚為4種亞型.

圖1　聚類數(shù)目k所對應(yīng)的Gap統(tǒng)計量

圖2　556個樣本的通路活性的K-means聚類圖

對于訓(xùn)練集樣本，基于通路活性特征的K-means聚類結(jié)果如圖2所示，樣本聚為4種亞型，分別含有75，134，201和146個樣本.通過Jaccard系數(shù)評價其聚類穩(wěn)定性.首先，對樣本進行10 000次bootstrap抽樣；然后，對每一次抽樣后利用K-means方法重新聚類(保持聚類數(shù)目為4)；最后，計算Jaccard系數(shù)均值.經(jīng)過計算，這4個聚類的Jaccard系數(shù)的均值分別為0.94，0.91，0.92，0.95.4種亞型的穩(wěn)定性指標均超過了0.85，說明這些亞型是穩(wěn)定非隨機的.驗證集樣本通過K-means聚類得到4種亞型，這4種亞型中樣本數(shù)目分別為50，79，29和66，對應(yīng)的平均Jaccard系數(shù)分別為0.74，0.75，0.82和0.67，可以認為它們是驗證集數(shù)據(jù)的有效聚類模式.

2.3亞型的生存特征及獨立數(shù)據(jù)集驗證

對于訓(xùn)練集識別的4種亞型，計算得到亞型內(nèi)樣本的平均生存時間(單位為月)分別為33.0，38.5，48.3和52.7.大多數(shù)的樣本(約75.9%)能分配到高風(fēng)險和低風(fēng)險的類當(dāng)中.其中有75個樣本屬于亞型1，其平均生存時間小于3a，具有較差的預(yù)后.而347個樣本屬于亞型3和亞型4，其平均生存時間均大于4a，具有較好的預(yù)后.預(yù)后最好的亞類4和預(yù)后最差的亞類1的平均生存時間具有明顯差異，前者平均生存時間約是后者的1.6倍，兩者的生存曲線如圖3(a)所示，兩條曲線具有顯著差異(log-rank檢驗，p=0.000 2).

對于驗證集中識別的4個亞型，其平均生存時間(單位為月)分別為52.1，34.0，21.2和34.6.預(yù)后最好的亞型1和預(yù)后最差的亞型3的平均生存時間也具有明顯差異，前者平均生存時間約是后者的2.5倍，兩者的生存曲線如圖3(b)所示，也具有顯著差異(log-rank檢驗，p=0.014).同樣的，近一半樣本也可被分配到高風(fēng)險和低風(fēng)險的類當(dāng)中.以上分析表明，通過整合多維基因組學(xué)改變所富集的通路進行定義的亞型，按其通路活性模式不同指示了不同的生存時間.

2.4亞型的通路活性特征及獨立數(shù)據(jù)集驗證

為了驗證通路活性的程度能否在一定程度上對亞型的預(yù)后情況進行判斷，本文對預(yù)后最差的亞型提取通路活性特征并驗證其在驗證集中的一致性.在訓(xùn)練集中，采用RankProduct方法計算預(yù)后最差的亞型1與其他亞型之間存在顯著差異的通路.在PFDR≤0.05的顯著性水平下，亞型1在通路p4，p2，p6和p14上的活性是明顯上調(diào)的.通路p4，p2，p6和p14涉及到了許多有關(guān)腫瘤發(fā)生和轉(zhuǎn)移的功能，如細胞粘附，血管生成，轉(zhuǎn)化生長因子β(TGF-β)的結(jié)合，和調(diào)節(jié)細胞增殖等.這也能較好解釋亞型1的預(yù)后最差的原因.在驗證集上，同樣地應(yīng)用RankProduct方法，預(yù)后最差的亞型3相對于其它亞型在通路p4，p2和p6上的活性也是顯著上調(diào)(p14通路在驗證集不顯著，但調(diào)控方向同樣是上調(diào)的).這說明通路活性和生存時間之間的對應(yīng)關(guān)系是可以重復(fù)的，通過少數(shù)通路的活性可以高重復(fù)性地定義亞型，可以很好地識別出生存不同的亞型其背后較為明確的生物學(xué)機制，這一結(jié)果對有針對性地治療患者有積極意義.

進一步比較臨床上常用的stage分期亞型和基于通路活性的亞型.根據(jù)臨床分期信息，得出亞型1中的樣本中處于StageIIIC和StageIV的占比分別為78.67%，17.33%，具有比其他亞型樣本更高階段的分期(Wilcoxon秩和檢驗，p=0.02).StageIIIC和StageIV在臨床上認為是已發(fā)生腫瘤轉(zhuǎn)移，尤其是StageIV意味著遠端轉(zhuǎn)移.通過比較亞型1和數(shù)據(jù)集中StageIIIC，StageIV的樣本平均生存時間，發(fā)現(xiàn)亞型1與StageIIIC(平均生存時間為42.8月)內(nèi)的樣本的生存時間存在著顯著區(qū)別(log-rank檢驗，p=0.013).亞型1與StageIV(平均生存時間為35.2月)內(nèi)的樣本沒有顯著差異，然而亞型1內(nèi)的樣本5年生存率仍是低于StageIV的.這表明在分子水平特征上識別的亞型1與臨床分期存在著差異，該差異也可作為臨床分型的補充.

3結(jié)束語

本文利用基因表達譜，結(jié)合拷貝數(shù)變異、DNA甲基化和miRNA表達等數(shù)據(jù)信息，得到了卵巢癌生存相關(guān)的通路并以通路活性作為特征來進行卵巢癌亞型的分析，最終得到了預(yù)后有著差異的亞型.同時，亞型預(yù)后與通路活性的對應(yīng)關(guān)系在獨立數(shù)據(jù)驗證中呈現(xiàn)一致性.其中，預(yù)后最差的亞型表現(xiàn)為TGF-β等通路活性顯著上調(diào)，文獻[3-4]也有相關(guān)報道.從預(yù)后異質(zhì)性的卵巢癌患者中識別出具備共同的通路活性特征的卵巢癌預(yù)后相關(guān)的亞型，對預(yù)測患者的預(yù)后及治療方案的選擇具有一定的參考價值.

參考文獻

[1]SIEGEL R L, MA J, ZOU Z, et al. Cancer statistics, 2014[J]. Ca A Cancer Journal for Clinicians, 2014, 64(1): 9-29.

[2]The Cancer Genome Atlas Research Network. Integrated genomic analyses of ovarian carcinoma[J]. Nature, 2011, 474(7353): 609-615.

[3]YANG D, SUN Y, HU L, et al. Integrated analyses identify a master microRNA regulatory network for the mesenchymal subtype in serous ovarian cancer[J]. Cancer cell, 2013, 23(2): 186-199.

[4]ZHANG W, LIU Y, SUN N, et al. Integrating genomic, epigenomic, and transcriptomic features reveals modular signatures underlying poor prognosis in ovarian cancer[J]. Cell reports, 2013, 4(3): 542-553.

[5]CRIJNS A P G, FEHRMANN R S N, dE JONG S, et al. Survival-related profile, pathways, and transcription factors in ovarian cancer[J]. PLoS medicine, 2009, 6(2): 181-193.

[6]DRESSMAN H K, BERCHUCK A, CHAN G, et al. An integrated genomic-based approach to individualized treatment of patients with advanced-stage ovarian cancer[J]. Journal of Clinical Oncology, 2007, 25(5): 517-525.

[7]SPENTZOS D, LEVINE D A, RAMONI M F, et al. Gene expression signature with independent prognostic significance in epithelial ovarian cancer[J]. Journal of Clinical Oncology, 2004, 22(23): 4700-4710.

[8]GUO Z, ZHANG T, LI X, et al. Towards precise classification of cancers based on robust gene functional expression profiles[J]. BMC bioinformatics, 2005, 6(1): 1-12.

[9]SU J, YOON B J, DOUGHERTY E R. Accurate and reliable cancer classification based on probabilistic inference of pathway activity[J]. PloS one, 2009, 4(12): 3875-3884.

[10]TOTHILL R W, TINKER A V, GEORGE J, et al. Novel molecular subtypes of serous and endometrioid ovarian cancer linked to clinical outcome[J]. Clinical Cancer Research, 2008, 14(16): 5198-5208.

[11]KAMBUROV A, PENTCHEV K, GALICKA H, et al. ConsensusPathDB: toward a more complete picture of cell biology[J]. Nucleic acids research, 2011, 39(Database issue): D712-D717.

[12]TIBSHIRANI R, WALTHER G, HASTIE T. Estimating the number of clusters in a data set via the gap statistic[J]. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 2001, 63(2): 411-423.

[13]HENNIG C. Cluster-wise assessment of cluster stability[J]. Computational Statistics & Data Analysis, 2007, 52(1): 258-271.

[14]HONG F, BREITLING R, MCENTEE C W, et al. RankProd: a bioconductor package for detecting differentially expressed genes in meta-analysis[J]. Bioinformatics, 2006, 22(22): 2825-2827.

Ovarian Cancer Subtype Analysis Based on Multi-dimensional Genomics and Pathway Activity

MENG Linghao, ZHANG Lin, LI Lihua

(SchoolofLifeInformationScience&InstrumentEngineering,HangzhouDianziUniversity,HangzhouZhejiang310018,China)

Abstract：The outcome of ovarian cancer was poor and varied greatly among patients. To reveal the molecular mechanisms underlying the differences in subtype prognosis, it’s necessary to understand the complicated carcinogenic mechanisms of ovarian cancer from the perspective of multi-dimensional genomics. Considered prognosis related markers identified from one-dimensional genomics often lack of consistency, a subtype identification method based on pathway activity which obtained from multi-dimensional genomics was proposed in this article. Firstly, by integrating such factors which regulate gene expression as DNA copy number, DNA methylation and microRNA with expression changes, 14 prognosis-related pathways were identified for ovarian cancer. Then, 4 stable subtypes with significant differences in survival were obtained and significant pathways such as the TGF-β pathway between the poorest prognosis subtype and others were selected by the knowledge that activity of prognosis-related pathway shows a different pattern among different subtypes. Finally, it showed that the poor prognosis subgroups exhibited a reproducible relationship between pathway activity pattern and their prognosis in independent validation. Identification of subtypes with significant differences of prognosis and consistent molecular mechanism underlying it intrinsically may provide a reference to prognosis prediction and treatment of ovarian cancer.

Key words：multi-dimensional genomics; survival analysis; pathway activity; cluster analysis; ovarian cancer subtypes

DOI：10.13954/j.cnki.hdu.2016.04.007

收稿日期：2015-11-26

基金項目：國家自然科學(xué)基金資助項目(61271063)；國家重點基礎(chǔ)研究發(fā)展計劃(973計劃)資助項目(2013CB329502)；浙江省自然科學(xué)基金資助項目(L215F010001)

作者簡介：孟令豪(1990-)，男，山東臨沂人，碩士研究生，生物醫(yī)學(xué)工程.通信作者：厲力華教授，E-mail：lilh@hdu.edu.cn.

中圖分類號：Q811.4

文獻標識碼：A

文章編號：1001-9146(2016)04-0029-06