胶质母细胞瘤中差异甲基化增强子区域调控的蛋白编码基因识别研究

赵潇潇 于秋红 嵇江淮 王世佳王仁东李冬果*

(1.首都医科大学生物医学工程学院，北京 100069；2.首都医科大学临床生物力学基础研究北京市重点实验室，北京 100069；3.首都医科大学附属北京天坛医院高压氧科，北京 100070；4.浙江肿瘤医院放射物理科，杭州 310022；5.浙江省放射肿瘤学重点实验室，杭州 310022)

多形性胶质母细胞瘤(glioblastoma,GBM)是最常见且最致命的肿瘤之一，被归类为Ⅳ级胶质瘤[1]。它是一种高度侵袭性肿瘤其特征是脑血管的改变和周围组织沿血管空间的逐渐侵袭[2]。GBM细胞通常会侵入距肿瘤块几厘米的地方，甚至可以跨入对侧半球[3]。目前GBM的治疗标准[4]包括手术切除，然后辅以放射治疗(以下简称放疗)和辅助治疗。但由于这些肿瘤的弥散性浸润性质，很少能实现完全切除，GBM患者在治疗结束后容易复发，中位生存时间约15个月，且预后较差[5-6]。因此，迫切需要在GBM准确的分子机制和可靠的治疗靶点的研究中取得新进展。

癌症是异常遗传和表观遗传事件的结果。表观遗传机制是无法用DNA序列的变化来解释的稳定遗传特征，与组蛋白修饰一样，DNA甲基化并不影响基因组DNA序列本身[7]，而是在CG二核苷酸胞嘧啶上增加一个甲基(CH3)基团。在哺乳动物中，5 mC参与了长期的沉默过程，如X染色体失活、基因组印记、生殖系基因的体细胞沉默和重复DNA元素的沉默[8]。根据DNA甲基化的异常变化对肿瘤的诊疗和预测生物标志物被认为具有广泛的前景。

近年来，多是对GBM中启动子异常甲基化动力学进行研究[9]，也有少量对增强子甲基化方面的研究[10]。异常甲基化增强子会导致多种疾病，其中包括多种癌症的基因表达异常[11]。研究[12]显示，与启动子相比，增强子区域的甲基化状态与靶基因表达的相关性更强。然而，只有少量研究[13-14]关注增强子区域且增强子甲基化模式在GBM中仍不清楚。

1 资料与方法

1.1 数据来源及预处理

本研究从基因表达综合数据库(gene expression omnibus,GEO)数据库中下载GBM肿瘤样本(GSE36278)和正常样本(GSE42861)的甲基化数据，包括136个肿瘤样本和58个正常样本。肿瘤样本和正常样本的DNA甲基化数据是通过Infinium HM450k平台用同样方法获得的。从基因组数据共享数据库(genomic data commons,GDC)数据门户下载患者临床数据和136例GBM样本的表达数据。人类全基因组的注释数据(V19)收集于GENCODE数据库[15]。

对于甲基化数据，每个探针的甲基化水平用β值表示。β值=Imeth/(Imeth+Iunmeth)，其中Imeth为甲基化强度，Iunmeth为未甲基化强度。为确保甲基化水平的准确性，移除在所有样本中的甲基化水平超过30%为缺失值的探针。然后，本研究使用 “DMwR” R包中的knnImputation函数[16]对甲基化数据缺失值进行处理。将下载的RNA-seq(FPKM)数据使用log2转换进行标准化，表达值为0的用数据集中的最小正数代替，最终获得蛋白编码基因(protein-coding genes，PCGs)的表达谱数据。

1.2 构建增强子区域和启动子区域

本研究整合已文献[17]和GPL13534注释文件，从而获得了位于增强子的探针。本研究将这两部分探针合并，并删除重复的探针。结果保留161 708个增强子探针用于后续分析。已有研究[18-19]显示，典型增强子区域的中位间隔大小约为1 000 bp。因此，本研究以一个增强子探针坐标为基准，分别向其上游和下游扩展500 bp的区间来构建增强子区域。重叠的增强子区域被连接，并扩展成较大的区间，然后计算增强子区域中CpG探针的平均值，作为增强子区域的DNA甲基化水平[20]。

由于基因不仅受增强子甲基化的调控，也受启动子甲基化的调控。为了获得只受差异甲基化增强子区域(differential methylation enhancer regions, DMERs)调控的基因，本研究需要重新注释启动子区域，鉴别出可能受启动子甲基化调控的基因，然后移除这些基因。本研究定义将基因转录起始位点(transcription start site, TSS)上游2 kb区域作为启动子区，基因的DNA甲基化水平为映射到其启动子区域的探针平均值[21-22]。

1.3 识别受差异甲基化增强子区域调控的靶基因

本研究利用基于线性模型设计的R包“limma”[23]来识别肿瘤与正常样本之间的DMERs和差异启动子甲基化基因(differential promoter methylation genes, DPMGs)。采用Benjamini-Hochberg(BH)方法[24]对P进行校正。同样，本研究也鉴定了DMERs以进行进一步分析。本研究选择log2(FC)≥0.2和adjustedP≤0.05 的增强子区域作为高甲基化增强子区域，将log2(FC)≤-0.2 和adjustedP≤0.05 的增强子区域作为低甲基化增强子区域。为了获得仅受DMERs调控的基因，本研究需要重新注释启动子区域，移除可能受启动子甲基化调控的基因。研究[25]显示，增强子和基因之间已知的最大距离约为1 Mbp。因此，本研究选择了位于同一染色体上的DMER-PCGs对，基因的TSS与增强子区域之间的最大线性距离为1 Mbp。本研究用皮尔森相关系数(Pearson correlation coefficient, PCC)计算DMER与基因表达之间的相关性。基因表达被增强子甲基化负调控，高甲基化增强子区域可下调甚至沉默基因表达，而低甲基化增强子区域则倾向于激活基因表达[26]，因此本研究只保留了PCC<-0.2的DMER-PCGs对，且P≤0.01。

1.4 功能富集分析

为了预测DMERs调控的PCGs的功能，采用富集分析方法，对受低甲基化增强子区域调控的PCGs进行功能和通路的显著性分析，使用R包“clusterProfiler”[27]来预测PCGs的功能。通过BH方法校正P，如果校正后的P≤0.05，该生物过程(biological process, BP)和京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Genomes, KEGG)就认为是显著的。通过R包‘ggplot2 ’进行显示。

1.5 生存分析

为了鉴别挑选出的PCGs是否具有良好的预后效果，本研究基于152个GBM患者的表达谱信息，依据中值将患者分为两组。所有的研究均使用R 4.0.2完成。Kaplan-Meier生存分析和log-rank检测被用来评估两组患者的生存差异，以P<0.05为差异有统计学意义。

2 结果

2.1 GBM中增强子区域甲基化谱的构建

为了研究在GBM中增强子的DNA甲基化模式，本研究应用计算策略将增强子探针注释到Infinium 450K阵列，从而构建GBM增强子区域甲基化谱。在本研究中，增强子探针的收集来自先前发表的文献的[17]补充文件和GPL13534注释文件。利用上述方法构造了113 178个非重叠增强子区域。其中增强子区域长度为1 000 bp占比79.99%,长度为>1 000 bp和<2 000 bp占比17.79%，其他占比2.22%(图1A)。同样，本研究将Infinium 450K阵列重新注释到基因的启动子区域。最终54 477个探针位于20 386个基因启动子区域。虽然每个基因都有几个探针定位到相应的启动子区域，但本研究只计算该基因启动子中DNA甲基化探针的平均值作为该基因的DNA甲基化水平。

2.2 识别受差异增强子区域调控的基因

在对增强子区域甲基化谱进行预处理后，本研究从136个GBM样本和58个正常组织样本中识别出16 287个DMERs，其中包括2 271个高甲基化增强子区域和14 016个低甲基化的增强子区域。类似地，本研究基于基因启动子甲基化谱共鉴别出4 889个DPMGs。为了获得表达严格受差异启动子甲基化调控的基因，本研究进一步计算每个DPMG甲基化值与其对应基因表达值之间的PCC。结果显示，共识别出167个可能受启动子甲基化调控的基因，其中有144个低甲基化基因和23个高甲基化基因，这些基因将会被排除在研究之外。为了确定哪些基因受DMERs调控，本研究结合多组学数据建立增强子与其靶基因(DMERs-相关基因)的关联模型(如上方法描述)。

由于已知增强子甲基化水平与染色质活性呈负相关[28]，本研究只保留了PCC<-0.2的DMER-基因对。最后，本研究得到了795对DMER-PCGs，其中包含有593个低甲基化增强子区域，82个高甲基化增强子区域和642个PCGs。整体来看，GBM中存在明显的低甲基化模式(图1B)。

2.3 靶基因的功能富集分析

为了评估DMERs调控的靶基因的生物学特性，本研究对642个PCGs进行功能富集分析。结果显示，PCGs富集在许多与肿瘤进程相关的生物过程，比如细胞凋亡、突触信号传送与调节、组蛋白H3-K79甲基化的调控(图1C)。对于KEGG通路富集分析，本研究发现 PCGs在 KEGG中富集出多条与 GBM 相关的通路，比如ErbB(EGFR)、神经营养因子、黏着斑黏附(图1D)。

图1 GBM中增强子区域的DNA甲基化模式

2.4 GBM中预后良好的PCGs的鉴别

为了评估这些表达受DMERs影响的PCGs是否可以作为GBM的预后因子，本研究将表达数据与患者临床信息相结合，基于生存分析的结果，获得了45个与GBM整体存活相关的PCGs(P<0.05，表1)。这些基因的异常表达可能是由于增强子异常甲基化调控的。它们可能是GBM潜在的诊断和治疗的生物标志物。在本研究中，DTX1和SLC2A3 在GBM中均属于表达上调的基因。DTX1和SLC2A3 基因高表达组患者的生存期显著缩短且预后差，可作为GBM患者的独立预后因素(图2)。

图2 DTX1和SLC2A3基因生存分析

表1 预后良好的PCGs

3 讨论

研究[12, 29]表明增强子通过甲基化介导表观遗传调控，对靶基因的表达进行异常调节，从而导致人类疾病。表观遗传调控已成为DNA甲基化等生物医学研究的热点[30]。癌症中DNA甲基化的变化被认为是开发有效的诊断、预后和预测标志物的靶点。最近的研究[31-33]也报道了一些基因启动子在GBM中的异常甲基化。然而，到目前为止，使用异常甲基化增强子识别GBM生物标志物的研究很少。

本研究通过构建增强子区域，对其进行差异甲基化分析并鉴别出可能受DMERs调控的PCGs，结果发现DMERs在GBM中多呈现一种低甲基化的模式。为了评估PCGs的生物学功能，对这些基因进行功能富集分析。本研究结果显示PCGs参与了许多与肿瘤发生及进程相关的生物过程。由此推断这些受DMERs调控的PCGs更有可能成为GBM潜在的诊断和治疗的生物标志物。其中EGFR 在GBM 和胶质瘤中所发挥的作用已被广泛研究，其可以通过促进血管生成和细胞侵袭从而促进肿瘤的恶化[34]。Neurotrophin 为 GBM 发展提供了微环境，并且促进肿瘤浸润细胞的增生[35]。Huber等[36]发现DTX1的过表达增加了细胞的迁移和侵袭，这与ERK激活、miR-21水平和内源性Notch水平相关。DTX1的Notch通路似乎是胶质母细胞瘤的一种致癌因子，这些发现提供了新的潜在的治疗靶点。研究[37-38]表明SLC2A3 在GBM中高表达,沉默期的表达可以抑制GBM细胞的增生和侵袭能力。本研究结果表明DTX1和SLC2A3 在GBM中表达的上调可能是受增强子区域的低甲基化调控，从而导致疾病的发生与发展。

然而，这项研究有一些需要强调的局限性。GBM的样本采集是一个问题，这一缺陷可能会在一定程度上影响最终的结果。此外，本研究结果基于生物信息学分析。由于技术和时间限制，本研究没有验证在GBM动物模型和GBM患者脑组织中的结果。未来的功能研究和分子实验仍需进一步探索其作用机制。

综上，本研究通过构建GBM中较为全面的增强子区域的甲基化谱，鉴定出GBM中表达可能受DMERs调控的PCGs，挖掘出与GBM预后相关的PCGs并加深对GBM中差异甲基化增强子调控模式的理解。同时，本研究为在表观遗传方面识别潜在的GBM生物标志物和治疗靶点提供了新的思路。