基于加权基因共表达网络分析建立结肠腺癌预后风险模型以及关键基因药物敏感性分析

2023-04-01黎蕴琪蒋利和陈闯广西医科大学附属肿瘤医院南宁53002右江民族医学院广西百色533000广西大学南宁530004

中南药学 2023年2期

黎蕴琪，蒋利和，陈闯*（. 广西医科大学附属肿瘤医院，南宁 53002；2. 右江民族医学院，广西百色533000；3. 广西大学，南宁 530004）

结直肠癌（colorectal cancer，CRC）是世界上最常见的消化系统恶性肿瘤之一。据统计，在全球范围内，CRC发病率排第三，病死率排第二[1]。结肠腺癌（colorectal adenocarcinoma，COAD）是CRC中最常见的病理分型。目前，根治性手术治疗是早期COAD最主要的治疗方法，而对中晚期患者来说，化疗是最常用且有效的选择，可限制肿瘤细胞的扩散和转移，进而控制病情。但化疗药物耐药性的产生和肿瘤的高复发率仍是COAD临床实践中的主要问题[2-3]。

近年来，生物信息学在挖掘癌症生物标志物[4]和抗肿瘤药物分子靶点[5]等研究方面发挥重要作用。加权基因共表达网络分析（WGCNA）主要运用于分析大样本基因表达数据，将基因根据表达模式相似性分为不同模块[6]，这种方法跟普通聚类方法的不同在于其将基因表达值的相关系数取了N次幂，从而使得相关系数分布更加符合无尺度网络分析，这更符合生物学规律[7]。相比于只关注差异表达的基因，WGCNA充分利用信息，把基因与表型的关联转换为基因集与表型的关联，免去了多重假设检验校正的问题。本研究基于WGCNA结合差异表达基因（DEGs）的方法，构建高性能的预后预测模型，同时对这些靶点进行化疗药物敏感性分析，旨在筛选出高风险肿瘤患者并对其进行精准的靶向治疗，以达到个体COAD治疗疗效最大化。

1 材料和方法

1.1 数据采集

微阵列数据集GSE39582和GSE41258以及这些数据集的相应临床数据从GEO数据库（https：//www.ncbi.nlm.nih.gov/geo）中下载，通过R软件进行数据归一化处理，利用limma包去除批次效应，通过ggord包绘制PCA图。GEO数据用于筛选DEGs，构建预后模型。从癌症基因组图谱（TCGA）数据库（http：//cancergenome.nih.gov/）下载COAD的RNA-seq数据和相应临床信息。该数据集包括424个COAD样本，用于验证预后模型和生存分析。

1.2 差异表达基因分析

采用limma软件包寻找mRNA的差异表达，并通过pheatmap包绘制热图，“AdjustedP＜0.05且|log2FC|＞1”被定义为阈值，得到的DEGs被用于后续分析。

1.3 WGCNA

在评估GSE39582和GSE41258的基因表达谱后，使用WGCNA包构建无标度的共表达网络，计算基因的皮尔逊相关系数以获得相关矩阵，选择适当的软阈值来测量基因之间的连接性，将邻接矩阵转换为拓扑重叠矩阵（TOM），并计算相应的差异（1-TOM）。具有相似表达模式的基因被分为同一个共表达模块。此外，我们计算了模块特征基因的差异性，对模块进行了层次聚类，并合并了相似的模块（abline＝0.25）。

1.4 预后风险模型的构建

将DEGs和关键模块基因取交集，得到共同基因用于进一步分析。将这些基因的表达谱和临床信息合并，再随机分为训练集（70%）和测试集（30%），使用训练集建立模型并在测试集进行验证。在训练集中，通过单变量Cox比例风险回归分析筛选与预后相关的基因，纳入Lasso回归分析，防止模型过度拟合，然后纳入多因素Cox回归分析，得出COAD预后风险模型，计算每个患者的风险评分，依据训练集风险评分的中位值，分别将训练集和测试集分为高风险组和低风险组。运用Kaplan-Meier生存分析和ROC曲线评估该预后风险模型的准确性。利用生存时间和基因风险模型分别绘制散点图和高低风险热图，利用风险评分结合临床信息绘制列线图。

1.5 基因集富集分析（GSEA）

为了进一步确认预后关键基因的潜在功能和识别高、低风险组的相关通路，使用GSEA软件（4.0.1）进行基因组富集分析，以P＜0.05为阈值，展示KEGG相关通路的富集结果。

1.6 药物敏感相关性分析

CellMiner是以美国国家癌症研究所癌症研究中心（NCI）所列出的60种癌细胞为基础建立的数据库[8]。从CellMiner数据库下载NCI-60中基因mRNA表达数据和药物活性相关数据，探究5个关键基因的表达与药物敏感性之间的关系。此外，通过数据库GSE81005分析5-氟尿嘧啶（5-FU）耐药的结直肠癌细胞系和野生型结直肠癌细胞系中HSD17B2mRNA的表达高低。

1.7 统计学分析

使用Wilcoxon检验分析连续变量。使用Kaplan-Meier 生存分析和 Cox 单因素分析进行预后分析。Pearson相关性分析用于评估5个关键基因的表达与药物敏感性之间的相关性。使用R软件（3.6.3版本）进行统计分析，P＜0.05为差异有统计学意义。

2 结果

2.1 数据处理

GSE39582和GSE41258数据集共包含73个正常样本和752个COAD组织样本。从这些样本中鉴定出490个DEGs，其中上调基因313个，下调基因177个（见图1A、B）。数据批次效应情况通过批次去除前后的PCA图进行评估（见图1C、D）。

图1 基因表达差异分析结果Fig 1 Gene expression differential analysis

2.2 WGCNA和关键模块的识别

为了找出与COAD相关的关键基因，笔者使用正常样本和COAD样本的基因表达矩阵进行了WGCNA分析。将软阈值β设置为5，平均连通性接近0（见图2C、D）。将具有相似表达模式的DEGs聚集到相同的模块中，最终得到11个共表达模块（见图2B），灰色模块被认为是无法被分配给任何模块的基因集合。绿松石色模块与正常样本表型相关系数为0.64，绿色模块与肿瘤样本表型相关系数为0.53，通过分析基因与模块的相关性（MM）和基因与临床特征（正常和肿瘤）的相关性（GS）之间的关联，发现绿色和绿松石色模块中MM和GS成正相关（见图2E、F），说明这些与性状高度相关的基因，在关键模块中也扮演着举足轻重的角色。因此选择绿色和绿松色为目的模块，并将DEGs和目的模块基因取交集筛选出247个基因用于进一步分析。

图2 WGCNA分析Fig 2 WGCNA analysis

2.3 预后模型的构建

将表达和生存数据合并后的740个GEO数据集样本分为训练集（70%）和测试集（30%）。首先使用训练集的生存数据对247个DEGs进行单因素Cox比例风险回归分析，鉴定出7个对预后有显著影响的基因（P＜0.01），分别为GZMB、HSD17B2、PDK4、PIGR、PXMP2、GINS1、EMP1。将上述基因纳入Lasso回归分析降维（见图3A、B），进一步使用多变量Cox比例风险回归分析，最终得到5个与COAD预后相关的风险基因（见图3C）。基于这5个关键基因进行风险预后模型的构建：Risk score＝-0.162GZMB+0.210HSD17B2+0.184PDK4-0.141PIGR-0.516PXMP2。生存分析结果均显示低风险组的患者生存状况明显优于高风险组（P＜0.001，见图4A、B）。训练集1年、2年和3年的随访预测AUC分别为0.690、0.709、0.699（见图4C），测试集的随访预测AUC，见图4D。从生存时间和风险评分绘制的热图和散点图中看出，随着风险得分的增加，死亡的患者也增加，存活时间相对减少，由此可见模型有相对较好的预测能力（见图5）。

图3 LASSO及多变量Cox回归分析进一步筛选预后基因Fig 3 Further screening of prognostic genes by the LASSO and multivariate Cox regression analysis

图4 训练集（A、C）、测试集（B、D）Kaplan-Meier曲线和时间依赖性ROC曲线Fig 4 Distribution of time-dependent ROC curves and Kaplan-Meier curves of the training set（A and C）and the test set（B and D）

图5 训练集（A）、测试集（B）风险评分的分布、患者的生存状态和5个基因的表达热图Fig 5 Distribution of risk scores，patient survival status，and the five-gene expression heat map of the training set（A）and the test set（B）

2.4 预后风险模型的临床应用

进一步通过单因素和多因素Cox回归分析评估COAD患者风险评分和临床病理信息与总生存期的关系。训练集和测试集的结果均显示风险评分可以作为独立危险因素来预测COAD患者的生存预后（P＜0.01），多因素Cox回归分析结果显示，风险评分在训练集（见图6B）和测试集（见图6D）中的HR分别为1.581和1.616。

图6 训练集（A、B）、测试集（C、D）关于临床特征因素和风险评分的独立预后分析Fig 6 Independent prognostic analysis of clinical characteristics and risk score of both the training set（A and B）and the test set（C and D）

2.5 模型验证

接下来，在TCGA验证集（424例样本）中使用上述模型公式计算每位患者的风险评分，将患者分为高风险和低风险组，并通过KM曲线比较预后。结果表明，该模型可以区分患者的预后（P＝0.013，见图7）。

图7 TCGA验证集的生存分析Fig 7 Survival on TCGA validation set

2.6 预后列线图的建立

笔者使用临床病理特征（年龄、性别、肿瘤分级、T分期、N分期）和5个基因的表达在GEO测试集具有完整临床资料的503例样本中建立了建立列线图来预测患者的生存率。结果显示列线图可以预测1年、2年和3年的生存率，接近理想模型（见图8A）。校准图显示出较佳的预测准确性（见图8B），并且预测的生存率与实际生存率基本重合。

图8 预后风险模型列线图（A）及校准曲线图（B）Fig 8 Nomogram（A）and calibration（B）curve of the prognostic risk model

2.7 高风险和低风险组的DEGs功能富集分析

为了探索5个关键基因与COAD的潜在生物学机制，笔者使用GSEA在训练集和测试集中对高、低风险组之间的DEGs进行了功能富集分析。训练集和测试集结果均显示低风险组的基因集在代谢以及细胞周期等相关途径有明显的富集作用（见图9），提示低风险组可能通过细胞周期通路影响肿瘤细胞凋亡进而影响患者生存。

图9 训练集（A）、测试集（B）GSEA富集分析Fig 9 GSEA enrichment analysis of the training set（A）and the test set（B）

2.8 药物敏感性分析

最后，笔者研究了5个关键基因在60种人类癌细胞系（NCI-60）中的表达水平及其与化疗药物之间的相关性。这些基因和多种肿瘤的化疗药物的耐药性增加有关（P＜0.05），其中HSD17B2mRNA的表达增加与多种治疗COAD的化疗药物的耐药性增加有关（见图10A～E），包括伊立替康（cor＝-0.262，P＝0.043，中晚期CRC的一线用药[9]）、5-氟脱氧尿苷（cor＝-0.318，P＝0.013）、氟尿苷（cor＝-0.334，P＝0.009，治疗结肠癌的肝转移[10]）、依托泊苷（cor＝-0.308，P＝0.017）、雷替曲塞（cor＝-0.298，P＝0.021）。另外，通过数据库GSE81005对5-FU耐药的CRC细胞系和野生型CRC细胞系进行分析，发现5-FU耐药的CRC细胞系中HSD17B2mRNA的表达明显高于野生型CRC细胞系（P＜0.001，见图10F）。这些结果表明HSD17B2在CRC的发生发展和耐药性中发挥着重要生物学作用。

图10 HSD17B2的表达与伊立替康（A）、5-氟脱氧尿苷（B）、氟尿苷（C）、依托泊苷（D）和雷替曲塞（E）的药物敏感性的相关性分析以及HSD17B2在CRC野生型细胞及其5-FU诱导耐药细胞系中的表达（F）Fig 10 Correlation between HSD17B2 expression and drug sensitivity of irinotecan（A），5-fluorodeoxyuridine（B），fluorouridine（C），etoposide（D）and raltitrexed（E）； expression of HSD17B2 in CRC wild type cells and its 5-FU-induced resistant cell lines（F）

3 讨论

COAD是全球范围内最常见的癌症类型之一，恶性程度和病死率较高[1，11]。目前COAD的治疗主要采取手术切除及化疗等传统治疗方法，对于早期肿瘤效果较好，但对于中、晚期肿瘤只能延长患者生命，缓解临床症状，不能从根本上提高患者的生活质量，患者预后仍较差。因此，寻找COAD中用于生存评估和靶向治疗的生物标志物仍至关重要。

本研究通过分析两个GEO数据集来鉴定DEGs，构建了WGCNA共表达网络，从中找出5个预后关键基因（GZMB、HSD17B2、PDK4、PIGR和PXMP2），构建COAD预后预测模型，并对其进行了验证，列线图和ROC曲线显示该模型具有较好的临床实用潜力。药物敏感性的分析显示，这些关键基因与多种化疗药物的耐药性相关。先前的研究表明，颗粒酶（GZMs）是细胞毒性T淋巴细胞和自然杀伤细胞的关键效应分子，在控制细胞内病原体和癌症方面发挥着关键作用[12-13]。颗粒酶B（GZMB）作为抗肿瘤和抗感染因子，在多种肿瘤中显著表达并与良好预后相关，如非小细胞肺癌、CRC、骨肉瘤等[14-16]。17-β羟基类固醇脱氢酶2型（HSD17B2）是一种参与雌激素和雄激素调节的酶[17]，其可催化睾酮和雌二醇（E2）分别转化为雄烯二酮和雌酮（E1）[18]。而雌二醇被报道可减少CRC细胞的增殖和促进凋亡[19]。本研究构建的风险模型也提示，随着HSD17B2的表达增加，COAD患者的预后更差，这可能与雌二醇的失活相关。Lee等[20]发现HSD17B2可以作为独立预测因子，其过表达与直肠癌患者的不良预后相关。丙酮酸脱氢酶激酶4（PDK4）是细胞能量代谢的关键调节剂，可以参与m6A调节的糖酵解和ATP生成[21]。研究表明PDK4的高表达可以促进肿瘤细胞的增殖、迁移和侵袭，例如CRC、宫颈癌、肝癌、膀胱癌细胞等[21-23]。Woolbright等[23]发现，抑制PDK4表达可能会加剧顺铂诱导的细胞死亡。聚合免疫球蛋白受体（PIGR）是黏膜免疫系统的关键组成部分，在非小细胞肺癌、鼻咽癌、上皮性卵巢癌等的表达降低[24-27]。PIGR被证实了和消化道肿瘤密切相关。研究发现，PIGR可以启动Yes-Dap12-Syk-Rac1/CDC42-MEK/ERK信号级联以促进肿瘤细胞的增殖和转化[28]，富含PIGR的细胞外囊泡可以促进肝癌细胞的侵袭性[29]。Ohkuma等[30]发现，PIGR的mRNA和蛋白水平是胰腺癌的独立预后因素，其高表达与患者预后不良有关。PXMP2是高等真核生物中最丰富的一种同源三聚体过氧化物酶体膜蛋白（PMP），具有通道形成活性[31-32]，参与广泛的代谢途径，包括脂质、氨基酸和羟基酸、嘌呤和活性氧物质的转化[33]。但PXMP2基因与人类疾病的关系尚未被充分了解。

对高风险组和低风险组进行GSEA功能富集分析，发现低风险组在细胞周期、DNA复制、RNA聚合酶和代谢相关通路显著富集。细胞周期失调和细胞周期蛋白依赖性激酶（CDK）激活导致的持续细胞增殖是癌症的标志[34]。目前，已经有3种选择性CDK4/6抑制剂获得美国食品药品监督管理局的批准[35]，多种CDK4/6抑制剂正处于治疗癌症的临床试验中[36]。此外，通过药物敏感性分析发现这5个关键基因与多种化疗药物相关，特别是HSD17B2与结肠癌化疗药物（伊立替康、5-氟脱氧尿苷、依托泊苷、雷替曲塞等）存在显著负相关。此外，5-FU耐药的细胞系中的HSD17B2表达显著高于野生组。5-FU是CRC治疗中使用最广泛的化疗药物之一。尽管近年越来越多的化疗药物出现，但以5-FU为基础的化疗（如 XELOX 和 FOLFOX）仍然是目前中晚期 CRC 的一线化疗药物。但由于化疗耐药性，其总体反应率仅为 10%～15%[37]，化疗耐药和严重的不良反应已成为影响化疗治疗效果的主要因素[38]。为了克服化学耐药性和减少不良反应，许多药物在临床试验中进行了测试，但迄今为止尚未取得重大进展，需要更多努力寻找能够克服5-FU耐药性且毒性较小的新药以改善CRC患者。本研究提示HSD17B2可能是5-FU化疗耐药的潜在生物标志物，针对HSD17B2的靶向治疗或许能在提高化疗敏感性的同时改善患者的预后。总之，预后模型中的5个靶基因可能通过影响细胞周期和细胞对化疗药物的敏感性，进而影响COAD的进展和预后。

综上，本研究筛选的5个靶基因是潜在COAD的预后生物标志物，HSD17B2更是有望作为降低COAD化疗药物耐药性的靶点。然而，本研究是基于多个公共数据库的数据分析，具有局限性，这些关键基因在结肠腺癌预后以及化疗耐药中的分子机制仍需要进一步的基础实验来证实。