加權重基因共表達網絡分析識別慢性腎臟病足細胞損傷關鍵節點基因MAGI2*

王萬鵬，唐伯玉，沈劍簫，李永剛，蒯巧林，高甄典，李娟，梁森林，陳皓瑜

（1.江蘇省漣水縣人民醫院，江蘇 淮安 225400；2.上海交通大學醫學院附屬仁濟醫院，上海 200127）

慢性腎臟?。╟hronic kidney disease,CKD）與足細胞損傷常?；橐蚬?，足細胞損傷機制極其復雜，多種理化、生物因素均可影響足細胞的功能。本研究擬采用加權基因共表達網絡分析（weighted gene coexpression network analysis,WGCNA）挖掘從美國國立生物技術信息中心（National Center for Biotechnology Information,NCBI）的基因表達綜合數據庫（gene expression omnibus,GEO）下載到的足細胞損傷體外模型以及CKD患者腎小球表達譜數據集，分析CKD狀態下腎小球足細胞損傷后加權基因共表達網絡關系，識別CKD狀態下足細胞損傷的相關基因模塊及模塊內樞紐（hub）基因，為進一步探索CKD足細胞損傷的發生、發展機制提供一種新思路。

1 材料與方法

1.1 數據來源

下載數據集GSE66107、GSE93798，GSE30528、GSE32591和GSE99339，其中GSE66107為嘌呤霉素氨基核苷（purinomycin amino nucleoside,PAN）誘導的足細胞損傷體外模型；GSE93798、GSE30528、GSE32591分別為通過激光微切割獲得的IgA腎病、糖尿病腎病（diabetic kidney disease,DKD）、狼瘡性腎炎（lupus nephritis,LN）患者腎小球mRNA表達數據，以上4個數據集用于差異表達基因（differentially expressed genes,DEG）分析；數據集GSE99339由多種CKD患者腎小球表達譜數據構成，疾病包括：LN、IgA腎病、膜性腎小球腎炎（membranous glomerulonephritis,MGN）、高血壓性腎?。╤ypertensive nephropathy,HT）、局灶節段性腎小球硬化（focal segmental glomerulosclerosis,FSGS）、腎腫瘤切除（tumor nephrectomy,TN）、微小病變腎?。╩inimal change disease,MCD）、DKD、 薄 基 膜 病（thin membrane disease,TMD），共133例患者用于WGCNA分析。見表1。芯片的探針注釋信息來自Affymetrix公司下載的原始文件。

1.2 數據預處理

采用R語言（3.4.0版）中的Affy包對原始數據進行預處理，包括ReadAffy函數讀取原始文件（后綴名為.CEL）；全局歸一化矯正（Robust Muti-array Average,RMA）函數歸一化及對數化；mas5calls函數檢測基因表達情況，函數結果表現為“表達（P）”或“缺失（A）”，本研究選擇在各數據集樣本中表達超過50%的探針用于后續研究；KNN法補充缺失值；limma包進行差異表達值計算，并用貝葉斯方法進行多重檢驗校正，最終選取閾值錯誤發現率（false discovery rate,FDR）＜0.05及倍數變化（fold change,FC）≥1.5得到DEG。

1.3 WGCNA

1.3.1 模塊構建及可視化 GSE99339數據集包含5個批次（bacth），使用sva包中ComBat函數去除批次效應后運行WGCNA包進行基因共表達網絡的構建與模塊鑒定。共表達網絡模塊內各基因利用DAVID 6.7（https://david.ncifcrf.gov/）在線工具進行基因本體（gene ontology,GO）分析。利用Cytoscape繪制WGCNA導出的加權共表達基因網絡，并通過插件Cluego（v2.3.5）與 Cluepedia（v1.3.5）關聯 String 數據庫，進行模塊內基因GO（功能）和京都基因組百科全書（Kyoto Encyclopedia of Genes and Genomes,KEGG）通路富集分析與蛋白間互作網絡（protein-protein interactions,PPI）可視化分析。Gluego插件參數設置為pV≤0.05，其余為默認參數。

1.3.2 足細胞特異性陽性標準基因 JU等[1]在過去的研究中使用高通量數據以機器學習為基礎發展處一種迭代算法，通過該算法分別定義統計學意義上的50種足細胞特異性標準基因（podocyte standard gene,PSG）。PSG包括：ACTN4、NPHS1及NPHS2等。

表1 數據集納入及分析結果

2 結果

2.1 足細胞損傷體外模型及CKD患者腎小球差異表達基因分析

經過預處理并刪除低表達探針（表達率＜50%數據集樣本數）后剩余探針分別為26 357、11 118、11 296及7 194個，合并重復探針和刪除無注釋探針后，最終分別得到6 154、3 882、1 422、840個DEG（FDR＜0.05，FC≥1.5）（見表1）。使用R語言ggplot包繪制DEG熱圖，可見所得DEG可以區分出大部分CKD患者或者損傷的足細胞（見圖1A～D）。本研究4個數據集共得到9 217個DEG，共同DEG（coDEG）15個，其中上調基因4個（見圖1E），下調基因11個（見圖1F），分別為CDKN1C，MYO5C，C1orf21，KBTBD11，NR3C2，FZD2，MAGI2，LPL，ZFPM2，USP46，PRKAR2B，TRIM38，DESI1，CTSS，ITGB2。

2.2 WGCNA構建基因共表達模塊

9 217 個DEG中有1 260個DEG表達趨勢在2個或2個以上數據集中存在相反（與對照組比較），刪除后剩余7 957個DEG，其中4 031個DEG存在于GSE99339預處理后的表達譜中，被用于WGCNA分析（見表1）。網絡模塊分析結果易受離群樣本的影響，因此在構建網絡之前去除離群的樣本數據顯得尤為重要。相關研究常采用芯片間相關度（inter.array correlation,IAC）來評估芯片數據的分布情況。IAC可以通過計算任意一對芯片上所有探針表達水平的相關系數獲得，芯片之間的關系可以用樹形圖來展示，平均IAC值較低的樣本，或者在樹形圖上無法聚類的樣本即為離群樣本，應被去除。本研究使用WGCNA包中的函數dist（）計算樣本間的IAC，并通過hclust（）函數繪制133個樣本4 031個DEG構成的樹形圖，結果未發現明顯離群值和疾病聚集（見圖2A）；根據動態分層剪切樹算法最終得到12個基因模塊（見圖2B），以顏色的英文命名，其中grey模塊包含518個DEG，表示未納入任何模塊的基因集合。

2.3 足細胞損傷相關基因模塊鑒定及分析

12個基因模塊與未歸類的grey模塊（見表2），本研究中共包含PSG 22個，其中尤以green和red模塊最多，分別為10和4個。同時通過DAVID進行模塊內基因功能富集發現green和red模塊所含基因與腎小球發育、細胞外基質解聚、細胞機械刺激等足細胞常見的病理狀態密切相關，因此筆者認為這兩個模塊與足細胞損傷關系較為密切。另一方面，與兩個模塊比較，發現green模塊中包含更多的coDEG，因此可以推測green模塊中的基因可能參與多種CKD發展中 的足細胞損傷。

圖1 各數據集差異基因表達熱圖及韋恩圖

圖2 WGCNA分析構建基因共表達模塊

2.4 足細胞損傷相關模塊加權網絡分析

hub基因是指在模塊中連接度最高的一系列基因，與全局網絡中的樞紐基因比較，模塊中的樞紐基因往往更有生物學意義。本文定義的WGCNAhub基因為模塊內平均連接度（Kwithin）前5%的基因[2]。結果顯示，green模塊中MAGI2基因可以同時作為即是green模塊的hub基因也是本研究中的coDEG和PSG（見表2），結果高度提示MAGI2不僅具有CKD足細胞損傷的分子標志物潛在價值，也可能在早期足細胞損傷中起到重要生物學作用。選擇關聯權重閾值≥0.05構建mRNA共表達網絡（見圖4A），并使用Cytoscape及ClueGO插件并聯合String數據庫構建green模塊內所有基因的KEGG/GO/PPI互作網絡，進一步挖掘模塊內基因參與的生物學進程及蛋白相互作用，結果顯示與MAGI2表達或蛋白互作高度相關的多種基因已經被報道參與足細胞損傷，包括NPHS2、FGF1、BMP7等（見圖4B）。

表2 構建出的模塊資料及所含基因總結

圖3 各模塊功能富集結果（取P值最小前3項）

圖4 green模塊內基因功能及相互作用網絡可視化

3 討論

CKD時足細胞損害機制極其復雜，涉及多種通路及細胞因子[3]。但傳統的生物學研究以單基因或蛋白為出發點，因此其僅能對生物系統的局部作出解釋，難以對系統的整體進行全面的探索。相比之下，WGCNA是一種從高通量數據中挖掘模塊（module）信息的算法。在本方法中module被定義為一組具有類似表達譜的基因，如果某些基因在一個生理過程或不同組織中總是具有相類似的表達變化，那么有理由認為該基因在功能上是相關的，可以把其定義為一個模塊（module）。這似乎有點類似于進行聚類分析所得結果，但不同的是，WGCNA的聚類準則具有生物學意義，而非常規的聚類方法（如利用數據間的幾何距離），因此該方法所得結果具有更高的可信度，更有利于從生物功能整體入手考慮基因功能以及內在關聯[4]。

本研究所納入的GSE66107數據集是PAN誘導的足細胞損傷的經典體外模型，PAN是嘌呤霉素的衍生物，具有選擇性的腎損害作用，PAN所致的腎臟損傷與人類腎病損傷一致，早期表現為足細胞足突融合和消失的微小病變（MCD），隨著PAN作用時間的延長和蓄積量的增加，足細胞凋亡增加，進一步發展為局灶節段性腎小球硬化（FSGS），而PAN作用的靶細胞即為足細胞。因此通過足細胞體外損傷模型的基因表達變化與另3個CKD患者腎小球基因變化取交集，有利于排除腎臟其他實質細胞（如系膜細胞和內皮細胞）表達變化的干擾，可以更明確探討CKD時足細胞內可能存在的基因表達異常。最終將4個數據集所得到的DEG結合含有133個CKD樣本表達譜數據進行WGCNA分析，并通過GO富集分析現red和green模塊在GO分類上與多種足細胞損傷相關進程較為密切：如機械刺激反應（cell response to mechanical stimulus）[5]、內質網應激反應（response to endoplasmic reticulum stress）[6-7]、腎小球發育（glomerulus development）、細胞外基質分解[8]（extracellular matrix disassembly）等。更有價值的是，筆者發現green模塊中的hub基因MAGI2也同為本研究中的coDEG和PSG，同時，網絡可視化結果顯示MAGI2與多種文獻已報道的在足細胞損傷或腎臟病中起到重要“著名”基因相關量，如FGF1[9]、BMP7等[10]。這意味著MAGI2可能具有作為足細胞損傷早期分子標志物的潛在價值，并且在足細胞損傷中起到重要的調控作用，同時通過共表達網絡的構建識別其高度相關的節點基因，對后續MAGI2的調控機制研究提供重要線索。

本研究不足之處在于從GEO數據庫中下載的基因表達譜無法得到相應的臨床資料，如腎小球濾過率（eGFR）、血肌酐（Scr）、尿蛋白等，如果可以加入該臨床參數作為參考，那么WGCNA分析時更可以將所得模塊與該臨床信息進行相關性分析，挖掘模塊可能所具有的生物學意義，發揮WGCNA的更大作用。

綜上所述，本研究通過運用一系列合理的生物信息學手段，分析CKD患者腎小球內與足細胞損傷相關度的15個差異變化基因，同時通過共表達分析構建出足細胞損傷相關基因模塊green，研究結果對于CKD，尤其是足細胞損傷的早期診斷提供線索，同時為足細胞保護相關研究提供依據及參考，是WGCNA在CKD足細胞損傷中的一次新型嘗試。