基于單細胞RNA測序技術篩選頸動脈粥樣硬化中巨噬細胞特征基因的研究

王建茹，李曉輝

0 引 言

動脈粥樣硬化(atherosclerosis，AS)是一種由免疫系統(tǒng)介導的以動脈壁炎癥和斑塊形成為特征的慢性動脈疾病，是多種心腦血管疾病的共同病理基礎[1]。頸動脈粥樣硬化(carotid atherosclerosis，CAS)是AS的一種臨床類型，因其解剖位置表淺，探查方便，常被作為全身AS的“窗口”，來判定AS的整體情況[2]。CAS發(fā)病機制與炎癥免疫反應、脂質代謝及非編碼RNA等關系密切，但其機制仍未完全闡釋清楚。眾所周知，巨噬細胞作為免疫細胞中重要的成員之一，在AS病理過程中的致病作用一直是研究的熱點和重點[3]。因此，深入探索巨噬細胞在CAS中的病理機制，挖掘關鍵基因，對于更好地理解CAS的發(fā)病機制以及防治靶點具有積極的作用。

單細胞RNA測序(single-cell RNAsequencing，scRNA-Seq)作為新一代的高通量測序技術，可進行單細胞的轉錄組分析，識別新的或罕見的細胞亞型，為深入探析不同細胞類型中單細胞的行為、機制及其與機體的關系提供了新的方法和思路，在心血管疾病中方面產生了革命性的影響[4-8]。近年來，已有多項研究利用scRNA-Seq技術在深化理解AS病理機制，促進潛在藥物靶點研發(fā)和臨床診斷標志物的篩選等方面做出了重要貢獻[9-15]。

現(xiàn)階段，高通量技術聯(lián)合生物信息學的方法已被廣泛用于CAS的研究中，為其病理機制、發(fā)展轉歸和診療等提供了新見解。本研究結合scRNA-Seq技術和生物信息學分析探索了人CAS組織中巨噬細胞的特征基因，以期更好的理解CAS的病理過程及巨噬細胞介導CAS的潛在作用機制，為其診斷、治療提供新的參考。

1 材料與方法

1.1 數(shù)據來源及預處理從GEO數(shù)據庫(http://www.ncbi. nlm.nih.gov/geo/)中下載CAS的scRNA-Seq數(shù)據集GSE159677。該數(shù)據集的樣本來自3例CAS患者，包括3個動脈粥樣硬化斑塊(scRNA-Seq-AS組)和3個匹配的頸動脈近端相鄰組織(scRNA-Seq-對照組)。利用hdf5r包和Seurat包分別對每個樣本的數(shù)據進行讀取，然后再分別對兩組的樣本數(shù)據進行合并后保存。

同時，下載基因芯片數(shù)據集GSE43292及其平臺文件GPL6244-17930。該數(shù)據集的組織來自32例CAS患者，包括32個動脈粥樣硬化斑塊組織(AS組)和32個匹配的遠端完整的頸動脈組織(對照組)。依據平臺文件對表達數(shù)據進行注釋，將探針矩陣轉化為基因矩陣，剔除基因表達量低的基因，基因對應多個探針則利用avereps函數(shù)取均值。然后，再利用normalizeBetweenArrays函數(shù)對數(shù)據進行矯正，并保存以進行后續(xù)的特征基因驗證和免疫細胞浸潤分析。

1.2scRNA-Seq數(shù)據的質控和降維聚類利用R語言讀取GSE159677中的scRNA-Seq-AS組和scRNA-Seq-對照組樣本數(shù)據，并對重復基因取均值。利用Seurat包對樣本數(shù)據進行過濾。然后，使用NormalizeData函數(shù)對數(shù)據進行標準化，再利用FindVariableFeatures函數(shù)計算細胞間基因的變異系數(shù)(標準差)。使用RunPCA函數(shù)對scRNA-Seq-AS組的樣本數(shù)據進行主成分分析( principal component analysis，PCA)，選取P< 0.05的主成分，使用RunTSNE函數(shù)進行t分布隨機鄰接嵌入(t-distributed stochastic neighbor embedding，t-SNE )聚類。

1.3細胞類型注釋對scRNA-Seq-AS組的樣本數(shù)據降維聚類后，使用Seurat包中的FindAllMarker函數(shù)找出每個細胞亞群(Cluster)相對其它Cluster差異表達的基因，然后利用SingleR包對不同的細胞亞群進行注釋。

1.4差異基因篩選利用FindAllMarker函數(shù)，以|log2(fold change, FC)|>1和矯正后的P<0.05為閾值，對scRNA-Seq-AS組和scRNA-Seq-對照組間的DEGs，以及scRNA-Seq-AS組不同Cluster間的DEGs進行篩選[16-17]。然后，將組間的DEGs和scRNA-Seq-AS組巨噬細胞的特異性DEGs取交集(即共同DEGs)用于后續(xù)分析。

1.5基因功能富集分析利用clusterProfiler包對共同DEGs進行GO和KEGG通路富集分析，并利用ggplot2、GOplot包進行可視化。

1.6蛋白互作網絡的構建及特征hub基因的篩選將共同DEGs上傳至STRING數(shù)據庫(https://string-db.org/)，選擇“Multiple proteins”模式，互作得分設置為中等置信度≥0.400，進行PPI分析，并將用Cytoscape3.7.2軟件構建PPI網絡。利用CytoHubba插件的MCC、DMNC、MNC、Degree、EPC算法來篩選PPI網絡中的巨噬細胞特征hub基因(即5種算法中Top10基因的交集基因)[17-18]。

1.7巨噬細胞特征hub基因的驗證從預處理后的GSE43292數(shù)據集中提取1.6中篩選出來的各樣本巨噬細胞特征hub基因的表達量，繪制箱體圖對特征基因的組間表達差異情況進行可視化，并利用pROC包繪制特征基因的ROC曲線，計算AUC值。

1.8免疫細胞浸潤分析CIBERSORT和ssGSEA是廣泛應用的量化免疫細胞浸潤狀態(tài)的方法[19]。利用這兩種方法對預處理后GSE43292的表達譜矩陣進行模擬分析，進而獲得所有樣本的免疫細胞的浸潤模式，利用Wilcoxon檢驗以P< 0.05為閾值篩選兩組間的差異免疫細胞；然后，提取AS樣本中巨噬細胞含量的相對比例。

1.9特征hub基因與AS中巨噬細胞浸潤狀態(tài)的相關性分析采用Pearson法對GSE43292數(shù)據集中特征hub基因的表達量與1.8中獲取的AS樣本巨噬細胞含量的相對比例進行相關性分析，并利用cor.test函數(shù)計算相關性系數(shù)。相關性系數(shù)大于0為正相關，相關性系數(shù)小于0為負相關，相關系數(shù)的絕對值代表強、弱或無相關，以P≤0.05為具有統(tǒng)計學意義。

2 結 果

2.1 scRNA-Seq數(shù)據質控結果預處理GSE159677數(shù)據集中的樣本后，scRNA-Seq-對照組和scRNA-Seq-AS組樣本中共獲得2756個細胞，scRNA-Seq-AS組樣本中獲得2095個細胞。利用Seurat包計算所有樣本中細胞的基因數(shù)(nFeature)、轉錄本測序Count數(shù)(nCount)、線粒體基因百分比(percent.mt)，過濾掉每個細胞測到的基因數(shù)目小于50，線粒體基因組占比大于5%的細胞。對篩選細胞間標準化方差較大的1500個基因進行后續(xù)分析，這些基因代表了細胞間的異質性。見圖1。

a、b、c、d、：分別為AS樣本中細胞的基因數(shù)、轉錄本測序Count數(shù)、所有細胞的線粒體比例、細胞間標準方差較大的前1500個變異基因；e、f、g、h：分別為scRNA-Seq-對照樣本和scRNA-Seq-AS樣本中細胞的基因數(shù)、轉錄本測序Count數(shù)、所有細胞的線粒體比例、細胞間標準方差較大的前1500個變異基因

2.2降維聚類及注釋結果利用PCA法對過濾后的樣本進行了線性降維，發(fā)現(xiàn)有20個主成分的P<0.05。然后，選取這20個主成分進行t-SNE聚類分析，結果發(fā)現(xiàn)1927個細胞被分為了12個亞群，不同細胞亞群被標記為了相應的顏色；SingleR包注釋后12個細胞亞群被分為了6個亞群，包括內皮細胞(3個細胞亞群)、巨噬細胞(3個細胞亞群)、單核細胞(3個細胞亞群)、T細胞(1個細胞亞群)、平滑肌細胞(1個細胞亞群)、組織干細胞(1個細胞亞群)。見圖2。

a：每個主成分的P值；b：各細胞亞群的聚類分布情況；c：注釋后各細胞亞群分布的t-SNE圖

2.3共同DEGs的篩選及其富集分析差異分析的結果顯示，巨噬細胞亞群和其它細胞亞群之間的有260個DEGs，其中上調88個，下調172個；組間有109個DEGs，其中上調35個，下調74個；兩者取交集后共獲得47個DEGs，其中上調27個，下調20個。利用clusterProfiler包、org.Hs.eg.db包等將篩選出的共同DEGs進行功能富集分析和可視化，見圖3、圖4。GO功能富集分析結果顯示，依據P<0.05確定了434個GO條目；生物學過程(biological process，BP)條目為355條，主要涉及中性粒細胞激活、中性粒細胞脫顆粒、膠原蛋白分解代謝過程等；細胞組分條目(cellular component，CC)為43條，主要涉及含膠原蛋白的細胞外基質、溶酶體腔等；分子功能(molecular function，MF)條目為36條，主要涉及鐵結合、鐵氧化酶活性、氧化還原酶活性等。KEGG通路富集分析結果顯示，依據P<0.05共映射出了6條信號通路，包括溶酶體、鐵死亡、凋亡、PPAR信號通路、自噬、補體與凝血級聯(lián)反應。

圖 3 共同DEGs的GO富集分析

2.4特征hub基因的篩選及驗證經STRING數(shù)據庫對共同DEGs進行PPI分析后，利用Cytoscape軟件構建了PPI網絡。該網絡包括84條邊和38個節(jié)點，其中上調基因23個，下調基因15個。對cytoHubba插件中MNC、DMNC、MCC、EPC、Degree算法排序前10的基因取交集共獲得5個交集基因，即凝聚素(CLU)、組織蛋白酶D(CTSD)、組織蛋白酶B(CTSB)、組織蛋白酶Z(CTSZ)、組織蛋白酶L(CTSL)，其中CLU為下調基因，其余為上調基因。如圖5所示，與其他細胞亞群相比，CLU在巨噬細胞的表達和分布最小，而CTSD、CTSB、CTSZ、CTSL則相反。隨后，該研究利用芯片數(shù)據集GSE43292驗證了5個特征hub基因在AS病理過程中的表達情況，結果與 scRNA-Seq數(shù)據分析的結果一致，見圖6；同時，5個特征hub基因表現(xiàn)出了良好的診斷AS的效能。見圖7。

圖 4 共同DEGs的KEGG通路富集分析

a：各細胞亞群中5個特征hub基因表達的氣泡圖；b、c、d、e、f：分別為CLU、CTSD、CTSB、CTSL、CTSZ在每個細胞分布的t-SNE圖

*P<0.001

2.5免疫細胞浸潤分析CIBERSORT結果所示，對照和AS組織樣本之間存在8種差異的免疫細胞，其中M0巨噬細胞在AS中的比例顯著增加。ssGSEA結果所示，對照和AS組織樣本之間存在27種差異的免疫細胞，其中巨噬細胞在AS中的比例顯著增加。兩種算法的結果均提示，巨噬細胞在AS組織中的浸潤比例較對照樣本高。見圖8。

圖 7 特征hub基因在GSE43292數(shù)據集中ROC分析的結果

a：CIBERSORT算法分析的組織樣本中22種免疫細胞浸潤差異分析的小提琴圖；b：ssGSEA算法分析的組織樣本中29種免疫細胞浸潤差異分析的小提琴圖

2.6特征hub基因與巨噬細胞的相關性分析CIBERSORT算法結果顯示，巨噬細胞含量的相對比例與CLU(r=-0.39,P=0.029)呈負相關性；與CTSD(r=0.83,P<0.001)、CTSB(r=0.76,P<0.001)、CTSL(r=0.85,P<0.001)、CTSZ(r=0.82,P<0.001)均表現(xiàn)出了明顯的強正相關性。ssGSEA算法結果顯示，巨噬細胞含量的相對比例與CLU(r=-0.51,P=0.002)呈一定的負相關性；與CTSD(r=0.87,P<0.001)、CTSB(r=0.91,P<0.001)、CTSL(r=0.86,P<0.001)、CTSZ(r=0.88,P<0.001)均表現(xiàn)出了明顯的強正相關性。

3 討 論

CAS作為AS最常見的臨床類型，可導致腦梗、心梗等不良事件[20-21]。研究顯示，巨噬細胞作為AS中一種豐富的功能異質性免疫細胞，在CAS的病理過程中也同樣發(fā)揮著重要作用[20]。巨噬細胞通過調控脂質代謝、炎癥反應、LC3相關吞噬作用、免疫代謝、胞葬等多種途徑在AS中發(fā)揮重要作用，被認為是治療AS的一個有潛力的作用靶點[1,3,22-24]。因此，深入挖掘巨噬細胞在CAS中潛在的分子機制和臨床診斷標志物，對于有效安全的管理CAS具有重要的意義。

本研究對數(shù)據集GSE159677進行分析，發(fā)現(xiàn)CAS樣本中巨噬細胞亞群有47個DEGs。富集分析顯示，這些DEGs可通過調控中性粒細胞的活化及相關免疫應答、膠原蛋白分解及代謝、鐵結合及鐵氧化酶活性、細胞凋亡、鐵死亡、PPAR通路等生物學功能和通路，來介導CAS的病理過程。既往研究報道，中性粒細胞明膠酶相關脂質沉積蛋白高表達于活化的中性粒細胞，可誘導巨噬細胞分泌促炎介質加劇癥狀性CAS局部及全身性的炎癥反應[25]。在CAS斑塊組織中，浸潤的巨噬細胞通過促進膠原蛋白分解、鐵死亡等，加速斑塊的不穩(wěn)定性[26-27]。

本文對共同DEGs進行PPI分析，篩選出5個巨噬細胞特征hub基因即CLU、CTSD、CTSB、CTSZ、CTSL。同時，利用芯片數(shù)據集對hub基因進行分析，結果發(fā)現(xiàn)hub基因在芯片數(shù)據和scRNA-Seq數(shù)據中的差異表達趨勢一致，并在診斷CAS方面表現(xiàn)出了良好的診斷效能。最后，進行免疫細胞浸潤分析，發(fā)現(xiàn)巨噬細胞在CAS組織中處于高浸潤狀態(tài)。相關性分析顯示，在CAS環(huán)境中CLU的表達水平與巨噬細胞的相對含量存在一定的負相關性；而CTSD、CTSB、CTSL、CTSZ分別與巨噬細胞具有明顯的強正相關性。CLU又稱載脂蛋白J，通過調節(jié)巨噬細胞參與CAS病程變化的機制可能包含兩個方面，一是抑制巨噬細胞浸潤和促炎性M1巨噬細胞的極化[28]；二是促進巨噬細胞源泡沫細胞的膽固醇逆流[29]。巨噬細胞分泌的CTSB參與了CAS斑塊的進展和破裂[30]。CTSL的表達隨著CAS的嚴重程度而增加，并與斑塊不穩(wěn)定和帽狀結構破裂有關，而巨噬細胞是易損和破裂的頸動脈斑塊中表達CTSL的主導細胞[31]。也就是說，CTSL可能通過調控巨噬細胞凋亡，參與壞死核的形成和不穩(wěn)定斑塊的發(fā)展。研究顯示，CTSD可以促進巨噬細胞泡沫化進而加重AS病變[32]?，F(xiàn)階段，尚未有CTSD、CTSZ調控巨噬細胞的功能介導CAS病理過程的報道。綜上所述，這5個hub基因的高表達在CAS中發(fā)揮著重要作用，將有助于闡釋CAS可能的分子機制，并為CAS的診斷和防治提供一定的思路。

本研究尚有一定的局限性：①由于GEO等公共數(shù)據庫里人CAS樣本較少，未能對多數(shù)據集、多批次的數(shù)據進行整合分析，導致納入分析的樣本有限。②本研究雖然利用外部基因芯片數(shù)據集對挖掘出的5個巨噬細胞特征基因進行了驗證，但結果仍需在后續(xù)的基礎和臨床研究中進行佐證。

總之，本研究利用生物信息分析和scRNA-Seq數(shù)據相結合的方法從免疫學的角度探索了CAS可能的病理機制，挖掘出5個巨噬細胞特征基因(即CLU、CTSD、CTSB、CTSL、CTSZ)。本研究的核心意義不僅有助于更好地了解CAS的病理機制，為后續(xù)以巨噬細胞為切入點來探索CAS的免疫學機制提供了方向，還為CAS潛在臨床診斷標志物的開發(fā)提供了一定的借鑒。