通過生物信息學方法構建系統性紅斑狼瘡患者的基因表達亞群

馬江磊，陳華秋，王光明，2*

（1.大理大學臨床醫學院，云南 大理 671000；2.大理大學第一附屬醫院基因檢測中心，云南 大理 671000）

系統性紅斑狼瘡（systemic lupus erythematosus，SLE）是一種以免疫調節異常、自身免疫復合物沉積為特征的慢性自身免疫性疾病，常累及多臟器損害，甚至危及患者的生命。女性患病率高于男性，比例約為10∶1〔1〕。在過去的半個多世紀，隨著人們對SLE認識逐漸增多，SLE患者的生存率和治療效果有所提高〔2〕，但SLE臨床癥狀表現多樣，仍有部分患者因診斷延誤或誤診，未得到及時治療，從而產生嚴重后果。因此，建立準確、科學的分類和有針對性的治療，對改善SLE患者的治療效果有十分重要的意義。

SLE發病機制復雜，目前已知除了與紫外線照射、氣候和地理分布等外界因素有關外，還發現多種與SLE有關的易感基因〔3〕。Wang等〔4〕通過對17 389名SLE患者和29 273名健康人群數據進行薈萃分析后發現，在歐洲及亞洲的總體人群中，信號轉導和轉錄激活因子（signal transducer and activator of transcription，STAT）4基因rs7574865、rs10168266多態性與SLE發病風險密切相關。Sandling等〔5〕采用以信號通路為中心的方法，對瑞典958名SLE患者和1 026名健康個體的1 832個候選基因進行了免疫通路基因測序，確定了T淋巴細胞分化和先天免疫2個涉及SLE易感性的主要獨立途徑，通路風險評分有可能將SLE患者臨床表現以外的分子進行分類，有助于臨床隨訪和治療的選擇。

目前，針對SLE的基因治療還不完善，為促進SLE治療新策略的開發，本研究對1 254例SLE患者的轉錄組數據進行匯總分析，建立共識聚類表達分析，分成相應的亞群，并比較各亞群之間在臨床特點、分子構成及細胞學功能方面的差異，為SLE的基因分型提供相關依據。

1 材料與方法

1.1 數據集的富集分析和主成分分析（principal componentanalysis，PCA）在NCBI的GeneExpression Omnibus（GEO）數據庫（https：//www.ncbi.nlm.nih.gov/geo/）中搜索“systemic lupus erythematosus”，下載GSE121239、GSE65391和GSE154851的微陣列數據集，各芯片數據集中只包含SLE患者及健康對照者的組織。GSE121239由GPL13158平臺文件提供，GSE65391和GSE154851分別來自GPL10558、GPL16699平臺文件。使用R 4.1.3軟件中的“limma”和“sva”數據包消除批次效應，進行歸一化〔6-7〕。將數據合并后得到1 254例SLE患者（以下稱“轉錄組”）和124例健康對照者（以下稱“對照組”）的基因表達數據集。其中臨床數據包括年齡、性別、SLE活動指數（systemic lupus erythematosus disease activity index，SLEDAI）等。

1.2 亞群的構建和亞群間臨床特點的比較1 254例SLE患者的基因表達矩陣利用R軟件中的“ggplot2”“rstatix”“reshape2”和“Consensus Cluster Plus”等數據包對矯正后的樣本數據集進行分析，得到具有相似差異表達基因（differentially expressed genes，DEGs）的SLE患者分型。每個亞群的聚類得分越高，說明各亞群基因在基因表達上的相似性越高，然后把每個亞群視為一個聚類進行分析。

使用R軟件中“ggplot2”和“ggpurb”數據包對各亞群的臨床特點進行分析。以患者性別作為離散變量，分析男性和女性在各亞群中所占的比例；以年齡和SLEDAI作為連續變量，進行差異富集分析得到相關箱式圖。

1.3 篩選亞群特異性基因及構建蛋白質相互作用（protein-protein interaction，PPI）網絡將各亞群與對照組以及各亞群之間進行比較后，從中篩選出均值差異>0.2、調整后P<0.05的基因，得到各亞群對應的DEGs，將所得到的DEGs取交集，進一步篩選出只在某個亞群中表達上調的基因，即亞群的特異性基因。取各群前100個特異性基因上傳到STRING（https：//stringdb.org/）網站，設置置信系數為0.4，生成特異性基因表達蛋白的PPI網絡圖，再將PPI網絡圖傳輸到Cytoscape軟件中，以Degree算法計算得到節點數最多的前10個蛋白并對其進行分析。

1.4 各亞群基因集的富集分析對基因集富集分析（gene set enrichment analysis，GSEA）。GSEA前對基因列表文件和基因集文件進行對應轉換，將每個組得到的列表和數據集文件傳輸到GSEA軟件（version 4.1.0）中，設置最大運行值為5 000，得到各亞群中特異性基因在對照組中的表達分析結果。

1.5 GO功能富集分析和KEGG信號通路分析選擇矯正后的基因數據集、各亞群中特異性上調基因文件以及含有連續變量的臨床特征文件，使用R軟件的“WGCNA”數據包進行分析。對樣本臨床特征與DEGs和模塊之間的相關性分析后通過動態樹切割法〔8〕，將基因數據集分成顏色不同的基因模塊。應用R軟件對分成的基因顏色模塊進行GO功能富集分析和KEGG信號通路分析。GO功能富集分析主要對生物學過程、細胞成分和分子功能進行富集分析。KEGG信號通路分析主要對分子生物學途徑進行分析。對KEGG信號通路分析的結果，篩選出每個模塊中差異最顯著的信號通路，分析各亞群與生物學信息之間更密切的聯系。

1.6 統計分析用SPSS 26.0軟件進行數據分析，數據以（±s）表示。組間比較采用單因素方差分析，不同組間的比較使用獨立樣本t檢驗，P＜0.05為差異有統計學意義。

2 結果

2.1 轉錄組數據的預處理及分組對“GSE121239”“GSE65391”和“GSE154851”3個數據集進行預處理及批次矯正后，得到總數據集的PCA圖。3個數據集在批處理前，各組間無明顯相關性。見圖1A。消除批次效應后，總體基因分布均勻集中。見圖1B。根據每個樣本中的基因與SLE相關情況進行聚類分析，設置2~10個聚類，根據一致性評分結果，可將樣本分為3個亞群，分別為亞群Ⅰ、亞群Ⅱ及亞群Ⅲ。表達差異富集模式圖的分析結果也表明分為3個亞群時，組間表達模式顯著不同。見圖2。

2.2 亞群間臨床特點的比較結果本研究納入“GSE65391”和“GSE154851”數據集中（“GSE121239”數據集中未標注性別信息）患者1 066例，男124例，女942例，男女性別比為1.00∶7.60。把性別作為每個亞群的離散變量進行亞群間分析比較后發現，性別在各亞群間無顯著差異，說明性別不是區分各亞群的要素。在3個亞群中，亞群Ⅰ的平均年齡為（15.24±5.10）歲，亞群Ⅱ的平均年齡為（16.09±5.76）歲，亞群Ⅲ的平均年齡為（15.09±6.00）歲，富集比較結果顯示，亞群Ⅰ和亞群Ⅱ、亞群Ⅱ和亞群Ⅲ之間年齡差異有統計學意義（P<0.05）。在3個亞群中，亞群Ⅰ的SLEDAI為（4.03±3.11），亞群Ⅱ的SLEDAI為（1.88±2.80），亞群Ⅲ的SLEDAI為（2.68±3.22）。SLEDAI在各亞群中差異均有統計學意義（P<0.05）。見圖3。

2.3 各亞群特異性基因與PPI網絡圖分析3個亞群的特異性上調基因數量分別為：亞群Ⅰ（2 134個）、亞群Ⅱ（105個）、亞群Ⅲ（2 289個），進一步篩選出各亞群中前10個表達差異最顯著的DEGs。見表1。選擇各亞群中上調前100個的特異性基因，對其表達的蛋白質進行PPI網絡分析，使用Degree算法得到節點>10的蛋白樞紐圖，其中，排名前10的蛋 白 分 別 是STAT3、TLR4、BRIX1、TLR2、PTEN、HSPA8、NHP2、NOP58、CAT、DDX18。見圖4。

表1 各亞群中前10個上調的差異表達基因

2.4 各亞群的GSEA結果通過對每個亞群中的DEGs與對照組進行比較分析，發現各組中黑色線及灰色線均靠左聚集，GSEA結果顯示各亞群的P值和FDR值均小于0.01，說明各亞群中的特異差異基因在對照組也存在明顯的表達差異。見圖5。

2.5 GO功能富集和KEGG信號通路分析結果在測試1~20閾值功率的結果中得到最佳功率值為9，在功率值為9的情況之下，基因之間的連通性符合無標度網絡分布，可分為8個顏色模塊。各基因顏色模塊的分析結果見表2～5。

表2 生物學過程富集分析結果

GO功能富集分析結果顯示，在生物學過程中，粉紅色模塊中的基因與血液凝固、止血途徑以及血小板聚集等聯系密切。在細胞成分中，黏著斑、細胞-基質連接以及血小板α顆粒等與粉紅色模塊密切相關。分子功能結果顯示，粉紅色模塊主要與結合肌動蛋白、結合整合素、結合膠原等有關。從KEGG信號通路富集結果中，篩選出每個顏色模塊中表達差異最顯著的信號通路，分別為：自然殺傷細胞介導的細胞毒作用、核糖體、線粒體自噬、細胞凋亡、血小板活化、造血細胞譜系、破骨細胞分化、甲型流感，以進一步分析與亞群之間的關系。

將上述信號通路與對照組及3個亞群中的富集表達情況進行分析，結果顯示，自然殺傷細胞介導的細胞毒作用、核糖體和造血細胞譜系在亞群Ⅲ中高表達，而在亞群Ⅰ中低表達；亞群Ⅱ在各信號路徑中表達無明顯差別，說明從生物學功能層面可以區別不同的亞群。

表3 細胞成分富集分析結果

表4 分子功能富集分析結果

續表4

表5 KEGG信號通路富集分析結果

3 討論

本研究通過對SLE患者的轉錄組基因進行分析，根據共識聚類表達，將所有基因分為3個亞群，對各亞群的臨床特點進行對比分析，探討了年齡、性別以及SLEDAI與各亞群間的聯系。利用GO功能和KEGG信號通路富集，對各亞群的生物學過程、細胞成分、分子功能和相關信號通路的關聯情況進行了分析，為SLE基因亞型的認識提供相應的依據。在所納入的1 254例SLE患者中，性別在各亞群中差異無統計學意義，而年齡和SLEDAI在各亞群中差異有統計學意義，說明年齡和SLEDAI在SLE的病理生理過程及預后中起到重要作用。有研究〔9〕發現，年齡越小，其臨床表現越嚴重，且發生治療副作用的風險更高，病死率也會增加。SLEDAI是評估SLE臨床狀態和治療效果的關鍵指標，多年來，在臨床工作中被廣泛應用〔10〕。SLEDAI-2K是對SLEDAI的修改，當SLEDAI-2K>4時，提示存在狼瘡疾病活動〔11〕。

利用STRING網站分析，得到了網絡節點最多的蛋白，通過表達蛋白研究其對應的基因可以看出，排名靠前的分別是STAT3、TLR4、BRIX1、TLR2。其中，STAT3是節點最多的基因，但它并不在每個亞群的前10個特異性上調基因中。STAT3是轉錄信號傳導子與激活子通路的重要成員，在調節炎癥和免疫應答反應中起到關鍵作用〔12〕。當T細胞功能缺陷時，導致STAT3信號通路中白細胞介素（interleukin，IL）-17和IL-21的異常表達，IL-17細胞因子可以介導中性粒細胞和其他促炎細胞因子的釋放，從而引起炎癥反應的發生〔13-14〕。Chen等〔15〕發現在SLE患者中，IL-17表達上調且與STAT3強度呈正相關，說明STAT3和SLE的發生有著密切聯系，這與樞紐圖分析的結果相同。

先天免疫系統是宿主抵御病原體侵襲的第一道防線，Toll樣受體（Toll-like receptor，TLR）是固有免疫受體的重要組成部分，其在B細胞參與自身免疫性疾病的機制中起到產生抗體、T細胞的抗原呈遞和促進細胞因子生成等重要作用〔16-17〕。目前已知存在于人類的TLR有11種（TLR1~11），TLR4是一種跨膜受體，屬于富含亮氨酸重復序列的超家族蛋白，在哺乳動物中是第一個被定義的TLR成員，它既能在質膜上發揮作用，也能在核內體上發揮作用〔18-19〕。TLR2是另一種TLR受體，當TLR2被識別并與配體結合后，會激活細胞內信號轉導通路，引起免疫應答，可介導多種生物學過程，包括多發性硬化癥、天皰瘡、特發性皮炎、白癜風和1型糖尿病等〔16，20〕。Ma等〔21〕發現TLR4在狼瘡腎炎的發生發展過程中具有重要的致病作用，可能為狼瘡的治療提供新的治療策略。然而，也有研究發現，TLR7和TLR9是SLE患者中B細胞的效應功能的主要作用因子，而并不是TLR4和TLR2〔16-17，22〕。綜上所述，對于TLR4和TLR2在SLE的發病發展中的作用還需要進一步的證實。

在過去的十幾年中，隨著精準醫療的迅速發展，越來越多治療手段被應用到疑難疾病的診療當中，特別是與基因相關的研究中〔23〕。2021年，我國科研工作者對SLE患者的轉錄組數據進行分析，建立了一個關于SLE患者動脈粥樣硬化風險預測模型〔24〕。在本研究中，通過利用SLE患者在轉錄組數據聚集高表達分析得到的3個基因亞群，并對各亞群及對照組進行差異基因表達、信號通路以及臨床特點等分析比較，闡述了SLE可能的亞型，為SLE的精準治療提供依據。當然，本研究仍存在幾點不足之處，從GEO數據庫中獲取的數據集樣本缺乏更加詳細的臨床信息；研究中的樣本年齡偏小，還需要結合更多數量的成年SLE樣本一起分析；基因的亞型是多方面的，需要更多的樣本量以及其他組學的相互結合分析作為依據。