基于多维基因组学的卵巢癌亚型分析

孟令豪，章　琳，厉力华

(杭州电子科技大学生命信息与仪器工程学院，浙江 杭州 310018)



基于多维基因组学的卵巢癌亚型分析

孟令豪，章琳，厉力华

(杭州电子科技大学生命信息与仪器工程学院，浙江 杭州 310018)

摘要：卵巢癌预后较差且个体差异很大，有必要从多维基因组学的角度来理解卵巢癌复杂的致癌机制，以期获得导致卵巢癌亚型间预后差异的分子机制.鉴于从基因组学数据中所识别的预后特征往往不具备较好的一致性，提出一种多维基因组学中应用通路活性对卵巢癌亚型进行区分的方法.首先，通过DNA拷贝数、DNA甲基化及miRNA等基因表达调控因素与基因表达变化相结合的方法来识别卵巢癌预后相关的14个通路.其次，基于卵巢癌预后相关的通路活性在不同亚型中有不同的模式，得到4种存在着生存差别的稳定亚型，并提取了其中预后最差的亚型区别于其它亚型的TGF-β等显著差异通路.最后，独立数据验证显示，该组显著差异通路的活性在验证集预后最差的亚型中也发生了一致地改变.识别预后显著差别的亚型及其内在一致的分子机制可以对卵巢癌的预后预测及治疗提供一定的参考.

关键词：多维基因组学；生存分析；通路活性；聚类分析；卵巢癌亚型

0引言

卵巢癌是较为常见的妇科恶性肿瘤之一，同时也是致死率最高的妇科肿瘤[1].卵巢癌预后较差并且不同患者的预后差别较大，目前影响卵巢癌预后的分子机制仍不明确.研究表明，卵巢癌可能在分子层面上存在着多种亚型.文献[2]基于基因表达谱发现了包括间叶型(mesenchymal)在内的4种转录组亚型，然而这4种亚型间并没有显著的生存差异.文献[3]在4种转录亚型的间叶型基础上，通过整合基因表达谱、拷贝数变异、DNA甲基化和miRNA表达等数据信息，进一步发现了具有显著生存差异的2种亚型.因此，对于癌症这种复杂疾病的研究应考虑遗传组和表观遗传组上的改变，而不是仅仅基于基因表达谱数据来理解癌症亚型[4].同时，大规模的癌症基因组计划已经积累了大量癌症样本的全基因组数据，为全面解读癌症发生发展的分子机制提供了前所未有的机会.

目前，已有一些研究人员利用高通量基因表达数据，得到了与生存显著相关的基因标志并用其对卵巢癌预后亚型进行区分[5-7]，然而不同研究得到的基因标志往往不具有一致性.比如，文献[5]利用功能类评分方法得到了88个与预后相关的基因标志，然而这些基因和文献[7]利用SPLASH方法得到的115个与预后相关的基因标志没有交叠.事实上，这类分析忽略了基因间的相互关系，基因通常是通过网络或通路上的相互关联来发挥作用的.为了解决这个问题，许多研究利用通路信息来获取与癌症预后相关的分子标志.比如，文献[5]对卵巢癌分析得到的17个通路与文献[6]得到的48个通路有16个是重叠的.进一步研究表明，基于通路信息的分类器通常能达到或超过基于基因的分类器的分类效果，并且能更好地解释卵巢癌亚型相关的生物学机制[8-9].因此，为了提取卵巢癌预后的一致特征并由此鉴别出预后有显著差别的亚型，本文提出了一种基于多维基因组学和通路活性的亚型划分方法.通过对卵巢癌数据的独立验证表明，本方法可以识别出存在着预后显著差异的卵巢癌亚型及可能导致预后差异的内在通路的活性差异.同时，通路活性变化与亚型预后的对应关系在独立验证中呈现一致性，这对有针对性地治疗卵巢癌患者有积极作用.

1数据与方法

1.1数据的采集与预处理

本文采用的卵巢癌的多维基因组和临床数据来源于TCGA数据库.选取了同时记录生存信息和多维基因组(基因表达、miRNA表达、DNA甲基化、拷贝数)数据的卵巢癌样本556例.其中，基因表达谱经过lowess标准化处理，原始的拷贝数改变由CBS分割后经GISTIC2.0软件提取获得，DNA甲基化数据采用的是映射到基因组后的beta值.选用的验证集来自GEO数据库的GSE9891[10]，共有可用卵巢癌样本224例，采集了临床数据和基因表达谱数据.

1.2预后相关的特征选择及打分

Cox回归模型对于样本的生存时间无分布要求，且可以应用于删失数据.首先，利用Cox回归模型，针对mRNA表达变化、miRNA表达变化、DNA甲基化异常和DNA拷贝数变异，分别筛选出其中对生存时间有显著影响的风险因子(wald检验，p≤0.05).对于上述筛选得到的风险因子，依据以下条件进一步筛选出候选的基因列表：1)mRNA表达变化与生存显著相关；2)影响mRNA表达的miRNA表达、DNA甲基化或DNA拷贝数变异也与生存显著相关.

其次，针对筛选出的基因列表，利用ConsensusPathDB在线工具[11]进行富集分析，得到相应的KEGG通路.最后，采用生存相关的候选基因所富集到的通路的活性作为区分卵巢癌亚型的特征，通路的活性值打分用通路内基因表达值的综合性度量[8](本文采用平均值)来表示.

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

采用非监督聚类分析中经典的K-means方法来区分卵巢癌亚型.其中，K-means聚类中存在着聚类数目的估计问题，采用常用的Gap统计量的方法来估计[12].Gap统计量方法中的关键问题是如何选择合适的参考分布.在多维情形下，选取参考分布的方法可以利用奇异值分解方法得到原数据的主成分，在服从均匀分布的主成分数据中生成参考分布，因为其产生的参考数据集往往更接近原始数据集.

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

(1)

1.4亚型稳定性检验

为了评价聚类分析得到的类簇是否具有聚类稳定性，本文采用集合间相似性的度量Jaccard系数来进行衡量.Jaccard系数定义为两个集合交集元素的数目和该两个集合并集元素的数目之比，对于两个类簇A和B进行相似性度量，表示为：

(2)

通常，Jaccard系数的值是在0～1之间.本文中，用Jaccard系数的bootstrap分布来估计聚类类簇的稳定性[13]，当Jaccard系数的均值小于0.60时，认为该类簇是相当不稳定的；当处于0.60～0.75之间时，认为该类簇可以作为从数据发掘的有效模式；当大于0.85时，认为该类簇具有高稳定性且是有意义的.

1.5提取亚型间显著差异的特征通路

RankProduct方法是一种基于倍数变化的秩的简单无参数的差异特征提取方法[14].此方法不要求数据集满足特定的分布，采用了简单而严谨的统计方法来确定每个特征的显著性水平，并允许灵活控制错误发现率.

(3)

2结果与分析

2.1卵巢癌预后特征及其相关通路的选取

为了筛选mRNA表达、miRNA表达、DNA甲基化和拷贝数变异数据中与卵巢癌预后显著相关的特征，使用单变量Cox回归模型来进行回归分析.在5%的显著性水平下，选取了1 710个mRNA表达、131个miRNA表达、1 449个DNA甲基化位点和2 901个基因拷贝数变异特征.对于上述选取的特征，进一步利用前述的规则筛选出候选的基因列表，得到570个预后相关的基因.使用ConsensusPathDB在线工具对所选的570个基因进行富集分析，得到14个预后相关的KEGG通路，如表1所示.

表1　预后相关的KEGG通路

2.2基于通路活性的卵巢癌亚型及其稳定性验证

利用上述得到的14个通路，以通路内基因表达值的均值作为通路的活性打分，对卵巢癌样本进行聚类以区分亚型.在确定样本的最优聚类数时，利用Gap统计量结合K-means聚类的方法，生成100个参考样本集，取k=1,2,…,10，分别计算G(k)，聚类数目k所对应的Gap统计量如图1所示.图1(a)中，当训练集中聚类数k=4时，Gap统计量满足最优聚类数目条件，可以确定训练集样本聚为4种亚型.图1(b)中，当验证集中聚类数k=4时，Gap统计量也满足最优聚类数目条件，可以确定验证集样本也聚为4种亚型.

图1　聚类数目k所对应的Gap统计量

图2　556个样本的通路活性的K-means聚类图

对于训练集样本，基于通路活性特征的K-means聚类结果如图2所示，样本聚为4种亚型，分别含有75，134，201和146个样本.通过Jaccard系数评价其聚类稳定性.首先，对样本进行10 000次bootstrap抽样；然后，对每一次抽样后利用K-means方法重新聚类(保持聚类数目为4)；最后，计算Jaccard系数均值.经过计算，这4个聚类的Jaccard系数的均值分别为0.94，0.91，0.92，0.95.4种亚型的稳定性指标均超过了0.85，说明这些亚型是稳定非随机的.验证集样本通过K-means聚类得到4种亚型，这4种亚型中样本数目分别为50，79，29和66，对应的平均Jaccard系数分别为0.74，0.75，0.82和0.67，可以认为它们是验证集数据的有效聚类模式.

2.3亚型的生存特征及独立数据集验证

对于训练集识别的4种亚型，计算得到亚型内样本的平均生存时间(单位为月)分别为33.0，38.5，48.3和52.7.大多数的样本(约75.9%)能分配到高风险和低风险的类当中.其中有75个样本属于亚型1，其平均生存时间小于3a，具有较差的预后.而347个样本属于亚型3和亚型4，其平均生存时间均大于4a，具有较好的预后.预后最好的亚类4和预后最差的亚类1的平均生存时间具有明显差异，前者平均生存时间约是后者的1.6倍，两者的生存曲线如图3(a)所示，两条曲线具有显著差异(log-rank检验，p=0.000 2).

对于验证集中识别的4个亚型，其平均生存时间(单位为月)分别为52.1，34.0，21.2和34.6.预后最好的亚型1和预后最差的亚型3的平均生存时间也具有明显差异，前者平均生存时间约是后者的2.5倍，两者的生存曲线如图3(b)所示，也具有显著差异(log-rank检验，p=0.014).同样的，近一半样本也可被分配到高风险和低风险的类当中.以上分析表明，通过整合多维基因组学改变所富集的通路进行定义的亚型，按其通路活性模式不同指示了不同的生存时间.

2.4亚型的通路活性特征及独立数据集验证

为了验证通路活性的程度能否在一定程度上对亚型的预后情况进行判断，本文对预后最差的亚型提取通路活性特征并验证其在验证集中的一致性.在训练集中，采用RankProduct方法计算预后最差的亚型1与其他亚型之间存在显著差异的通路.在PFDR≤0.05的显著性水平下，亚型1在通路p4，p2，p6和p14上的活性是明显上调的.通路p4，p2，p6和p14涉及到了许多有关肿瘤发生和转移的功能，如细胞粘附，血管生成，转化生长因子β(TGF-β)的结合，和调节细胞增殖等.这也能较好解释亚型1的预后最差的原因.在验证集上，同样地应用RankProduct方法，预后最差的亚型3相对于其它亚型在通路p4，p2和p6上的活性也是显著上调(p14通路在验证集不显著，但调控方向同样是上调的).这说明通路活性和生存时间之间的对应关系是可以重复的，通过少数通路的活性可以高重复性地定义亚型，可以很好地识别出生存不同的亚型其背后较为明确的生物学机制，这一结果对有针对性地治疗患者有积极意义.

进一步比较临床上常用的stage分期亚型和基于通路活性的亚型.根据临床分期信息，得出亚型1中的样本中处于StageIIIC和StageIV的占比分别为78.67%，17.33%，具有比其他亚型样本更高阶段的分期(Wilcoxon秩和检验，p=0.02).StageIIIC和StageIV在临床上认为是已发生肿瘤转移，尤其是StageIV意味着远端转移.通过比较亚型1和数据集中StageIIIC，StageIV的样本平均生存时间，发现亚型1与StageIIIC(平均生存时间为42.8月)内的样本的生存时间存在着显著区别(log-rank检验，p=0.013).亚型1与StageIV(平均生存时间为35.2月)内的样本没有显著差异，然而亚型1内的样本5年生存率仍是低于StageIV的.这表明在分子水平特征上识别的亚型1与临床分期存在着差异，该差异也可作为临床分型的补充.

3结束语

本文利用基因表达谱，结合拷贝数变异、DNA甲基化和miRNA表达等数据信息，得到了卵巢癌生存相关的通路并以通路活性作为特征来进行卵巢癌亚型的分析，最终得到了预后有着差异的亚型.同时，亚型预后与通路活性的对应关系在独立数据验证中呈现一致性.其中，预后最差的亚型表现为TGF-β等通路活性显著上调，文献[3-4]也有相关报道.从预后异质性的卵巢癌患者中识别出具备共同的通路活性特征的卵巢癌预后相关的亚型，对预测患者的预后及治疗方案的选择具有一定的参考价值.

参考文献

[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

基金项目：国家自然科学基金资助项目(61271063)；国家重点基础研究发展计划(973计划)资助项目(2013CB329502)；浙江省自然科学基金资助项目(L215F010001)

作者简介：孟令豪(1990-)，男，山东临沂人，硕士研究生，生物医学工程.通信作者：厉力华教授，E-mail：lilh@hdu.edu.cn.

中图分类号：Q811.4

文献标识码：A

文章编号：1001-9146(2016)04-0029-06