基于WGCNA联合LASSO筛选宫颈癌预后枢纽基因

郭依琳 王 璐 徐 臻 赵 虎 韩少聪 王武亮

宫颈癌是严重威胁女性健康的恶性肿瘤,在全球范围内发生率和病死率均居女性恶性肿瘤第4位[1]。在我国宫颈癌的发生率和病死率分别位居第6位和第8位,尽管近年发生率有所下降,但发病年龄逐渐呈年轻化[2, 3]。HPV是宫颈癌发病的主要原因[4]。早期宫颈癌患者以手术治疗为主,预后较好。中晚期及复发的宫颈癌患者以放化疗治疗为主,预后较差[5]。目前可以用来预测宫颈癌患者预后的分子标志物较少。因此,寻找合适的预后生物学标志物和治疗靶点可能有助有提高宫颈癌患者的总生存率。

近年来,随着基因组学和生物信息学的不断发展,利用高通量测序和基因芯片技术筛选宫颈癌预后相关标志物具有重要研究价值。加权基因共表达网络分析(weighted gene co-expression network analysis,WGCNA)是一种识别与表型性状相关的基因模块和关键枢纽基因的方法[6]。最小绝对值选择与收缩算子(least absolute selection and shrinkage operator,LASSO)是一种通过构造惩戒函数压缩部分低权重回归系数,明确两变量之间的关联程度的回归算法[7]。WGCNA联合LASSO算法可以更好的筛选出相关性更高的枢纽基因。因此,本研究基于基因表达汇编(Gene Expression Omnibus,GEO)数据库、肿瘤基因组图谱(The Cancer Genome Atlas,TCGA)数据库和基因型和基因表达量关联(Genotype Tissue Expression,GTEx)数据库,运用WGCNA联合LASSO算法筛选出宫颈癌预后相关的分子标志物,并初步分析其与宫颈癌肿瘤微环境(tumor microenvironment,TME)中免疫细胞浸润的相关性,为宫颈癌发生、发展中的分子机制研究提供依据。

资料与方法

1.数据的获取和处理:从GEO数据库选取GSE9750和GSE52903数据集。GSE9750数据集包括33例宫颈癌组织和24例正常宫颈组织;GSE52903数据集包括55例宫颈癌组织和17例正常宫颈组织。从TCGA数据库选取304例子宫颈癌组织和3例正常宫颈组织。从GTEx数据库选取10例正常宫颈组织。下载芯片及RNA测序转录组数据和对应的临床特征数据。在R软件(版本号4.1.2)中利用“sva”包中的“combat()”函数对GSE9750和GSE52903数据集进行标准化处理,合并和校正,主成分分析(Principal component analysis,PCA)检验校正的数据。利用“limma”包中的“normalizeBetweenArrays()”函数对TCGA和GTEx数据集进行合并和校正。

2.差异表达基因的筛选和功能富集分析:在R软件中使用“limma”包筛选宫颈癌中差异表达的基因(differentially expressed gene,DEG),筛选标准P<0.05且|log2差异倍数(fold change,FC)|>2。绘制DEG的火山图。在R软件中使用“cluserProfiler”包对宫颈癌和正常宫颈中DEG进行基因本体论(Gene Ontology,GO)功能学。错误发现率(false discovery rate,FDR)值<0.05被认为差异有统计学意义。

3.WGCNA模块构建及可视化:在R软件中使用“WGCNA”包,删除GEO数据集(GSE9750和GSE52903)中基因表达离散程度较小(sd≤0.7)的基因,使用“goodsamplegenes()”函数检查缺失值。使用“pickSoftThreshold()”函数计算共表达网络的软阈值。当软阈值等于4,共表达网络更逼近无尺度网络。构建加权邻接矩阵,基于拓扑重叠矩阵(TOM)的相异度(1-TOM)进行层次聚类,构建相关模块。将模块与临床特征数据关联后,绘制模块身份(module membership,MM)及基因显著性(gene significance,GS)的散点图,明确模块内基因的显著性。

4.枢纽基因的筛选:WGCNA中按照MM >0.8且GS >0.5筛选模块内的基因。将筛选到的模块基因与差异表达基因取交集,并在R软件中使用“glmnet”包,在α=1的条件下,选择合适的λ值,用十折交叉验证法实现内部验证,确定最佳枢纽基因。使用TCGA可视化分析工具GEPIA,进一步筛选出与预后相关的枢纽基因。

5.枢纽基因在宫颈癌和正常宫颈组织的表达:在GEO数据集(GSE9750和GSE52903)和TCGA联合GTEx数据集中,检测预后相关的枢纽基因在宫颈癌组织与正常宫颈组织之间的表达情况。人类蛋白图谱数据库(human protein atlas,HPA)分析与预后相关的枢纽基因在宫颈癌组织与正常宫颈组织之间的蛋白表达情况。

6.免疫浸润和免疫检查点基因相关性分析:利用聚类分析和单样本基因集富集分析(single sample gene set enrichment analysis,ssGSEA)计算宫颈癌组织和正常宫颈组织中免疫细胞浸润程度,绘制小提琴图。计算预后相关的枢纽基因与28种免疫细胞以及36个免疫检查点基因的相关性,以P<0.05为差异有统计学意义。

结 果

1.纳入患者临床特征分析:本研究共纳入392例宫颈癌组织和54例正常宫颈组织的转录组表达数据。GSE9750、GSE52903和TCGA数据集具有完整的临床病理数据。3组数据集在年龄分布和病理类型之间比较,差异无统计学意义(P<0.05),但在临床分期和生存状态之间比较,差异有统计学意义(P<0.05,表1)。

表1 392例宫颈癌组织的临床病理特征[n(%)]

2.宫颈癌中差异表达的基因和功能富集分析:将GSE9750和GSE52903数据集合并且校正后,利用PCA主成分分析检验校正后的数据,未发现存在批次效应(图1)。在88例宫颈癌样本和41例正常宫颈样本中,根据设置的筛选条件,共发现109个DEG,其中表达上调38个,表达下调71个,并绘制火山图(图2)。根据P值选择表达差异最显著的前10个基因(表2)。对109个DEG进行GO功能注释分析。

图1 校正后的GSE9750和GSE52903数据集PCA分析

图2 宫颈癌中差异表达基因的火山图

表2 宫颈癌中差异表达最显著前10个基因

GO功能结果注释发现,差异表达基因可能参与109条通路的功能,其中生物过程(biological process,BP)富集度最大的3条通路为皮肤发育、表皮发育和表皮细胞分化;细胞组分(cellular component,CC)富集度最大的3条通路为角化包膜、中心颗粒体和DNA复制前起始复合物;分子功能(molecular function,MF)富集度最大通路为丝氨酸酶活性(图3)。

图3 宫颈癌中差异基因的GO功能富集分析

3.WGCNA网络构建及显著性模块确认:对GSE9750和GSE52903数据集中所有基因构建共表达网络。当软阈值等于4时,共表达网络接近为无尺度网络。确定软阈值后,采用动态剪切法,获得9个模块(图4A)。进一步分析各模块的特征向量,发现turquoise模块(包含126个基因)与宫颈癌的相关性最高(r=0.87,P=2×10-40),并且GS和MS最高,确认为显著性块(图4B)。

图4 WGCNA网络构建及显著性模块确认A.基因共表达网络和共表达模块的层次聚类树;B.特征模块与临床特征之间的关系。红色为正相关,蓝色为负相关

4.枢纽基因的筛选:将turquoise模块中126个基因与差异表达109个基因取交集,得到27个基因。利用LASSO模型进行内部验证,最终确定了7个枢纽基因(图5A)。GEPIA数据库进一步分析得到3个预后相关的枢纽基因,分别为MCM2、APOD和RAD54L(图5中B～D)。

图5 宫颈癌中筛选枢纽基因A.LASSO模型建立;B.MCM2表达水平的生存曲线图;C.APOD表达水平的生存曲线图;D.RAD54L表达水平的生存曲线图

图6 宫颈癌中预后相关枢纽基因的表达情况A.GSE9750和GSE52903数据集;B.TCGA和GTEx数据集。*P<0.001

5.预后相关枢纽基因在宫颈癌和正常宫颈组织表达情况:在GSE9750和GSE52903数据集中,与正常宫颈组织比较,MCM2和RAD54L在宫颈癌组织中高表达,而APOD则低表达(图6A)。在TCGA联合GTEx数据集中,3个枢纽基因的表达结果与GEO数据集(GSE9750和GSE52903)结果一致(图6B)。为了在蛋白水平验证3个预后相关的枢纽基因的表达情况,笔者在HPA数据库中分析了其在宫颈癌组织和正常宫颈组织中免疫组化的表达情况。结果发现MCM2和RAD54L在宫颈癌组织中的表达高于正常组织,而APOD在宫颈癌组织和正常宫颈组织都属于低表达(图7)。

6.预后相关的枢纽基因与免疫相关性分析:ssGSEA方法计算宫颈癌组织和正常宫颈组织中免疫细胞细胞浸润程度,初步探究宫颈癌TME中免疫细胞浸润程度。结果发现激活的B淋巴细胞、激活的CD4+T细胞、激活的树突状细胞、调节性T细胞和辅助性T细胞在宫颈癌组织中表达上调;而肥大细胞、嗜酸性粒细胞和中性粒细胞中则在宫颈癌组织中表达下调(图8A)。进一步探究3个预后相关的枢纽基因与28种免疫细胞以及36个免疫检查点基因的相关性。结果发现,在28种免疫细胞中,MCM2基因与CD56bright.NK细胞呈明显正相关,与肥大细胞呈明显负相关;APOD基因与浆细胞样树突细胞和激活的B细胞呈明显正相关;RAD54L基因与大部分免疫细胞呈负相关(图8B)。在36个免疫检查点基因中,MCM2基因与免疫检查点LGALS9呈明显正相关,与NRP1呈明显负相关;APOD基因与CD70和TNFSF9呈明显正相关;RAD54L基因与CD44呈明显正相关(图8C)。

讨 论

宫颈癌传统的治疗方式主要为手术切除、放疗和化疗,不同患者对放化疗的敏感度差异很大,并且即使患者接受上述治疗后仍会出现复发和转移的可能[8, 9]。宫颈癌的复发和转移是目前造成宫颈癌患者死亡的主要原因,临床上迫切需要发现更多分子标志物用于宫颈癌的预后预测。因此,本研究基于GEO数据库中GSE9750和GSE52903数据集,利用WGCNA和LASSO筛选宫颈癌枢纽基因,GEPIA数据库进一步筛选预后相关的枢纽基因;然后在GEO数据集(GSE9750和GSE52903)和TCGA联合GTEx数据集中比较预后相关的枢纽基因在宫颈癌和正常宫颈组织中的表达情况,并在HPA数据库中验证;最后利用ssGSEA分析宫颈癌TME免疫细胞浸润情况,探究3个预后相关的枢纽基因与免疫细胞浸润和免疫检查点基因表达的相关性。

近年来,随着计算机技术及人工智能等在生物医学方面的快速发展,生物信息学已经成为基于大数据进行分子标志物筛选的有力手段之一。方萌等[10]研究利用WGCNA和机器深度算法确定了5个miR-141-3p的靶基因FOXA1、DMBX1、TMEM98、RHPN1、SRMS,并构建了预后模型。Liu等[11]研究利用WGCNA方法,在GSE26511数据集中,发现ACKR1基因可能与宫颈癌淋巴结转移和预后有关。同样,Wang等[12]也利用相同的方法发现RRM2对宫颈癌的诊断和预后具有一定的价值。然而,运用单一的生物信息学分析方法,可能造成干扰数据过多,对结果的准确性造成影响。因此本研究选了WGCNA和LASSO模型两种方法联合分析增加研究结果的准确性和真实性。

本研究在GSE9750和GSE52903数据集中共获得109个差异表达的基因和126个WGCNA筛选turquoise模块基因,取交集后得到27个枢纽基因,利用LASSO模型确定7个枢纽基因。利用GEPIA数据库进一步分析最终得到3个预后相关的枢纽基因,分别为MCM2、APOD和RAD54L。MCM2基因编码微小染色体维持蛋白2,是DNA复制起始的主要调控因子之一。有研究发现,MCM2的过表达与多种恶性肿瘤的发生、发展相关[13]。Das等[14]研究发现,MCM2在宫颈癌组织和宫颈癌细胞系中高表达,然而MCM2的表达水平与临床病理分期无关,这在Amaro等[15]研究中也被证实。APOD基因编码载脂蛋白D,有研究发现,APOD参与多种免疫反应、细胞凋亡和肿瘤的发生、发展[16]。然而,APOD在宫颈癌方面的相关研究较少,仅Wang等[17]通过生物信息学技术确定了包括APOD基因在内的10个免疫相关性基因。RAD54L基因在同源重组修复中发挥重要作用。有研究发现RAD54L在脉络丛癌的发生、发展中起到重要作用,但其在宫颈癌中的研究较少。因此,本研究在GEO数据集(GSE9750和GSE52903)和TCGA联合GTEx数据集中比较了3个预后相关的枢纽基因在宫颈癌组织和正常宫颈组织中的表达情况,发现与正常宫颈组织比较,MCM2和RAD54L在宫颈癌组织中高表达,而APOD则低表达,为后续3个枢纽基因在宫颈癌中的作用机制提供前期研究的依据。

有研究发现,肿瘤的进展不仅受到肿瘤本身生物学行为特征的影响,同时也有受到TME的影响。TME是肿瘤生长、转移和调节肿瘤免疫反应的关键因素[18]。本研究采用ssGSEA方法,计算了宫颈癌组织和正常宫颈组织中28种免疫细胞浸润程度,发现大部分免疫细胞,如激活的B淋巴细胞、激活的CD4T细胞、激活的树突状细胞等在宫颈癌组中表达高于正常宫颈组织。另外,笔者还比较了3个预后相关的枢纽基因与免疫细胞和免疫检查点基因表达的相关性,发现APOD基因与大多数免疫细胞的表达呈明显负相关。特别地,3个预后相关的枢纽基因与部分免疫检查点基因的表达呈正相关。这些结果提示筛选的3个预后分子标志物与宫颈癌TME的免疫浸润水平存在关联,可能对宫颈癌的免疫治疗可以提供一定的参考依据。

综上所述,本研究利用WGCNA联合LASSO方法筛选出了3个与预后相关的枢纽基因,MCM2、APOD和RAD54L,探究了其在宫颈癌和正常宫颈组织中的表达,并初步分析了其与宫颈癌TME的免疫浸润水平的关系,为宫颈癌的预后预测和免疫治疗提供指导意义。