加权基因共表达网络与机器学习算法确定宫颈癌中免疫浸润相关的基因签名与预后分析①

方萌 赵晗池 晴佳（江汉大学医学院，武汉430056）

宫颈癌（cervical cancer，CC）是世界上第四大最常见的女性癌症。近年来，其发病率呈现年轻化趋势［1］。尽管已经证实CC是由HR-HPV（HPV16和18）持续感染引起的［2］。然而，CC的高复发率和异质性导致其预后较差。因此，开发新的诊断和预后生物标记非常迫切。MicroRNA（miRNA）的异常表达被认为是各种癌症发展的关键因素［3］。作为miR‑NA家族的成员，miR-141-3p影响多种癌症进展，并且在CC的发生和发展过程中发挥着重要作用［4-5］。已有研究报道miR-141-3p通过靶向单个基因影响CC细胞的生长和侵袭［6］。但在CC的预后预测方面，基于预后模型的多基因签名比单基因的准确度更高。考虑到免疫微环境对肿瘤发生发展的关键作用，深入签名基因的免疫浸润有可能帮助改善CC的预后预测。先前研究尚未探索miR-141-3p的靶基因与免疫微环境之间的相关性。因此，全面分析miR-141-3p的靶基因与总体生存率之间的相关性，建立可靠的免疫相关的多基因签名以准确预测CC病人的生存至关重要。机器学习在快速、准确的确定基因标记物中起到了重要作用［7-9］。并且，诺谟图比普通的风险评分模型拥有更好的预后预测能力。在以往的研究中，已结合实验和生物信息学方法来研究肝癌与其他疾病的预后标记物与免疫调节分子机制［10-12］。本文通过TCGA数据库分析miR-141-3p在CC中的预后性能。使用WGCNA、COX回归、LASSO COX回归建立并验证了与免疫浸润相关的多基因签名与诺谟图，以预测患者的生存。

1 材料与方法

1.1 数据预处理CC（TCGA数据库中简写为CESC）的mRNA，miRNA表达数据和临床数据可从TCGA数据库（https：//tcga-data.nci.nih.gov/tcga/）下载。表达数据进行log2转换和归一化处理，排除无随访消息或<1 d的时间临床数据。

1.2 miR-141-3p靶基因的筛选使用miRWalk2.0（http：//zmf.umm.uni-heidelberg.de/apps/zmf/mir‑walk2/）中的12个数据库（Microt4，miRWalk，mirbridge，miRanda，miRDB，miRMap，Pictar2，PITA，MiRNAMap，RNAhybrid，RNA22和Targetscan）来确定miR-141-3p的靶基因。然后，使用R语言中的“limma”程序包分析CESC数据。以FDR<0.01，|log2FC|≥2为阈值筛选差异表达基因（DEGs）。最后，通过Venn Plot分析DEGs与预测基因之间的重叠基因（候选基因）。

1.3 加权基因共表达网络分析通过R软件中的“WGCNA”软件包分析候选基因的表达和临床文件。通过计算基因对之间的Pearson相关性，模块特征基因与性状（包括疾病状态）的相关性，确定与预后相关的基因模块。

1.4 基于LASSO COX回归的风险预测模型通过单变量Cox回归分析探讨了每个模块基因对总体生存的影响。将P<0.05和与生存相关的模块基因整合到机器学习算法（LASSO回归）中，以鉴定预后风险特征。使用R包“timeROC”绘制tROC曲线，预测1年、3年和5年总体生存。

1.5 基于支持向量机免疫细胞浸润基于支持向量机的CIBERSORT算法用于计算每个患者中22个免疫细胞的比例。通过wilcoxon检验评估了高危和低危患者之间免疫细胞浸润的差异（P<0.05）。TIMER用于分析签名基因与免疫细胞浸润之间的相关性。

1.6 基因集富集分析通过“clusterprofiler”软件包（Perm=1 000，minGSSize=100和P.value=0.01）分析了高危和低危患者之间的信号通路多样性。“limma”程序包用于计算基因的logFC值。

2 结果

2.1 miR-141-3p的临床价值miR-141-3p在肿瘤组织中的表达显著上调（P=0.003 2，图1A）。ROC曲线表明miR-141-3p的表达可以很好地将肿瘤与正常样品区分开（AUC=0.871 19，图1B）。以303个肿瘤组织中miR-141-3p的平均值作为临界值，将所有患者分为低表达和高表达患者。排除生存时间缺或生存时间为0的情况得到的患者样本为280个。生存分析结果表明，miR-141-3p的低表达可显著改善患者的预后，而高表达对应不良的预后（P=0.022 0，图1C）。tROC曲线分析表明，miR-141-3p对患者的预后具有很强的预测能力。1年、3年和5年的AUC分别为0.762、0.715和0.741（图1D）。

图1 miR-141-3p的临床价值Fig.1 Clinical value of miR-141-3p

2.2 鉴定预后相关基因模块miR-141-3p靶基因包含1 485个交集基因（图2A）。选择软阈值=5以符合无标度网络法则（图2B）。基于基因和样本数据的层次聚类建立的网络热图（图2C）。在9个模块和多个临床表型参数的相关分析结果中，发现黄色模块同时与T.stage（cor=0.28，P=1e-06），N.stage（cor=-0.23，P=5e-05）和AJCC阶段（cor=0.22，P=2e-04）具有很强的相关性（图2D）。因此，确定黄色模块为与CC预后相关的功能性miR-141-3p下游基因集。

图2 WGCNA分析Fig.2 WGCNA analysis

2.3 风险预测模型与诺谟图随机将280例CC患者的样本分为训练集（n=140）和验证集（n=140）。通过单变量COX回归确定了与训练集中的患者OS显著相关的31个基因（P<0.05）。利用LASSO回归进一步优化结果，最后有5个靶标显著影响患者的生存率（表1，图3A）。基于LASSO COX分析的结果，建立了5个靶标签名风险预后模型。模型公式为：RS=0.012 1×FOXA1+0.007 7×DMBX1+0.149 3×TMEM98+0.014 2×RHPN1+0.026 4×SRMS基 于 训练集的生存分析表明，与高风险相比，低风险可以显著改善患者的预后（P=0.000 54）。tROC曲线分析表明，该风险预后模型具有显著的预后预测效果，其1年、3年和5年AUC分别为0.810、0.821和0.679（图3B）。此外，基于验证集K-M曲线分析（P=0.000 86）和tROC分析也证实了上述结果（图3C）。

表1 与OS相关的5个预后靶标Tab.1 5 prognostic targets related to OS

图3 风险预测模型构建与验证Fig.3 Risk prediction model construction and verification

通过单变量和多变量COX回归评估了患者的风险评分和多种临床特征。分析显示，在训练和验证组中，风险评分和AJCC分期与患者的OS显著相关（表2）。通过组合风险评分和AJCC分期而建立的诺谟图（C-index=0.83，图4A）在预测患者生存率方面明显优于单独使用风险评分（C-index=0.77）和AJCC分期（C-index=0.70）。基于训练集的研究发现。1年、3年和5年OS预测表明，由诺谟图预测的生存率与最佳预测性能非常匹配（图4B）。与以上结果一致，诺谟图（C-index=0.68，图4C）与预后因素风险评分（C-index=0.67）和AJCC分期（C-index=0.56）相比，具有更好的预后和预测能力（验证集）。同时，校准曲线还验证了验证集中诺谟图的预后预测性能（图4D）。

表2 预后信息的单因素和多因素分析Tab.2 Univariate and multivariate analyses of prognostic informations

图4 建立预后模型的诺谟图Fig.4 Establishment of nomogram of prognostic model

2.4 低危和高危人群免疫细胞浸润和GSEA通路差异基于CIBERSORT的280个肿瘤样品中22个免疫细胞的浸润率（图5A、B）。CD4+T细胞，巨噬细胞M0和M2是肿瘤微环境中3个最丰富的免疫细胞。有趣的是，巨噬细胞和T细胞的浸润率最大，约占免疫细胞的65%。在TCGA训练集中，高危和低危患者之间存在5种类型（P<0.05）的免疫细胞浸润差异（图5C）。根据TCGA验证集，在不同风险组中存在6种类型（P<0.05）的免疫细胞浸润差异（图5D）。并且通过TIMER获得了预后标志介导的免疫浸润变化（图6）。结果表明，树突状细胞的浸润与5个基因在预后标志物中的表达显著相关。基于全基因组GSEA分析，在高危和低危患者中获得了6条重要的KEGG通路（图7）。

图5 基于CIBERSORT算法分析22种免疫细胞浸润Fig.5 Analysis of 22 immune cells infiltration according to CIBERSORT algorithm

图6 基于TIMER的5种预后基因表达与肿瘤浸润性免疫细胞的相关性分析Fig.6 Correlation analysis between 5 prognostic gene expression and tumor infiltrating immune cells based on TIMER

图7 GSEA获得了高危组和低危组间存在差异的6条通路Fig.7 GSEA obtained 6 pathways with differences be⁃tween high-risk group and low-risk group

3 讨论

首先通过TCGA数据集挖掘了miR-141-3p的临床价值。结果表明，miR-141-3p在CC组织中的表达明显上调，与LI等［6］的研究一致。miR-141-3p高表达水平表明CC患者生存期较差。tROC曲线表明，miR-141-3p的表达对患者的预后具有强大的预测能力。然后，对miR-141-3p的潜在靶标进行了全面分析。WGCNA指出目标基因模块与T.Stage，N.stage和AJCC阶段显著相关。

越来越多的研究表明，基因标记在预测肿瘤的预后中起着重要的作用［13-15］。迄今为止，还未有研究通过将WGCNA和LASSO Cox回归方法相结合以鉴定miR-141-3p的预后靶基因标志物来提高CC预后的预测。在本次研究中，通过对miR-141-3p的候选基因进行了WGCNA和LASSO Cox回归确定了5个与预后相关的靶基因。在以前的报道中，FOXA1、DMBX1、TMEM98、RHPN1和SRMS与CC细胞迁移发展显著相关［16-20］。

作为一种机器学习算法，LASSO COX已有大量研究［21-23］。对于生存数据，最常用的统计模型是Cox比例风险回归模型［24］。而其他传统机器学习算法如随机森林、人工神经网络、贝叶斯算法并没有专门开始针对生存数据的模块。本研究的重点是预测CC患者的总体生存。作为机器学习算法的LAS‑SO COX算法能够识别潜在的风险因素，评估预测、拟合优度，以及结果的临床相关性［25］。故本文采用LASSO COX回归作为机器学习算法。LASSO COX模型的优点在于可以直接解释患病风险和生存之间的关系。另一方面，机器学习技术是无假设和数据自适应，这意味着它们可以有效地用于对复杂数据进行建模。当然，LASSO COX算法也存在一定的局限。例如，该算法要求输入数据进行适当的处理，并且计算过程中会对参数进行适当的调整，以避免算法性能下降［26］。

研究表明，多基因签名在预测肿瘤的预后中起着至关重要的作用。LIANG等［27］提出了一种3-miRNA签名，用于CC的独立预后预测。XIE等［28］还确定了1个8基因签名来预测CC患者的预后。此外，LI等［29］使用组蛋白家族基因签名来预测CC患者的生存。同样，DING等［30］提出了1个3基因签名，以通过Cox回归改善CC的生存预测。先前的研究使用单变量Cox回归分析方法来预测CC患者的生存。但当前对CC预后签名的研究用到的算法不够深入，对签名基因的免疫浸润角色也极少关注。本文通过结合WGCNA、LASSO Cox回归以及支持向量机算法得到的5个基因签名和诺谟图，有效改善了CC预后的预测。WGCNA使用无监督聚类来识别基因模块，默认方法是分层聚类。分层聚类是一种用于确定多维空间中相似数据点聚类的常用方法，它可以使用动态树切割方法来确定模块［31-32］。动态树切割方法成功地识别了使用静态切割方法无法识别的分支，产生预后相关的基因簇。为了结果的准确性，本实验进一步通过单变量COX回归分析基因簇以获得与总体生存相关的基因，利用LASSO COX算法则解决了肝癌从高维数据获得的预测模型可能具有过拟合的风险，构建了预后的特征。基于风险评分和临床特征，通过多变量COX回归建立了诺谟图。而研究结果也证实，诺谟图在预测患者生存率方面明显优于单独使用风险评分（C-index=0.77）和AJCC分期（C-index=0.70）。因而，基于预后签名的风险评分显示出生存预测具备潜在优势。

确定5个基因的特征后，构建5个基因的预后模型并研究其预后价值。发现训练组和验证组的高风险和低风险水平的预后显著不同。所有高危人群的OS均低于低危人群。tROC分析表明5个基因签名具有很强的预测能力。肿瘤的进展不仅受肿瘤本身特征的影响，而且还受肿瘤微环境的影响。值得注意的是本实验的预后标志与肿瘤微环境的免疫浸润水平显著相关。根据CIBERSORT分析，发现在不同的风险组中，B细胞，浆细胞和巨噬细胞M1、M2具有较大的丰度差异。此外，GSEA还揭示了高风险组与低风险组相比有6种不同的免疫状态。因此，研究结果可为CC患者更好的治疗提供依据。

总之，综合运用加权基因共表达网络、COX回归、LASSO回归算法，开发并验证了基于miR-141-3p的预后靶基因的新型5个基因签名（FOXA1、DMBX1、TMEM98、RHPN1和SRMS），并进行了预后分析。结果显示，与免疫浸润相关的5个基因签名可以显著改善CC患者的预后。因此，本次研究提出的5个基因签名与诺谟图有潜力为CC患者提供预后评估工具。