基于TCGA数据库构建雌激素受体阳性乳腺癌患者的预后风险模型

杨 燕 贺 鑫 龚鹏举 宋文静 魏 蕾 张京伟

1武汉大学中南医院甲状腺乳腺外科 湖北 武汉 430071；2武汉大学基础医学院病理与病理生理学教研室 湖北 武汉 430071

乳腺癌是女性中最常见的恶性肿瘤，也是女性癌症死亡的主要原因之一［1］。先前的研究表明，雌激素受体（estrogen receptor,ER）阳性乳腺癌约占乳腺癌的70%，雌激素受体α（ERα）在乳腺癌的发生和发展过程中起重要作用［2-4］。ERα是一种核受体，由几个结构域组成：募集转录调节因子的N端激活功能域1（AF1）和C端激活功能域2（AF2），与靶基因的雌激素反应元件结合的DNA结合域（DBD）［5］。ERα的活性主要受雌激素调节，ERα的激活促进了许多调控致癌过程的靶基因的转录，包括细胞周期、细胞增殖和上皮/管腔细胞分化等过程［6，7］。内分泌治疗通过拮抗雌激素与ERα的结合，下调ERα或抑制雌激素的产生来抑制ERα信号通路的激活，结合手术和放射疗法，显著改善了ER阳性乳腺癌患者的临床结局［8］。然而，仍有部分ER阳性的乳腺癌患者在内分泌治疗过程中产生获得性耐药，出现复发或转移［9-11］。因此，迫切需要可以预测ER阳性乳腺癌患者预后的生物标志物。

微小RNA（miRNA）是一类调节基因表达的短非编码单链RNA，长度约为20～24个核苷酸，可通过与靶基因的3'-非翻译区（3'-UTR）碱基配对抑制信使RNA（mRNA）的翻译或降解mRNA从而调节靶基因的表达［12］。许多证据表明，miRNA在细胞增殖分化、血管生成和代谢等各种生物过程发挥重要作用［13-15］。多种恶性肿瘤发生发展也与miRNA的表达失调相关［16，17］，例如结肠癌［18，19］、肺癌［20］、胃癌［21］、鼻咽癌［22，23］等。有研究显示，乳腺癌中也存在多种miRNA表达失调，且部分与ERα状态有关［24］。例如，miRNA-519a的高表达与ER阳性乳腺癌患者较差的生存率相关，且miRNA-519a上调可导致ER阳性乳腺癌患者对他莫昔芬的耐药。

因此，本研究旨在通过分析TCGA数据库的miRNA测序数据构建miRNA-临床的预测ER阳性乳腺癌患者预后的模型，用于有效选择高危患者和预测ER阳性乳腺癌患者的3年和5年生存率，从而有利于制定更有效的个体化治疗决策。

1 资料与方法

1.1 下载与处理原始数据2019年11月15日从TCGA官方网站获取乳腺癌患者miRNA表达谱和临床数据，miRNA表达谱数据中总共包括1 096例乳腺癌组织和104例正常乳腺组织，临床数据中包括1 174例乳腺癌患者临床病理信息。纳入标准：（1）病理学确诊为乳腺癌。（2）既有miRNA的样品测序数据又有完整的临床预后资料（生存状态及生存时间）。排除标准：（1）雌激素受体表达情况不明确。（2）生存时间未知或生存时间≤0 d。最后共纳入730例乳腺癌患者，其中555例雌激素受体阳性患者，175例雌激素受体阴性患者。从GEO中下载数据集GSE37405（GPL 13703），该数据集包含60例ER阳性乳腺癌患者的miRNA测序数据，患者均接受他莫昔芬治疗，其中乳腺癌复发患者和未复发患者各30例。

1.2 筛选ER阳性乳腺癌特异性差异miRNA以|log2FC|≥2，FDR＜0.05为筛选条件，使用R3.6.0软件中的edgeR包分别筛选出ER阳性乳腺癌差异表达的miRNA以及ER阴性乳腺癌中差异表达的miRNA，并使用ggplot2和pheatmap包绘制火山图和热图。使用VennDiagram包绘制韦恩图筛选出ER阳性乳腺癌特异性上调和下调的差异miRNA。

1.3 构建风险评分公式和miRNA-临床预后模型筛选得到的ER阳性乳腺癌特异性差异miRNA进行Lasso回归分析，筛选出与生存相关miRNA，根据每个miRNA的回归系数，构建基于miRNA表达的风险评分公式。绘制受试者工作曲线（ROC）评估miRNA风险评分的预测性能，根据ROC曲线选择最佳临界风险评分将ER阳性乳腺癌患者分为低风险组和高风险组，比较两组患者之间生存的差异及miRNA的表达情况。进一步通过单因素和多因素Cox回归分析结合临床病理因素和miRNA风险评分构建预测ER阳性乳腺癌患者预后的预后模型，并使用rms包绘制列线图。

1.4 评估风险评分公式和miRNA-临床预后模型的预测性能绘制风险评分公式和miRNA-临床预后模型的ROC曲线并计算曲线下面积（AUC）评估模型的预测能力。绘制校准曲线评估模型的准确性。

1.5 靶基因预测和功能富集分析通过Starbase预测miRNA的候选靶基因。选择至少3个数据库重叠的基因作为miRNA的靶基因。使用Cytoscape可视化miRNA与靶基因之间的相互作用网络。最后，使用DAVID 6.8进行靶基因的GO功能富集和KEGG通路富集分析。

1.6 统计学分析正态分布的数值变量采用均数±标准差描述，t检验用于组间比较。偏态分布的数值变量采用中位数（四分位数间距）描述，Mann-WhitneyU检验用于组间比较。分类变量采用频数（构成比）描述，χ2检验或Fisher精确检验用于组间比较。Lasso回归分析用于筛选与预后相关miRNA并构建风险评分公式。使用单因素Cox回归和多因素Cox回归分析构建miRNA-临床预测模型，并使用R中rms包绘制列线图。survival ROC包绘制ROC曲线评估了列线图预测能力，校准曲线用于评估模型的准确性。使用Kaplan-Meier方法绘制生存曲线，并使用Log-rank检验比较生存率。P＜0.05被认为具有统计学意义。应用R3.6.0进行数据分析。

2 结果

2.1 患者的基线特征共555名ER阳性乳腺癌患者纳入研究，用随机抽样法按1∶1比例分为训练集278例和验证集277例，训练集用于构建预测模型，验证集用于验证预测模型效能。两组基本情况比较差异均无统计学意义（P＞0.05），两组具有可比性。所有患者的5年总生存率为91.7%。

2.2 筛选ER阳性乳腺癌与生存相关的特异性差异miRNA从TCGA数据库中获得了乳腺癌的miRNA测序数据，包括555例ER阳性乳腺癌，175例ER阴性乳腺癌和104正常组织样本。使用R3.6.0软件中的edgeR包分别筛选出ER阳性乳腺癌差异表达的miRNA和ER阴性乳腺癌差异表达的miRNA（|log2FC|≥2，错误发现率FDR＜0.05）。根据两组上调的miRNA和下调的miRNA分别绘制韦恩图，选择ER阳性乳腺癌非交集部分miRNA为ER阳性乳腺癌特异性miRNA，其中5个miRNA在ER阳性乳腺癌特异性上调，17个miRNA在ER阳性乳腺癌特异性下调。为了进一步筛选与生存相关的miRNA，在训练集中对22个miRNA进行了Lasso回归，使用交叉验证以建立风险模型，最终筛选出4个与生存相关的miRNA，分别为miR-331、miR-615、miR-653、miR-887（见图1）。

图1 ER阳性与ER阴性乳腺癌下调的miRNA的韦恩图（A）和上调的miRNA的韦恩图（B）以及Lasso回归分析筛选miRNA(C)和交叉验证结果(D)

2.3 风险评分公式的建立和miRNA-临床预后模型的构建根据Lasso回归得到相应的系数，建立立风险评分公式：风险评分=0.098×miR-331exp+0.122×miR-615exp+0.102×miR-653exp+0.113×miR-887exp。根据此方程，计算每个患者的风险评分，并根据ROC曲线的最佳临界风险评分（风险得分为1.757）将训练集患者分为低风险组（n=189）的和高风险组（n=89），验证集也分为低风险组（n=73）和高风险组（n=204）。Kaplan-Meier生存分析表明，在训练集（P=0.000 26）和验证集（P＜0.000 1）中，高风险组的预后比低风险组差（图2）。随后，将风险评分作为变量，结合其他临床病理特征进行单变量和多变量Cox回归分析来识别与预后相关的危险因素，结果显示TNM分期（HR=5.068，95%CI：1.799～14.278，P=0.002）、风险评分（HR=1.191，95%CI：1.030～1.377，P=0.018）和术后是否放疗（HR=0.411，95%CI：0.176～0.965，P=0.041）是预后的独立危险因素（表1）。将上述危险因素用于构建miRNA-临床预后模型并绘制列线图，在列线图中，每个变量的得分相加计算总分，通过总分评估患者的3年和5年生存率（图3）。

图2 训练集（A）和验证集（B）中高低风险评分组患者的生存曲线

表1 影响ER阳性乳腺癌预后的临床病理因素的单因素和多因素Cox回归分析

图3 基于4个miRNA表达的预测ER阳性乳腺癌患者3年和5年生存率的列线图

2.4 评估风险评分和miRNA-临床预后模型的预测效能绘制风险评分和基于miRNA-临床预后模型的ROC曲线并计算AUC值以评估模型的预测能力。在训练集和验证集中，miRNA-临床预后模型预测3年生存率的AUC分别为0.768和0.909，预测5年生存率的AUC分别为0.849和0.860（图4）。

图4 训练集（A）和验证集（B）评估miRNA-临床预后模型3年和5年预测能力的ROC曲线

2.5 miRNA在他莫昔芬治疗患者中的表达从基因公共表达数据库（GEO）中下载数据集GSE37405（GPL 13703），该数据集包含30例他莫昔芬治疗后未复发患者及30例他莫昔芬治疗后复发患者的miRNA测序数据，分析4个miRNA在两组间的表达水平，结果显示，与未复发组相比，复发患者miR-331的表达水平显著升高，而两组之间的miR-615、miR-653和miR-887的水平没有明显差异（图5）。miR-331可能和ER阳性乳腺癌他莫昔芬耐药相关。

图5 他莫昔芬治疗后未复发组和复发组中4个miRNA的表达水平

2.6 mi-331的靶基因的GO功能富集分析和KEGG通路富集分析GO功能富集分析表明，miR-331的靶基因细胞组成（图6A）主要富集在细胞核和核质，生物过程（图6B）主要参与在RNA聚合酶Ⅱ转录的终止、核转录mRNA poly（A）尾巴的缩短和内皮细胞增殖的调控，分子功能（图6C）主要参与蛋白质结合。KEGG通路分析（图6D）结果表明，miR-331的靶基因主要参与mRNA监视途径，Rap1信号通路和MAPK信号通路相关。

图6 miR-331的GO功能富集和KEGG通路富集

3 讨论

高通量测序结合生物信息学对基因组学数据进行分析，可以从分子角度探索与恶性肿瘤诊断，治疗与预后密切相关的标志物。通过分子标志物对患者生存率的预测有助于制定个体化的临床决策。许多研究［25-29］表明，miRNA与ER阳性乳腺癌的发生和进展有关。Mulrane等［30］研究发现miR-187在乳腺癌中的高表达可能与更具侵略性的表型相关，且miR-187与乳腺癌预后不良相关，是预测患者生存率的独立预测因子。Okuda等［31］研究表明，miR-7可以通过下调KLF4基因的表达进而抑制乳腺癌细胞的远处转移。有研究者［32］发现ER阳性患者的肿瘤中miRNA-375水平显著高于ER阴性乳腺癌患者，且miRNA-375在他莫昔芬耐药细胞中下调，其重新表达可以抑制MTDH基因表达从而恢复肿瘤细胞对他莫昔芬敏感性，同时也抑制肿瘤细胞侵袭和逆转EMT样特性。此外，Gong等［33］建立了一个基于miRNA的模型来预测激素受体阳性HER2阴性乳腺癌患者的复发风险。因此，在多种癌症中基于miRNA的预测模型对预测患者的生存或复发具有重要意义。传统的生物学方法验证miRNA和患者的预后之间的关联既耗时又昂贵，随着生物学数据的积累，我们可以利用大量数据分析miRNA和肿瘤的相关性并开发预测模型来预测肿瘤患者的预后及治疗反应，通过另一数据集的验证检验模型的可靠性，建立基于miRNA的预测模型既对个体化治疗方案的选择具有重要的参考价值。也可为miRNA与肿瘤的相关性进一步的实验验证提供理论依据。

在本研究中，我们从TCGA数据库下载miRNA表达数据，通过差异表达分析和Lasso回归发现了4种ER阳性乳腺癌特异性表达上调的miRNA（miR-331、miR-615、miR-653、miR-887）并建立了miRNA表达的风险评分公式。根据ROC曲线的最佳临界风险评分将患者分为低风险组和高风险组，风险评分与ER阳性乳腺癌患者的总生存率显著相关。通过Cox回归分析风险评分和临床病理因素与预后的关系，构建了miRNA-临床预测模型用于预测ER阳性乳腺癌患者的预后。该模型预测3年和5年生存率的AUC值分别为0.768和0.849，表明其具有良好的预测性能。模型中包含的4个miRNA中，部分miRNA已有试验证实与乳腺癌相关。Jiang等［34］发现miR-331在乳腺癌细胞中上调，从而促进乳腺癌细胞的增殖、迁移和侵袭，miR-331的过表达与淋巴结转移、TNM分期和预后不良有关。Lei等［35］发现miR-615-3p在乳腺癌细胞和组织中特别是在转移的乳腺癌细胞和组织中显著升高，miR-615-3p通过抑制PICK 1基因表达促进了肿瘤细胞的EMT过程。因此，我们构建的miRNA-临床的列线图可能是ER阳性乳腺癌患者生存预测的重要工具，有助于制定个性化治疗策略。随后我们利用GEO数据库中GSE37405数据集分析了4个miRNA在他莫昔芬治疗后复发和未复发两组之间的表达，发现miR-331在他莫昔芬治疗后复发组中表达明显高于未复发组，提示miR-331可能与ER阳性乳腺癌患者发生他莫昔芬耐药相关。进一步进行了miR-331靶基因的GO和KEGG途径的分析，结果显示miR-331可能在蛋白质结合、转录、mRNA监视通路、MAPK信号通路中起关键作用。

综上所述，我们构建的预测ER阳性乳腺癌患者生存的miRNA-临床风险模型，该模型结合了4个miRNA的表达情况和临床危险因素，可预测ER阳性乳腺癌患者的3年和5年生存率，且建立的风险评分公式可有效地识别ER阳性乳腺癌中的高危患者。因此，该预测模型对制定个性化治疗决策有重要意义。但在该模型广泛应用于临床实践之前，仍需要进行多中心、大规模、前瞻性研究以验证该预测工具。