基于加权基因共表达网络挖掘卵巢癌相关基因△

李康梅，陈泯锜，戴明明，陈秀榕，黎明星，何国珍#

广西中医药大学1基础医学院，2赛恩斯新医药学院，3研究生院，南宁 530200

卵巢癌是妇女生殖系统致死率最高的恶性肿瘤，发病率和病死率均较高，并呈逐年上升趋势[1-2]。数据表明，中国每年有14 080人死于卵巢癌[3]。尽管肿瘤的治疗不断发展，但仍有2/3以上的卵巢癌患者未能早期发现和确诊，生存率仍无明显改善，且其发病机制和起源尚未清楚[4]。因此，研究卵巢癌的发生发展机制，寻找相关的枢纽基因对卵巢癌的早期诊治及预后判断具有一定意义。

加权基因共表达网络分析（weighted gene coexpression network analysis，WGCNA）是一种基于高通量基因表达谱的算法，广泛应用于各种疾病的基因共表达网络识别中，以揭示基因间的相关性，寻找显著相关的基因模块。借助WGCNA算法构建基因模块，可挖掘基因模块信息，对相关基因模块中的枢纽基因进行研究。本研究基于WGCNA对美国国立生物技术信息中心（National Center for Biotechnology Information，NCBI）的基因表达综合数据库（gene expression omnibus，GEO）下载的卵巢癌基因表达谱数据集GSE18520进行挖掘相关基因，构建卵巢癌基因共表达网络，识别卵巢癌相关基因模块及模块内枢纽（Hub）基因，为进一步探索卵巢癌的发生发展机制提供一种新思路。

1 材料与方法

1.1 数据下载和数据预处理

从GEO数据库（https://www.ncbi.nlm.nih.gov/geo/）下载基因表达谱数据集GSE18520，共53个卵巢癌样本和10个正常样本，以筛选卵巢癌的差异表达基因（differentially expressed gene，DEG）。对下载的数据进行预处理，包括提取样本信息、构建基因表达矩阵以及将探针名转化为基因名等。

1.2 对卵巢癌进行DEG筛选

使用“limma”R包筛选卵巢癌样本和正常样本之间的DEG。将调整后的P＜0.05和|log2 fold change（FC）|＞1作为筛选条件。

1.3 卵巢癌中DEG的PPI网络

使用 STRING（http://www.string-db.org/）数据库评估基因列表背后的相互关系。同时，对网络进行MCODE分析，以degree cut-off=5，node score cut-off=0.2，k-core=2，maximum depth=100为参考值选择PPI网络的核心模块。

1.4 卵巢癌中DEG的共表达分析

安装R软件的WGCNA包，为DEG构建一个无标度的共表达网络，选取合适的加权系数β，将相关矩阵转化为邻接矩阵。建树的方法采用动态混合剪切法，将相异度作为距离测度，设定最小模块尺寸为30，进行模块识别并绘制基因树状图。

1.5 核心模块及核心基因的识别

WGCNA中的每一个模块都在STRING数据库中进行分析，通过模块内统计分析，根据每个基因与该模块的相关系数即模块隶属度（module membership，MM）≥0.8且定义基因的显著性（gene significance，GS）P＜0.05筛选出核心基因。

1.6 功能注释和通路分析

为了研究大量基因的功能注释，使用“Clusterpro-filer”软件包对枢纽模块进行了基因本体（gene ontology，GO）功能富集分析和京都基因以及基因组百科全书（Kyoto encyclopedia of genes and genomes，KEGG）通路富集分析。

1.7 枢纽基因的验证和生存分析

利用基因表达谱交互分析数据库（gene expression profiles interactive analysis database，GEPIA，http://GEPIA.cancer-pku.cn/index.html）对核心基因进行进一步验证和生存分析。以P＜0.01在GEPIA数据库中绘制枢纽基因的生存曲线。

1.8 统计学方法

采用R语言平台Rstuio对卵巢癌临床特征和基因过表达模块相关系数进行计算。

2 结果

2.1 卵巢癌的DEG筛选

在R软件的“limma”软件包中，对53个卵巢癌样本和10个正常样本的GSE18520基因表达谱进行分析。共检测到2474个DEG（上调864个，下调1610个）。（图1）

图1 GSE18520基因表达谱中所有DEG的火山图

2.2 枢纽模块的鉴定

根据在MCODE中评估的分数得出的前3个簇作为核心模块。其中MCODE 1有40个节点和592条边，在这些集群中得分最高（图2）。

2.3 筛选软阈值

共表达网络符合无尺度网络，即出现连接度为k的节点的对数lgk与该节点出现的概率的对数lg[P（k）]呈负相关，且相关系数应＞0.8。使用R软件WGCNA包构建权重共表达网络，使用分析包自动选择的软阈值计算得到软阈值β=20（scale freeR^2=0.89）。

2.4 WGCNA 的构建和关键模块的识别

使用WGCNA R包，可以通过平均连锁聚类将具有高度相关表达模式的DEG放入模块中，并挖掘出8个模块（图3A）。设置阈值为0.25以合并相似模块（图3B），并相应生成6个模块。无法纳入任何模块的基因放入灰色模块，在后续分析中剔除。其他5个模块分别用蓝色、棕色、黑色、粉色和蓝绿色表示（图3A）。GS和MM分析在蓝色和蓝绿色模块中显示出非常显著的相关性（图3C）。

图2 在MOCDE中评估的PPI网络中得分最高的枢纽模块

2.5 枢纽模块的GO功能富集分析和KEGG通路富集分析

对蓝色和蓝绿色枢纽模块进行GO功能富集分析以及KEGG通路富集分析。GO功能富集分析分别有生物过程（biological process，BP）涉及的调节神经递质水平和细胞组成（cellular component，CC）涉及的细胞皮层、核外膜、细胞器外膜、有机外膜（表1）。KEGG通路富集分析主要有癌症蛋白聚糖、视黄醇的新陈代谢、黏着斑、WNT信号通路、肌动蛋白骨架的调节、代谢途径（表2）。

2.6 枢纽基因的验证和生存分析

图3 WGCNA的构建和关键模块的识别

表1 蓝色和蓝绿色枢纽模块GO功能富集分析

表2 蓝色和蓝绿色枢纽模块KEGG富集通路分析

通过GEPIA数据库对蓝色模块和蓝绿色模块的482个核心基因进行验证和生存分析。根据模块相关性大于0.9进行筛选，得到以下8个基因：AC104667.3、BTD、DDX26B、K+通道子族B成员1（K+channel subfamily B member 1，KCNB1）、前列腺素F2α受体（prostaglandin F2 receptor，PTGFR）、脑糖原磷酸化酶（brain glycogen phosphorylase，PYGB）、runt相关转录因子1（runt-related transcriptionfactor1，RUNX1T1）、TMEM91。其中AC104667.3、PYGB在卵巢癌中表达上调，BTD、DDX26B、KCNB1、PTGFR、RUNX1T1、TMEM91在卵巢癌中表达下调。AC104667.3、DDX26B的高表达与更好的总生存率有关，BTD、KCNB1、PTGFR、PYGB、RUNX1T1、TMEM91的高表达与预后不良有关（图4）。

3 讨论

本研究主要通过生物信息学分析GSE18520（53个卵巢癌样本和10个正常样本）的基因表达谱中卵巢癌发生及发展过程中的枢纽基因。研究结果显示，2474个基因在肿瘤组织中有差异表达。在这2474个基因中，有864个上调基因和1610个下调基因。利用Cytoscape软件的MCODE插件选择PPI网络中最重要的前3个模块作为枢纽模块，并对最重要的枢纽模块（MCODE1）进行GO功能富集分析以及KEGG通路富集分析的进一步研究。

图4 不同基因表达水平的卵巢癌患者的生存曲线

为了在卵巢癌中发现具有潜在价值的基因，在蓝色模块和蓝绿色模块的482个枢纽基因中根据相关系数大于0.9筛选出8个枢纽基因（AC104667.3、BTD、DDX26B、KCNB1、PTGFR、PYGB、RUNX1T1、TMEM91）。Tan 等[5]研究表明，BTD在宫颈癌中表达上调，作为宫颈癌的潜在基因，其意义还需进一步研究。DDX26B是一种RNA解旋酶，是不同的RNA进行代谢所必需的[6]。刘志荣等[7]研究表明DDX26B基因在早期大肠癌和腺瘤中不表达，而在正常大肠组织中表达。KCNB1通道的氧化是一种神经脆弱的机制，这种机制在脊椎动物的大脑中普遍存在[8]。如小鼠衰老过程中KCNB1通道的中度氧化可以通过影响突触功能影响神经元网络[9]。Yu等[10]通过使用药物直接撞击KCNB1通道氧化激活的毒性路径，改善小鼠颅脑损伤（traumatic brain injury，TBI）的结果，从而阐明了脊椎动物神经毒性的机制，为人类创伤性脑损害提供一种新的治疗办法。研究认为PTGFR和神经生长因子（nerve growth factor，NGF）蛋白在子宫腺肌病痛经的发生和发展中起作用，血清NGF蛋白水平及血浆PTGFR蛋白水平与子宫腺肌病痛经程度呈正相关，通过抑制其表达可能为子宫腺肌病痛经提供新的治疗线索[11-12]。研究表明，PTGFR激活可能诱发子宫内膜修复移植PTGS-2基因表达和刺激血管内皮生长因子（vascular endothelial growth factor，VEGF）、结缔组织生长因子（connective tissue growth factor，CTGF）、转化生长因子β1（transforming growth factor β1，TGF-β1），引发蛋白激酶C（protein kinase C，PKC）信号通路激活和基因表达[13]。PYGB目前已被报道参与多种肿瘤的发生，糖原代谢与许多其他关键代谢途径交叉，在肿瘤微环境中发挥着关键作用，抑制糖原代谢可能是一种有前途的抗癌治疗策略[14-15]。Zhou等[16]研究表明，PYGB在卵巢癌组织中表达上调，高水平的PYGB表达与卵巢癌患者预后不良显著相关，这与本研究结果一致。另外敲除PYGB后可显著抑制卵巢癌细胞的增殖、侵袭和迁移，并发现miRNA-133a-3p/PYGB轴在卵巢癌治疗中是一个重要的诊断指标和治疗靶点，可能成为未来治疗卵巢癌的靶标[16]。Sun等[4]研究发现，EPB41L4A-AS2通过与miRNA-103a联合诱导RUNX1T1表达可导致卵巢癌细胞的迁移和侵袭，证明了EPB41L4AAS2、miRNA-103a和RUNX1T1的相互作用在卵巢癌的发生和发展中起关键作用，但关于miRNA-103a和RUNX1T1的表达如何影响卵巢癌的发生和发展还需要进一步的研究。Yeh等[17]研究也表明了RUNX1T1在卵巢癌细胞系中表达下调，可以抑制卵巢癌细胞的生长，通过卵巢癌中异常的TGF-β/SMAD家族成员4（SMAD family member 4，SMAD4）信号建立表观遗传沉默，但RUNX1T1在卵巢癌中的作用机制尚不明确，还需进一步研究。另外，RUNX1T1蛋白也作为诊断肝转移的一个有前途的生物标志物[4]。而AC104667.3、TMEM91目前尚未有相关研究。

综上所述，本研究通过运用一系列合理的生物信息学手段，结合文献报道以及本研究结果，发现8个基因均与卵巢癌的发生和发展相关。而PYGB和RUNX1T1基因与卵巢癌的生存及预后密切相关，可能通过细胞周期相关通路影响卵巢癌的发生、进展及预后，为卵巢癌的进一步研究提供依据及参考。