基于加权基因共表达网络分析筛选SCAP基因预测脓毒症发生

郝 谦，陈 勇

脓毒症是一种异质性、复杂的病理生理综合征，可引起多器官功能障碍，死亡率高。据估计，全世界每年约有800万人死于脓毒症，主要原因为脓毒症休克，而循环、细胞和代谢异常可显著增加死亡率[1，2]。这些不良的结果可能是由于缺乏对脓毒症分子机制的理解。近年来发展的高通量测序技术可以快速明确脓毒症发病机制，以及脓毒症相应的标记物，但是传统的筛选差异基因表达极易遗漏调控过程中的核心分子，同时仅基于差异基因筛选难以对生物系统整体探索研究。通过系统的绘制个体生物网络互作图可以精准地寻找出与预后相关的核心分子[3，4]，加权基因共表达网络分析（weighted gene co-expression network analysis，WGCNA）可以有效地解决上述问题[5，6]。笔者拟利用该方法探索出与脓毒症发生、预后相关的关键基因，为今后的研究提供线索。

1 材料与方法

1.1 实验材料

选择Gene Expression Omnibus（GEO）数据库（https://www.ncbi.nlm.nih.gov/geo/）。GEO数据库为目前最常用的非肿瘤研究数据库，其内含有比较齐全及数据量最多的数据库。

1.2 方法

1.2.1 数据检索方法与要求

通过GEO数据库检索脓毒症，下载时间：2020年12月28日。

纳入分析数据集要求为：①数据量大于100；②含有脓毒症患者和健康对照组；③脓毒症患者有预后情况。

1.2.2 加权基因共表达网络构建和模块筛选

WGCNA是一种常用的基因共表达网络筛选技术，已被广泛应用于识别和筛选各种疾病生物标记物[7]。第一步，对每个样本基因名称进行转换，对基因表达值进行质检，检测是否有缺失值，同时对每个样本表达值进行层次聚类分析，筛选是否存在离群样本，基于上述筛选，对缺失基因表达值及离群样本予以剔除。第二步，通过R软件中的“WGCNA”分析包构建基因共表达网络[4，8]。第三步，构建基因之间的相关矩阵，确定软阈值大小。第四步，构建拓扑重叠矩阵（topological overlap matrix，TOM），用于度量基因的网络连通性[3，9]。第五步，基于TOM数据将表达谱相似的基因归为同一个基因模块，并进行连锁层次聚类[4，8]。第六步，计算每个模块在每个数据集中所占权重，计算出权重最大模块，用于后续研究。

1.2.3 核心基因筛选

在对基因模块分析时，通过对上述模块权重计算获得占比最大模块TURQUOISE模块，为进一步识别模块中核心靶基因，首先通过对模块基因进行皮尔森相关分析，并进一步通过STRING（search tool for recurring instances of neighbouring genes）网站进行了蛋白与蛋白之间的相互作用分析（protein protein interaction，PPI）（选择连接数最大10个），并利用Cytoscape软件中的MCODE插件对模块中基因进一步评分筛选。综合上述的3个结果筛选真正的核心靶基因，SCAP。

1.2.4 GO功能学、KEGG通路分析和基因富集分析

为了解TURQUOISE模块主要涉及的功能学和通路。采用标准富集计算方法进行基因本体（gene ontology，GO）功能分析和京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路分析用以筛选与其相关的功能和通路。为了解SCAP基因在脓毒症患者中的潜在功能，利用实验所纳入的2个数据集进行基因富集分析（gene set enrichment analysis，GSEA）。其注释基因集c2.cp.kegg.v5.2，选择基因通路列表（gene lists per pathway，GMT）数据库作为参考基因集。将P＜0.05和富集评分（enrichment score，ES）＞0.3作为截止标准，并选择2个数据集中一致的通路认为是其可能影响通路。

1.3 统计学方法

统计分析采用的是SPSS 24.0，GraphPad Prism 7.0和R 3.4.1，同时用上述软件进行图像生成处理。t检验用来分析两个组别之间的平均数的差异。绘制受试者工作特性 （receiver operating characteristic，ROC）曲线，评估核心靶基因的预测能力，利用曲线下面积（area under curve，AUC）评估灵敏度和特异度。P＜0.05为差异有统计学意义。

2 结果

2.1 筛选样本结果

检索GEO数据库并按筛选要求得到了2个数据集GSE26378和GSE54514。GSE54514含有163例患者全血芯片数据，其中健康对照组36例，脓毒症患者127例，脓毒症患者中其中死亡31例；含有男性64例，女性99例；平均年龄55.56岁（标准差17.21岁）；使用基因芯片平台为GPL6947。用于实验中进行外部验证数据集GSE26378有103例患者全血芯片数据，其中健康对照组21例，脓毒症患者82例，脓毒症患者中其中死亡12例；平均年龄53.75岁（标准差3.21岁）；使用基因芯片平台为GPL570。

2.2 构建共表达模块

基于21 085个基因在163例脓毒症患者的表达数据，利用WGCNA方法对每个样本基因表达值及样本表达值进行层次聚类分析，163例样本无剔除离群样本。进一步利用基于无序列网络的WGCNA方法，首先对软阈值进行筛选（图1），当软阈值选择为6时即可获得较好的连接关系。将基因进行模块化富集分析（图2），将基因依据其各相关表达量进一步进行分类。共筛选获得了22个相应的基因模块（图3），依据各模块在各个样本上皮尔森系数绝对值相加为最高者认定为响应系数最高模块，最后筛选获取了TURQUOISE模块，TURQUOISE模块共含有1 377个基因。

2.3 TURQUOISE模块相关功能分析及核心分子的筛选

为寻找TURQUOISE模块中主要富集的细胞功能和通路，基于GO功能富集及KEGG数据库中信号通路的上下游关系，依据TURQUOISE模块中所含有的1 377个相关基因筛选表达响应基因参与的信号通路的关联图（图4A、B）。研究结果提示，发现多条极为相关的信号通路，即RNA调控表达异常、RNA剪切组成等。为进一步寻找TURQUOISE模块中真正的核心分子，同时为避免遗漏真正的核心分子，将TURQUOISE模块中的基因上传到STRING网站中，进行PPI，并利用Cystacape软件中复杂的分子检测（molecular complex detection，MCODE）插件对PPI进行评分筛选位于调控核心位置的基因（图5），并结合PPI网络所获得的基因连接数（纳入连接数前10位的基因）用于筛选真正的核心分子。将上述3种方法所筛选的结果综合分析，最后筛选出了固醇调节元件结合蛋白裂解激活蛋白（srebpcieavage activating protein，SCAP）基因为核心靶基因。

2.4 SCAP基因核心靶基因验证

为了进一步验证筛选的基因其在健康对照组和脓毒症组中，以及在死亡组和存活组中表达情况，实验利用建模数据集GSE54514，对SCAP基因进一步验证（图6）。SCAP基因在脓毒症患者存活组中表达较死亡组显著更高（9.036±0.032vs8.857±0.067。P＜0.01），SCAP基因在健康对照组表达较脓毒症组显著更高（9.190±0.073 vs 8.99±0.030。P＜0.01）。同时在ROC曲线下分析发现，SCAP基因对于发生脓毒症和脓毒症患者中死亡情况也有较好的预测能力（图7）。其在预测脓毒症患者死亡上的ROC AUC为0.64，P=0.02。在预测是否发生脓毒症上其AUC=0.62，P=0.03。上述结果进一步阐明实验所筛选的SCAP基因具有较好的预测脓毒症发生和死亡能力。

2.5 SCAP基因核心靶基因外部验证

为了明确实验所筛选的SCAP基因是否真正可以有效地预测脓毒症的发生和预后情况，进一步利用外部数据集GSE26378中对所筛选的SCAP基因表达进行验证其在健康对照组和脓毒症组中，以及在死亡组和存活组中表达情况（图8），SCAP基因在脓毒症患者存活组中表达较死亡组显著更高（0.886±0.030 vs 0.730±0.061。P=0.045），SCAP基因在健康对照组表达较脓毒症组显著更高（1.154±0.088 vs 0.863±0.078。P＜0.01）。同时在ROC曲线下分析发现，SCAP基因对于发生脓毒症和脓毒症患者中死亡情况也有较好的预测能力（图9）。其在预测脓毒症患者死亡上的AUC=0.69，P=0.04。在预测是否发生脓毒症上其AUC=0.71，P＜0.01。上述结果进一步阐明实验筛选的SCAP基因在外部数据集中也有较好的预测脓毒症发生和死亡能力，进一步验证了SCAP基因的可靠性。

2.6 GSEA功能学分析

为了解SCAP基因主要涉及何种功能调控，实验进一步利用GEO数据库中GSE54514和GSE26378数据集并采用GSEA其所涉及的下游通路。发现SCAP基因主要涉及了T细胞受体及糖代谢。见图10、11。

3 讨论

脓毒症是由于炎性细胞因子的失调导致体内细胞因子风暴最终造成多器官衰竭和死亡[10，11]，是目前临床中所面临的一个主要的公共卫生问题；在过去的10年里，重症监护医学取得了巨大的进步，但是脓毒症的死亡率却没有实质性的改善[12]。就目前研究结果而言，由于脓毒症发病过程复杂，可能涉及到复杂的全身炎症效应、免疫功能障碍、凝血功能异常、组织损伤及宿主对不同感染病原微生物和其毒素的异常反应等多个方面，与机体多系统、多器官病理生理改变密切相关，脓毒症的根本发病机制尚未明确[13]。笔者研究为寻找脓毒症发生的机制，通过精准的WGCNA算法识别出与脓毒症发生和预后显著的SCAP基因。SCAP基因在内部和外部数据集进行了相应的验证，均获得了良好的结果，提示笔者研究所筛选的基因具有较强的重复性和可靠性。

SCAP是细胞脂肪合成酶的表达调控因子，在调节体内甘油三酯和胆固醇水平中起着重要作用[14]。目前已有较多研究显示，SCAP通过结合胆固醇调节元件结合蛋白 （sterol-regulatory element binding proteins，SREBP）并将它们从内质网运输到高尔基体进行蛋白水解处理，蛋白水解释放的SREBP转录因子进入细胞核，从而促进胆固醇的合成和摄取[15]。当内质网中的胆固醇堆积阻止了SCAP/SREBP复合物的运输，随即中止SREBP的蛋白水解过程，导致靶基因转录水平下降[16]。因此，SREBP不再被加工，胆固醇合成和摄取受到抑制，胆固醇的稳态得到恢复。目前已有较多的研究揭示了SCAP在脂质代谢中所发挥的功能，其可能影响多种代谢通路等[17，1,8]。但目前尚无SCAP在脓毒症中的作用及相关机制的报道。

笔者通过分析GEO中含有脓毒症发生和预后相关的数据集GSE54514，首先利用了基于无序基因表达网络的WGCNA方法，识别出对每个样本具有最大影响作用的TURQUOISE模块，为筛选出TURQUOISE模块中最为核心的靶基因，即可能为数据集最为核心的靶基因，运用了多种网络识别算法，最终获得了SCAP基因。为进一步验证笔者研究上述算法是否真正可行，利用内部数据集GSE54514通过验证SCAP基因在脓毒症患者和健康对照组中SCAP基因表达显著在健康对照组更高，提示SCAP基因高表达可能预示患者发生脓毒症概率更低；同时在脓毒症患者中死亡组和存活组患者的SCAP基因表达差异，结果提示存活组患者SCAP基因表达更高。上述结果进一步提示SCAP基因在脓毒症发生和发展过程中可能作为一个保护因素。为进一步验证笔者研究通过内部数据集筛选所获得数据是否可靠，再次利用GEO数据集筛选获得了含脓毒症和健康对照组数据集GSE26738。并通过差异表达分析和ROC曲线验证了SCAP基因在脓毒症患者中较健康对照组低表达，同时在死亡患者中表达也较存活患者低。进一步说明了SCAP基因作为脓毒症中保护因子的可靠性。为进一步探索SCAP基因所涉及的通路和功能，实验通过GSEA分析筛选了SCAP基因下游通路。结果提示SCAP基因可能通过T细胞受体及糖代谢两个功能产生作用，也有研究显示其与脓毒症发生、发展相关。进一步佐证了笔者研究的结果。因此，笔者研究结果对未来探索或者应用SCAP基因表达与脓毒症患者发生和预后奠定了一定的基础，具有较强的指导意义。但不可否认的是，笔者研究的结果还需要进一步在临床实际样本中获得更为可靠的验证，以及通过一系列的体内/体外实验验证SCAP基因的临床和生物学意义。

总之，笔者采用WGCNA等生物信息学方法利用GEO数据库研究了脓毒症患者发生和预后相关基因，并发现了SCAP基因的低表达与脓毒症发生和较差的预后相关。同时，笔者研究利用外部数据集对上述结果进一步验证。上述的SCAP基因有可能成为新的预后指标，有助于脓毒症患者个性化治疗及临床预后判断。