候圣祥 ,张敬芳,胡 俊,刘 璐,陈天禹,庞 燊
雄激素性脱发(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的潜在靶标。
本项目的AGA基因表达矩阵(GSE90594)来自GEO数据库。包含14份AGA患者头皮样本,和14份正常人样本,平台为GPL17077。使用R语言程序包(V4.0.1)对芯片的数据加以分析,对样本进行聚类分析后去除离群样本,利用WGCNA程序包分析芯片数据,选取与AGA正、负相关性最高的模块,使用clusterprofiler软件包完成基因组百科全书(KEGG)及基因本体论(GO)、limma软件包分析基因的表达。
利用WGCNA程序包分析芯片数据,选取β值为5(图1)用来估算网络拓扑重叠TOM,通过层次聚类分析链接值的差异获得的基因聚类树,基于动态分层剪切数法运算来获得模块。计算基因模块的特征值(ME),联系临床信息,对ME进行分段聚类并重新排列树状图,设置高度值0.7为分割线,合并相似程度较高的基因模块,再用剪切后的模块替换新的聚类树和模块图。
图1 WGCNA中的软阈值筛选图
选取与AGA正、负相关性最高的模块各一个,通过对模块内基因分析基因临床特征相关性(GS)和基因模块表达水平(MM)的计算和取绝对值,筛选出高度正相关或负相关的基因,从而识别关键枢纽基因。
本项目选取与临床显著相关的基因模块,计算模块内基因的MM和GS值,将|MM|>0.8且|GS|>0.6的基因导出,它们对具有模块主要的调控作用。用Cytoscape软件绘制模块中共表达基因网络关系[5]通过基因之间加权的共表达关系使用Cytoscape绘制网络图,通过MCODE[6]插件寻找Hub基因。
选取与临床表型显著相关的基因模块,使用clusterProfiler软件包进行GO分析和KEGG富集分析[7]。
LASSO具有较强的预测值和相关性高,并且可以用于高维数据等功能[6,7]。为了区分AGA与对照组,通过glmnet包提取了基因的表达谱以构建LASSO模型。结果使用来自LASSO分析的回归系数为每个样本创建模型索引,以下列公式加权所选基因的表达值:索引= ExpGene1×Coef1 + ExpGene2 × Coef2 +ExpGene3×Coef3。“Coef”是基因的回归系数,它是从LASSO Cox回归得出的,“Exp”表示基因的表达值。为了评估LASSO模型识别AGA的能力,使用pROC软件包[8]在自身进行ROC曲线分析。
GSE90594芯片数据及临床信息的下载及预处理,过滤其中探针信息注释不全和重复的基因,最终获得29 499个基因对应27个样本的表达矩阵,绘制样本的分层聚类图与对应的临床信息的热图(图2)。
图2 样本对应性状热图
通过WGCNA包的算法,根据无尺度网络分布拟合,选取5作为软阈值,并计算基因间的相关性矩阵和TOM,使用TOM构建基因间分层聚类树,同时使用动态剪切树的方法把基因分成20个模块(图3)
图3 模块的动态切割图
通过计算各个模块与临床表型之间的关系,绘制基因模块与临床表型热图,从基因模块与临床表型的热图中可以看出brown模块与turquoise模块与临床分型显著相关。能体现出该评分高低在正常人和AGA患者间基因表达的差异有非常强的联系。
通 过 |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模块基因共表达网络图
选取与临床高度相关的brown模块和turquoise模块进行GO与KEGG富集分析(图6)。GO分析发现brown模块的基因主要参与淋巴细胞活化的调节、免疫效应过程的调节、T淋巴细胞活化、细胞粘附的正向调节、细胞激活的正调节、B淋巴细胞活化、B淋巴细胞活化的调节等生物功能;turquoise模块的基因主要参与表皮发育、皮肤发育、细胞器分裂、分化形成角质形成细胞、染色体分离、角化、减数分裂细胞周期、翻译起始等生物功能。KEGG分析发现brown模块的基因参与内吞作用、造血细胞系、弓形虫病、T淋巴细胞受体信号通路、补体与凝血级联、利什曼病、原发性免疫缺陷等重要通路;turquoise模块的基因参与河马(hippo)信号通路、Wnt信号通路、细胞周期、黑素生成、基底细胞癌等通路。
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曲线
本项目对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选取的对象为发病样本,缺少发病前期样本,无法分析基因在疾病各个时间及空间的转录水平;且为头皮样本,无法进一步分析在不同种类细胞中转录组的表达。今后可做不同发病时间段的单细胞转录组数据处理,可以更加清晰的表现各个基因在不同时间段,不同种类细胞中表达的差异。