基于条件关联互补基因的乳腺癌预后分析

张屹 闫双双 张成 王冠方 黄海云 韩珊珊

摘 要：为了提高乳腺癌患者的生存率，改善病人的临床治疗效果，从分子机制上研究了乳腺癌的致病基因。首先对113个正常组织和1 109个癌症组织的表达量进行差异分析，然后对差异表达的基因采用条件联合分析方式对互补基因进行分组，并用逐步Cox回归挑选出一组基因拟合预后模型。研究结果显示：VWCE，SPDYC，CRYBG3，DEFB1，SEL1L2，NMNAT2 6个基因对患者生存率是有害的，AMZ1，GJB2，CXCL2，ALDOC 4个基因对患者生存率是有利的，最终确定10个基因的预后模型能够显著地将样本分为高风险组和低风险组，并且对乳腺癌患者5年和10年的生存率进行了预测，依赖时间的AUC值均可达0.7以上。所提方法能够利用基因与基因之间的关联性，很好地对高维数据进行降维，消除基因与基因之间的共线性问题，10个基因的预后模型可以对患者的临床预测提供帮助。

关键词：生物数学;乳腺癌;条件关联基因;预后模型;临床预测

中图分类号：Q786 文献标识码：A

文章编号：1008-1542（2020）04-0349-07

doi：10.7535/hbkd.2020yx04008

乳腺癌是女性常见的一种癌症，对女性身心健康有非常严重的影响。世界癌症研究基金会数据显示，2018年，全球有超过200万新增乳腺癌患者[1]。其中有研究表明被诊断为Ⅰ期和Ⅱ期（肿瘤比较小，或只见局部扩散）的乳腺癌患者，5年的生存率大约为80%以上;但对于Ⅲ和Ⅳ期的晚期乳腺癌患者，5年的存活率仅有20%左右[2-3]。可见，对乳腺癌患者的治疗效果还有待提高。手术和药物治疗对乳腺癌患者的治疗有所改善，但是近一半的乳腺癌患者仍会复发，并因复发而导致死亡[4-5]。传统治疗是通过肿瘤大小、分期等风险因素作为乳腺癌患者的风险评估和确定治疗计划的依据。但是这些临床病理风险因素无法将高风险和低风险的乳腺癌患者区分开，亦无法预测哪类患者更能受益于化疗的治疗方案。

隨着高通量技术的不断提高，开始采用测序和基因芯片的数据来查找致癌基因或确立疾病的生物标志物。这种数据的格式为在n个样本中测得m个基因的基因表达量，结果应为n×m维的矩阵，xij表示第i个样本中测的第j个基因的值。针对同一基因在n个样本中的测量数据所组成的向量称为该基因的基因表达向量，可以表示为xi=[x1i，x2i，…，xni]。基因表达数据最明显的特性就是样本少而基因（特征）多，这种高维度的基因表达数据增加了鉴定出与癌症相关基因的挑战难度，严重影响了对疾病的正确预测。而且基因与基因之间也有着复杂的联系，不同基因之间的关联有可能会影响病人的生存期。国内外对如何从众多基因中挑选出与癌症预后相关的基因一直都有研究。LIU等[6]对差异表达的基因，先使用单变量Cox比例回归模型选出与预后显著相关的，再使用逐步多元Cox回归确立了最终的7个基因的预后模型。BIERMANN等[7]在单变量Cox回归显著基因的基础上，将贝叶斯迭代模型平均应用于多变量Cox回归模型，确定了最终18个基因的预后模型。SHI等[8]通过构建基因共表达网络和挖掘与生存相关的模块来鉴定可作为乳腺癌预后标志物的基因。关于最后的预后模型的建立，通常采用基于基本的统计方法Cox回归来选择最后的基因，采用一致性指数（C-Index）来评估模型拟合，用风险评分（线性预测因子）将不同风险组的患者分开，最后通过ROC曲线函数下的时间依赖区域AUC对临床预测能力进行验证[9]。

目前已发现与乳腺癌相关的某些分子标记，多个基因生物标记物对于癌症的预后比单个生物标记物更为准确[10-13]，但当前临床预测和治疗的结果仍不能令人满意。主要的难点在于发现基因与基因之间复杂的联系并有效地减少基因个数，为建立良好的预后模型奠定基础。因此，本文借鉴了YANG 等[14]研究SNP对复杂性状影响中采用的条件联合多SNP分析的部分方法，先寻找互补基因，再对互补基因进行Cox回归，选择与乳腺癌预后相关的基因，确定模型能够很好地将样本分为高风险组和低风险组，通过了对秩数的检验，并对乳腺

癌患者5年和10年的生存率进行了预测，检验了此预后模型的准确性。这种从分子机制上确定疾病的分子标志物对病人治疗的临床结果有了很大改善，研发分子诊断标记物来监测和预测各种癌症的疗效是有实际意义的[15]。

通过确立出的生物标志物能够很好地将不同风险的个体区分开，采取不同的治疗方案，改善治疗效果，减少病痛，延长患者寿命。

1 材料及方法

针对基因表达量的数据特点，提出了基于条件筛选互补基因并进行Cox回归拟合的预后模型，对确定的预后模型进行测试验证，算法过程如图1所示。

1.1 数据的下载与预处理

从TCGA（the cancer genome atlas）数据库中下载乳腺癌的基因表达量数据和对应的临床数据，由于乳腺癌的基因表达量数据包含1 222个样本（其中正常组织样本113个，癌症组织样本1 109个），样本量较大，不能从网页直接下载基因表达量的数据，需下载官方下载工具gdc-client并在终端命令行中下载。将下载后的每个样本中的各个基因的表达量数据进行合并，合并后的基因表达量数据为矩阵形式，行代表的是1 222个样本，列代表的是5万多个基因，基因表达量数据有3种格式，本文中用的是counts值格式的数据。对ENSEMBL格式的基因名字进行转换，转换成只对人类基因命名的Official Gene Symbol格式，并对样本中80%以上表达量为0的低质量基因进行剔除，减少基因个数。整理好基因表达量数据后，对113个正常组织和1 109个癌症组织的表达量进行差异分析。首先对数据进行标准化，然后创建分组矩阵和设计矩阵，查看上调基因和下调基因，筛选出与癌症显著相关的基因，降低维数，以利于下一步模型的计算与构建。其中差异表达基因的筛选标准为log FC>2和FDR<0.05，将满足条件的基因与临床数据合并，其中的临床数据包括生存时间、生存状态、年龄、组织学亚型等。

1.2 寻找多组关联互补基因

对已有的高维基因表达量数据进行差异分析，虽然很大程度降低了维数，但是并不能直接对此时的差异基因进行预后模型估计，此时的基因不能满足Cox回归模型的前提假设，存在多重共线性问题。因此提出了一种新的方式寻找多组互补关联基因。

假设有多个基因对生存时间有线性影响，采用式（1）来表示此模型：

1.3 Cox回归

对筛选出来的每一组基因分别进行多元Cox回归时，并不是一组里的每个基因都是有用的。因此，将数据分为训练集和测试集，在训练集中对每一组数据进行逐步Cox回归，逐步Cox回归采取AIC信息准则作为衡量模型拟合的优劣，AIC越小，说明模型拟合得越好，多个基因的Cox回归模型如式（18）所示：

其中h（t，X）表示的是在时刻t时，m个基因影响下的危险率，h0（t）称为基准危险率。w1是自变量的偏回归系数，一般取RRj=exp（wj）為相对危险度，通过相对危险度来直观地解释Xj的取值是否对h（t，X）的取值产生影响。如果RRj>1，说明第j个基因对癌症样本有着危害影响;如果RRj<1，说明第j个基因对癌症样本有着保护影响;如果RRj=1，说明第j个基因基本不会对癌症样本产生影响。因此，可以通过拟合好的模型提取对乳腺癌病人生存时间有着主要影响的基因。Cox比例风险模型一般采用一致性指数即C-Index来评价模型的预测能力，一致性指数越接近于1，说明模型的预测能力越强，预测的准确率越高。

1.4 风险评分及预后分析

对训练集和测试集的基因表达数据利用确定好的多个基因的Cox回归模型对风险评分进行预测，风险评分表达式如式（19）所示：

预测出的风险评分的中位数为阈值，将每个样本的风险评分与阈值进行比较，大于阈值的样本为高风险组，小于样本阈值的为低风险组。对训练集、测试集和整个数据集分别绘制K-M生存曲线，并用对数秩检验，检验两组患者的风险是否有显著区别，如果检验的p值<0.05，说明通过了显著性检验，选出来的基因确定的Cox回归预测模型能够很好地预测患者的风险。然后再对患者5年和10年的临床预测能力绘制ROC，并计算出AUC值，越接近于1，说明临床预测能力越强，准确率越高。

2 实验结果及分析

2.1 差异分析结果

下载后的基因表达量完成数据的预处理后为1 222个样本和21 487个基因，对113个正常组织和1 109个癌症组织进行差异分析，使用R语言的limma包，将阈值设定为log FC>2和FDR<0.05，对2个条件同时满足的差异表达基因绘制火山图，如图2所示，差异基因中有278个上调基因和868个下调基因。

2.2 预后基因的确定

1）将差异分析后得到的显著基因与临床数据合并，分别进行单基因Cox回归和单基因的一元线性回归分析，选出与乳腺癌患者总生存期之间显著的基因，检验每个基因回归系数值，2个回归中以p<0.05作为阈值，筛选显著基因并进行合并，其中线性回归中显著基因有211个，Cox回归中显著基因有438个，对2种单基因回归分析中显著基因取并集后的612个基因按关联互补方式选择开始进行分组。

2）提取单基因的一元线性回归系数及显著性检验的p值， p值最小的一个作为第一组的第一个入选基因，剩下的基因叫做备选组，计算在入选基因组的基础上再从备选组中选择一个基因加入的条件阈值，如果有p<0.05，选择最小的p值加入，以此类推，加入下一个基因，直至没有p值小于0.05，这样就选出了第一组基因，共4个。

3）重复步骤2），至所有的基因都分了组，共分47个组。通过验证证明了组中基因个数小于15，进行拟合Cox回归，C-Index均小于0.7，所以将组中基因个数小于15的组删掉，剩下15个组。

4）完成数据分组后，对935个样本数据集采用计算机生成随机数的方式划分训练集735个样本和测试集200个样本，对每一组基因分别进行逐步Cox回归，最终选出了第8组的10个基因作为乳腺癌的预后基因。 通过这10个基因确定了乳腺癌的预后模型，Cox回归结果如表1所示。

最终多变量Cox回归模型的10个基因中有6个基因（VWCE，SPDYC，CRYBG3，DEFB1，SEL1L2，NMNAT2）的风险比值均大于1，说明对乳腺癌患者的生存率有着危害的影响，剩下的4个基因（AMZ1，GJB2，CXCL2，ALDOC）的风险比值均小于1，说明对乳腺癌患者的生存率有着有利的影响。对10个基因的回归系数进行显著性检验，8个p值均小于0.05，剩下2个p值小于0.1，由于对多变量模型的预测能力强而被保留在模型中。

2.3 预后模型的测试

通过10个基因在训练集中估算乳腺癌患者的生存风险评分：

最终的计算结果表明，风险评分的中位数为1.013，平均数为1.473，最小值为0.052，最大值为26.788。采用训练集中风险评分的中位数将训练集、测试集和整个数据集中的乳腺癌患者分为高风险组和低风险组，绘制生存风险曲线。训练集中高风险组和低风险组中样本个数分别为367个和368个，测试集中高风险组和低风险组中样本个数分别为88个和112个。K-M分析显示这2组的生存曲线显著不同，与低风险组相比，高风险组的患者表现出更短的生存时间，且通过了对秩数检验，p值都分别小于0.05，如图3所示。

为了评估这10个基因预后模型的准确率，分别对训练集、测试集和整个数据集绘制了时间依赖性的ROC曲线，

如图4所示，在训练集中，5年和10年的AUC值分别为 0.75和0.81。

如图5所示，在测试集中，5年和10年的AUC值分别为0.77和0.78。以上这些结果表明，风险评分是乳腺癌患者临床结果的有力预测指标。

3 结 论

研究并设计出能够区分基因与基因之间复杂关系、更好地降低维数、更为简单高效的癌症预后模型，具有十分重要的理论意义和临床应用价值。本研究使用了一种新的方式来寻找关联互补基因，对每组基因分别进行Cox回归，得到了一组稳健的预后模型，并分别通过了训练集和测试集的验证，通过10个预后基因，可以将高风险组和低风险组显著分开。

通过在GeneCard上查找已经选择出来的基因的分子功能，发现VWCE是β-catenin信号传导途径中的调控元件，也是化学预防细胞癌的靶标。SPDYC通过结合和激活CDK1与CDK2，能够促进整个细胞周期的进展。GJB2表达于耳蜗，在表皮的上基底层和乳腺及子宫内膜的上皮细胞中弱表达。防御素β1（DEFB1）基因与COPD相關，在COPD患者中，人DEFB1的表达上调会影响肺功能的下降[16]。CXCL2是一种趋化因子，炎症介质，趋化因子在协调免疫细胞向炎症或损伤部位的协调募集中起作用[17]。比如其他趋化因子CCL2已被报道能够负面调节管腔B乳腺癌细胞的过程，包括自噬和坏死[18]。有文献证明CXCL2通过抑制ERK1/2信号传导通路，减弱成骨细胞分化[19]。NMNAT2是维持健康轴突的重要生存因素[20]。

总之，确定的10个乳腺癌的预后生物标记物基因被证明具有高预测能力，且会随着时间的推移保持稳定研究结果有助于开发更有效的预后工具，最终改善患者预后。得到的这些标志性基因在其他文献中虽然没有显示出与乳腺癌直接关联，但是这些基因与疾病相关，说明有致病影响，还需要作进一步的实验研究，了解这些基因之间的相互作用，对临床结果进行预测，制定治疗方案，从而为高危乳腺癌患者提供更好的治疗选择，对低危患者减少过度治疗。此方法也可作为筛选其他癌症预后标志物的方法，结果有待进一步测试验证。

参考文献/References：

[1] BRAY F， FERLAY J， SOERJOMATARAM I， et al. Global cancer statistics 2018： GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries[J]. CA： A Cancer Journal for Clinicians， 2018， 68（6）：394-424.

[2] CHEN L， LINDEN H M， ANDERSON B O， et al. Trends in 5-year survival rates among breast cancer patients by hormone receptor status and stage[J]. Breast Cancer Research and Treatment， 2014， 147（3）：609-616.

[3] SAADATMAND S， BRETVELD R， SIESLING S， et al. Influence of tumour stage at breast cancer detection on survival in modern times： Population based study in 173，797 patients[J]. Nederlands Tijdschrift Voor Geneeskunde， 2016，160： A9800.

[4] BREWSTER A M， HORTOBAGYI G N， BROGLIO K R， et al. Residual risk of breast cancer recurrence 5 years after adjuvant therapy[J]. Journal of the National Cancer Institute， 2008， 100（16）：1179-1183.

[5] ADES F， TRYFONIDIS K， ZARDAVAS D. The past and future of breast cancer treatment from the papyrus to individualised treatment approaches[J]. Ecancer Medical Science， 2017.doi：10.3332/ecancer.2017.746.

[6] LIU L， CHEN Z， SHI W， et al. Breast cancer survival prediction using seven prognostic biomarker genes[J]. Oncology Letters， 2019， 18（3）：2907-2916.

[7] BIERMANN J， NEMES S， PARRIS T Z， et al. A novel 18-marker panel predicting clinical outcome in breast cancer[J].Cancer Epidemiology and Prevention Biomarkers， 2017， 26（11）： 1619-1628.

[8] MALLETT S， ROYSTON P， WATERS R， et al. Reporting performance of prognostic models in cancer： A review[J]. BMC Medicine， 2010， 8（1）：21.

[9] HEAGERTY P J， LUMLEY T， PEPE M S. Time-dependent ROC curves for censored survival data and a diagnostic marker[J]. Biometrics， 2000， 56（2）： 337-344.

[10]閆丽娜， 覃婷， 王彤. LASSO 方法在 Cox 回归模型中的应用[J]. 中国卫生统计， 2012， 29（1）：58-60.

YAN Lina， QIN Ting， WANG Tong. The application of LASSO in the Cox model[J]. Chinese Journal of Health Stats， 2012， 29（1）：58-60.

[11]线云开， 孙明立， 于兆进， 等. 基于TCGA数据库筛选乳腺癌不良预后相关mi RNAs及风险评估[J]. 解剖科学进展， 2019， 25（1）：38-40.

[12]孙景波， 陈嘉炜， 王植治， 等. NUF2 基因在乳腺癌中的表达及临床意义[J]. 南方医科大学学报， 2019， 39（5）： 591-597.

SUN Jingbo， CHEN Jiawei， WANG Zhizhi， et al. Expression of NUF2 in breast cancer and its clinical significance[J]. Journal of Southern Medical University， 2019，39（5）： 591-597.

[13]SHI H J， ZHANG L， QU Y J， et al. Prognostic genes of breast cancer revealed by gene coexpression networkanalysis[J]. Oncology Letters， 2017，14（4）： 4535-4542.

[14]YANG J， FERREIRA T， MORRIS A P， et al. Conditional and joint multiple-SNP analysis of GWAS summary statistics identifies additional variants influencing complex traits[J]. Nature Genetics， 2012， 44（4）：369-375.

[15]CH O， WILLIAM C S. Molecular diagnostics for monitoring and predicting therapeutic effect in cancer[J]. Expert Review of Molecular Diagnostics， 2011， 11（1）：9-12.

[16]ELLEN A， GNTHER G， BULLWINKEL J， et al. Increased expression of beta-defensin 1（DEFB1） in chronic obstructive pulmonary disease[J]. Plos One， 2011. doi：10.1371/journal.pone.0021898.

[17]BORO M， BALAJI K N. CXCL1 and CXCL2 regulate NLRP3 inflammasome activation via G-protein-coupled receptor CXCR2[J]. The Journal of Immunology， 2017， 199（5）：1660-1671.

[18]FANG W B， YAO M， JOKAR I， et al. The CCL2 chemokine is a negative regulator of autophagy and necrosis in luminal B breast cancer cells[J]. Breast Cancer Research and Treatment， 2015， 150（2）：309-320.

[19]YANG Y， ZHOU X Y， LI Y J， et al. CXCL2 attenuates osteoblast differentiation by inhibiting the ERK1/2 signaling pathway[J]. Journal of Cell Science， 2019， 132（16）： jcs230490.

[20]GILLEY J， COLEMAN M P， BARRES B A. Endogenous Nmnat2 is an essential survival factor for maintenance of healthy axons[J]. Plos Biology， 2010， 8（1）： e1000300.