冠狀動脈粥樣硬化性心臟病患者心外膜脂肪組織的生物信息學研究

柴晏，趙玉青，郭旭男，王東英，邊云飛

動脈粥樣硬化（atherosclerosis，AS）導致的心血管疾病（cardiovascular disease，CVD）是全球死亡的主要原因，世界衛生組織預測2030年將有2 360萬人死于CVD［1］。AS形成機制復雜，主要涉及血管內皮損傷、單核巨噬細胞黏附、脂質沉積等［2］，但這些因素如何影響AS進程，其機制尚未完全清楚。因此找尋AS新的診斷靶點，預防AS的發生，延緩其進展尤為重要。利用基因芯片技術對臨床患者標本進行檢測分析，篩選出一些具有價值的基因，進而深入研究冠狀動脈粥樣硬化性心臟病（coronary artery disease，CAD）的發病機制具有重要的臨床意義。肥胖是CVD的獨立危險因素，但肥胖的主要評價指標體質指數（BMI）與AS患者死亡率呈“U”型關聯，這種現象被稱為肥胖悖論，因此脂肪組織不再被認為是單純的“能量倉庫”，而是一種代謝活躍的內分泌和旁分泌器官［3］。心外膜脂肪組織（epicardial adipose tissue，EAT）在解剖上與冠狀動脈緊密相連，同AS、心房顫動、心力衰竭等CVD密切相關［4-6］。多項研究發現功能失調的EAT通過分泌外泌體（exosome，EXO）和生物活性物質促進AS疾病進展［7］，但其作用機制仍需進一步研究。本研究通過對美國國家生物技術信息中心（NCBI）的基因表達綜合數據庫（Gene Expression Omnibus database，GEO）中 CAD患者EAT和EXO的基因測序數據進行生物信息學分析，探討EAT參與AS的分子機制及其與EXO的聯系，旨在為CAD的診斷和治療提供新的理論基礎。

1 材料與方法

1.1 數據集的下載與預處理 以“epicardial adipose tissue OR EAT”為關鍵詞，選擇“Homo sapiens”為研究對象檢索GEO數據庫，根據數據集提供的相關信息，獲得所需芯片GSE64554、GSE120774，根據臨床信息將EAT的測序數據分為CAD組和健康對照組。下載series matrix的數據，采用“sva”包“ComBat”方法去除批次效應，合并基因表達矩陣。通過exoRBase 2.0數據庫獲取CAD組與健康對照組血液EXO基因表達譜。

1.2 差異基因的篩選 通過“limma”R語言包［8］對合并后數據集中CAD組與健康對照組EAT間差異表達基因（differentially expressed genes，DEGs）和CAD組與健康對照組EXO間DEGs進行篩選，以P＜0.05為篩選條件，并將結果可視化為火山圖和熱圖。

1.3 加權基因共表達網絡（WGCNA）的構建與CAD相關模塊的識別 為了防止合并數據集對基因關聯度的影響，使用“WGCNA”R語言包構建GSE64554數據集基因共表達網絡并識別出模塊。選擇滿足無尺度網絡的標準軟閾值β構建基因共表達網絡，轉化拓撲重疊矩陣（topological overlap matrix，TOM），采用動態剪枝法對模塊進行初步劃分。根據GSE64554數據集樣本臨床信息，決定納入的臨床信息包括疾病和年齡，得到相關性最高的模塊進行后續分析。計算基因與模塊的相關系數（module membership，MM）來評估基因與模塊之間的相關性，以|MM|＞0.8篩選模塊內樞紐基因（hub gene）。

1.4 基因功能的富集分析 通過“clusterProfiler”R語言包［9］將所獲基因進行富集分析，以FDR＜0.05為篩選標準，進行通路富集分析。通過GO富集分析，得到所選基因所富集的生物過程（biological process，BP）、細胞組分（cellular components，CC）及分子功能（molecular function，MF）。通過KEGG通路富集分析，得到所選基因所富集的信號通路。將所獲基因導入Metascape在線數據庫［10］進行基因的富集分析以及轉錄因子預測。

1.5 蛋白質-蛋白質相互作用（protein-protein interaction，PPI）網絡的構建與hub 基因篩選 將所選基因導入STRING數據庫（https://string-db.org）構建PPI網絡，將PPI網絡導入Cytoscape（Version 3.9.0）軟件［11］通過CytoHubba插件MCC算法獲得連接度最高的10個基因作為關鍵基因，并進行繪圖。

1.6 免疫細胞浸潤分析 在Cibersort網站（https://cibersort.stanford.edu/）下載22種免疫細胞的基因表達矩陣，采用Cibersort反卷積算法，計算GSE64554芯片中46個樣本的免疫細胞浸潤情況，并作圖對組織的免疫細胞浸潤情況進行比較。

1.7 臨床樣本的采集及實時熒光定量PCR（quantitative real-time PCR，qRT-PCR） 選取山西醫科大學第二醫院2021年6—10月行冠狀動脈造影患者20例，其中CAD患者及健康者各10例，留取外周血儲存在-80 ℃冰箱中備用。引物由武漢賽維爾生物科技有限公司合成，見表1。按照總RNA快速提取試劑盒（RP4001，百泰克公司）的說明提取血漿中的總RNA，cDNA的合成按照HiScriptⅢ RT SuperMix for qPCR（+gDNA wiper）試劑盒的步驟進行，第一步去除總RNA中殘余的DNA，第二步反轉錄cDNA，得到的cDNA保存于-20 ℃冰箱。以cDNA為模板，進行qRT-PCR，反應程序為：95 ℃3 min；95 ℃ 5 s；60 ℃ 10 s，72 ℃ 30 s，40 個循環；95 ℃15 s，65 ℃ 5 s，95 ℃ 50 s。結果用2-ΔΔCT進行計算分析。

表1 PCR引物Table 1 Primer sequences for qRT-PCR

1.8 統計學方法 采用Excel對錄入整理數據，應用SPSS 26.0軟件對數據進行統計學分析，符合正態分布的連續變量采用（±s）表示，兩組間比較采用獨立樣本t檢驗，以P＜0.05為差異有統計學意義。

2 結果

2.1 CAD組與健康對照組EAT間的DEGs

2.1.1 DEGs的獲取 采用“limma”包獲取CAD組與健康對照組EAT間DEGs，以P＜0.05作為篩選差異基因的閾值，得到差異表達的1 511個mRNA，其中表達上調956個，表達下調555個，見圖1、2。

圖1 CAD組與健康對照組EAT間DEGs火山圖Figure 1 Volcano plot of differentially expressed genes in epicardial adipose tissue between CAD group and healthy control group

圖2 CAD組與健康對照組EAT間DEGs熱圖Figure 2 Heatmap of differentially expressed genes in epicardial adipose tissue between CAD group and healthy control group

2.1.2 EAT間DEGs的GO/KEGG富集分析與PPI網絡對合并后CAD組較健康對照組EAT間DEGs進行GO富集分析，共獲得符合篩選條件的GO術語共800條，其中BP 567條，主要涉及細胞壓力反應、細胞DNA損傷刺激反應等。獲得CC 127條，主要涉及細胞質基質、MHC蛋白復合物、內質網腔側的組成部分等。獲得MF 106條，主要涉及MHCⅡ類蛋白復合物、特殊核糖核酸酶激活等，見圖3。共富集到KEGG通路20個，主要涉及RNA降解、病毒性心肌炎、移植物抗宿主病、1型糖尿病等信號通路，見圖4。

圖3 EAT間DEGs的GO富集分析Figure 3 GO enrichment analysis of differentially expressed genes in epicardial adipose tissue

圖4 EAT間DEGs的KEGG富集分析Figure 4 KEGG enrichment analysis of differentially expressed genes in epicardial adipose tissue

將所獲EAT間DEGs輸入STRING數據庫，構建PPI網絡，其中實際上有互作關系的節點有1 507個，有1 258條邊，每個節點的平均Degree為1.67分。將PPI網絡圖導入Cytoscape軟件，使用CytoHubba插件根據MCC算法得到連接度前10的基因并作圖，分別為 RPS27A、PSME4、PSMA1、PSMA5、GLI3、GLI1、PSME2、PSMB9、SUFU和SHH，見圖5。

圖5 EAT間DEGs的PPI網絡與關鍵基因Figure 5 Differentially expressed genes in epicardial adipose tissue and the key genes identified in the PPI network

2.1.3 EAT間DEGs的Metascape富集分析 對合并后CAD組與健康對照組EAT間DEGs進行Metascape富集分析，發現DEGs主要富集于嗅覺轉導、細胞對DNA損傷刺激反應、單純皰疹病毒1型感染、RNA代謝、調節細胞對壓力的反應、適應性免疫系統等，TRRUST數據庫預測MHCⅡ反式激活蛋白（CIITA）轉錄因子可能參與了EAT間DEGs的調控，見圖6、7。

圖6 EAT間DEGs的Metascape富集分析Figure 6 Metascape enrichment analysis of differentially expressed genes in epicardial adipose tissue

圖7 EAT間DEGs的TRRUST數據庫預測轉錄因子Figure 7 TRRUST database-based analysis of differentially expressed genes in epicardial adipose tissue to predict transcription factors

2.2 EAT組織的WGCNA分析

2.2.1 EAT組織WGCNA軟閾值的選擇與模塊的初步劃分 本研究最終以β=5作為軟閾值來構建共表達網絡，設置每個基因模塊中基因數目最小為30，設定基因的聚類高度上限為0.995。最終得到156個網絡模塊，模塊中的基因個數為30～2 169。根據拓撲重疊程度，制作模塊相關性圖，見圖8、9。

圖8 軟閾值參數網絡拓撲結構的分析Figure 8 Network topology for soft thresholding powers

圖9 基因聚類樹狀圖Figure 9 Dendrogram of gene clustering

2.2.2 臨床信息關聯與hub基因的獲取 將各個劃分出的模塊與臨床信息進行關聯分析，得到與AS患者的EAT樣本臨床信息最為相關的模塊進行后續研究。結果發現，AS患者EAT與green4（r=0.50，P=0.02）、slateblue（r=0.54，P=8.0×10-3）、darkseagreen4（r=0.55，P=6.4×10-3）、darkolivegreen1（r=0.55，P=7.0×10-3）、sienna4（r=0.56，P=5.6×10-3）、brown3（r=0.57，P=4.2×10-3）模塊呈正相關，與 deeppink（r=0.50，P=0.01）、lavender（r=0.60，P=2.6×10-3）、lavenderblush3（r=0.56，P=5.9×10-3）呈負相關，以|MM|＞0.8篩選所選模塊內hub基因，見圖10。

圖10 基因模塊與臨床信息關聯性Figure 10 Correlation between gene modules and clinical information

2.3 關鍵基因的獲得 將DEGs同模塊內hub基因取交集，獲得差異表達的模塊內hub基因，將其作為EAT組織參與AS病理生理過程的關鍵基因，其中上調hub基因為DDX47、FEM1C、NOL11、SRP54、ABI1、PATL1、BNIP2，下調hub基因為C1orf159、CHCHD4，見圖11。

圖11 差異性表達基因與模塊內hub基因Venn圖Figure 11 Venn diagram of differentially expressed genes and hub genes in the module

2.4 免疫細胞浸潤分析 通過Cibersort反卷積算法，計算GSE64554芯片中CAD組和健康對照組EAT的免疫細胞浸潤情況。研究發現，CAD組EAT組織較健康對照組中初始CD4+T細胞表達豐度升高（P=0.048），靜息樹突細胞表達豐度減低（P=0.044），見圖12、13。

圖12 GSE64554數據集EAT組織免疫細胞浸潤柱狀圖Figure 12 Bar graph of immune cell infiltration in epicardial adipose tissue between CAD patients and controls included in dataset GSE64554

圖13 GSE64554數據集EAT組織免疫細胞浸潤小提琴圖Figure 13 Violin plot of immune cell infiltration in epicardial adipose tissue between CAD patients and controls included in dataset GSE64554

2.5 CAD組與健康對照組EXO間DEGs 采用“limma”包選獲取CAD組與健康對照組EXO間DEGs，以P＜0.05，|logFoldchange|＞1作為篩選差異基因的閾值，共得到DEGs 1 658個，其中上調278個，下調1 380個，見圖14、15。

圖14 CAD組與健康對照組EXO間火山圖Figure 14 Volcano plot of differentially expressed genes in exosomes between CAD patients and heathy controls

圖15 CAD組與健康對照組EXO間DEGs熱圖Figure 15 Heatmap of differentially expressed genes in exosomes between CAD patients and heathy controls

2.6 EAT來源EXO 將CAD組與健康對照組EAT、EXO間DEGs取交集并做Venn圖，共獲得129個基因，其中上調16個，下調113個，經篩選，將BPI、BIRC5、CXCL12、RNASE1、F2R作為關鍵基因并通過qRT-PCR進行驗證，見圖16。

圖16 CAD患者EAT組織、EXO間DEGs Venn圖Figure 16 Venn diagram of differentially expressed genes in epicardial adipose tissue and exosomes from CAD patients

2.7 EAT來源EXO間DEGs的驗證 收集CAD組和健康對照組外周血，運用qRT-PCR驗證生物信息學分析所得的EAT來源EXO基因BPI、BIRC5、CXCL12、RNASE1、F2R的mRNA水平。結果顯示，與健康對照組相比，CAD患者的5個所選基因的表達水平差異有統計學意義（P＜0.05），其中BPI、BIRC5、CXCL12和RNASE1的mRNA水平升高，F2R基因的mRNA水平下降（P＜0.05），見表2。

表2 CAD組與健康對照組EXO差異性基因表達（±s）Table 2 Differentially expressed genes in exosomes between CAD patients and healthy controls

表2 CAD組與健康對照組EXO差異性基因表達（±s）Table 2 Differentially expressed genes in exosomes between CAD patients and healthy controls

注：CAD=冠狀動脈粥樣硬化性心臟病；BPI、BIRC5、CXCL12、RNASE1、F2R為心外膜、EXO間DEGs

組別 樣本數 BPI BIRC5 CXCL12 RNASE1 F2R CAD 組 10 2.14±0.32 1.98±0.23 2.91±0.89 4.47±1.56 0.45±0.06健康對照組 10 0.99±0.17 1.02±0.11 1.14±0.44 1.28±0.56 1.00±0.12 t值 -9.945 -12.070 -5.630 -6.001 12.637 P值 0.011 0.010 ＜0.001 0.009 ＜0.001

3 討論

EAT作為一種特殊的脂肪組織，解剖結構上緊貼冠狀動脈，在病理環境中其向炎癥表型轉化，通過分泌促炎細胞因子和EXO參與血管內皮損傷、血管重塑、免疫細胞黏附等AS病理生理過程。EAT的炎性反應是AS的全新診斷治療靶點［5］，OIKONOMOU等［12］通過影像組學技術將心外膜脂肪衰減系數（fat attenuation index，FAI）應用于CAD患者的早期篩查，而目前已廣泛應用于臨床的鈉-葡萄糖協同轉運蛋白2（sodiumdependent glucose transporters 2，SGLT-2）抑制劑達格列凈則可以通過減輕EAT炎癥而減少糖尿病患者心血管事件的發生率［13］，但EAT在CAD中的作用機制仍需進一步研究。

此次通過生物信息學分析的方法從數據集GSE64554和GSE120774中篩選出CAD組和健康對照組EAT間DEGs 1 511個，通過富集分析發現DEGs主要富集于細胞質基質、MHC蛋白復合物、MHCⅡ蛋白復合物、RNA降解、抗原加工和呈遞等。目前AS被認為是一種慢性炎癥性疾病，MCHⅡ作為T細胞的呈遞抗原，通過抗原識別、細胞活化、產生細胞因子等方式在炎癥過程和CAD進展中發揮作用［14］。通過Meatscape富集分析中的TRRUST數據庫，筆者預測到CIITA轉錄因子可能在EAT中發揮了重要的調節作用。研究發現，CIITA通過其轉錄后修飾甲基化［15］、去乙酰化［16-17］調節MHCⅡ的轉錄而參與AS。通過免疫細胞浸潤分析本研究發現CAD患者EAT中幼稚CD4+T細胞豐度升高，早在2015年IBORRA等［18］便指出可以通過觀察幼稚CD4+T細胞的體外分化程度評估AS程度，GADDIS等［19］指出AS通過影響糖酵解而損害幼稚CD4+T細胞功能。綜上所述，MHCⅡ在AS中發揮了重要作用，但其在CAD患者EAT中的作用和調控機制仍需進一步研究。

通過篩選CAD患者與健康對照組EAT間DEGs，構建EAT組織WGCNA，獲得了DDX47、FEM1C、NOL11、SRP54、ABI1、PATL1、BNIP2、C1orf159、CHCHD4、RPS27A共10個基因作為EAT參與CAD病理生理過程的關鍵基因。目前僅有關鍵基因被發現在CVD中發揮作用，XIA等［20］通過生物信息學技術發現免疫相關基因RPS27A可能參與AS的病理過程，HUANG等［21］指出RPS27A基因在糖尿病視網膜病變的關鍵基因且同幼稚CD4+T細胞密切相關。LI等［22］發現MBNL1通過對ABI1的剪切調控AS中平滑肌細胞的表型轉換。因此，關鍵基因仍需進一步進行細胞和動物實驗證明其在EAT中的差異表達以及相關功能。

內皮細胞、平滑肌細胞、心肌細胞、脂肪細胞和血小板來源的EXO通過轉運蛋白質、RNA、DNA和其他生物活性物質在缺血再灌注損傷、血管鈣化、AS和心臟重塑中發揮了重要作用，已成為CAD的全新診斷治療靶點［23］。脂肪組織是循環中EXO的來源之一［24］，冠狀動脈血管周圍脂肪組織是一種特殊的EAT［25］，LIU等［26］發現血管周圍脂肪組織來源EXO中miR-382-5p通過ABCA1/ABCG1信號通路調節泡沫巨噬細胞形成。LI等［27］發現血管周圍脂肪組織來源EXO中miR-221-3p通過介導血管重塑參與AS。ZHAO等［7］發現芒果苷誘導血管周圍脂肪組織衍生的EXO可以改善AS中的內皮功能障礙，因此EAT來源EXO在AS中的作用具有重要的研究價值。

為進一步探討EAT在CAD中的作用機制，本研究將CAD患者與健康對照組的EAT間DEGs、EXO間DEGs取交集，共獲得129個基因，并篩選表達程度較高、有意義的基因RNASE1、BPI、BIRC5、CXCL12、F2R。其中前4個基因為CAD患者中上調基因。RNASE1是一種可分泌的耐熱蛋白，具有強大的心臟保護功能［28］，其通過保護內皮細胞在炎癥中的損傷而成為血管穩態的關鍵調控者［29］，在CAD中，EAT可能通過EXO釋放RNASE1進而調節冠狀動脈內皮細胞的炎性損傷。血管生成趨化因子介導細胞趨化、白細胞脫粒和血管生成，CXCL-12同AS密切相關［30］，其通過CXCL-12-CXCR4軸參與平滑肌細胞表型轉化保護血管穩態［31］，特異性過表達巨噬細胞IGF-1通過減少CXCL-12介導的單核細胞募集和增加ABCA1依賴性巨噬細胞脂質外流增加AS斑塊的穩定性［32］，但目前尚無關于EXO內CXCL-12的相關研究。BIRC5是凋亡抑制蛋白（IAP）家族的成員，LI等［33］通過對小鼠頸動脈單細胞測序發現BIRC5可能與血流紊亂所致血管巨噬細胞聚集相關。早在2002年便有報道指出BPI的單核苷酸多態性同心肌梗死密切相關［34］，BLEIJERVELD等［35］通過循環中粒細胞的深度蛋白質組分析發現了BPI可以作為嚴重AS性冠狀動脈狹窄的生物標志物，這表明BPI可以作為晚期CAD乃至預測急性心血管事件的標志物。F2R是血管壁炎癥和血栓形成的關鍵調節劑，其多態性同氯吡格雷療效密切相關［36］，但目前少有F2R在AS中的研究。

綜上所述，本研究共確定了DDX47、FEM1C、NOL11、SRP54、ABI1、PATL1、BNIP2、C1orf159、CHCHD4、RPS27A共10個關鍵基因，RNASE1、BPI、BIRC5、CXCL12、F2R共5個可能源自EAT EXO的基因進行了驗證，證明其在CAD中的診療意義，此外還探討了CAD患者EAT免疫浸潤情況，為CAD的基礎研究提供理論支持。本研究有其局限性，本研究主要針對mRNA，但蛋白質才是生命活動的執行者，mRNA層面的改變并不能等同蛋白質層面的改變，因此本研究所得結論仍需動物、細胞實驗進一步驗證。

作者貢獻：柴晏、趙玉青、郭旭男負責課題設計、文章撰寫；柴晏、郭旭男負責數據庫檢索及數據整理；趙玉青、王東英負責實驗實施及實驗數據處理；邊云飛負責文章的質量控制及審校。

本文無利益沖突。