基于加權基因共表達網絡分析識別雄激素性脫發(fā)的樞紐基因

候圣祥 ，張敬芳，胡 俊，劉 璐，陳天禹，龐 燊

雄激素性脫發(fā)（AGA）是導致臨床成年人脫發(fā)的最主要的疾病。在我國，男性患病率約為21.3%，女性患病率約為6.0%[1]。近年來隨著對AGA的流行病學調查的深入研究，發(fā)現遺傳因素在發(fā)病的時間及程度過程中起著重要作用。盡管AGA不是威脅生命的疾病，但它會嚴重影響患者的精神狀態(tài)和生活質量。目前，2%～5%米諾地爾溶液可通過局部應用促進頭發(fā)生長[2]，而5α-R抑制劑非那雄胺和度他雄胺是經美國聯邦食品藥品監(jiān)督管理局（FDA）批準可用于治療AGA的口服合成藥物[3]。但是這些藥物均具有不良反應、功效有限，且必須持續(xù)用藥的缺點[4]。因此需要更多方法來發(fā)現潛在的關鍵（樞紐、Hub）基因，進而揭示其分子機制。加權基因共表達網絡分析（WGCNA）又稱基因權重共表達網絡分析，常用于分析樣品基因與表型相關聯系。因此，本研究基于基因表達數據庫（GEO）中AGA數據，利用WGCNA挖掘AGA相關的基因網絡，發(fā)現脫發(fā)中特異表達的關鍵基因，進一步分析和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，通過層次聚類分析鏈接值的差異獲得的基因聚類樹，基于動態(tài)分層剪切數法運算來獲得模塊。計算基因模塊的特征值（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]。為了區(qū)分AGA與對照組，通過glmnet包提取了基因的表達譜以構建LASSO模型。結果使用來自LASSO分析的回歸系數為每個樣本創(chuàng)建模型索引，以下列公式加權所選基因的表達值：索引= ExpGene1×Coef1 + ExpGene2 × Coef2 +ExpGene3×Coef3?！癈oef”是基因的回歸系數，它是從LASSO Cox回歸得出的，“Exp”表示基因的表達值。為了評估LASSO模型識別AGA的能力，使用pROC軟件包[8]在自身進行ROC曲線分析。

2 結果

2.1 芯片數據處理及加權基因共表達模塊的構建

GSE90594芯片數據及臨床信息的下載及預處理,過濾其中探針信息注釋不全和重復的基因，最終獲得29 499個基因對應27個樣本的表達矩陣，繪制樣本的分層聚類圖與對應的臨床信息的熱圖（圖2）。

圖2 樣本對應性狀熱圖

通過WGCNA包的算法,根據無尺度網絡分布擬合，選取5作為軟閾值，并計算基因間的相關性矩陣和TOM，使用TOM構建基因間分層聚類樹，同時使用動態(tài)剪切樹的方法把基因分成20個模塊（圖3）

圖3 模塊的動態(tài)切割圖

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分析發(fā)現brown模塊的基因主要參與淋巴細胞活化的調節(jié)、免疫效應過程的調節(jié)、T淋巴細胞活化、細胞粘附的正向調節(jié)、細胞激活的正調節(jié)、B淋巴細胞活化、B淋巴細胞活化的調節(jié)等生物功能；turquoise模塊的基因主要參與表皮發(fā)育、皮膚發(fā)育、細胞器分裂、分化形成角質形成細胞、染色體分離、角化、減數分裂細胞周期、翻譯起始等生物功能。KEGG分析發(fā)現brown模塊的基因參與內吞作用、造血細胞系、弓形蟲病、T淋巴細胞受體信號通路、補體與凝血級聯、利什曼病、原發(fā)性免疫缺陷等重要通路；turquoise模塊的基因參與河馬（hippo）信號通路、Wnt信號通路、細胞周期、黑素生成、基底細胞癌等通路。

2.5 構建LASSO模型

LASSO模型是AGA的潛在預測指標。本文提取了中樞基因的表達譜以構建LASSO模型（圖7）。使用LASSO方法，鑒定了5個具有非零回歸系數的基因，并且lambda.min的值= 0.02903?；诨虻哪Ｐ退饕匆韵鹿絼?chuàng)建：索引=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信號通路對毛囊干細胞的生長發(fā)育有著密切相關性，且在AGA中Wnt/β-catenin信號通路活性常常出現異常。Wnt/β-catenin信號通路在毛發(fā)周期期間對毛囊形態(tài)、發(fā)育和再生中起關鍵作用[8]。在毛發(fā)周期中，毛囊干細胞的端粒期-生長期轉變是由Wnt/β-catenin信號通路與BMP途徑之間的平衡相互作用驅動的[9]。在毛囊隆起和繼發(fā)性發(fā)芽中，表皮干細胞和黑素細胞干細胞共同經歷Wnt/β-catenin信號通路激活參與生發(fā)[10]?？紤]到 Wnt/β-catenin 信號通路的改變會導致AGA中毛囊的狀態(tài)的改變，因此可以通過以下方式調節(jié)Wnt/β-catenin信號通路，促進頭發(fā)再生及發(fā)育[11]：①調節(jié)Wnt信號通路配體的分泌；②改變配體-受體的結合；③促進β-catenin向細胞核移位。對于AGA，由于β-catenin表達不穩(wěn)定，導致處于生長期的毛囊干細胞減少。通過上調Wnt/β-catenin信號通路中β-catenin蛋白的活性可以延長發(fā)囊生長期。AGA患者的頭皮在很大程度上保留了靜止的毛囊干細胞，但缺乏祖細胞[12]。顯然，毛囊干細胞向祖細胞轉化的減少歸因于雄激素介導的Wnt/β-catenin信號通路活性的異常[13]。調節(jié)雄激素受體（AR）反式激活并穩(wěn)定Wnt/β-catenin信號通路可以治愈，逆轉并預防AGA中的脫發(fā)，并將不良反應降至最低。除新療法外，靶向激活Wnt/β-catenin信號通路途徑對頭發(fā)的生長也顯示出促進的效果。

綜上所述，PDGFRA、PMP22、ZCCHC24、COL6A1、ISLR、PRRX1的 表 達 與AGA呈 正 相 關，PKP1、CALN1、PNMAL2、PPP5D1、GJB6、DSC2、GJA3的表達與AGA呈負相關。本研究不足之處在于GSE90594選取的對象為發(fā)病樣本，缺少發(fā)病前期樣本，無法分析基因在疾病各個時間及空間的轉錄水平；且為頭皮樣本，無法進一步分析在不同種類細胞中轉錄組的表達。今后可做不同發(fā)病時間段的單細胞轉錄組數據處理，可以更加清晰的表現各個基因在不同時間段，不同種類細胞中表達的差異。