共表达网络分析筛选肾透明细胞癌进展和预后相关基因

肖观发 石胜军 胡万里

肾癌是西方国家10种最常见的癌症之一，占所有癌症的3%～5%[1]。2015年中国约有6.68万肾癌新发病例和2.34万死亡病例[2]。肾细胞癌(renal cell carcinoma， RCC)约占所有肾癌的90%，其中大部分(80%～90%)为肾透明细胞癌(clear cell renal cell carcinoma， ccRCC)[3]。RCC早期诊断和随访的生物标志物尚缺乏。据统计，超过50%的RCC是偶然发现的，约30%的患者诊断时已发生转移；此外，30%～50%的RCC患者在随访过程中会发生转移[4-5]。对于局限性RCC，手术是唯一有高质量的治疗方法，而对于转移性RCC患者则需要全身治疗[3]。然而，ccRCC通常对放、化疗不敏感，靶向治疗因其靶向特异性和低毒性可能是非手术治疗的最佳选择[3,6]，但是ccRCC肿瘤内的分子异质性可能会影响靶向治疗的应答，对靶向治疗的抵抗也是一个主要问题[7-8]；尽管行全身治疗，转移性RCC患者的预后仍很差。因此，早期诊断、个体化治疗和随访成为治疗成功的关键，寻找与ccRCC诊断、治疗和预后相关的生物靶点至关重要。

随着全基因组测序技术的快速发展，越来越多的肿瘤相关基因被快速检测出，肿瘤研究取得突破性进展[9]。虽然具有相似表达模式的基因可能在功能上是相关的，但在许多研究中，基因之间的高度相关性往往被忽略。加权基因共表达网络分析(weighted gene co-expression network analysis, WGCNA)是一种系统生物学算法，用于描述基因芯片样本中基因之间的相关模式，以及高共表达基因或模块簇与样本外部特征之间的关系[10]。简而言之，WGCNA可以识别高共表达基因簇，也可以筛选与临床特征相关的枢纽基因和模块。相关网络促进了基于网络的基因筛选方法，这些方法可用于识别候选生物标志物或治疗靶点。

蛋白相互作用(protein-protein interaction, PPI)，通常被理解为发生在细胞或活体生物中的蛋白之间分子对接的物理接触[11]，是一种新型分子治疗靶点[12]。预测PPI的生物数据库STRING可以构建PPI网络，以蛋白为节点，相互作用为边缘[13]。在PPI网络中“高度连接”的基因更有可能是重要蛋白质[14]。

本研究尝试构建加权共表达网络和PPI网络，以识别与ccRCC进展和预后显著相关的生物标志物。

对象与方法

一、数据采集与预处理

从基因表达谱(gene expression omnibus, GEO)数据库(http://www.ncbi.nlm.nih.gov/geo/)下载所需基因芯片集。GSE68417作为训练集构建共表达网络，包括14例正常肾脏组织样本和29例ccRCC组织样本。另外，GSE53757作为测试集来验证得到的结果，包括72例正常肾脏组织样本和72例ccRCC组织样本。此外，本研究还从癌症基因组图谱(cancer genome atlas, TCGA)数据库(https://genome-cancer.ucsc.edu/)下载了RNA测序数据集和相关临床数据，并将其作为一个测试集来验证研究结果。

二、差异表达基因筛选

采用R软件中的limma包[15]对正常肾组织标本和ccRCC组织标本之间的差异表达基因(differentially expressed genes， DEGs)进行筛选。错误发现率(false discovery rate， FDR)<0.05和│log2 FC│≥1.0为筛选条件，其中FC为两组间差异表达倍数。

三、加权共表达网络构建与枢纽模块筛选

利用R软件中的WGCNA包构建DEGs共表达网络。首先检测DEGs的表达数据图谱，检验其是否为良好的样本和基因，并在后续分析中剔除离群样本。步骤如下：①计算各基因间Pearson相关系数。②构建加权邻接矩阵aij=│sij│β(sij代表基因i和j的Pearson相关系数；aij代表i基因与j基因的邻接系数；β是一个软阈值)。③将邻接矩阵转化为拓扑重叠矩阵(TOM)，并计算基因间相异度矩阵dissTOM=1-TOM，然后对dissTOM进行层次聚类。④设置基因树状图的最小模量为30，合并相似度高于0.8的模块。

四、枢纽模块和枢纽基因的筛选

将各模块与临床性状关联后，计算各模块显著性(module significance， MS)。MS值越高，模块就越重要。通过比较各MS，与某一临床特征高度相关的模块被视为枢纽模块。然后计算基因显著性(gene significance， GS)和模块身份(module membership， MM)。枢纽模块中│MM│>0.8、│GS│>0.2的基因被视为备选枢纽基因。另外将该模块中的基因上传至STRING数据库[13]构建PPI网络，选取点度中心性≥14的基因作为备选枢纽基因。最终，筛选两个网络中共有的基因即为枢纽基因。

五、枢纽基因的验证

利用数据集GSE53757、RNA测序数据和人类蛋白质图谱数据库[16](http://www.proteinatlas.org)，探索正常肾组织和ccRCC组织中枢纽基因的mRNA和蛋白表达情况。在测试集GSE53757和RNA测序数据中进行线性回归分析，验证枢纽基因表达与ccRCC进展的关系。绘制枢纽基因的受试者工作曲线图，评价其对局限性(病理Ⅰ/Ⅱ期)和非局限性(病理Ⅲ/Ⅳ期)ccRCC的鉴别能力。利用RNA测序数据对枢纽基因进行生存分析。

六、GO及KEGG富集分析

利用R软件中的clusterProfiler包[17]对枢纽模块中的基因进行GO富集和KEGG通路分析。FDR<0.05被认为具有统计学意义。

七、基因集富集分析

数据集GSE53757用于枢纽基因的基因集富集分析(gene set enrichment analysis, GSEA)[18]。本数据中72例ccRCC样本根据枢纽基因的中位表达(高、低)分为两组。选择c2.cp.kegg.v6.2.symbols.gmt作为参考基因集，FDR<0.05被认为具有统计学意义。

结 果

一、DEGs筛选、共表达网络构建及枢纽模块确定

经过数据预处理，得到训练集GSE68417的表达式矩阵。选取共1 678个DEGs (634个上调基因，1 044个下调基因)构建共表达网络。选择β=4(R2=0.89)为软阈值(图1A),确定7个模块(图1B)。各模块与ccRCC的进展(分级)相关，通过比较各模块间的MS，确定棕色模块为枢纽模块(图1C)。

A:软阈值;B:基因模块;C:枢纽模块筛选图1 WGCNA软阈值、临床性状相关模块的确定

二、枢纽基因确定

在棕色模块中，通过│MM│>0.8、│GS│>0.2筛选得到27个备选枢纽基因。在PPI网络中，9个基因连接度≥14。取两者交集得到VWF、TEK和FLT1这3个基因，其被认为是与ccRCC进展密切相关的枢纽基因。

三、枢纽基因验证

基于数据集GSE53757和RNA测序数据,发现枢纽基因在正常肾组织和ccRCC组织中的表达差异(图2)。此外,免疫组织化学染色从人类获得蛋白质图谱数据库也证明了枢纽基因的表达(图3)。在检测集GSE53757和RNA测序数据中，线性回归分析证实所有枢纽基因与ccRCC进展(病理分期)呈负相关(P<0.05)(图4)。绘制受试者工作曲线评估枢纽基因对局限性(病理分期Ⅰ/Ⅱ)和非局限性(病理分期Ⅲ/Ⅳ) ccRCC的鉴别能力，曲线下面积值均>0.5(图5)。此外，生存分析显示，枢纽基因表达越低，整体生存状况越差，说明枢纽基因可能是ccRCC预后的生物标志物(图6)。

A:基于GSE53757数据集(VWF、TEK、FLT1);B:基于RNA测序数据(VWF、TEK、FLT1)图2 正常肾组织与ccRCC组织间枢纽基因转录水平的验证

A:VWF;B:TEK;C:FLT1图3 基于人类蛋白质图谱数据库正常肾组织和ccRCC组织间枢纽基因的翻译水平

A:基于GSE53757数据集(VWF、TEK、FLT1);B:基于RNA测序数据(VWF、TEK、FLT1)图4 各枢纽基因表达水平与ccRCC进展(病理分期)的相关性

A:基于GSE53757数据集(VWF、TEK、FLT1); B:基于RNA测序数据(VWF、TEK、FLT1)图5 各枢纽基因的受试者工作曲线

A:VWF;B:TEK;C:FLT1图6 RNA测序数据库中ccRCC各枢纽基因的总体生存分析

四、GO和KEGG富集分析

为进一步了解DEGs在枢纽模块中的作用，进行GO和KEGG富集分析。GO分析结果显示，枢纽模块的DEGs主要富集于10个生物过程，包括血管发育调控、肾小球发育、肾小球血管发育、内皮细胞迁移、肾脏系统血管发育、肾脏血管发育、上皮细胞迁移、上皮组织迁移、内皮细胞迁移调控和组织迁移。此外，KEGG通路富集分析结果显示，DEGs主要富集于HIF-1信号通路、Rap1信号通路、胆固醇代谢、碳代谢和类风湿关节炎(图7)。

A:GO分析;B:KEGG通路富集分析图7 枢纽模块中DEGs的GO和KEGG富集分析

为了解枢纽基因的潜在影响机制，进行GSEA研究。FDR<0.05被认为具有统计学意义，18个基因集普遍富集于枢纽基因低表达组，且大多数集中在与免疫相关的通路上。

讨 论

RCC是世界范围内最常见的肿瘤之一。当前迫切需要更有效的ccRCC诊断、治疗和预后的生物标志物。本研究通过构建共表达网络和PPI网络，以识别与ccRCC进展高度相关的枢纽基因，最终确定3个枢纽基因(VWF、TEK和FLT1)。基于数据集GSE53757、RNA测序数据和人类蛋白质图谱数据库的分析，这些基因被确认为是与ccRCC进展和预后高度相关的候选生物标志物，可有效区分局限性和非局限性ccRCC。高共表达基因通过WGCNA连接并分组成模块，不同的模块包含不同的功能相关基因。为了挖掘与ccRCC进展相关的潜在机制，本研究对枢纽模块的基因进行了GO富集和KEGG通路分析，结果显示血管生成相关通路在ccRCC进展中发挥了重要作用。

VWF编码一种参与止血、免疫应答、炎症、血管生成和肿瘤转移的糖蛋白[19]。关于VWF在肿瘤转移中的作用说法大相径庭[20]。一些研究发现VWF是肿瘤性疾病的一种消极预后标志物，可促进肿瘤的进展和转移[21-23]；然而，另有研究表明VWF可以作为肿瘤细胞存活的抑制剂，减少转移[24-25]。TEK编码血管生成素1和2的受体，属于酪氨酸激酶Tie2家族。近期研究表明，TEK高表达与乳腺癌生存率低相关[26-27]。Park等[28]发现TEK活化可诱导肿瘤血管正常化，促进血液灌注和化疗药物递送，显著减少乳酸酸中毒，降低肿瘤生长和转移。Ha等[29]研究表明TEK是一种潜在的ccRCC预后基因。TEK在肿瘤中的作用仍在探索中，它可能是一个极具研究价值的癌症治疗靶点。FLT1编码血管内皮生长因子受体家族成员，在血管生成中发挥重要作用。Qian等[30]证明，FLT1调控一组炎症反应基因，影响乳腺肿瘤转移。FLT1似乎与癌症进展相关，可能作为一种预后指标[31-33]。以上枢纽基因与肿瘤进展及预后密切相关，可能成为肾癌预后及治疗靶点。

GSEA研究枢纽基因的潜在影响机制，结果显示18个基因集普遍富集于枢纽基因低表达组，且大部分集中在与免疫相关的通路上。因此，我们可以假设这些枢纽基因可能通过调节免疫相关通路来影响ccRCC的预后。

目前已有关于ccRCC进展相关基因的研究报道，但其可靠性有待进一步验证。因此，本研究选取更多数据集来进行验证，大大提升了研究结果的可信度；另外，通过对数据库的分析，不仅验证枢纽基因的转录水平，更验证其翻译水平。本研究最大的局限性是缺乏外部验证，需进一步研究加以证实。

综上所述，本研究确定并验证了与ccRCC的进展和预后密切相关的3个枢纽基因(VWF、TEK和FLT1)，其可能是ccRCC的候选生物标志物。