GWGS芯片整合分析識別外周血細胞中與骨質疏松癥相關的樞紐基因

黎海霞,許 陽,鄧凱麗,李 然,肖蘇妹

(中山大學公共衛生學院,中國廣東廣州510030)

骨質疏松癥是一種由于骨密度(bone mineral density,BMD)降低和骨微結構退化,引起骨強度降低,使得骨折風險增加的常見疾病[1]。目前世界上約2億人患有該疾病,其中80%是年齡大于60歲的女性[2]。在中國,骨質疏松癥在過去10年中的發病率顯著升高,從2008年之前的14.9%增加至2012-2015年間的27.9%,而且相同年齡組中女性的發病率顯著高于男性[3～4]。

骨質疏松癥是一種由遺傳因素、環境因素以及它們的相互作用共同決定的復雜疾病。遺傳因素對骨質疏松癥的發生發展有著重要影響。BMD是骨質疏松癥臨床診斷的重要指標,也是最常用的研究表型。研究證實BMD具有較高的遺傳率(0.5～0.85)[5]。低BMD的病理機制包括破骨細胞的壽命延長和成骨細胞、骨細胞的過早凋亡[6]。外周血單核細胞(peripheral blood monocytes,PBMs)和B淋巴細胞是骨質疏松癥相關研究的重要細胞類型,均參與形成破骨細胞。體外研究顯示,PBMs在合適的微環境中可以分化為有功能的破骨細胞[7～8]；B 淋巴細胞能通過 RANKL(receptor activator of NF-κB ligand)通路促進破骨細胞的形成[9]。PBMs及B淋巴細胞與破骨細胞關系密切,通過參與破骨細胞生成,最終影響骨吸收和骨形成之間的平衡。

現在已有一些基于表達譜芯片鑒定骨質疏松癥相關基因的研究報道[10～12]。但是這些研究的結論并不一致,原因可能是單個研究樣本量太少使得統計學功效不足,而對芯片數據展開整合分析能增加研究樣本量,提高結果的穩健性。由于不同研究間的P值差異可能較大,其中P值最小的研究會很大程度上決定整合P值的結果,所以傳統的基于P值的meta分析對于整合小樣本量的研究并不合適。基于基因組的總體顯著性(genomewide global significance,GWGS)的整合分析方法則是首先將不同數據集差異表達基因(differentially expressed genes,DEGs)的表達倍數絕對值由大到小進行排序,再根據整合后的序號對基因進行排序篩選,這種方式消除了因不同研究的P值差異很大對結果產生的影響,因此更適合用于整合小樣本量的研究[13]。本研究應用GWGS方法來整合分析公共數據庫中可獲得的PBMs及B淋巴細胞的基因表達譜芯片數據,以鑒定外周血細胞中與骨質疏松癥相關的樞紐基因。

1 材料與方法

1.1 基因表達數據集的篩選

從 GEO(Gene Expression Omnibus,http：//www.ncbi.nlm.nih.gov/geowebcite)和 ArrayExpress(http：//www.ebi.ac.uk/arrayexpress/)數據庫中尋找骨質疏松癥相關的表達譜芯片數據。分別用關鍵詞“骨質疏松癥”“骨密度”“BMD”進行檢索,然后根據以下納入排除標準進行篩選。納入標準：1)研究對象為成年女性；2)研究類型為病例-對照研究,且樣本量大于10；3)研究樣本來自全血組織。排除標準：1)研究對象為繼發性骨質疏松癥患者,或患有影響骨代謝的疾病；2)公共數據庫上無法獲得基因表達譜芯片原始數據。

1.2 數據預處理及差異表達基因的獲得

下載納入研究的原始文件后,用R語言(v3.3.3)中的Bioconductor包進行數據處理。原始數據用Affy包進行讀取[14],用RMA(Robust Multichip Average)進行背景校正[15],并在R中用層次聚類和主成分分析進行質量控制。探針編號轉為對應的基因后,如果多個探針對應同一個基因,則用這些探針的平均值代表這個基因的表達值。納入的研究用limma包識別DEGs[16]。

1.3 GWGS整合分析

1.4 基因功能注釋和蛋白質相互作用網絡的構建

對篩選出的200個基因通過DAVID數據庫(Database for Annotation,Visualization and Integrated Discovery,v6.8,https：//david.ncifcrf.gov/)進行GO(gene ontology)富集分析和KEGG(kyoto encyclopedia of genes and genomes)通路分析。用STRING(Search Tool for the Retrieval of Interacting Genes/Proteins,v11,https：//string-db.org/)進行蛋白質相互作用(protein-protein interaction,PPI)網絡的構建[17],并用Cytoscape進行可視化處理[18]。PPI中的節點(note)代表相應的蛋白質,邊(edge)代表兩個蛋白質之間存在相互作用,度(degree)代表某蛋白質與網絡中其他蛋白質存在直接相互作用的個數。度的數值越大代表該基因越可能是骨質疏松癥相關的樞紐基因(hub gene)。

2 結果

2.1 納入研究的基本情況

根據篩選流程和納入排除標準,本研究共納入3個涉及骨質疏松癥研究的表達譜基因芯片數據集(GSE56815、GSE7429 和 GSE7158)。具體信息如表1所示,其中2個研究來自PBMs,相應的芯片平臺是Affymetrix Human Genome U133A Array；1個研究來自B淋巴細胞,相應的平臺是Affymetrix Human Genome U133 Plus 2.0(Affymetrix,Santa Clara,CA,USA)。納入的3個研究的研究對象均為成年女性,年齡范圍為20～60歲。經過質量控制,把層次聚類中與總體距離最遠以及在主成分分析結果中明顯偏離總體的樣本剔除,最終的高低BMD樣本數分別為60例和59例。根據世界衛生組織對BMD的定義,這60例高BMD樣本可歸類為骨量正常人群,59例低BMD樣本可歸類為骨質疏松癥或骨量減少人群[19]。

2.2 GWGS整合分析

進行數據預處理后,共獲得12 420個共同基因。針對每一個數據集,采用limma包篩選出DEGs,剔除重復DEGs后,共獲得4 069個DEGs。基于GWGS方法計算獲得DEGs的整合GWGS值,最終選出整合GWGS值最大的200個基因。表2列出了GWGS值排序為前20的基因的簡稱、全稱及其在染色體中的位置。

2.3 功能注釋和PPI分析

對GWGS值排序前200的DEGs進行GO富集分析和KEGG通路分析。GO分析結果顯示,富集最顯著的前3個條目分別為脂多糖的細胞反應(cellular response to lipopolysaccharide,GO：0071-222,P=3.30E-05)、凋亡過程(apoptotic process,GO：0006915,P=1.60E-04)和炎癥反應(inflammatory response,GO：0006954,P=2.69E-04)。表 3展示了排序前10的GO富集條目的具體信息。KEGG通路分析顯示,篩選出的基因富集于骨質疏松癥相關信號通路。其中,最顯著的通路為破骨細胞分化通路(osteoclast differentiation,hsa04380,P=2.54E-03)。表4展示了排序前10的KEGG富集通路條目的具體信息。由200個DEGs構建的PPI網絡結果如圖1所示,該網絡由196個節點和511條邊組成。圖中深色的節點代表樞紐基因,節點越大,代表連接度越大。結果顯示其中10個基因的連接度大于20,提示為與骨質疏松癥相關的樞紐基因(表5)。樞紐基因中具有最大度值的基因是位于14號染色體上的AKT1基因(度=44)。

表1 納入研究的基本特征(n=3)Table 1 Basic characteristics of the included study(n=3)

表2 GWGS值最大的前20個基因Table 2 The top 20 genes with the highest GWGS values

3 討論

本研究應用GWGS方法整合分析了3個骨質疏松癥相關的表達譜芯片數據,對獲得的4 069個DEGs展開分析,選出GWGS值最大的前200個基因,通過構建PPI網絡,篩選出10個外周血細胞中與骨質疏松癥相關的樞紐基因。

對于篩選到的這10個樞紐基因,既往研究已發現其中9個基因參與骨代謝過程。AKT1基因是度最大的樞紐基因,影響著破骨細胞和成骨細胞的功能。小鼠模型的研究顯示Akt1可以控制破骨細胞和成骨細胞的分化[20]。在小鼠研究中,降低Akt1的表達可同時影響破骨細胞和成骨細胞的功能,使得骨轉換降低,從而引起骨量減少[21]。另外有研究表明,Akt1缺乏導致小鼠股骨BMD、股骨皮質厚度和體積以及小梁厚度的降低[22]。以上證據顯示,AKT1可能在骨質疏松癥中起重要作用。但是,現階段針對AKT1在人體內的功能研究相對較少。因此,需要進一步的研究來證實AKT1在人類骨質疏松癥中的潛在作用。STAT1基因在骨代謝中扮演著重要的角色。研究發現在骨質疏松癥的小鼠模型中,Stat1的表達顯著上調[23～24]。Kwak等[25]發現,RANKL可通過p38 MAPK通路促進Stat1編碼的蛋白質的磷酸化,影響PBMs的遷移和黏附。小鼠研究發現,Cxcl10基因在T細胞和巨噬細胞向炎癥關節浸潤并導致骨破壞中起關鍵作用,它能部分逆轉PRMT5(protein arginine methyltransferase 5)抑制劑對破骨細胞生成的抑制作用[26]。Lee等[27]的研究發現CXCL10能促進破骨細胞的分化,在成骨細胞和骨髓細胞的共培養中CXCL10以劑量依賴的方式誘導破骨細胞的形成。Arg1基因在破骨細胞分化過程中下調,可能通過調節一氧化氮的產生,負性調節破骨細胞分化[28]。Ptpn6(也叫作Shp-1)基因參與調節破骨細胞的存活,抑制Ptpn6可上調破骨細胞存活的基礎水平,降低活性氧對破骨細胞存活的影響[29]。Zhang等[30]的研究發現Ptpn6通過向含有Traf6(TNF receptor associated factor 6)的復合物聚集來調節RANKL誘導的破骨細胞生成。Chodynicki等[31]的研究發現,CTSD基因能激活破骨細胞并在膽脂瘤引起的骨組織破壞中起主要作用。CAMP基因編碼的蛋白質是一種抗菌肽,具有多種生物學功能和促進成骨的潛力。有研究表明Camp通過P2X7R(purinergic ligand-gated ion channel 7 receptor)和MAPK(mitogen-activated protein kinase)信號途徑抑制炎癥及破骨細胞骨吸收功能[32]。EGR1基因的表達下降可導致前列腺癌病人的骨轉移,溶骨性骨損傷面積減少,以及骨和腫瘤界面的破骨細胞數量減少[33]。LCN2基因可抑制破骨細胞的增殖和分化,導致破骨細胞的形成受損。在骨髓源性巨噬細胞中,如果Lcn2過表達或加入重組Lcn2蛋白,會抑制多核破骨細胞的形成[34]。

表3 GO分析排序前10的條目Table 3 Top 10 items of GO analysis

表4 KEGG富集通路排序前10的條目Table 4 Top 10 items of KEGG enrichment pathway

圖1 PPI網絡分析結果節點越大,代表度越大；深色節點代表樞紐基因,淺色節點代表其他基因。Fig.1 PPI network analysis resultsA larger node means a greater degree；the dark nodes are the hub genes,and the light nodes are the other genes.

表5 樞紐基因基本信息(n=10)Table 5 Basic information of the hub genes(n=10)

此外,有1個樞紐基因ELANE(neutrophil elastase),目前并未有研究報道其參與骨質疏松癥的發生發展。基于79個人類組織和細胞類型的芯片數據(Genealas U133A,http：//biogps.org/)的研究結果顯示,ELANE基因在人體的骨髓組織中高表達[35]。基于96個小鼠組織和細胞類型的芯片數據(Genealas MOE430)的研究結果顯示,Elane基因同時在小鼠骨髓和骨中有高表達[36]。這些研究結果提示,ELANE基因很可能與骨質疏松癥有潛在的聯系,對其可以進一步展開相關研究。

當然,本研究也存在一些局限性。首先,研究組織的樣本來源不完全一致,雖然3個研究的樣本均來自外周血細胞,但有2個為PBMs,1個為B淋巴細胞；其次,數據集研究對象均為女性,而且研究對象年齡跨度較大,從20歲到60歲。此外,本研究也僅進行了表達譜芯片數據的分析,相關基因與骨質疏松癥的確切關系還需要實驗證據的支持。

4 結論

本研究通過對3個骨質疏松癥相關表達譜芯片數據的GWGS整合分析,檢測到10個外周血細胞中與骨質疏松癥相關的樞紐基因。其中9個基因已有研究報道與骨質疏松癥相關,另外1個基因ELANE和骨質疏松癥的關系有待進一步探索。本研究結果將有助于進一步理解骨質疏松癥的分子致病機理。