基于加权基因共表达网络分析识别雄激素性脱发的枢纽基因

候圣祥 ，张敬芳，胡 俊，刘 璐，陈天禹，庞 燊

雄激素性脱发（AGA）是导致临床成年人脱发的最主要的疾病。在我国，男性患病率约为21.3%，女性患病率约为6.0%[1]。近年来随着对AGA的流行病学调查的深入研究，发现遗传因素在发病的时间及程度过程中起着重要作用。尽管AGA不是威胁生命的疾病，但它会严重影响患者的精神状态和生活质量。目前，2%～5%米诺地尔溶液可通过局部应用促进头发生长[2]，而5α-R抑制剂非那雄胺和度他雄胺是经美国联邦食品药品监督管理局（FDA）批准可用于治疗AGA的口服合成药物[3]。但是这些药物均具有不良反应、功效有限，且必须持续用药的缺点[4]。因此需要更多方法来发现潜在的关键（枢纽、Hub）基因，进而揭示其分子机制。加权基因共表达网络分析（WGCNA）又称基因权重共表达网络分析，常用于分析样品基因与表型相关联系。因此，本研究基于基因表达数据库（GEO）中AGA数据，利用WGCNA挖掘AGA相关的基因网络，发现脱发中特异表达的关键基因，进一步分析和AGA相关的蛋白功能及通路，为AGA建立相关的基因功能和信号通路分析，以期获得治疗AGA的潜在靶标。

1 材料与方法

1.1 数据收集与数据预处理

本项目的AGA基因表达矩阵（GSE90594）来自GEO数据库。包含14份AGA患者头皮样本，和14份正常人样本，平台为GPL17077。使用R语言程序包（V4.0.1）对芯片的数据加以分析，对样本进行聚类分析后去除离群样本，利用WGCNA程序包分析芯片数据，选取与AGA正、负相关性最高的模块，使用clusterprofiler软件包完成基因组百科全书(KEGG)及基因本体论(GO)、limma软件包分析基因的表达。

1.2 加权共表达模块的计算

利用WGCNA程序包分析芯片数据，选取β值为5（图1）用来估算网络拓扑重叠TOM，通过层次聚类分析链接值的差异获得的基因聚类树，基于动态分层剪切数法运算来获得模块。计算基因模块的特征值（ME），联系临床信息，对ME进行分段聚类并重新排列树状图，设置高度值0.7为分割线，合并相似程度较高的基因模块，再用剪切后的模块替换新的聚类树和模块图。

图1 WGCNA中的软阈值筛选图

1.3 共表达模块与临床表现相关性的分析

选取与AGA正、负相关性最高的模块各一个，通过对模块内基因分析基因临床特征相关性(GS)和基因模块表达水平(MM)的计算和取绝对值,筛选出高度正相关或负相关的基因,从而识别关键枢纽基因。

1.4 构建基因网络与寻找Hub基因

本项目选取与临床显著相关的基因模块，计算模块内基因的MM和GS值，将|MM|＞0.8且|GS|＞0.6的基因导出，它们对具有模块主要的调控作用。用Cytoscape软件绘制模块中共表达基因网络关系[5]通过基因之间加权的共表达关系使用Cytoscape绘制网络图，通过MCODE[6]插件寻找Hub基因。

1.5 GO和KEGG富集分析

选取与临床表型显著相关的基因模块,使用clusterProfiler软件包进行GO分析和KEGG富集分析[7]。

1.6 LASSO模型的构建和ROC曲线分析

LASSO具有较强的预测值和相关性高,并且可以用于高维数据等功能[6,7]。为了区分AGA与对照组，通过glmnet包提取了基因的表达谱以构建LASSO模型。结果使用来自LASSO分析的回归系数为每个样本创建模型索引，以下列公式加权所选基因的表达值：索引= ExpGene1×Coef1 + ExpGene2 × Coef2 +ExpGene3×Coef3。“Coef”是基因的回归系数，它是从LASSO Cox回归得出的，“Exp”表示基因的表达值。为了评估LASSO模型识别AGA的能力，使用pROC软件包[8]在自身进行ROC曲线分析。

2 结果

2.1 芯片数据处理及加权基因共表达模块的构建

GSE90594芯片数据及临床信息的下载及预处理,过滤其中探针信息注释不全和重复的基因，最终获得29 499个基因对应27个样本的表达矩阵，绘制样本的分层聚类图与对应的临床信息的热图（图2）。

图2 样本对应性状热图

通过WGCNA包的算法,根据无尺度网络分布拟合，选取5作为软阈值，并计算基因间的相关性矩阵和TOM，使用TOM构建基因间分层聚类树，同时使用动态剪切树的方法把基因分成20个模块（图3）

图3 模块的动态切割图

2.2 共表达模块与临床表型的相关性分析

通过计算各个模块与临床表型之间的关系，绘制基因模块与临床表型热图，从基因模块与临床表型的热图中可以看出brown模块与turquoise模块与临床分型显著相关。能体现出该评分高低在正常人和AGA患者间基因表达的差异有非常强的联系。

2.3 共表达网络的可视化

通 过 |MM|＞ 0.8且 |GS|＞ 0.6的标准筛选brown模块和turquoise模块，分别在brown模块和turquoise模块中筛选出73、138个符合标准的基因。分别将筛选出的基因上传至STRING数据库，并通过Cytoscape对共表达网络基因间的相互作用关系的可视化，使用Cytoscape绘制网络图（图4），通过MCODE筛选枢纽基因(hubgenes)。在brown模块寻找出PDGFRA、PMP22、ZCCHC24、COL6A1、ISLR、PRRX1 和turquoise模 块 寻找 出PKP1、CALN1、PNMAL2、PPP5D1、GJB6、DSC2、GJA3共计13个枢纽基因（图5）。

图5 brown和turquoise模块基因共表达网络图

2.4 GO和KEGG富集分析

选取与临床高度相关的brown模块和turquoise模块进行GO与KEGG富集分析（图6）。GO分析发现brown模块的基因主要参与淋巴细胞活化的调节、免疫效应过程的调节、T淋巴细胞活化、细胞粘附的正向调节、细胞激活的正调节、B淋巴细胞活化、B淋巴细胞活化的调节等生物功能；turquoise模块的基因主要参与表皮发育、皮肤发育、细胞器分裂、分化形成角质形成细胞、染色体分离、角化、减数分裂细胞周期、翻译起始等生物功能。KEGG分析发现brown模块的基因参与内吞作用、造血细胞系、弓形虫病、T淋巴细胞受体信号通路、补体与凝血级联、利什曼病、原发性免疫缺陷等重要通路；turquoise模块的基因参与河马（hippo）信号通路、Wnt信号通路、细胞周期、黑素生成、基底细胞癌等通路。

2.5 构建LASSO模型

LASSO模型是AGA的潜在预测指标。本文提取了中枢基因的表达谱以构建LASSO模型（图7）。使用LASSO方法，鉴定了5个具有非零回归系数的基因，并且lambda.min的值= 0.02903。基于基因的模型索引按以下公式创建：索引=PRRX1 *(注:表示相乘)（1.9166078）+GJA3 * （-1.0668399）+DSC2*（-2.9185647）+COL6A1*（4.0211617）+CALN1*（0.3254326）。ROC曲线分析表明基于5基因的模型的训练集的AUC在训练集中为0.98，这表明LASSO模型可以用作AGA的生物标志物。

图7 5基因lasso模型预测自身ROC曲线

3 讨论

本项目对GSE90594数据集转进行权重基因共表达网络分析，成功构建20个模块，其中brown模块和turquoise模块与AGA密切相关的新基因共表达网络模块。笔者进一步探索与AGA相关的生物学过程。

功能富集分析表明与AGA高度相关的模块参与了细胞分化、表皮形成等生物过程，并显著参与WNT等通路。Wnt/β-catenin信号通路对毛囊干细胞的生长发育有着密切相关性，且在AGA中Wnt/β-catenin信号通路活性常常出现异常。Wnt/β-catenin信号通路在毛发周期期间对毛囊形态、发育和再生中起关键作用[8]。在毛发周期中，毛囊干细胞的端粒期-生长期转变是由Wnt/β-catenin信号通路与BMP途径之间的平衡相互作用驱动的[9]。在毛囊隆起和继发性发芽中，表皮干细胞和黑素细胞干细胞共同经历Wnt/β-catenin信号通路激活参与生发[10]。考虑到 Wnt/β-catenin 信号通路的改变会导致AGA中毛囊的状态的改变，因此可以通过以下方式调节Wnt/β-catenin信号通路，促进头发再生及发育[11]：①调节Wnt信号通路配体的分泌；②改变配体-受体的结合；③促进β-catenin向细胞核移位。对于AGA，由于β-catenin表达不稳定，导致处于生长期的毛囊干细胞减少。通过上调Wnt/β-catenin信号通路中β-catenin蛋白的活性可以延长发囊生长期。AGA患者的头皮在很大程度上保留了静止的毛囊干细胞，但缺乏祖细胞[12]。显然，毛囊干细胞向祖细胞转化的减少归因于雄激素介导的Wnt/β-catenin信号通路活性的异常[13]。调节雄激素受体（AR）反式激活并稳定Wnt/β-catenin信号通路可以治愈，逆转并预防AGA中的脱发，并将不良反应降至最低。除新疗法外，靶向激活Wnt/β-catenin信号通路途径对头发的生长也显示出促进的效果。

综上所述，PDGFRA、PMP22、ZCCHC24、COL6A1、ISLR、PRRX1的 表 达 与AGA呈 正 相 关，PKP1、CALN1、PNMAL2、PPP5D1、GJB6、DSC2、GJA3的表达与AGA呈负相关。本研究不足之处在于GSE90594选取的对象为发病样本，缺少发病前期样本，无法分析基因在疾病各个时间及空间的转录水平；且为头皮样本，无法进一步分析在不同种类细胞中转录组的表达。今后可做不同发病时间段的单细胞转录组数据处理，可以更加清晰的表现各个基因在不同时间段，不同种类细胞中表达的差异。