通过加权基因共表达网络分析法探索胰腺癌特异表达的关键基因及表达网络

陈立材 成 雨

1 滨州医学院第二临床医学院 山东 烟台 264003；2 滨州医学院烟台附属医院 山东 烟台 264100

胰腺癌是一种进展迅速，恶性程度极高的消化道肿瘤，早期诊断困难，5年生存率极低，其中90%以上为胰腺导管腺癌(pancreatic ductal adrenocarcinoma，PDAC)[1-2]。随着近年来PDAC的发病率的上升，研究其发生发展的机制，寻找潜在治疗靶点，开发新型的治疗手段显得尤为重要；而针对全基因组的芯片及测序技术的成熟为探索肿瘤标志物提供了技术支持与数据支撑。本研究基于肿瘤基因组图谱(the cancer genome atlas，TCGA)数据库，采用加权基因共表达网络分析(weighted gene co-expression network analysis，WGCNA)法[3]，探索胰腺癌特异表达的关键基因及表达网络，为进一步发现胰腺癌的生物标志物，确立新的诊断及治疗靶点提供思路。

1 材料与方法

1.1 数据的获取与差异表达基因的确定 在TCGA数据库(https://portal.gdc.cancer.gov/)下载并处理胰腺癌数据库的mRNA表达RNA seq数据，同时下载临床资料数据(Clinical)。本研究从TCGA数据库中下载了182例转录组数据，其中178例胰腺癌患者，4例健康对照组。对TCGA数据库的基因矩阵信息进行了预处理，并转化为基因名。在R语言环境下运行limma包，将logFC>1，矫正后的P<0.05确定为差异表达的基因(differently expressed genes，DEGs)，进行后续的共表达网络的构建。

1.2 WGCNA构建基因模块流程 基于DEGs，用R语言的WGCNA法构建胰腺癌权重共表达网络[3]。首先根据R2>0.9，根据真实生物网络状态的无尺度网络确定加权系数β(软阈值)，在确定邻接函数参数β后，构建不同分支和颜色表达的不同基因模块的分层聚类树最后，根据基因间Pearson相关系数将相关矩阵转换为邻接矩阵，进一步转化为拓扑重叠矩阵(topological overlap matrix，TOM)。在下一步分析中，基因表达模块被归类为不同的模块(Module)

1.3 模块与临床特征关联分析 将模块相关的网络矩阵(module eigengene，ME)与临床性状的Pearson相关系数进行计算。P<0.05确定为显著差异。定义显著性P值的以10为底的对数为基因显著性(gene significance，GS)，再将每一个模块显著(module significance，MS)定义为模块中所包含基因的GS的平均值。通过分析MS与GS，取相关系数最高的模块用于后续分析。

1.4 蛋白质相互作用(protein-protein interaction，PPI)网络分析及核心基因的确定 通过在线分析网站 STRING[4](http://www.string-db.org/)得到DEGs 的蛋白质相互作用网络，以 TSV格式导出，所得源文件导入到Cytoscape[5](http://www.cytoscape.org/, version 3.2.0)软件中，使用Cytoscape将基因信息进行可视构建，并生成相关的网络结构图，之后用插件cytoHubba 进行核心基因分析, 同时采用 MCC 算法, 选取排名前10位的基因为核心基因。

2 结果

2.1 胰腺癌组织和正常组织的DEGs 通过对数据的标准化及预处理，共有14 869个基因，通过设定的阈值(logFC>1，P<0.05，FC: Fold Change)筛选之后，得到106个差异表达的基因，其中表达上调的基因70个，表达下调的基因36个，见图1。

A.热图；B.火山图。

2.2 WGCNA分析

2.2.1 网络构建及模块识别 为了尽量满足无尺度网络分布前提条件，需要探索邻接矩阵权重参数β取值。通过设置网络构建参数选择范围，计算无尺度分布拓扑矩阵。计算相应的模型选择统计量绘制图形，见图2，图中的横轴代表权重参数power，纵轴代表对应的网络中log(k)与log(p(k))相关系数的平方。相关系数的平方取值越高说明该网络越逼近无网络尺度的分布。最上面的树形图表示基于TOM的系统聚类，Dynamic tree cut表示网络模块前后合并模块。本研究选取相关系数平方值首次达到0.9时的power的取值。根据TOM进行层次聚类得到基因的系统聚类树。

A.邻接矩阵权重参数power选择图，红色线表示相关系数的平方值达到0.9的标准线；B.基于TOM的基因系统聚类树的识别结果，图中不同的颜色代表不同多基因模块。

2.2.2 关键模块的确定 对各个模块和样本临床信息进行关联分析，从模块和性状热图中可以发现紫色(MEpurple)模块与胰腺癌相关程度最高，见图3A。计算紫色模块中基因GS和MM相关系数(cor=0.58)进一步验证此结果的可信度，见图3B。

A.热图；B.关键模块。

2.3 PPI蛋白网络分析及核心基因确定 通过在线分析网站 STRING，对关键模块中DEGs进行分析，得到PPI蛋白网络相互作用图，见图4，进一步使用Cytoscape将基因信息进行可视化及网络构建，并用插件cytoHubba 进行核心基因分析, 确定排名前10位的基因为核心基因，分别为PKP3，EPCAM，RAB25，CBLC，AP1M2，PRP15L，B3GNT3，ESRP1，AGR2，ARHGEF16，见图5。

图4 PPI蛋白网络构建

图5 Top10核心基因确定

3 讨论

随着发病率的不断增加，胰腺癌逐渐成为世界范围内最致命的恶性肿瘤之一[6-7]。虽然胰腺癌在治疗方面取得了一定的进展，但是由于其缺乏典型的临床表现及特异性的肿瘤标志物，并常伴有血管神经浸润及早期远处转移，其预后往往不佳[8]。目前探索胰腺癌发生发展机制，构建关键基因表达网络从而发现其早期生物标志物已成为研究热点。近年来，基因芯片及测序技术的进步，为肿瘤疾病的深入研究提供了可能性。基于这些技术，癌症基因组研究项目将人类全部癌症的基因组变异图谱进行绘制，收录于TCGA数据库。目前，TCGA数据库已收录30多种癌症数据及临床信息，总计超过一万例患者的基因组序列，供科研人员免费下载使用[9]。本研究基于TCGA数据库，共下载了182例患者的转录组基因信息，其中包括178例胰腺癌患者，4例对照组。经过数据的预处理，共发现106个差异表达的基因，其中表达上调的基因70个，下调的基因36个。

WGCNA算法是一种构建基因共表达网络的经典算法。WGCNA基于高通量mRNA表达芯片数据，假定基因网络服从无尺度网络，通过定义共表达矩阵和邻接函数，并将其转换为拓扑矩阵，从而识别与疾病关联的基因集合模块，从生物功能整体考虑基因功能及其联系，弥补了传统方法的缺陷。通过将临床信息与模块相关联，还可进一步获得与临床特征相关的基因，有助于基于疾病模型的临床特征构建相关基因的表达网络[3, 10]。目前有很多研究运用了WGCNA算法对肿瘤疾病的基因表达谱进行分析，并取得了有意义的进展[10-13]。本研究中，运用WGCNA算法，我们发现可以发现紫色(MEpurple)模块与胰腺癌相关程度最高(cor=0.58)。通过PPI网络分析，并将结果导入Cytoscape软件中进行基因可视化，确定排名前10位的核心基因为确定排名靠前10的基因为核心基因，分别为PKP3，EPCAM，RAB25，CBLC，AP1M2，PRP15L，B3GNT3，ESRP1，AGR2，ARHGEF16。这10个基因可能是胰腺癌的发生发展的关键基因。

PKP3是桥粒斑菲素蛋白家族中的一员[14]，一般位于所有含有桥粒的复层上皮以及单层上皮组织，可以与FXR1、PABPC1等RNA结合蛋白紧密结合，并且在细胞受到氧化应激等外源性刺激时出现应激颗粒，提示PKP3与RNA代谢、基因转录后调节密切相关，并且参与肿瘤的调控[15-16]。EPCAM是一种跨膜糖蛋白, 参与多种细胞活动包括增殖、迁移、分化等。另外，EPCAM也可介导细胞粘附，参与细胞内信号转到等。其表达与肿瘤的恶性程度相关[17-18]。在本研究中，PKP3和EPCAM是与胰腺癌相关度最高的基因，可能是影响胰腺癌生存和预后的关键基因，也可能作为胰腺癌诊断和治疗的潜在靶点，需要动物及临床试验进行进一步的验证。