生物信息學與機器學習識別診斷絕經后骨質疏松癥患者的生物標志物

李樹棟 梁學振 李剛

山東中醫藥大學第一臨床醫學院,山東 濟南 250014

絕經后骨質疏松癥(postmenopausal osteoporosis, PMOP)是世界范圍內最常見的疾病之一,其特征是雌激素缺乏以及隨著年齡的增長而持續的鈣流失,PMOP能夠導致骨脆性增加從而增加骨折發生風險[1]。研究表明目前全球有近30 %的絕經后婦女患有PMOP[2]。由PMOP引起的腕部及髖部骨折正逐年增加[3]。目前診斷PMOP最常見的方式是雙能X線檢查,但其診斷與應用存在局限性[4],近年來學者也嘗試從炎癥指數、基因易感性等方面尋找新的診斷手段[5-7]。研究表明,雌激素能夠在促進成骨分化和成骨細胞成熟、抑制破骨細胞形成的同時誘導破骨細胞凋亡,使得骨的形成和吸收處于相對穩定狀態,絕經后由于雌激素水平下降使得這種平衡狀態被打破,從而導致PMOP發生[8]。免疫學研究表明當雌激素水平下降時可誘發淋巴T細胞活化從而作用于骨髓間充質細胞及成骨細胞,以起到調節骨代謝平衡的作用,這被認為在骨質疏松癥發病機制中具有關鍵作用[9]。研究表明,免疫細胞分泌的細胞因子如IL-6、IL-17、IL-1和 TNF等都可以直接或間接作用于破骨細胞,促使破骨細胞分化發育,加速骨吸收,直接導致骨質疏松[10-12]。由于免疫細胞可以在不同的水平上與骨細胞相互作用,免疫細胞的改變被證明也是PMOP的發病機制之一。本研究旨在通過生物信息學和機器學習,尋找PMOP新的潛在生物標志物,為臨床診斷和防治PMOP提供新的思路。

1 材料和方法

1.1 資料

從美國國家生物技術信息中心(National Center for Biotechnology Information, NCBI)的基因表達數據庫(Gene Expression Omnibus, GEO)數據庫(https://www.ncbi.nlm.nih.gov/geo/)中檢索與PMOP相關的數據集,篩選標準為PMOP患者和健康人群的樣本,排除動物系模型。最終下載符合要求的數據集GSE56815(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE56815)和數據集GSE7429(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE7429),數據集基于GPL96平臺([HG-U133A] Affymetrix Human Genome U133A Array),GSE56815數據集包含40個樣本,其中20例PMOP患者的樣本,20例絕經后的健康人群樣本。GSE7429數據集包含20個樣本,其中10例PMOP患者樣本,10例健康人群樣本。

1.2 方法

1.2.1數據差異表達分析:使用R軟件(https://www.r-project.org/)中的Limma包(https://bioconductor.org/packages/release/bioc/html/limma.html)對樣本數據集GSE5681和數據集GSE7429進行差異表達,篩選標準設定為 adj.P<0.05,∣log2 fold change (FC)∣> 1.2差異表達數據用火山圖及聚類圖進行可視化展示。

1.2.2GO分析:使用R軟件ClusterProfiler包(https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html)對得出的差異基因進行基因本體分析。其中GO 分析主要包括生物學過程 (biological process, BP)、細胞組分 (cellular component, CC) 和分子功能 (molecular function, MF)三大部分。

1.2.3蛋白互作網絡構建:為進一步識別明顯差異表達基因,利用 STRING(search tool for the retrieval of interacting genes) 數據庫對獲得的差異基因進行相互作用分析 (protein-protein interaction,PPI),參考值設置confidence評分為>0.4,通過 Cytoscape軟件使PPI 網絡可視化。通過設置Degree數值對 PPI 網絡中的每個節點進行相關性評估,篩選出意義最大的前10個基因。

1.2.4篩選生物標志物:采用兩種機器學習的方式,通過R中的glmnet包(https://cran.r-project.org/web/packages/glmnet/index.html)構建LASSO模型,以確定與PMOP和健康標本顯著相關的基因。在此基礎上,通過支持向量機遞歸特征消除從元數據隊列中篩選出最佳基因,建立SVM-REF分類模型,進一步確定生物標志物在PMOP中的診斷價值。基于LASSO回歸分析,SVM-REF分析和PPI蛋白互作網絡篩選出候選生物標志物。

1.2.5CIBERSORT分析:通過 CIBERSORT進行反卷積分析,評估22個免疫細胞的免疫浸潤情況,并評估這些免疫細胞在正常樣本與PMOP樣本中和關鍵基因表達之間的關聯。

2 結果

2.1 差異基因

對數據集GSE56815和數據集GSE7429進行數據差異表達分析,共找出30個差異基因,其中9個為表達上調基因(OSTF1、RAB8B、SERPINF1、FUCA1、EIF5A、TCL1A、IGLL3P、HIST1H2AC、EIF2S3),21個為表達下調基因(C1D、KCNJ2、PMAIP1、HIRA、ADRB2、ZEB2、ELANE、CEACAM8、LCN2、LTF、CAMP、GBP1、S100A12、EGR2、SCD5、CRISP3、HP、MMP8、DEFA4、S100P、ARG1)。

2.2 GO富集分析

通過 R 軟件 ClusterProfiler 包進行 GO 功能分析共得到84個富集條目,得到 BP條目62個、CC條目15個、MF條目7個,BP主要富集在殺死其他組織細胞、對真菌的反應、對真菌的防御反應、細胞殺傷等方面。CC條目主要富集在分泌顆粒腔、細胞質囊腔、囊泡腔。MF主要富集在絲氨酸型內肽酶活性、絲氨酸水解酶活性、脂多糖、鐵離子結合。

2.3 蛋白互作網絡構建

通過 STRING數據庫預測并構建30個差異基因的PPI網絡,取Degree值最大的前十個基因即ELANE、CAMP、LTF、S100A12、LCN2、ARG1、MMP8、DEFA4、CEACAM8、HP。網絡由 10個節點和 34條邊組成(圖1)。

圖1 構建PP網絡后取Degree前10個節點

2.4 篩選生物標志物

通過LASSO分析共確定12個可能成為PMOP診斷標志物的候選基因,通過SVM-REF分析共確定29個候選基因。基于LASSO分析、SVM-ERF分析、PPI蛋白互作網絡得到三者有且僅有一個共同基因S100A12,該基因可能是診斷PMOP的關鍵基因。

2.5 S100A12在PMOP中的表達及診斷意義

研究發現與健康樣本相比較,S100A12在PMOP樣本中的表達顯著下調(圖2A)。對S100A12進行進一步ROC分析得到其AUG=0.658(圖2B)。

注:A :S100A12在PMOP樣本中的表達明顯下調;B:S100A12的ROC分析。

2.6 S100A12與免疫浸潤的關系

通過研究PMOP患者標本和健康人群標本中的S100A12表達水平以及免疫細胞浸潤狀態,以確定兩者之間的關系。通過CiberSort算法檢測在PMOP樣本中免疫細胞浸潤情況。此外在健康樣本和PMOP樣本之間的Macrophages M0表達異常。免疫細胞浸潤分析顯示S100A12與免疫細胞T cells CD8、T cells CD4 memory resting、T cells CD4 memory activated、Plasma cells、Monocytes、Mast cells resting、Macrophages M0、Macrophages M1 、Macrophages M2、Eosinophils、Dendritic cells resting等相關。提示S100A12可能通過調節多個免疫細胞參與了PMOP的進展,表明S100A12可能成為診斷PMOP的生物標志物。

3 討論

絕經后骨質疏松癥常出現在女性絕經期后,其主要危險因素是雌激素缺乏、吸煙、鈣的缺失[13]。PMOP會增加脊柱、髖部、前臂遠端及肱骨近端等部位的骨折發病風險,研究表明,絕經期女性在這些部位發生骨折的風險高達40 %以上[14]。

絕經后骨質疏松癥具體的發病機制尚未明確,研究表明作為多系統共同參與的代謝性疾病,PMOP的發病機制已經不再是孤立的“雌激素-骨代謝紊亂-骨質流失”,其與雌激素缺乏誘導的氧化應激、腸道菌群、鐵過載等多因素參與有關[15],采用多靶點、多系統聯合用藥也將是PMOP新的治療思路。

及早發現PMOP的高危患者有助于預防骨丟失,及時診斷PMOP對于其治療顯得尤為重要,但臨床實踐中仍然缺乏早期診斷手段。隨著轉錄組學的進步、生物信息學的發展及高通量數據庫的建立,為研究諸如此類疾病的潛在生物標志物提供了有效的數據支撐,高通量數據庫及相關研究表明能夠通過調控相關基因的差異表達從而影響相關疾病的發生和發展[16-22]。而近些年對PMOP發病機制的研究得出其與免疫細胞息息相關[9,11,15]。因此本研究旨在通過挖掘GEO數據庫中的數據信息分析差異表達基因,并通過機器學習證實其通過調節免疫細胞參與PMOP的發生,以期為絕經后骨質疏松癥的防治提供理論基礎和數據支撐。

研究中通過對GEO數據庫進行分析,得出PMOP標本與健康標本共存在30個差異基因,GO分析顯示,這30個差異基因主要參與殺死其他組織細胞、對真菌的反應、對真菌的防御反應、細胞殺傷、細胞質囊腔、囊泡腔、絲氨酸型內肽酶活性、絲氨酸水解酶活性、鐵離子結合等。這些結果表明,這些差異基因可能對于PMOP的發展起著重要作用。

為了篩選出PMOP的潛在生物標志物,本研究對30個差異基因構建PPI蛋白互作網絡取出Degree前10位的基因,同時對30個差異基因進行LASSO模型分析和SVM-REF,篩選出PMOP患者中潛在的關鍵基因,最終識別出S100A12。S100A12是S100家族成員,是一種位于粒細胞胞質中的低分子量的鈣結合蛋白。其表達過程多由中性粒細胞和巨噬細胞完成[23]。研究表明S100A12增強了成熟破骨細胞相關分子NFATc1、ACP5、CALCR和ITGβ3的表達,通過RAGE和TLR4通路直接促進破骨細胞形成和骨吸收。同時S100A12給藥顯著增強了破骨細胞的骨吸收能力[24]。然而,S100A12在PMOP中的作用機制仍不明確。本研究首次發現在PMOP樣本與健康樣本中S100A12的表達水平異常。ROC分析進一步證實可以通過S100A12對PMOP樣本和健康樣本進行鑒別篩選。研究表明S100A12的表達水平在PMOP樣本中明顯下調,因而S100A12可能成為診斷PMOP的潛在生物標志物。

近年來,越來越多的研究表明,免疫細胞浸潤在PMOP的發生和發展過程中起著至關重要的作用。從免疫系統的角度,評估免疫細胞浸潤,識別浸潤免疫細胞成分的多樣性,對于揭示PMOP的分子水平因果關系和尋找新的免疫治療靶點至關重要[25-26]。本研究發現,S100A12與免疫細胞T cells CD8、T cells CD4 memory resting、T cells CD4 memory activated、Plasma cells、Monocytes、Mast cells resting、Macrophages M0、Macrophages M1 、Macrophages M2、Eosinophils、Dendritic cells resting等相關,因而筆者推測S100A12通過調節多種免疫細胞參與PMOP的發生發展。

S100A12在骨代謝疾病中的研究相對較少,其研究主要集中在心腦血管疾病和肺部疾病[27-31],本研究表明S100A12參與免疫細胞浸潤,其可能是骨免疫環境中的關鍵分子,進一步實驗可能為研究S100A12找到新的方向。盡管研究認為S100A12可能成為診斷POMP的潛在生物學標志物,但該類研究仍存在部分限制。首先,數據集及包含樣本的數量較少,缺乏大樣本研究和動物實驗從外部驗證生物信息學分析;其次,研究僅能支持與S100A12免疫細胞間的相關性分析,并不能揭示其因果關系,其作用機制有待研究,仍需要更多研究來揭示基因和免疫細胞之間復雜的相互作用關系。

綜上所述,基于GEO高通量數據庫,借助生物信息學和機器學習,通過對PMOP樣本與健康樣本比較篩選得到S100A12是PMOP發生發展過程中的關鍵基因。本研究結果為揭示PMOP的潛在分子水平因果關系提供了思路,并為進一步探討和治療PMOP提供了潛在靶點和參考方向。