基于CIBERSORT反卷積算法研究強直性脊柱炎差異基因的免疫浸潤機制及相關中藥預測的生物信息學分析*

丁 燦，沈 浮，劉 鑫，楊 雷

（1.湖南中醫藥大學第一附屬醫院，湖南 長沙 410007；2.永州市中醫醫院，湖南 永州 425000）

強直性脊柱炎（ankylosing spondylitis，AS）是一種常見的血清陰性型脊柱炎，以慢性炎癥和新骨形成為主要特征，常累及脊柱及骶髂關節[1]。流行病學提示其全球發病率為0.1%～1.4%，男性發病率高于女性，晚期AS可導致活動受限和關節僵硬，甚者致殘而嚴重影響患者生活質量[2-3]。AS的確切原因尚不清楚，目前一般認為其與遺傳及環境因素導致的自身免疫機制紊亂所致的慢性炎癥相關[4-5]。近年來，越來越多學者認為，免疫細胞浸潤在AS的發生、發展和治療中起著至關重要的作用。如樹突狀細胞、巨噬細胞、自然殺傷細胞、輔助性T1細胞和輔助性T17細胞都與AS的發展有重要關系[6]。因此，評估免疫細胞浸潤和確定免疫細胞成分的差異對闡明AS進展的分子機制和開發新的免疫治療靶點具有重要價值。在臨床實踐中，AS的治療主要依賴于非甾體抗炎藥及生物制劑[7]，而此類藥物長期使用具有心血管、胃腸道及腎臟不良反應的潛在風險[8]，而中藥的較好療效、低副作用和低成本，以及對于免疫失衡的調節干預，使得其成為治療AS較為優越的替代療法。

隨著高通量遺傳分析的出現，基因表達譜分析成為了鑒定各種疾病差異表達基因的有效方法，有助于探索發病機制和開發生物標志物。迄今為止，已有多項研究使用微陣列表達來鑒定AS發病機制中的關鍵基因[9-10]。CIBERSORT是一種分析工具，它使用微陣列數據或RNA測序數據來評估樣本中免疫細胞的表達，并獲得各種免疫細胞比率[11]。本研究應用基因芯片數據庫（gene expression omnibus，GEO）中AS患者與正常人群的全血差異樞紐基因進行相關生物信息學分析，采用廣泛用于評估鑒定22種免疫細胞相對含量的CIBERSORT反卷積算法進行差異基因的免疫浸潤分析研究，并篩選與免疫浸潤相關的關鍵生物標志物，從而預測免疫相關靶點及中藥復方，旨在為AS的診斷和治療提供依據。

1 資料與方法

1.1 數據收集 以“Ankylosing Spondylitis”為關鍵詞檢索GEO基因表達綜合數據庫，篩選得到“GSE25101”芯片為本次研究分析的對象。該芯片主要包含16例AS患者及16名年齡性別對應的健康人全血基因表達譜，通過GPL6947 Illumina HumanHT-12 V3.0表達珠芯片平臺將探針轉換成基因名稱，利用R語言程序中“Limmar”安裝包對每個樣本的微陣列數據進行預處理，剔除缺失值，當基因對應于多個探針時，取最大值[12]。

1.2 關鍵差異樞紐基因的蛋白網絡互作（PPI）及功能注釋 在STRING數據庫（https://string-db.org/）獲取蛋白互作網絡模型（protein-protein Interaction，PPI）網絡，在Cytoscape軟件中利用CentiScape插件計算每個靶基因的DC（度）值，以平均度值的2倍進行核心靶基因的篩選并進行可視化[13]。運用R語言程序中的“clusterProfiler”插件包進行差異基因的基因本體論（Gene Ontology，GO）與京都基因和基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）富集分析。

1.3 免疫浸潤矩陣整理 本研究利用CIBERSORT反卷積算法進行免疫細胞浸潤分析，包括T細胞、B細胞、NK細胞等22種免疫細胞在內的轉錄特征矩陣的模擬計算，設定次數為100次，對P＜0.05的數據進行后續分析。

1.4 免疫細胞間相關性分析及組間差異分析 為了分析各類免疫細胞之間的相關性，本研究在CIBERSORT反卷積算法以P＜0.05篩選可信樣本數據，然后計算可信樣本數據中不同免疫細胞間的Pearson相關系數，并運用秩和檢驗分析AS患者與健康人群之間的差異。

1.5 關鍵差異樞紐基因與免疫細胞浸潤水平關系的相關性分析 為進一步研究關鍵差異樞紐基因與疾病相關性程度，本研究將篩選所得的差異樞紐基因與差異性較為顯著的免疫細胞進行相關性分析，分析靶基因與各免疫細胞表達是否存在相關性并計算相關性系數r值。|r|＜0.4為弱相關，0.4＜|r|＜0.7為中等強度相關，|r|＞0.7為高度相關。

1.6 AS免疫浸潤過程中潛在中藥及復方的篩選預測 將AS差異核心靶基因及GO富集的與免疫浸潤相關生物學過程（biological process，BP）導入Coremine Medical數據庫（https://coremine.com/）及本草組鑒數據庫（http://herb.ac.cn/）中，篩選具有潛在效應機制的中藥。Coremine Medical數據庫中的藥物篩選標準為P＜0.05；本草組鑒數據庫中篩選標準為P-value＜0.05，且FDR_BH＜0.05。通過閱讀文獻及古籍，基于上述所預測的藥物進行復方預測，使得復方盡可能的包含較多的預測藥物。

2 結果

2.1 AS差異表達基因獲取結果 GEO數據庫檢索及篩選結果顯示，AS組與正常組的差異基因為135個，其中上調基因為99個，下調基因為36個。圖1A熱圖展示了前20個上調及下調基因的情況，圖1B為差異基因的分布情況。

圖1 AS 相關基因探針差異表達的熱圖及火山圖

2.2 AS關鍵樞紐基因獲取 利用STRING數據庫對差異基因進行蛋白網絡互作分析（PPI），將分析結果的TSV文件導入Cytoscape中進行可視化展示。圖2A中節點顏色越深代表該節點的度值越大，在PPI模型中的重要程度越高。圖2B可視化展示了度值大于2倍平均度值（15.8）的關鍵樞紐基因，其中度值大于15.8的有14個，分別是SNRPG、LSM3、PSMA6、MRPL22、RPL26、COX7C、PSMA4、RPS17、RPL31、SLIRP、NDUFS5、RPS24、NDUFB3、SRP14。

圖2 AS 差異基因的PPI 網絡圖及核心靶基因

2.3 GO及KEGG富集分析結果 運用R語言程序中的“ClusterProfiler”插件包分析AS差異基因進行GO及KEGG富集分析。差異基因的GO分析中共獲取生物過程（biological process，BP）77條，細胞組成（cellular component，CC）28條，分子功能（molecular fuction，MF）17條，對其中基因富集數目最多的10條進行條形圖展示。根據基因的富集數目，表1展示了與免疫相關的生物學過程，主要涉及中性粒細胞脫粒（Neutrophil degranulation）、免疫反應中的中性粒細胞激活（Neutrophil activation involved in immune response）、先天免疫反應的調節（Regulation of innate immune response）、核糖體生物發生（Ribosome biogenesis）、氧化磷酸化（Oxidative phosphorylation）、NF-κB轉錄因子活性的正調控（Positive regulation of NF-κB transcription factor activity）、中性粒細胞趨化性（Neutrophil chemotaxis）、細胞防御反應（Cellular defense response）、白細胞介素-1介導的信號通路（Interleukin-1-mediated signaling pathway）、T細胞介導細胞毒性（T cell mediated cytotoxicity）。（見圖3）在KEGG富集分析中，差異基因主要涉及核糖體合成、氧化磷酸化等相關通路，與代謝、炎性反應相關通路聯系密切。

表1 AS 與對照組差異基因免疫浸潤相關GO 富集分析列表

圖3 AS 差異基因的KEGG 及相關基因可視化

2.4 AS患者與正常組的免疫細胞分布特征 通過CIBERSORT反卷積算法以P＜0.05樣本進行篩選，最終獲取32個可信樣本，如圖4A所示。熱圖左側16個樣本（GSM616684～GSM616699）為正常組，右側16個樣本（GSM616668-GSM616683）為AS組，結果顯示，單核細胞（Monocytes）、CD8 T細胞（T cells CD8）、中性粒細胞（Neutrophils）、初始CD4+T細胞（T cells CD4 naive）在兩組中含量均較高，且AS組中單核細胞（Monocytes）含量高于正常組。圖4B顯示，免疫細胞比例箱圖中進一步展示了22種免疫細胞的浸潤比例。

圖4 AS 組與正常組中免疫相關細胞的浸潤組成

2.5 免疫細胞間相關性分析 未活化的CD4+記憶性T細胞（T cells CD4 memory resting）與γδT細胞（T cell gamma delta）呈中度正相關（r=0.52），CD8+T細胞（T cells CD8）與活化的自然殺傷細胞（NK cell activated）呈中度正相關（r=0.51），而CD8+T細胞（T cells CD8）與中性粒細胞（Neutrophils）呈中度負相關（r=-0.65），初始CD4+T細胞（T cells CD4 naive）與未活化的CD4+記憶性T細胞（T cells CD4 memory resting）呈中度負相關（r=-0.69）。（見圖5）

圖5 AS 組中免疫細胞間的相關性分析

2.6 AS組與正常組的免疫細胞浸潤差異分析 通過小提琴圖對不同組織樣本免疫細胞浸潤差異分析進行可視化，以P＜0.05為顯著差異性。結果顯示，AS組血清中單核細胞（Monocytes）明顯多于正常組（P=0.002），而正常組血清中調節性T細胞（T cells regulatoy，Tregs）明顯少于AS組（P=0.03）。（見圖6）

圖6 AS 組與正常組間的免疫細胞浸潤差異分析小提琴圖

2.7 關鍵差異樞紐基因與免疫細胞浸潤水平關系的相關性分析 將AS關鍵樞紐基因與上述所得差異性顯著的免疫細胞（單核細胞、調節性T細胞）進行相關性分析，對相關性較為顯著的基因進行展示。結果顯示，SNRPG（r=-0.26）、LSM3（r=-0.14）、PSMA6（r=-0.27）、MRPL22（r=-0.18）與單核細胞呈負相關（P＜0.05），但相關性較弱；SNRPG（r=-0.30）、LSM3（r=-0.16）、PSMA6（r=-0.20）、MRPL22（r=-0.16）與調節性T細胞呈負相關（P＜0.05），但相關性較弱。（見圖7）

圖7 AS 關鍵樞紐基因與免疫細胞浸潤水平相關性分析

2.8 干預AS免疫細胞浸潤相關生物途徑的中藥及復方預測通過Coremine Medical及本草組鑒預測具有潛在治療作用的中藥。圖8顯示了AS核心差異靶基因及相關免疫反應生物學過程相關的中藥，其中雷公藤、蜈蚣、冬蟲夏草、紅花、人參、銀杏葉等度值較高，與AS相關免疫細胞浸潤的關系最密切。藥物歸經多歸為肝、腎、肺經，如圖9A所示；功效多與清熱補虛及祛風濕相關，如圖9B所示。表2為AS免疫細胞浸潤相關生物途徑的中藥復方預測結果，涉及補腎活血湯、化痰逐瘀散結湯、寄生湯等多種方劑，功效多為補腎健脾、活血通絡，兼祛風散寒化痰。

表2 AS 免疫細胞浸潤相關的中藥復方預測

圖8 AS 免疫細胞浸潤相關的中藥預測網絡圖

圖9 AS 免疫細胞浸潤相關的中藥歸經功效統計

3 討論

AS是一種多種復雜因素誘發的疾病，隨著現代醫學的進步，AS的診斷與治療已獲得了飛速的發展，但由于其病因的不確切性，目前尚無治愈的療法[14]。因此深入研究AS的發病機制，為AS的治療提供潛在的靶點對于AS的診治具有重要意義。本研究通過GEO數據庫篩選合適的基因芯片，比較正常人群與AS患者血清樣本的基因表達譜，得到差異基因并通過蛋白網絡互作獲得其核心靶基因。核心靶基因分別是SNRPG、LSM3、PSMA6、MRPL22、RPL26、COX7C、PSMA4、RPS17、RPL31、SLIRP、NDUFS5、RPS24、NDUFB3、SRP14。有研究表明SNRPG與腸道中巨噬細胞活化相關[15]，而產生巨噬細胞遷移抑制因子（MIF）的巨噬細胞在AS患者腸道中富集，MIF在AS患者的關節液中處于高水平，由此推測SNRPG干預AS的機制可能與激活腸道中MIF相關。LSM3基因在AS患者中表達下降，LSM3基因作為預催化的前體參與前體mRNA剪接，而前體mRNA在加工和剪接過程中可能會影響AS患者其他基因的表達[16]。PSMA4、PSMA6基因是蛋白酶體功能所需的20S蛋白酶體核心，蛋白酶體復合物是一種蛋白水解系統，通過主要組織相容性復合物在AS適應性免疫系統的功能中發揮關鍵作用[17]。MRPL22在細胞質中生成，屬于基因編碼線粒體核糖體蛋白的成員，對線粒體氧化磷酸化至關重要，且在AS中調節凋亡誘導因子方面發揮關鍵的作用[18]。RPL26通常構成參與核糖體亞基翻譯的細胞過程，而敲低該基因的表達，能使成纖維細胞樣滑膜細胞中的NF-κB信號通路去磷酸化，進一步降低AS相關風險基因的轉錄[19]。COX7C是細胞色素氧化酶（COX）的組成部分[20]。這對能量的生產至關重要，與本研究中的KEGG富集分析中的核糖體生物發生、氧化磷酸化密切相關，并且MRPL22、RPS17、RPL31、NDUFS5、RPS24等核糖體相關蛋白可能與AS發病的能量代謝有關。而其他基因尚缺乏與AS相關的研究，可能是AS的潛在治療靶點，需進一步發掘驗證。本研究中差異基因蛋白互作中的部分樞紐基因有AS的相關基礎研究的結果支持，進一步說明了本研究預測結果的可靠性，但其他基因的表達或功能在AS的發病基礎研究中鮮有提及。這些基因可能為AS診斷和治療靶點的研究提供新的思路。

本研究進一步對AS相關差異基因進行GO及KEGG富集分析。結果顯示，差異基因的生物學過程不僅有T細胞介導細胞毒性反應參與，還有更為豐富的中性粒細胞活化介導的相關免疫反應，包括中性粒細胞趨化性、中性粒細胞脫粒等生物過程。目前現代醫學認為，AS是一種由T淋巴細胞介導的慢性炎癥性疾病，尤其是CD4+T細胞的紊亂在其病程的發展中發揮關鍵性作用，最近的研究表明，先天免疫細胞在炎癥和新骨形成過程中也很重要，而中性粒細胞越來越被認為是自身炎癥和自身免疫性疾病的介質[21]，且中性粒細胞通過諸如TLR受體、吞噬、脫顆粒和活性酶的釋放來實現上述過程，即細胞以網狀物的形式釋放其染色質[22]。NF-κB介導的免疫調節在AS中發揮著樞紐作用，Th1/Th2的平衡在免疫調節中發揮重要作用。I-κB激酶的激活使得NF-κB的抑制作用解除，從而激發不同Th因子的合成和釋放，進而調節免疫紊亂[23]。白細胞介素-1（IL-1）介導的信號通路在AS差異基因的富集上也表現出良好的聚集性。IL-1由炎癥滑膜中活化的巨噬細胞分泌，能啟動免疫細胞的募集和炎癥反應。有研究顯示，IL-1多態性與AS易感性相關，IL-1基因簇中的IL-1A、IL-1B和IL-1RN可觸發AS先天性免疫反應[24]。除此之外，AS的差異基因KEGG通路富集中涉及核糖體生物發生、氧化磷酸化等相關通路，與代謝、炎癥反應相關通路聯系密切。

本研究利用CIBERSORT反卷積算法探索免疫細胞浸潤在AS發病機制的作用，結果顯示，AS組患者血清中單核細胞較正常組明顯增多，而調節性T細胞（Treg）則明顯減少。Treg主要分為胸腺來源的Treg細胞（tTreg）及誘導調節性T細胞（iTreg），其是CD4+T細胞的一小部分，在維持免疫耐受和自身免疫方面起著關鍵作用[25]。Treg在AS中的含量變化是有爭議的，但多數學者認為Treg比例下降是誘發AS重要因素，其機制可能為AS患者Treg細胞中STAT5磷酸化相對較少，FOXP3的平均熒光強度降低，導致Treg細胞不能有效抑制幼稚T細胞的增殖，進而導致疾病的進展[26]。單核細胞是專業抗原提呈細胞的前體。有研究顯示，AS患者的血單核細胞的微陣列分化標記物上調，與正常單核細胞相比，AS單核細胞參與了多種炎癥途徑的蛋白質水平升高[27]。免疫浸潤相關性分析結果中顯示，CD8+T細胞與活化的自然殺傷細胞呈中度正相關性。在一項AS患者外周血淋巴細胞亞群變化規律的研究中，其NK細胞含量伴隨T細胞顯著升高，表明AS的免疫亢進伴有非特異性[28]。γδT細胞介于固有免疫與獲得性免疫之間，在免疫反應中發揮重要作用。本研究中，CD4+記憶性T細胞與γδT細胞呈中度正相關，說明兩者在AS的發病中發揮協同作用。而CD8+T細胞與中性粒細胞呈中度負相關性的基礎研究尚缺乏大規模的驗證，故本預測具有一定程度的參考意義。

AS在中醫學屬“大僂”“僂痹”等范疇。病機為五臟內傷、外感六淫，又以肝腎虧虛、督脈虛空為要，病性多為本虛標實[29]。中醫藥治療AS療效穩固且持久，對于調節患者免疫失衡具有獨特的優勢。本課題組通過Coremine Medical及本草組鑒預測具有潛在治療作用的中藥，其中雷公藤、蜈蚣、冬蟲夏草等出現的頻數較高且度值較高，提示與AS相關免疫浸潤的關系最密切。藥物歸經多歸肝、腎、肺經，功效多與清熱補虛及祛風濕相關，此與AS的病機相符。其中雷公藤中的主要成分雷公藤甲素在體外可減少膠原形成、堿性磷酸酶活性和鈣沉積，其對AS的免疫調節作用已經得到證實。雷公藤甲素可能通過上調外周血中的CD4+CD25+CD127+Treg和下調IL-17來實現免疫調節作用[30-31]。蜈蚣在骨痹的治療中使用頻率較高，因其走竄力強，臟腑經絡及氣血不通之處皆能通之，并且蜈蚣在關節炎的局部炎癥細胞浸潤中表現為顯著的抑炎反應及軟骨的保護作用[32]。冬蟲夏草素目前尚缺乏在AS免疫調節方面的基礎研究，但其在類風濕關節炎小鼠模型中有著顯著的平衡免疫的功能，具體機制為調控MMP-3、TIMP-1而抑制炎癥細胞因子IL-1β、IL-10，進而平衡Ⅱ型膠原免疫紊亂[33]。在方劑預測中，補腎活血湯、寄生湯中多數藥物與預測中藥相符合。有臨床研究表明，補腎活血湯能顯著改善AS患者的疼痛、晨僵不適及活動性炎性指標[34]，且有報道稱補腎活血湯含藥血清對CD4+CD25+Treg增殖具有明顯促進作用[35]。而在本研究中，樞紐差異基因SNRPG、LSM3、PSMA6、MRPL22與Tregs呈患者負相關，但相關性較弱。故筆者推測，上述靶基因可能作為補腎活血湯干預AS免疫浸潤相關機制的靶點，但需進一步實驗驗證。

本研究局限之處在于數據分析來源依賴于GEO數據庫中所納入的AS探針樣本，需要進一步實驗驗證，而單一芯片納入樣本的數量局限可能造成偏倚。總之本課題組利用生物信息學方法對AS的免疫浸潤機制進行分析并預測相關生物途徑所涉及的中藥及復方，對于后續研究及臨床用藥具有參考意義。