慢性腎臟病中Cyr61加權基因共表達網絡分析

(青島大學附屬醫院腎內科,山東 青島 266003)

目前,世界范圍內慢性腎臟病(Chroic Kidney Disease，CKD)的患病率逐年增高,已經成為全球重要的公共衛生問題之一，美國CKD的患病率約為13%[1-2],而在我國CKD的患病率約為10.8%[3]。CKD會持續進展為終末期腎病，并增加病人心血管事件及死亡風險。對CKD的發病機制進行更深入的研究，將有助于對CKD進行預防和治療。富含半胱氨酸蛋白61(Cysteine rich angiogenic inducer 61，Cyr61)是一種具有肝素結合活性的分泌蛋白，其作為基質相關信號分子參與細胞增殖、分化、轉化、凋亡及應激等生物過程。研究顯示，缺血再灌注腎損傷早期，Cyr61蛋白在組織及尿液中均明顯升高[4]。我們的前期研究顯示，在缺血性急性腎損傷(Acute kidney injury，AKI)早期，缺氧狀態下過表達Cyr61能夠促進腎小管上皮細胞增殖，抑制細胞凋亡。本文采用生物信息學數據分析方法，利用網絡公共數據庫基因芯片篩選CKD表達差異基因，構建基因共表達網絡，探討Cyr61基因在CKD腎小球及腎間質中的表達情況，并探討Cyr61基因是否為CKD樞紐基因。

1 材料與方法

1.1 基因芯片數據的獲取

對公共的基因芯片數據庫NCBI GEO[5]和Array Express[6]進行CKD相關樣本檢索，篩選條件：①物種設定人類；②轉錄組表達譜芯片；③包含正常腎臟對照組；④不重復的數據。得到GEO數據庫中GSE47185[7]、GSE37463[8]、GSE35489[9]、GSE3-2592[8]共4個數據集373個樣本(平臺GPL11670、平臺GLP14663)，樣本的臨床信息已被貢獻者上傳至nephroseq(www.nephroseq.org, 11 2016, University of Michigan, Ann Arbor, MI)。GEO提供原始數據文件(TAR of CEL)和經過預處理的矩陣文件(TXT)，由于不同實驗對數據的預處理方法不一致，本研究下載并使用原始數據及平臺注釋進行統一的數據處理。

1.2 方法

1.2.1數據預處理與差異基因篩選 Affy芯片預處理一般分為背景處理、標準化處理及總匯3步。本文應用R軟件(v.3.3.1)對所有樣本數據進行數據分析，使用Affy包(v.1.50.0)[10]的RMA(Robust Multiarray Average)方法進行三合一預處理與表達基因篩選，根據疾病類型將373個樣本分為7個疾病組和正常組(living kidney donors, LD)，并按照樣本的來源分為腎小球組與腎間質組，使用線性回歸模型limma包[11]對各組進行表達差異計算，設定校正后P值(FDR)與對數化表達變化倍數(log2 fold change，logFC)作為篩選差異基因閾值。整理各疾病組基因表達數據，從而得到Cyr61基因的表達情況。

1.2.2加權基因共表達網絡分析(Weighted Gene Co-expression Network Analysis，WGCNA) 本實驗應用WGCNA對CKD進行基因共表達分析。WGCNA是構建基因共表達網絡的典型系統生物學算法，WGCNA首先將在不同的樣本中存在高度相關性的基因劃分在同一個共表達網絡，或稱之為屬于同一模塊，隨后將相關模塊與外界信息(臨床數據、單核苷酸多態性、蛋白組學等)進行關聯，獲得模塊與各外界信息的相關程度，最后在與疾病相關程度高的模塊中尋找樞紐基因、預測調控通路等。參考WGCNA使用手冊[12]，將表達差異基因變化倍數對數化，使用R語言中WGCNA包(v.1.51)進行WGCNA。

1.2.3樞紐基因的篩選及Cyr61相關基因分析基因或蛋白在疾病發生過程中不是單獨發揮作用而是共同發揮作用的，有多種方法可以用來描述基因與基因之間的相關性，WCGNA使用拓撲重疊法(Topological overlap measure，TOM)來計算基因之間的相關程度，使用TOM計算的結果更加具有生物學意義[13-14]。本研究使用WGCNA包提供的“networkScreening”函數在與慢性腎臟疾病發生相關性最高的模塊中尋找樞紐基因，獲得基因與疾病相關性加權P值、矯正后的q值[15]、加權后相關系數cor及FishZ檢驗結果，其中q數值越小說明該基因與CKD的發生相關性越強。使用“exportNetworkToCytoscape”函數獲得Cyr61基因與其他樞紐基因的相關性數據，將數據導入Cytoscape(v.3.4.0)[16]，并按照拓撲重疊權重構圖，得到與Cyr61基因高度相關的樞紐基因。蛋白互作網絡分析使用STRING數據庫(http://www.string-db.org)[17]。

1.3 統計學分析

采用Limma軟件包以moderated t-statistic[18]進行兩組間比較；WGCNA采用Pearson法或者TOM法。以P<0.05為差異有統計學意義。

2 結 果

2.1 Cyr61基因在CKD中的表達

7種CKD及正常對照共373個樣本原始數據均來源于NCBI GEO(表1)。為了研究Cyr61基因在腎臟疾病不同部位的表達情況，本研究將樣本分為腎小球組與腎間質組，并使用R軟件獲取未經對數化處理的基因表達數據，結果顯示Cyr61基因在CKD腎小球中表達普遍降低(圖1)，而在腎間質中部分表達降低(圖2，表2)。應用線性回歸模型limma包對各組進行表達差異計算，獲取Cyr61基因在各組中的變化倍數，結果顯示Cyr61基因在高血壓腎病(hypertensive nephropathy, HTN)、微小病變性腎病(minimal change nephropathy, MCD)、膜性腎病(membranous nephropathy, MGN)、狼瘡性腎炎(Lupus nephritis，LN)、局灶階段性腎小球硬化(lupus nephritis and focal segmental glomerulosclerosis, FSGS)的腎小球及間質中降低,差異具有統計學意義(FDR<0.05，logFC<-0.8)，在IgA腎病(IgA nephropathy, IgAN)腎間質中降低，差異具有統計學意義(FDR<0.01，logFC<-1)。

表1 GEO數據庫中實驗編號及樣本的數量

圖1 Cyr61基因在CKD腎小球中的相對表達

圖2 Cyr61基因在CKD腎間質中的相對表達

2.2 Cyr61基因與CKD的關系

由于Cyr61基因在各慢性腎臟疾病間質中普遍表達降低，本研究對GPL14663平臺上各腎間質組的2 134個差異基因進行基因共表達分析。按照無尺度網絡的標準，本研究選擇相關系數等于0.87的鄰接矩陣權重參數(軟閾值)β=8構建基因模塊(圖3),并根據動態混合剪切法，共得到9個基因模塊(圖4，表3)，通過計算模塊內各基因表達量與樣本特征向量的皮爾森相關系數，尋找與CKD發生顯著相關的基因模塊，其中Brown模塊基因的顯著性(Gene significance，GS)高于其他模塊(圖5)，因此該模塊與CKD相關性最高。隨后，應用WGCNA的“networkScreening”函數判斷樞紐基因，研究結果顯示Cyr61(q<0.01)是CKD發生的樞紐基因(表4)。

表3 各模塊中基因的數量

上圖紅線顯示相關系數等于0.85，β=8最接近紅線，下圖顯示β=8時相關系數等于0.87。

圖3加權系數β以及無尺度網絡特性檢驗

圖4 基因層次聚類樹狀圖及對應模塊

柱狀圖高度表示整個模塊基因顯著性的平均值。

基因q值corZ值基因q值corZ值Cyr610-0．780-14．494DUSP10-0．908-25．173JUN0-0．819-16．613BHLHE400-0．704-11．504JUNB0-0．894-23．235BTG20-0．789-14．921NFKBIA0-0．774-14．194IER20-0．750-13．162ZFP360-0．926-28．403TRIB10-0．785-14．736IER30-0．793-15．131CDKN1A0-0．656-10．098EGR10-0．729-12．356NR4A10-0．867-20．206SGK10-0．871-20．599RGS20-0．717-11．964PPP1R15A0-0．755-13．373KLF100-0．827-17．072LDLR0-0．743-12．906SERPINE10-0．808-15．942

注：僅展示q最小的前20個基因。

2.3 Cyr61與CKD發生過程中相關蛋白的關系

圖6a、b展示了q值最小的30個基因之間的相互關系，Cyr61作為樞紐基因在CKD中與眾多樞紐基因之間有相互作用。使用Cytoscape按拓撲重疊權重構建Cyr61與其他樞紐基因網狀圖，圖7展示出Cyr61在與CKD發生相關性最高的Brown模塊中拓撲重疊最高的15個基因，其中MAFF、JUN、KLF6拓撲重疊系數大于0.10，ATF3、JUNB拓撲重疊系數大于0.09，EGR1、MYC、GDF15、HBEGF、GEM、CCL2、RHOB、TRIB1的拓撲重疊系數大于0.06。本文研究將這15個基因與STRING PPI網絡進行對比，其中MAFF、KLF6、GDF15、HBEGF、GEM、CCL2、RHOB以及TRIB1沒有與Cyr61基因相關實驗。

a、b圖分別基于拓撲重疊與相關系數，其中每一行及一列分別代表一個基因，顏色越深表示兩個基因之間的拓撲重疊或者相關系數越高，兩者之間的相關度越高。

圖6熱圖展示CKD中顯著性最高的前30個基因

基因之間的距離越近，提示兩基因之間的相關度越高。

3 討 論

已有研究表明，Cyr61在許多疾病及損傷中表達發生變化[19-20]，特別是在急性損傷中，Cyr61可以作為早期生物標志物反應疾病的損傷程度和活動性[21-23]。在腎臟疾病的研究中，Cyr61與缺血再灌注腎損傷[4]和腎小球腎炎[24]存在一定的關系。免疫組化結果顯示，Cyr61在IgAN和微小病變型腎病腎小球組織中表達降低[24]。本研究結果與其一致，但免疫組化結果顯示，Cyr61在各腎臟病腎間質中的表達穩定。而本研究同時顯示，在基因芯片上Cyr61基因在各種慢性腎臟疾病間質中的表達是顯著降低的，這可能是由于基因的豐富度不一定與翻譯產物蛋白質呈線性關系有關。

我們的前期研究顯示，Cyr61可抑制轉化生長因子-β(transforming growth factor-β,TGF-β)受體表達，對TGF-β誘導的腎小管上皮細胞凋亡有拮抗作用，同時Cyr61還可以參與調解細胞周期蛋白表達，促進細胞增殖[25-27]，并且在缺氧條件下，Cyr61表達可以保護腎小管上皮細胞免于凋亡[28]。最新的研究表明，在急性梗阻性腎纖維化模型中，Cyr61參與TGF-β介導的炎癥反應，拮抗Cyr61表達可暫時改善間質的炎癥反應[29]。以上研究結果證明，Cyr61在腎臟疾病發病過程中發揮一定的作用。在急性腎損傷、CKD中，不同的病理損害以及不同的腎組織中Cyr61表達有差異性，此為預防、診斷和治療腎臟疾病提供了新的研究方向。

本研究應用生物信息學的方法，整合公共數據庫373個樣本，從Cyr61的角度分析其在各種CKD中的表達情況，并首次證明Cyr61在HTN、MCD、MGN、LN、FSGS的腎小球及間質中表達降低，在IgAN只有腎間質中表達顯著降低。在此基礎上，本研究對CKD腎間質的2 134個差異基因進行基因共表達分析，結果顯示Cyr61所處的基因模塊與CKD的發生相關性最高，在此模塊中Cyr61與CKD相關性加權q值<0.01，從而證明Cyr61屬于CKD的樞紐基因之一，這說明Cyr61在CKD的發生過程中可能起到重要的作用。而在CKD腎小球中，Cyr61與CKD的發生相關性卻較低，這也提示Cyr61在CKD腎間質中參與的功能比在腎小球中多。本實驗進一步對Cyr61共表達基因進行研究分析，結果顯示MAFF、JUN[30]、KLF6、ATF3[31-32]、JUNB與Cyr61呈高度相關性，目前尚沒有相關實驗證實MAFF、KLF6、JUNB與Cyr61的相關性，這為進一步研究Cyr61在CKD中的作用提供了理論基礎。

目前生物信息學分析廣泛應用于醫學領域，為疾病的預防、診斷、治療及預后判斷提供了新的方向，使人們對于疾病發生機制有了更深的理解，生物信息學在醫學中的應用也是“精準醫療”的要求。隨著計算機硬件和軟件的不斷發展，算法的不斷改進，生物信息學將會在醫學領域承擔更加重要的角色。

[1] COLLINS A J, FOLEY R N, CHAVERS B, et al. United States Renal Data System 2011 Annual Data Report: Atlas of chronic kidney disease & end-stage renal disease in the United States[J]. Am J Kidney Dis, 2012,59(1 Suppl 1):A7,e1-420.

[2] COVIC A, KOTHAWALA P, BERNAL M, et al. Systematic review of the evidence underlying the association between mine-ral metabolism disturbances and risk of all-cause mortality, cardiovascular mortality and cardiovascular events in chronic kidney disease[J]. Nephrol Dial Transplant, 2009,24(5):1506-1523.

[3] ZHANG L, WANG F, WANG L, et al. Prevalence of chronic kidney disease in China: A cross-sectional survey[J]. Lancet, 2012,379(9818):815-822.

[4] MURAMATSU Y, TSUJIE M, KOHDA Y, et al. Early detection of cysteine rich protein 61 (CYR61, CCN1) in urine following renal ischemic reperfusion injury[J]. Kidney Int, 2002,62(5):1601-1610.

[5] BARRETT T, WILHITE S E, LEDOUX P, et al. NCBI GEO: Archive for functional genomics data sets-update[J]. Nucleic Acids Res, 2013,41(Database issue):D991-995.

[6] KOLESNIKOV N, HASTINGS E, KEAYS M, et al. ArrayExpress update-simplifying data submissions[J]. Nucleic Acids Res, 2015,43(Database issue):D1113-1116.

[7] JU W, GREENE C S, EICHINGER F, et al. Defining cell-type specificity at the transcriptional level in human disease[J]. Genome Res, 2013,23(11):1862-1873.

[8] BERTHIER C C, BETHUNAICKAN R, GONZALEZ-RIVERA T, et al. Cross-species transcriptional network analysis defines shared inflammatory responses in murine and human lupus nephritis[J]. J Immunol, 2012,189(2):988-1001.

[9] REICH H N, TRITCHLER D, CATTRAN D C, et al. A molecular signature of proteinuria in glomerulonephritis[J]. PLoS One, 2010,5(10):e13451.

[10] IRIZARRY R A, HOBBS B, COLLIN F, et al. Exploration, normalization, and summaries of high density oligonucleotide array probe level data[J]. Biostatistics, 2003,4(2):249-264.

[11] RITCHIE M E, PHIPSON B, WU D, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies[J]. Nucleic Acids Res, 2015,43(7):e47.

[12] LANGFELDER P, HORVATH S. WGCNA: An R package for weighted correlation network analysis[J]. BMC Bioinformatics, 2008,9:559.

[13] LI A, HORVATH S. Network neighborhood analysis with the multi-node topological overlap measure[J]. Bioinformatics, 2007,23(2):222-231.

[14] YIP A M, HORVATH S. Gene network interconnectedness and the generalized topological overlap measure[J]. BMC Bioinformatics, 2007,8:22.

[15] STOREY J D, TIBSHIRANI R. Statistical significance for genomewide studies[J]. Proc Natl Acad Sci U S A, 2003,100(16):9440-9445.

[16] SHANNON P, MARKIEL A, OZIER O, et al. Cytoscape: A software environment for integrated models of biomolecular interaction networks[J]. Genome Res, 2003,13(11):2498-2504.

[17] SZKLARCZYK D, FRANCESCHINI A, WYDER S, et al. STRING v10: Protein-protein interaction networks, integrated over the tree of life[J]. Nucleic Acids Res, 2015,43(Database issue):D447-452.

[18] CUI X, CHURCHILL G A. Statistical tests for differential expression in cDNA microarray experiments[J]. Genome Biol, 2003,4(4):210.

[19] GAO L, FAN Y, HAO Y, et al. Cysteine-rich 61 (Cyr61) upregulated in pulmonary arterial hypertension promotes the proliferation of pulmonary artery smooth muscle cells[J]. Int J Med Sci, 2017,14(9):820-828.

[20] REID S E, KAY E J, NEILSON L J, et al. Tumor matrix stiffness promotes metastatic cancer cell interaction with the endothelium[J]. EMBO J, 2017,36(16):2373-2389.

[21] KLINGENBERG R, AGHLMANDI S, LIEBETRAU C, et al. Cysteine-rich angiogenic inducer 61 (Cyr61): A novel soluble biomarker of acute myocardial injury improves risk stratification after acute coronary syndromes[J]. Eur Heart J, 2017,38(47):3493.

[22] WOO Y J, SEO Y, KIM J J, et al. Serum CYR61 is associa-ted with disease activity in graves’ orbitopathy[J]. Ocul Immunol Inflamm, 2017:1-7.

[23] LIN J, LI N, CHEN H, et al. Serum Cyr61 is associated with clinical disease activity and inflammation in patients with systemic lupus erythematosus[J]. Medicine (Baltimore), 2015,94(19):e834.

[24] SAWAI K, MUKOYAMA M, MORI K, et al. Expression of CCN1 (CYR61) in developing, normal, and diseased human kidney[J]. Am J Physiol Renal Physiol, 2007,293(4):F1363-1372.

[25] 徐巖,國敏,董輝. 含胱氨酸蛋白61對缺氧狀態下腎小管上皮細胞的保護作用[J]. 中華腎臟病雜志, 2015,31(6):451-455.

[26] 徐巖,沈學飛,宋年華. 富含半胱氨酸蛋白61對人腎小管上皮細胞增殖與細胞周期的影響[J]. 中華腎臟病雜志, 2013,29(4):273-276.

[27] 徐巖,沈學飛,宋年華. 富半胱氨酸蛋白61對轉化生長因子β1誘導腎小管上皮細胞凋亡的影響[J]. 中華腎臟病雜志, 2012,28(10):813-814.

[28] MA R, ZHANG J, LIU X, et al. 7,8-DHF Treatment induces Cyr61 expression to suppress hypoxia induced ER stress in HK-2 cells[J]. Biomed Res Int, 2016, 2016:5029797.

[29] LAI C F, CHEN Y M, CHIANG W C, et al. Cysteine-rich protein 61 plays a proinflammatory role in obstructive kidney fibrosis[J]. PLoS One, 2013,8(2):e56481.

[30] KEENAN S W, HILL C A, KANDOTH C, et al. Transcriptomic responses of the heart and brain to anoxia in the western painted turtle[J]. PLoS One, 2015,10(7):e0131669.

[31] KIM I, LEE S H, JEONG J, et al. Functional profiling of human MeCP2 by automated data comparison analysis and computerized expression pathway modeling[J]. Healthc Inform Res, 2016,22(2):120-128.

[32] JONGBLOETS B C, VAN GASSEN K L, KAN A A, et al. Expression profiling after prolonged experimental febrile seizures in mice suggests structural remodeling in the hippocampus[J]. PLoS One, 2015,10(12):e0145247.