基于m5C 相关基因构建三阴性乳腺癌预后预测模型及药物敏感性分析

安外尔·约麦尔阿卜拉，布尔兰·叶尔肯别克，孙莉莉，刘富中，迪丽娜尔·叶尔夏提，郭文佳

新疆医科大学附属肿瘤医院，乌鲁木齐 830011

三阴性乳腺癌（triple-negative breast cancer，TNBC）是一种缺乏雌激素受体（estrogen receptors，ER）、孕激素受体（progesterone receptors，PR）和人表皮生长因子受体2（human epidermal growth factor receptor 2，HER2）表达的乳腺癌亚型，占所有乳腺癌的10%～20%［1］。TNBC的预后较差，易于转移和复发，治疗难度大，常规内分泌治疗和化疗药物易产生耐药［2］。因此，全面了解TNBC的分子机制将有助于探索更有效的治疗方法。

免疫治疗是一种利用人体免疫系统来攻击癌细胞的治疗方法。在TNBC 中，免疫治疗的主要目标是增强免疫细胞对癌细胞的攻击能力［3］。免疫检查点是一种负向调节机制，如程序性死亡-1（procedural death-1，PD-1）、程序性死亡配体-1（programmed death ligand-1，PD-L1）、狐猴酪氨酸激酶3（lemur tyrosine kinase 3，LMTK3），可以抑制免疫细胞的活性，从而保护正常组织不受免疫攻击［4-5］。然而，肿瘤细胞可以利用免疫检查点来逃避免疫攻击。因此，免疫检查点抑制剂可以通过阻断免疫检查点来激活免疫细胞，从而攻击肿瘤细胞。在TNBC 中，PD-1 可通过与PD-L1 的结合来抑制T 细胞的活性［6］，而PD-1/PD-L1 抑制剂可以通过阻断PD-1/PD-L1 信号通路来激活T 细胞，进而攻击肿瘤细胞［7］。目前，PD-1/PD-L1 抑制剂已被证明对TNBC具有一定的疗效。

5-甲基胞嘧啶（5-methylcytosine，m5C）是一种重要的RNA 修饰，其在转录后的RNA 分子中广泛存在，并参与调控RNA 的稳定性、翻译和功能［8］。近年来，越来越多的研究表明，m5C 修饰与肿瘤的发生和发展密切相关［9］。m5C 甲基化的形成过程由甲基转移酶（也称为“写入器”）如NSUN和DNMT 家族成员催化，并可由去甲基化酶如TET 家族和结合蛋白如YBX1（称为“阅读器”）动态调节［10］，此外，m5C甲基化的异常表达与多种恶性肿瘤的发生发展有关［8］。最近的研究还表明，m5C相关基因表达与肺癌和胰腺癌患者的预后相关［11-12］，而尚未有研究报道m5C 相关基因在TNBC预后中的作用。本研究基于m5C 的相关基因，构建了TNBC 的预后预测模型，并对化疗药物的敏感性进行了分析。同时，进行了免疫浸润分析和高低风险组免疫检查点表达的分析，以期为TNBC的治疗提供新的策略。

1 材料与方法

1.1 数据来源

从GEO公共数据库（https：//www.ncbi.nlm.nih.gov/geo/）获得了GSE76275 和GSE38959 数据集的基因表达谱，用于鉴定差异表达基因以及富集分析。从TCGA（https：//portal.gdc.cancer.gov/）数据库下载了100 个TNBC 样本的RNA-seq 数据以及相应的临床信息作为训练集。选取GSE58812 数据集中107 个样本的基因表达谱以及相应的临床信息用作验证集。详细研究流程见图1。

图1 研究流程图Fig.1 Flowchart of the study

1.2 鉴定GEO数据库的差异表达基因

把GSE76275 数据集的265 个样本和GSE38959数据集的47个样本分为TNBC和非TNBC组，并进行主成分分析，然后用Limma R包分析TNBC和非TNBC组之间的基因表达差异（P＜0.05，|log2FC|＞1），随后用火山图展示了差异分析结果并标注了差异最显著的前10个基因，绘制热图显示了两组之间前30个基因的差异表达谱。

1.3 富集分析

我们对GSE76275 和GSE38959 两个数据集差异表达的196 个基因用Cluster Profiler 包进行了基因本体（gene ontology，GO）和京都基因与基因组百科全书（Kyoto encyclopedia of genes and genomes，KEGG）富集分析。利用提供的分子数值通过GO plot 包计算每个富集条目对应的z-score值，使用ggplot 2包对富集分析结果进行可视化。

1.4 筛选预后有关的m5C相关基因

从先前发表的文献中获得了m5C 甲基化调节基因，包括NSUN2、NSUN3、NSUN4、NSUN5、NSUN6、NSUN7、DNMT2、DNMT3A、DNMT3B、TET2和ALYREF共11 个，并用Pearson 方法对m5C 调节基因和196 个差异表达基因的相关性分析，以P＜0.05、相关系数|r|＞0.03 为筛选条件选出了99 个基因，随后对这99 个基因用Logrank-test方法鉴定了预后相关的5个m5C相关基因。

1.5 预后预测模型的构建

使用Predict 函数基于5 个基因的表达计算了风险评分，并按风险评分中位数把样本分为高风险组和低风险组构建预后预测模型，训练集和验证集模型的预测效果通过ROC 曲线、Kaplan-Meier 生存分析、列线图和校准曲线进行评估。

1.6 免疫浸润分析

根据Cibersort 的基因集和Charoentong 的研究，采用ssGSEA 来定量每个TME 细胞浸润的丰度。为了控制由肿瘤纯度引起的偏差，我们通过使用估计算法计算肿瘤纯度来调整每个肿瘤微环境（tumor microenvironment，TME）细胞亚型的富集分数，共评估了28种人类TME细胞亚型。

1.7 药物敏感性分析

基于癌症药物敏感性基因组学（cancer drug sensitivity genomics，GDSC）（https：//www.cancerRxgene.org）数据库中的岭回归模型，TCGA-TNBC 中某些药物的半最大抑制浓度（half maximal inhibitory concentration，IC50）值通过“pRRophetic”包计算，t检验比较高风险组和低风险组之间某些药物的IC50值差异，进行了药物敏感性预测。

2 结果与分析

2.1 鉴定GEO数据库的差异表达基因

对GSE76275数据集中的198个TNBC和67个非TNBC 样本进行了主成分分析，结果如图2A 所示；对样本进行差异分析发现有668个差异表达基因，结果用火山图（图2C）展示并在图上注释了差异最显著的10个基因，热图（图2E）展示了前30个基因的表达谱；同理，GSE38959数据集中30个TNBC 和17 个非TNBC 样本的主成分分析结果如图2B 所示，发现有1 628 个差异基因，结果用火山图（图2D）展示并在图上注释了差异最显著的10 个基因，热图（图2F）展示了前30个基因的表达谱。

图2 筛选差异表达基因Fig.2 Screening differentially expressed genes

2.2 差异表达基因功能富集分析

对上述2 个数据集差异表达基因使用韦恩图取交集，显示有196个共差异表达基因（图3）。为了探索这些基因在TNBC 中潜在的生物学功能和靶标通路，进行KEGG/GO 分析。气泡图（图4A）、弦图（图4B）、柱状图（图4C）、圈图（图4D）可视化富集分析结果发现，它们主要跟泌尿生殖系统发育、细胞器裂变、腺发育、微绒毛膜、浓缩的染色体、着丝粒区域、染色体、着丝粒区、DNA结合转录激活活性、RNA聚合酶Ⅱ特异性-转录辅激活子结合等生物学过程相关，也与细胞周期、卵母细胞减数分裂等通路相关（表1）。

表1 差异表达基因功能富集特征Table 1 Characters of functional enrichment of differentially expressed genes

图3 GSE76275和GSE38959差异表达基因韦恩图Fig.3 Venn diagram of differentially expressed genes of GSE76275 and GSE38959

图4 差异表达基因功能富集分析Fig.4 Functional enrichment analysis of differentially expressed genes

2.3 预后模型的构建

为了筛出在TNBC 中起到关键作用的m5C相关基因，进行了相关性分析和单因素Cox 分析，结果共筛选出99 个m5C 相关基因（图5），其中LMO4、BCL11A、UGT8、PSAT1、SLC6A14等5 个基因与预后相关（表2）。TCGA-TNBC 作为训练集，GSE58812 作为验证集基于以上5 个基因计算了风险评分，按风险评分中位数分成高风险和低风险2 组，结果发现高风险组患者比低风险组患者的预后更差（P=0.0056、P=0.028）（图6A、B）。风险因子图显示，随着风险值的增加，患者死亡率显著增加（图6C、D）。在验证集中，高风险组和低风险组相比，这5 个基因显著低表达；而在验证集，高风险组中的BCL11A、PSAT1、LMO4显著低表达，UGT8和SLC6A14表达差异没有统计学意义（图6E、F）。

表2 与预后有关的m5C基因Table 2 m5C genes related with prognosis

图5 m5C相关基因的鉴定Fig.5 Identification of m5C-related genes

图6 预后模型的构建Fig.6 Construction of prognostic model

2.4 模型具有良好的预测效能

为验证模型的预测效能，我们构建了列线图（图7A）及校正曲线（图7B），结果显示风险评分和肿瘤分期相关，模型预测效能良好。通过多因素分析和ROC 曲线分析进一步评估了模型的价值，风险评分结果显示是独立的TNBC 预后因子（图7C）。训练集1、3、5 年AUC 值分别为0.85、0.79、0.79（图7D），而验证集中生存时间在一年内的患者不足于计算AUC 值，因此图中显示了3、5年的生存ROC 曲线，AUC 值分别为0.68、0.72（图7E）。

图7 模型预测效能的验证Fig.7 Validation of model prediction efficacy

2.5 不同风险组与免疫浸润之间的相关性分析

为了探索高低风险组中的TME 浸润情况，研究分析了PD-1、LMTK3 等免疫检查点基因在高低风险组中的表达，发现与高风险组相比低风险组的PD-1、LMTK3 呈高表达（P=0.031、P=0.0044）（图8A、B）。基质评分和免疫评分显示，在高低风险组中基质评分没有显著差异（P=0.24）（图8C），免疫评分在高风险组中更高，差异有统计学意义（P=0.022）（图8D）。进一步分析28 种免疫细胞在不同风险组中的浸润情况，结果显示中性粒细胞、17型辅助性T细胞、明亮自然杀伤细胞、γδT 细胞、巨噬细胞、髓源性抑制细胞、单核细胞、T滤泡辅助细胞、活化的CD8T 细胞在高风险组中是高表达的，而2型辅助性T细胞在低风险组中呈现高表达（图8E）。

图8 高低风险组间免疫浸润分析Fig.8 Analysis of immune infiltration between high and low risk groups

2.6 高低风险组与化疗药物敏感性分析

为了探究模型与药物敏感性之间关系，研究用pRRophetic 包预测了常用的化疗药物IC50值在高低分险组中的差异，结果发现多西他赛（docetaxel）、多柔比星（doxorubicin，DOX）、紫杉醇（paclitaxel）等药物IC50值在高低风险组中的差异没有统计学意义（P＞0.05 图9），ARFGAP1 抑制剂（QS11）在低风险组中的IC50值更低，低风险组患者对QS11 更名为敏感（图9）。Akt1/2/3 抑制剂（mk2206）、组蛋白去乙酰化酶（HDACs）抑制剂（ms275）、PI3K 和 mTOR 抑制剂（NVP.BEZ235）、CDK4/CDK6 选择性抑制剂（PD.0332991）等药物在高风险组中的IC50值比低风险组低，说明高风险组患者对这些药物更为敏感（图9）。

图9 高低风险组间化疗药物IC50值差异分析Fig.9 Analysis of differences in IC50 values for chemotherapy drugs between high and low risk groups

3 讨论

乳腺癌是妇女中最常见的癌症类型之一，全球每年估计有150 万新病例，是全世界常见的死亡原因［13］。TNBC 占侵袭性乳腺癌的15%～20%，其主要表现为ER 和PR 的表达缺乏以及HER2的扩增不足［14］，TNBC 比非TNBC 更具侵袭性和增殖性，并且具有更差的预后和存活率［15］，然而TNBC 缺乏确定的治疗靶标（例如ER 和HER2），非特异性化疗仍然是TNBC 患者的主要治疗选择［16］。化疗耐药性是癌症治疗的主要障碍，DOX（一种蒽环类DNA 损伤剂）作为TNBC 的一线方案，在TNBC 患者中重复给药时常常导致耐药［17］。尽管已有报道表明RNA 修饰与疾病发病机制和癌症肿瘤发生有关［18-19］，但是TNBC 和m5C 相关基因之间的潜在关系尚不清楚。本研究通过对GEO 数据集的分析，鉴定了196 个差异表达基因，并对这些基因进行了KEGG/GO 富集分析，结果表明差异表达基因主要涉及到了一些与癌症相关的生物学过程和通路，如DNA 结合转录激活活性、RNA 聚合酶Ⅱ特异性-转录辅激活子结合、细胞周期、卵母细胞减数分裂等。通过相关性分析和单因素分析发现m5C 相关的5个基因在TNBC 中具有重要的生物学意义。进一步，我们使用了TCGA-TNBC 和GEO 数据集作为训练集和验证集，基于m5C 相关的5 个基因构建了TNBC 预后预测模型，此模型可以有效地预测TNBC 患者的预后。我们还进行了免疫浸润分析，发现高风险组和低风险组在免疫细胞浸润方面存在显著差异，高风险组中大部分免疫抑制相关免疫细胞的丰度较高，PD-1、LMTK3 等免疫检查点也在高风险组中呈现高表达，这也可能是其免疫逃逸的原因之一。此外，本研究还使用pRRophetic 包进行了药物敏感性分析，发现高风险组和低风险组对某些化疗药物的敏感性存在差异。这些结果为临床治疗提供了一些有价值的参考。

RNA 甲基化是RNA 转录后调控中最重要的表观遗传修饰之一，也是近年来的研究热点。RNA 甲基化包括m6A、m1A、m5C、m7G、Nm 等几种类型，其中RNA m5C 甲基化是指RNA 的第5个胞嘧啶发生甲基化修饰［20-21］。m5C 甲基化通过影响免疫细胞，从而重塑TME。例如，TET 家族可以影响许多免疫细胞表型的功能，包括B 细胞、浆细胞、树突状细胞和Tregs［22-23］。在这项研究中我们鉴定了5 个与m5C 相关的基因，包括SLC6A14、BCL11A、UGT8、LMO4、PSAT1，构建了预后模型。不同组间的免疫浸润分析结果表明，高风险组与低风险组相比，PD-1、LMTK3 呈高表达。在高风险组中免疫评分更高，中性粒细胞、17型辅助细胞、CD56 明亮的自然杀伤细胞、γδT 细胞、巨噬细胞、髓源性抑制细胞、单核细胞、T 滤泡辅助细胞、活化的CD8T 细胞在高风险组中高表达，而2型辅助细胞在低风险组高表达。

综上所述，本研究的主要贡献在于构建了一个基于m5C 相关基因的TNBC 预后预测模型，并对化疗药物敏感性、免疫浸润和免疫检查点表达等方面进行了深入研究。这些结果为TNBC 的个性化治疗提供了新的思路和方法。但是，本研究也存在一些局限，例如样本量较小、数据来源有限等。未来需要进一步扩大样本量、整合多个数据集进行验证，并结合临床实践进行进一步探索。