基于单细胞RNA测序技术筛选颈动脉粥样硬化中巨噬细胞特征基因的研究

王建茹，李晓辉

0 引 言

动脉粥样硬化(atherosclerosis，AS)是一种由免疫系统介导的以动脉壁炎症和斑块形成为特征的慢性动脉疾病，是多种心脑血管疾病的共同病理基础[1]。颈动脉粥样硬化(carotid atherosclerosis，CAS)是AS的一种临床类型，因其解剖位置表浅，探查方便，常被作为全身AS的“窗口”，来判定AS的整体情况[2]。CAS发病机制与炎症免疫反应、脂质代谢及非编码RNA等关系密切，但其机制仍未完全阐释清楚。众所周知，巨噬细胞作为免疫细胞中重要的成员之一，在AS病理过程中的致病作用一直是研究的热点和重点[3]。因此，深入探索巨噬细胞在CAS中的病理机制，挖掘关键基因，对于更好地理解CAS的发病机制以及防治靶点具有积极的作用。

单细胞RNA测序(single-cell RNAsequencing，scRNA-Seq)作为新一代的高通量测序技术，可进行单细胞的转录组分析，识别新的或罕见的细胞亚型，为深入探析不同细胞类型中单细胞的行为、机制及其与机体的关系提供了新的方法和思路，在心血管疾病中方面产生了革命性的影响[4-8]。近年来，已有多项研究利用scRNA-Seq技术在深化理解AS病理机制，促进潜在药物靶点研发和临床诊断标志物的筛选等方面做出了重要贡献[9-15]。

现阶段，高通量技术联合生物信息学的方法已被广泛用于CAS的研究中，为其病理机制、发展转归和诊疗等提供了新见解。本研究结合scRNA-Seq技术和生物信息学分析探索了人CAS组织中巨噬细胞的特征基因，以期更好的理解CAS的病理过程及巨噬细胞介导CAS的潜在作用机制，为其诊断、治疗提供新的参考。

1 材料与方法

1.1 数据来源及预处理从GEO数据库(http://www.ncbi. nlm.nih.gov/geo/)中下载CAS的scRNA-Seq数据集GSE159677。该数据集的样本来自3例CAS患者，包括3个动脉粥样硬化斑块(scRNA-Seq-AS组)和3个匹配的颈动脉近端相邻组织(scRNA-Seq-对照组)。利用hdf5r包和Seurat包分别对每个样本的数据进行读取，然后再分别对两组的样本数据进行合并后保存。

同时，下载基因芯片数据集GSE43292及其平台文件GPL6244-17930。该数据集的组织来自32例CAS患者，包括32个动脉粥样硬化斑块组织(AS组)和32个匹配的远端完整的颈动脉组织(对照组)。依据平台文件对表达数据进行注释，将探针矩阵转化为基因矩阵，剔除基因表达量低的基因，基因对应多个探针则利用avereps函数取均值。然后，再利用normalizeBetweenArrays函数对数据进行矫正，并保存以进行后续的特征基因验证和免疫细胞浸润分析。

1.2scRNA-Seq数据的质控和降维聚类利用R语言读取GSE159677中的scRNA-Seq-AS组和scRNA-Seq-对照组样本数据，并对重复基因取均值。利用Seurat包对样本数据进行过滤。然后，使用NormalizeData函数对数据进行标准化，再利用FindVariableFeatures函数计算细胞间基因的变异系数(标准差)。使用RunPCA函数对scRNA-Seq-AS组的样本数据进行主成分分析( principal component analysis，PCA)，选取P< 0.05的主成分，使用RunTSNE函数进行t分布随机邻接嵌入(t-distributed stochastic neighbor embedding，t-SNE )聚类。

1.3细胞类型注释对scRNA-Seq-AS组的样本数据降维聚类后，使用Seurat包中的FindAllMarker函数找出每个细胞亚群(Cluster)相对其它Cluster差异表达的基因，然后利用SingleR包对不同的细胞亚群进行注释。

1.4差异基因筛选利用FindAllMarker函数，以|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)用于后续分析。

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

1.6蛋白互作网络的构建及特征hub基因的筛选将共同DEGs上传至STRING数据库(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数据集中提取1.6中筛选出来的各样本巨噬细胞特征hub基因的表达量，绘制箱体图对特征基因的组间表达差异情况进行可视化，并利用pROC包绘制特征基因的ROC曲线，计算AUC值。

1.8免疫细胞浸润分析CIBERSORT和ssGSEA是广泛应用的量化免疫细胞浸润状态的方法[19]。利用这两种方法对预处理后GSE43292的表达谱矩阵进行模拟分析，进而获得所有样本的免疫细胞的浸润模式，利用Wilcoxon检验以P< 0.05为阈值筛选两组间的差异免疫细胞；然后，提取AS样本中巨噬细胞含量的相对比例。

1.9特征hub基因与AS中巨噬细胞浸润状态的相关性分析采用Pearson法对GSE43292数据集中特征hub基因的表达量与1.8中获取的AS样本巨噬细胞含量的相对比例进行相关性分析，并利用cor.test函数计算相关性系数。相关性系数大于0为正相关，相关性系数小于0为负相关，相关系数的绝对值代表强、弱或无相关，以P≤0.05为具有统计学意义。

2 结 果

2.1 scRNA-Seq数据质控结果预处理GSE159677数据集中的样本后，scRNA-Seq-对照组和scRNA-Seq-AS组样本中共获得2756个细胞，scRNA-Seq-AS组样本中获得2095个细胞。利用Seurat包计算所有样本中细胞的基因数(nFeature)、转录本测序Count数(nCount)、线粒体基因百分比(percent.mt)，过滤掉每个细胞测到的基因数目小于50，线粒体基因组占比大于5%的细胞。对筛选细胞间标准化方差较大的1500个基因进行后续分析，这些基因代表了细胞间的异质性。见图1。

a、b、c、d、：分别为AS样本中细胞的基因数、转录本测序Count数、所有细胞的线粒体比例、细胞间标准方差较大的前1500个变异基因；e、f、g、h：分别为scRNA-Seq-对照样本和scRNA-Seq-AS样本中细胞的基因数、转录本测序Count数、所有细胞的线粒体比例、细胞间标准方差较大的前1500个变异基因

2.2降维聚类及注释结果利用PCA法对过滤后的样本进行了线性降维，发现有20个主成分的P<0.05。然后，选取这20个主成分进行t-SNE聚类分析，结果发现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信号通路、自噬、补体与凝血级联反应。

图 3 共同DEGs的GO富集分析

2.4特征hub基因的筛选及验证经STRING数据库对共同DEGs进行PPI分析后，利用Cytoscape软件构建了PPI网络。该网络包括84条边和38个节点，其中上调基因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则相反。随后，该研究利用芯片数据集GSE43292验证了5个特征hub基因在AS病理过程中的表达情况，结果与 scRNA-Seq数据分析的结果一致，见图6；同时，5个特征hub基因表现出了良好的诊断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数据集中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)均表现出了明显的强正相关性。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)均表现出了明显的强正相关性。

3 讨 论

CAS作为AS最常见的临床类型，可导致脑梗、心梗等不良事件[20-21]。研究显示，巨噬细胞作为AS中一种丰富的功能异质性免疫细胞，在CAS的病理过程中也同样发挥着重要作用[20]。巨噬细胞通过调控脂质代谢、炎症反应、LC3相关吞噬作用、免疫代谢、胞葬等多种途径在AS中发挥重要作用，被认为是治疗AS的一个有潜力的作用靶点[1,3,22-24]。因此，深入挖掘巨噬细胞在CAS中潜在的分子机制和临床诊断标志物，对于有效安全的管理CAS具有重要的意义。

本研究对数据集GSE159677进行分析，发现CAS样本中巨噬细胞亚群有47个DEGs。富集分析显示，这些DEGs可通过调控中性粒细胞的活化及相关免疫应答、胶原蛋白分解及代谢、铁结合及铁氧化酶活性、细胞凋亡、铁死亡、PPAR通路等生物学功能和通路，来介导CAS的病理过程。既往研究报道，中性粒细胞明胶酶相关脂质沉积蛋白高表达于活化的中性粒细胞，可诱导巨噬细胞分泌促炎介质加剧症状性CAS局部及全身性的炎症反应[25]。在CAS斑块组织中，浸润的巨噬细胞通过促进胶原蛋白分解、铁死亡等，加速斑块的不稳定性[26-27]。

本文对共同DEGs进行PPI分析，筛选出5个巨噬细胞特征hub基因即CLU、CTSD、CTSB、CTSZ、CTSL。同时，利用芯片数据集对hub基因进行分析，结果发现hub基因在芯片数据和scRNA-Seq数据中的差异表达趋势一致，并在诊断CAS方面表现出了良好的诊断效能。最后，进行免疫细胞浸润分析，发现巨噬细胞在CAS组织中处于高浸润状态。相关性分析显示，在CAS环境中CLU的表达水平与巨噬细胞的相对含量存在一定的负相关性；而CTSD、CTSB、CTSL、CTSZ分别与巨噬细胞具有明显的强正相关性。CLU又称载脂蛋白J，通过调节巨噬细胞参与CAS病程变化的机制可能包含两个方面，一是抑制巨噬细胞浸润和促炎性M1巨噬细胞的极化[28]；二是促进巨噬细胞源泡沫细胞的胆固醇逆流[29]。巨噬细胞分泌的CTSB参与了CAS斑块的进展和破裂[30]。CTSL的表达随着CAS的严重程度而增加，并与斑块不稳定和帽状结构破裂有关，而巨噬细胞是易损和破裂的颈动脉斑块中表达CTSL的主导细胞[31]。也就是说，CTSL可能通过调控巨噬细胞凋亡，参与坏死核的形成和不稳定斑块的发展。研究显示，CTSD可以促进巨噬细胞泡沫化进而加重AS病变[32]。现阶段，尚未有CTSD、CTSZ调控巨噬细胞的功能介导CAS病理过程的报道。综上所述，这5个hub基因的高表达在CAS中发挥着重要作用，将有助于阐释CAS可能的分子机制，并为CAS的诊断和防治提供一定的思路。

本研究尚有一定的局限性：①由于GEO等公共数据库里人CAS样本较少，未能对多数据集、多批次的数据进行整合分析，导致纳入分析的样本有限。②本研究虽然利用外部基因芯片数据集对挖掘出的5个巨噬细胞特征基因进行了验证，但结果仍需在后续的基础和临床研究中进行佐证。

总之，本研究利用生物信息分析和scRNA-Seq数据相结合的方法从免疫学的角度探索了CAS可能的病理机制，挖掘出5个巨噬细胞特征基因(即CLU、CTSD、CTSB、CTSL、CTSZ)。本研究的核心意义不仅有助于更好地了解CAS的病理机制，为后续以巨噬细胞为切入点来探索CAS的免疫学机制提供了方向，还为CAS潜在临床诊断标志物的开发提供了一定的借鉴。