加权基因共表达网络挖掘并验证调控前列腺癌转移的枢纽基因

张河元，陈南辉2，，王晓红，高白云，凌木安，陈 果，吴志明，李宇同，钟伟枫5，，潘 斌

1广东省梅州市人民医院泌尿外科，广东 梅州 514021；2南方医科大学南方医院泌尿外科，广东 广州 510515；3南方医科大学第三附属医院肾内科，广东 广州 510630；4暨南大学附属第一医院泌尿外科，广东 广州 510630；5广州市第十二人民医院泌尿外科，广东 广州 510630；6中山大学肿瘤防治中心泌尿外科，广东 广州 510060

前列腺癌（PCa）是泌尿系统最常见的肿瘤之一，世界范围内，前列腺癌发病率在所有男性恶性肿瘤中位居第二；因环境恶化及人口老龄化等原因，我国前列腺癌的发病率也在逐年上升［1］。虽然手术、根治性放疗等治疗手段对于早期前列腺癌有很好的疗效，而一旦前列腺癌患者发生肿瘤转移，其预后差且会严重影响其生活质量［2］。目前为止，前列腺癌转移的分子机制仍未明确，且临床治疗缺乏有效的治疗手段。因此，探索新型、有效的治疗方法成为现今国内外学者研究的重点。阐明前列腺癌发生和转移的分子机制，寻找可治疗前列腺癌的关键靶标至关重要。

传统的“单疾病单基因”的研究模式不能从多基因系统角度阐述疾病的发生发展，近年来，随着基因芯片、RNA-seq等高通量测序技术的兴起，产生了大规模的组学数据，而利用系统生物学分析的新方法，加权基因共表达网络分析（WGCNA）可大规模分析组学数据，并发掘与疾病相关的枢纽基因［3］，取得了一系列研究成果［4,5］。WGCNA可分析很多疾病的关键基因，比如影响成年人胰岛素敏感性的关键基因［6］，动脉粥样硬化合并糖尿病的功能基因模块［7］、影响脓毒症预后的关键基因［8］、阿尔茨海默病发病有关的关键基因［9］等，以及各种癌症的关键致病基因的识别，比如乳腺癌循环肿瘤细胞转移的关键基因分析［10］、舌鳞状细胞癌密切相关的枢纽基因分析［11］、龙葵抗肺腺癌的潜在生物靶标的探讨［12］及肾透明细胞癌进展相关基因的筛选［13］。但目前缺少关于前列腺癌的关键基因及预后系统性分析，因此，本研究采用WGCNA方法，对大样本的前列腺癌研究队列进行分析，筛选与前列腺癌发生和转移的枢纽基因，并通过TCGA公共数据库和细胞模型进行验证，以了解前列腺癌发生和转移进程中的关键基因，为临床治疗前列腺癌提供新的思路。

1 材料和方法

1.1 差异基因筛选

从NCBI的GEO公共数据库（http://www.ncbi.nlm.nih.gov/geo/）下载编号为GSE6919的前列腺癌基因芯片数据，该数据集采用GPL92、GPL93 及GPL8300 平台。共171个组织样本，包含了前列腺癌旁组织81例、前列腺癌组织65例及前列腺转移组织25例。利用R软件中的Affy包［14］对原始数据进行读取，通过PCA进行样本主成分分析，通过R软件的limma包［15］对基因表达矩阵进行差异基因分析，设置筛选值为：校正后P 值（FDR）＜0.05和对数化表达变化倍数（|log2FC|）＞1.2，作为筛选差异基因的阈值，总共进行两次分析：癌原位组ⅤS 癌旁组，癌转移组ⅤS癌原位组，得到相应组间的差异表达基因（DEGs），并绘制火山图及热图。

1.2 加权基因共表达网络分析

使用R软件的WGCNA包［3］，通过基因注释选择方差最大的前7500个基因进行WGCNA。计算各基因间的Pearson相关系数，构建无尺度网络并选择适当的阈值β进行网络构建。采用两步法构建，将邻接矩阵转化为拓扑重叠矩阵TOM，通过层次聚类生成聚类树，通过dynamic cut进行聚类合并。通过计算基因显著性（GS）以及模块显著性（MS），计算基因及模块，基因与临床样本分组信息的显著性。

1.3 枢纽基因的筛选

计算各基因模块身份（MM）以衡量基因在各个模块内的重要性。设置参数|MM|＞0.8及|GS|＞0.2为标准，筛选与临床性状密切相关的模块的枢纽基因。通过与相应组间的DEGs进行交叉，取两者的交集即为枢纽基因。使用DAⅤID数据库（http://david.abcc.ncifcrf.gov/）对特定模块的基因进行GO 分析。若错误发现率（FDR）＜0.05，认为差异具有统计学意义。

1.4 TCGA数据库验证与GO富集分析

通过R软件下载TCGA数据库前列腺癌研究对列，共得到498例前列腺癌样本，对照样本为52例。分析枢纽基因在对照组及前列腺癌样本组的表达，R软件采用COX分析前列腺癌患者的预后生存曲线。

1.5 细胞转染与Transwell 细胞迁移和侵袭

人前列腺癌细胞系PC3及LNCap由暨南大学附属第一医院中心实验室保存。细胞与RPMI 1640培养液中培养，待细胞融合至50%，采用Lipofectamine 2000试剂盒进行转染，HNRNPA2B1 siRNA及对照片段由GenePharma 吉玛公司提供，序列为：5'-GCAACCUU CUAACUACGGUtt-3' 和 5'-ACCGUAGUUAGAAG GUUGCtt-3'。方法步骤参照之前文献进行［16］。细胞转染HNRNPA2B1 siRNA片段后，收集各组处理细胞，用100 μL无血清重悬转入Transwell小室，进行迁移和侵袭（铺Matrigel胶），在37 ℃，5%CO2的细胞培养箱孵育24 h候取出小室，4%多聚甲醛固定20 min，PBS洗涤1次，结晶紫染色10 min，PBS洗涤干净，并于显微镜下观察细胞是否穿透小孔，并进行拍照，200倍光镜选取5个视野计数穿膜细胞数，实验重复3次。

1.6 MTT 实验

PC3及LNCap细胞分别接种于96孔板并进行小分子干扰转染，48 h后收集细胞，每孔加入20 μL MTT溶液（5 mg/mL），继续培养4 h，每孔加入100 μL二甲基亚砜，置摇床上低速振荡10 min。在酶标仪测量各孔的吸光度值A490nm。

1.7 流式细胞术与细胞凋亡

PC3和LNCap细胞转染后48 h后，用0.25%胰酶-EDTA消化，并用预冷的PBS洗涤3次，2000 r/min离心5 min，收集细胞。在检测细胞凋亡的试管中先后加入1×Buffer 100 mL、Annexin Ⅴ-FITC 染液5 μL、PI 染液1 μL，避光15 min，再加入1×Buffer 400 μL，上流式细胞仪进行检测。

1.8 克隆形成实验

细胞转染后，然后将细胞悬液转移到1.5 mL EP管中混匀并稀释，按照20 000/孔的细胞浓度接种于6孔板，1次/3 d观察及换液，克隆10 d左右，观察克隆团形成。吸出所有的培养基并清洗，用4%的多聚甲醛固定15 min，然后用0.1%结晶紫染色15 min后清洗，干燥后进行拍照，最后对克隆细胞团计数并进行统计学分析。

1.9 统计学处理

采用SPSS 20.0 统计软件进行数据分析，计量资料均由均数±标准差表示，两组间差异分析使用Student'st-test进行统计，P＜0.05 表示差异有统计学意义。

2 结果

2.1 PCA主成分分析

该队列研究共收集了171个组织样本，包含了前列腺癌旁组织81例、前列腺癌组织65例及前列腺转移组织25例。使用R软件读取及预处理芯片数据后，共得到可注释的8585个基因。对3组样本进行了PCA主成分分析，原位癌组和癌旁组聚类分不开，而癌转移组的基因成分和前面两组有显著区别，可单独聚类（图1）。

图1 基于3个组别每个样本的基因表达情况，进行PCA主成分分析Fig.1 PCA principal component analysis based on gene expression profile of each sample in 3 groups.

2.2 差异基因筛选

通过R语言的limma包进行了原位癌组ⅤS癌旁组、癌转移组ⅤS 原位癌组的比较。以FDR＜0.05 和|log2FC|＞1.2 为阈值，进行差异基因的筛选。结果表明，原位癌与癌旁组的基因表达差异不大，只发现了9个上下调差异基因（图2A），转移癌与原位癌基因差别明显，共筛选到425个差异基因（图2B）。

图2 DEGs火山图Fig.2 Volcano map for the differentially expressed genes (DEGs).A:Volcano map for DEGs between in situ and adjacent tissues of PCa.B:Volcano map for DEGs between metastatic and in situ tissues of PCa.

将转移癌和原位癌的DEGs挑出来对其表达量进行热图分析。结果显示，DEGs差异显著，分组聚合明显（图3）。

图3 DEGs热图展示差异基因Fig.3 DEGs heat map showing the differential genes.The x-axis represents the genes,and the y-axis represents the sample groups;red represents up-regulated genes and green the down-regulated genes.

2.3 WGCNA鉴定基因模块

选取样本间MAD方差最大的前7500个基因进行WGCNA 分析。以相关系数0.9 作为标准，使用pickSoftThredhold函数，选择邻接矩阵权重参数（软阈值）β=6为标准构建基因模块（图4）。采用两步法，设置每个基因模块最小基因数目为30，设置0.25为砍枝和合并模块高度，最终得到9个模块，其中grey模块表示未入组的基因集。

图4 WGCNA分析Fig.4 WGCNA analysis.A:Topological structure analysis of soft threshold parameters.B:Hierarchical clustering tree based on the difference of adjacent values.

2.4 模块与临床性质关联性分析

转移组与Green模块及Purple模块密切相关，原位癌与Brown模块密切相关（图5）。

图5 模块与临床特征相关性Fig.5 Correlation between modules and clinical characteristics.Red represents positive correlation,the numbers represent correlation coefficients,and the numbers in parentheses represent P values.

2.5 关键模块的功能分析与信号通路分析

GO 富集分析前列腺癌转移组的Green 模块及Purple模块，及癌原位组的Brown模块，显示：Purple模块主要与干细胞分化、细胞分裂、免疫呈递、炎症反应等生物学过程相关，说明转移过程中干细胞相关信号通路在这过程中的重要作用；而Brown模块则主要是细胞生长、氨基酸代谢、细胞凋亡等生物学过程相关（图6）。

图6 各基因模块的GO富集分析相关生物学过程Fig.6 GO enrichment analysis of the biological processes related with each gene module.

2.6 TCGA数据库枢纽基因的验证

模块基因与DEGs取交集，得到了与前列腺癌发生相关的多个枢纽基因（如BDH1、PAK4、EXTL3）及与前列腺癌转移密切相关的枢纽基因（如NKTR、CTBP2、HNRNPA2B1）。TCGA数据库分析结果显示在两组之间均差异显著（图7）。以中位数为界分成基因表达量高和低两组进行生存分析，结果显示枢纽基因均与前列腺癌病人的预后密切相关（图8）。

图7 TCGA数据库中各个基因在前列腺癌及对照组中的表达量Fig.7 Expression of BDH1 (A), PAK4 (B), EXTL3 (C), NKTR (D), CTBP2 (E) and HNRNPA2B1 (F) in prostate cancer and control group in the TCGAdatabase.

图8 TCGA数据库中基因表达量与前列腺癌病人的生存分析Fig.8 Gene expression in TCGA database and survival probability of prostate cancer patients based on expression level of BDH1(A),PAK4(B),EXTL3(C),NKTR(C),CTBP2(E)and HNRNPA2B1(F).

2.7 转移相关枢纽基因的实验验证

小分子干扰片段成功的将前列腺癌细胞系的HNRNPA2B1基因进行了沉默（图9A、B）。基因沉默HNRNPA2B1后细胞的生长受抑制（图9C、D），显著诱导细胞凋亡（图9E、F），抑制前列腺癌细胞系的克隆形成能力（图10A、B），且细胞的迁移和侵袭能力显著下降（图10C～F）。

图9 小分子沉默HNRNPA2B1观察对前列腺癌细胞系（PC3及LNCap）生长和凋亡的影响Fig.9 Effect of HNRNPA2B1 siRNA on growth and apoptosis of prostate cancer cell lines (PC3 and LNCap).A,B:Western blotting showing the efficiency of HNRNPA2B1 siRNA (NC is the control group).C,D:MTT showing the growth of cells in each group at 24,48,72 and 96 h.E,F:Flow cytometry showing cell apoptosis in each group.*P＜0.05,**P＜0.01,***P＜0.01.

图10 小分子沉默HNRNPA2B1观察对前列腺癌克隆形成、细胞系迁移和侵袭的影响Fig.10 Effect of HNRNPA2B1 siRNA on colony formation,migration and invasion of the prostate cancer cell lines.A-B:Colony formation experiment.C-F:Transwell assay to observe cell migration and invasion(Scale bar:20 μm).*P＜0.05,**P＜0.01.

3 讨论

前列腺癌的潜伏期长、患病率高、发病率高、死亡率高，是一个亟待解决的医学问题［17］。由于前列腺癌发病的分子机制非常复杂，涉及多种基因改变［18］。因此深入阐明前列腺癌发生和转移的分子机制有助于为开发前列腺癌治疗新靶点提供依据。

在本研究中，我们通过GEO数据库下载了前列腺癌大样本对列研究全基因组芯片数据GSE6919，经处理注释后，进行PCA主成分分析，通过R语言limma包分析癌原位与癌旁、癌转移与癌原位组别的差异表达基因（DEGs），发现该数据来源的原位癌与癌旁基因分组不明显，而癌转移组则具有明显的聚类且基因表达有显著差异；WGCNA构建基因共表达网络，分析识别与前列腺原位癌及转移癌组相关模块，对密切相关模块进行功能注释（GO），并利用可视化筛选出关键基因，发现了与前列腺癌发生相关的多个枢纽基因，如BDH1、PAK4、EXTL3及与前列腺癌转移密切相关的枢纽基因，如NKTR、CTBP2、HNRNPA2B1。进一步利用TCGA公共数据对关键基因进行验证，发现BDH1、PAK4、EXTL3、NKTR、CTBP2及HNRNPA2B1基因在前列腺癌中表达有差异且与前列腺癌病人的预后密切相关。我们选择了与前列腺癌转移密切相关的高表达基因HNRNPA2B1进行实验验证。合成并验证了HNRNPA2B1 siRNA的小分子干扰片段，对前列腺癌转移灶细胞系LNCap和PC3细胞系转染后，发现干扰HNRNPA2B1后细胞的生长增殖能力变缓，细胞凋亡增加，克隆形成能力减弱，且迁移和侵袭能力显著下降。通过数据库分析挖掘与细胞实验的双重验证，显示HNRNPA2B1基因的高表达与前列腺癌转移密切相关。

HNRNPA2B1 在人体组织中发挥功能，如信使RNA 加工、参与端粒的维持和细胞增殖的调控［19,20］。HNRNPA2B1 在多种肿瘤细胞系中高表达，HNRNPA2B1在胃癌［21］、乳腺癌［22］、肝癌［23］、胶质母细胞瘤［24］、肺癌及胰腺癌中高表达，并参与细胞凋亡和上皮-间充质转化（EMT），其可作为癌症发展阶段的一种生物标志物［18,25］。HNRNPA2B1是ΚRAS激活AΚT-mTOR信号以及PDAC细胞存活和小鼠肿瘤形成所必需的，且HNRNPA2B1可能是胰腺癌的治疗靶点［26］。HNRNPA2B1可通过调节信号转导子和转录激活子3（STAT3）以及细胞外信号调节激酶1/2（ERΚ1/2）信号通路，从而增加细胞的致瘤潜能［22］。抑制HNRNPA2B1可抑制丙酮酸激酶同工酶M2（PΚM2）［27］和肿瘤蛋白P53 诱导核蛋白2（TP53INP2）的选择性剪接［28］，已知PΚM2是肿瘤中葡萄糖代谢的关键酶，而TP53INP2是侵袭性细胞运动的重要组成部分。HNRNPA2B1是低氧诱导肿瘤存活的重要分子靶点［29］，它参与低氧诱导肿瘤存活。近期研究还发现HNRNPA2B1 可识别病原DNA和放大IFN-α/β生产。当HNRNPA2B1易位到细胞质中，激活TBΚ1-IRF3通路，导致IFN-α/β生产，还可促进N6-甲基腺苷（m6A）的修饰和CGAS、IFI16 和STING mRNAs的核质转运［30］。提示HNRNPA2B1功能多样，是多种肿瘤的致癌因子和潜在的治疗靶点。

癌症基因组图谱(TCGA)对大量人类肿瘤进行了剖析和分析，发现DNA、RNA、蛋白质和表观遗传水平上的分子畸变［31］。由于国家癌症研究所和各研究所的严格控制，TCGA数据质量很高，癌症/数据集之间的可比性比其他研究要高得多。我们通过GEO数据库得到前列腺癌基因芯片数据后，进行加权基因共表达网络分析（WGCNA），找到差异基因得到BDH1、PAK4、EXTL3、NKTR,CTBP2、HNRNPA2B1共6个枢纽基因。再通过TCGA数据库对上述基因进行验证，发现了与前列腺癌转移密切相关的枢纽基因，如NKTR、CTBP2、HNRNPA2B1。为进一步验证其与前列腺癌转移的相关性，我们选取高表达基因HNRNPA2B1，进行了细胞功能的实验验证，发现沉默HNRNPA2B1基因后细胞生长、增殖、迁移和侵袭能力显著减弱。因此权基因共表达网络分析（WGCNA）是鉴别疾病发生发展枢纽基因的可靠工具。

综上，我们基于加权基因共表达网络分析GEO数据库中的前列腺癌基因芯片数据，加权基因共表达网络分析可计算各基因间的Pearson相关系数，通过选择适当的阈值β进行网络构建。通过计算基因及模块显著性以及模块与临床性状显著性，找到前列腺癌进展的枢纽基因；再通过TCGA数据库及细胞实验验证了枢纽基因与前列腺癌转移的关系。因此加权基因共表达网络分析可以发现前列腺癌发生和转移的枢纽基因，对研究前列腺癌发生发展的调控机制及临床防治靶点提供了新思路。