0代码复现：如何对关联疾病的Hub基因进行筛选和验证？

桑格助手 2023-05-29 18:18

近年来，特别是在亚洲，克罗恩病（CD）的患病率和发病率迅速上升。病因一般认为是多种因素引起的，包括遗传易感性、环境和肠道菌群导致的粘膜免疫反应异常和上皮屏障功能受损。关于通过Sangerbox工具复现文章内容的教程，数据和视频链接放在最后，不要忘记领取哦。

文章名：Neutrophil Infiltration Characterized by Upregulation of S100A8, S100A9, S100A12 and CXCR2 Is Associated With the Co-Occurrence of Crohn's Disease and Peripheral Artery Disease

发表期刊：Frontiers in Immunology

影响因子：8.786

发表日期：2022.6.20

研究结果：

1、共鉴定出54个常见的DEGs(2个下调，52个上调)。

2、经验证，S100A8/S100A9/S100A12/CXCR2为Hub基因（AUC>0.7）。

3、中性粒细胞浸润与Hub基因上调有关。中性粒细胞活化、中性粒细胞趋化、中性粒细胞迁移等免疫过程通路与S100A8、S100A9、S100A12和CXCR2在CD和PAD中的高表达均显著相关。

数据采集

GSE95095、GSE134431、GSE120642和GSE111889数据集从基因表达综合(GEO)数据库下载。以克罗恩病和外周动脉疾病为关键词搜索相关基因表达数据集。PAD组每组至少15个样本，CD组每组至少50个样本；样本均为成年人。GSE120642数据集包含15个健康对照和36个PAD样本的腓肠肌RNA测序结果。GSE111889包括50名健康对照组和126名CD患者的活检样本。GSE95095（CD）和GSE134431（PAD）作为验证集。

鉴定差异表达基因

用R包edgeR对GSE111889数据进行归一化并鉴定DEGs。对于GSE120642，健康和PAD的DEGs用limma包分析。计算差异表达基因的差异倍数(FCs，P<0.05、|logFC|>1.0），用维恩图获取两个数据集的共有DEGs 。

在CD和健康样本间筛选出850个DEGs（图1A）。与健康对照组相比，PAD患者中有433个DEGs(图1B)。绘制维恩图，识别出63个DEGs(图1C)。剔除表达趋势相反的基因后，GSE120642和GSE111889共发现54个表达趋势相同的DEGs，其中52个上调，2个下调。

图1 共有的差异表达基因

基因功能富集分析

为了进一步揭示DEGs功能，用R“cluster profiler”包进行基因本体论(GO)注释和京都基因和基因组百科全书(KEGG)通路富集分析。使用标注、可视化和集成发现数据库(DAVID) 在线工具对重要模块进行GO功能富集和KEGG途径富集分析（P<0.05）。

GO分析表明，从生物学过程来看，基因主要富集于中性粒细胞趋化、中性粒细胞迁移和粒细胞趋化。细胞组成上，基因富集在分泌颗粒腔、细胞质囊泡腔和囊泡腔。分子功能方面，基因主要富集于RAGE受体结合、G蛋白偶联受体结合和免疫受体活性(图2A)。

KEGG分析显示，DEGs主要富集在病毒蛋白与细胞因子和细胞因子受体、细胞因子-细胞因子相互作用受体和IL-17信号通路的相互作用中(图2B)。

图2 基因功能富集分析

构建DEGs的PPI网络

基于已识别的常见DEGs，使用互作用基因数据库检索工具(STRING) (https://cn.string-db.org/)构建PPI网络。用Cytoscape 3.9.0绘图(置信度>0.4)。

STRING技术筛选DEGs的PPI网络，包括50个节点和89条边（P<1.0e-16，图3A)。MCODE检测关键的聚类模块。从使用通用DEGs构建的PPI网络中检索了两个模块。模块1包含8个节点和27条边，聚类评分(密度*成员数)为7.714(图3B)。模块2有3个节点和3条边，得分为3(图3B)。对最重要的模块1进行GO分析和KEGG富集分析。模块的基因功能与RAGE受体结合和中性粒细胞趋化有关，主要富集在IL-17信号通路。

用插件CytoHubba识别Hub基因，用插件MCODE从网络复合体中筛选核心基因的重要模块（Degree Cutoff=2, Node Score Cutoff=0.2, K-Core=2, MaxDepth=100）。同时使用多种拓扑分析算法（MCC、MNC、Degree、EPC）预测和探索PPI网络中排名前15的Hub基因。四种算法交集显示了13个候选Hub基因:HK3、S100A12、FCGR1A、S100A9、SLC11A1、FPR1、SAA1、S100A8、CXCR2、FCN1、AQP9、CSF3R和HP(图3C)。

图3 PPI网络与关键模块

Hub基因验证

用GSE95095和GSE134431验证Hub基因的表达水平。GSE95095包含24个CD和12个对照样品。GSE134431包含8份PAD的糖尿病溃疡样本及13份对照样本。两组数据集的比较采用Wilcoxon检验（P<0.05）。用pROC包生成受试者工作特征(ROC)曲线，评估Hub基因预测的准确性。用GeneMANIA数据库分析以Hub基因为核心的DEGs的基因-基因相互作用网络。

CD在GSE95095中验证，PAD在GSE134431中验证。结果表明，与正常组织相比，CD和PAD中只有四个中心基因显著上调。Hub基因为S100A8、S100A9、S100A12和CXCR2(图4)。

图4 验证集中显著上调的Hub基因

Hub基因效能评价及PPI构建

对所有数据集中这四个Hub基因绘制ROC曲线，评估诊断特征的准确性(图5)。在所有数据集中，这四个基因的AUC值都大于0.7，说明预测能力优秀。表明S100A8、S100A9、S100A12和CXCR2是诊断PAD和CD的有前景的标志物。

图5 Hub基因的ROC曲线

通过GeneMANIA对四个Hub基因及其20个相互作用基因进行PPI分析，预测共定位、共享蛋白结构域、共表达、与通路间的相关性(图6)。预测基因位于外圈，Hub基因位于内圈。网络显示这些基因在粒细胞趋化、髓系白细胞迁移、粒细胞迁移、中性粒细胞迁移、白细胞趋化和白细胞迁移中富集（图6）。

图6 基因-基因相互作用网络

Hub基因与免疫浸润

ssGSEA(单样本基因集富集分析)算法表示每个样本中特定基因集的绝对富集程度，用于量化CD或PAD组织中免疫细胞的浸润，确定数据集样本的免疫浸润水平。基于GeneMANIA生物学功能分析结果，通过Spearman相关系数确定Hub基因与23个免疫细胞间的相关性。

分析4个Hub基因与23个免疫细胞的相关性。对于GSE120642的PAD样本，中性粒细胞浸润水平与S100A8、S100A9、S100A12和CXCR2显著相关（图7A）。此外，S100A8和S100A9与B细胞、MDSCs、T细胞和树突状细胞相关。所有Hub基因均与中性粒细胞呈显著相关(r>0.7, P<0.05，图7B)。

在GSE111889的CD样本中，除活化B细胞、浆细胞样树突细胞、CD56 dim NK细胞、CD56 bright NK细胞和Type17辅助T细胞外，S100A8、S100A9、S100A12和CXCR2与大多数细胞类型呈正相关(图7C)。所有Hub基因均与中性粒细胞显著相关(r≥0.7,P<0.05，图7D)。

图7 Hub基因与免疫浸润的相关性

Hub基因的GSEA结果

基因表达水平设置为群体表型，GSEA 评估两组间的相关通路和分子机制（P<0.05、|归一化富集评分(NES)|>1、假阳性率(FDR) q<0.25），以确定GSE120642、GSE111889数据集中，Hub基因高表达组和低表达组间的差异调控通路，确定在PAD和CD中激活的信号通路。

对“gobp_阳性调控中性粒细胞激活”、“gobp_中性粒细胞外渗”、“gobp_中性粒细胞激活参与免疫反应”、“gobp_中性粒细胞趋化”、“gobp_中性粒细胞迁移”和“biocarta_il17通路”基因集进行GSEA分析。中性粒细胞活化、中性粒细胞趋化、中性粒细胞迁移等途径的激活与GSE120642(图8A)和GSE111889(图8B)中S100A8、S100A9、S10012和CXCR2的高表达密切相关。