慢性腎臟病中Cyr61加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析

(青島大學(xué)附屬醫(yī)院腎內(nèi)科,山東 青島 266003)

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

1 材料與方法

1.1 基因芯片數(shù)據(jù)的獲取

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

1.2 方法

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

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

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

1.3 統(tǒng)計(jì)學(xué)分析

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

2 結(jié) 果

2.1 Cyr61基因在CKD中的表達(dá)

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

表1 GEO數(shù)據(jù)庫中實(shí)驗(yàn)編號(hào)及樣本的數(shù)量

圖1 Cyr61基因在CKD腎小球中的相對(duì)表達(dá)

圖2 Cyr61基因在CKD腎間質(zhì)中的相對(duì)表達(dá)

2.2 Cyr61基因與CKD的關(guān)系

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

表3 各模塊中基因的數(shù)量

上圖紅線顯示相關(guān)系數(shù)等于0.85，β=8最接近紅線，下圖顯示β=8時(shí)相關(guān)系數(shù)等于0.87。

圖3加權(quán)系數(shù)β以及無尺度網(wǎng)絡(luò)特性檢驗(yàn)

圖4 基因?qū)哟尉垲悩錉顖D及對(duì)應(yīng)模塊

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

基因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個(gè)基因。

2.3 Cyr61與CKD發(fā)生過程中相關(guān)蛋白的關(guān)系

圖6a、b展示了q值最小的30個(gè)基因之間的相互關(guān)系，Cyr61作為樞紐基因在CKD中與眾多樞紐基因之間有相互作用。使用Cytoscape按拓?fù)渲丿B權(quán)重構(gòu)建Cyr61與其他樞紐基因網(wǎng)狀圖，圖7展示出Cyr61在與CKD發(fā)生相關(guān)性最高的Brown模塊中拓?fù)渲丿B最高的15個(gè)基因，其中MAFF、JUN、KLF6拓?fù)渲丿B系數(shù)大于0.10，ATF3、JUNB拓?fù)渲丿B系數(shù)大于0.09，EGR1、MYC、GDF15、HBEGF、GEM、CCL2、RHOB、TRIB1的拓?fù)渲丿B系數(shù)大于0.06。本文研究將這15個(gè)基因與STRING PPI網(wǎng)絡(luò)進(jìn)行對(duì)比，其中MAFF、KLF6、GDF15、HBEGF、GEM、CCL2、RHOB以及TRIB1沒有與Cyr61基因相關(guān)實(shí)驗(yàn)。

a、b圖分別基于拓?fù)渲丿B與相關(guān)系數(shù)，其中每一行及一列分別代表一個(gè)基因，顏色越深表示兩個(gè)基因之間的拓?fù)渲丿B或者相關(guān)系數(shù)越高，兩者之間的相關(guān)度越高。

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

基因之間的距離越近，提示兩基因之間的相關(guān)度越高。

3 討 論

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

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

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

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

[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對(duì)缺氧狀態(tài)下腎小管上皮細(xì)胞的保護(hù)作用[J]. 中華腎臟病雜志, 2015,31(6):451-455.

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

[27] 徐巖,沈?qū)W飛,宋年華. 富半胱氨酸蛋白61對(duì)轉(zhuǎn)化生長因子β1誘導(dǎo)腎小管上皮細(xì)胞凋亡的影響[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.