基于机器学习识别与振荡剪切应力相关的特征基因

Identification of Signature Gene Relate to Oscillatory Shear Stress Based on Machine Learning

摘要: 背景：动脉粥样硬化(AS)是中风的根本原因，而振荡剪切应力或低剪切应力是动脉粥样硬化的原因之一。我们研究的目的是通过机器学习筛选出受剪切应力调控的基因，并分析其生物学功能，以帮助动脉粥样硬化的诊断和治疗。方法：从基因表达综合数据库(GEO)中检索了GSE20739、GSE16706和GSE1518数据集。首先使用limma软件包筛选差异基因，然后使用蛋白质–蛋白质相互作用(PPI)网络删除一些不相互作用的基因。利用三种机器学习算法——最小绝对收缩和选择算子(LASSO)、随机森林和支持向量机(SVM)——筛选特征基因。通过受试者工作特征(ROC)曲线评估特征基因的诊断性能。结果：筛选出四个受剪切应力调控的特征基因(ADM, RAMP2, GNA14, FABP5)，这些基因与动脉粥样硬化相关。ADM和GNA14与振荡剪切应力之间存在强相关性，ADM和FABP5之间存在基因共表达。结论：我们的研究识别了四个受剪切应力调控的基因，这些基因影响动脉粥样硬化的发生和发展，这可能有助于未来动脉粥样硬化的预防和治疗。

Abstract: Background: Atherosclerosis (AS) is the root cause of stroke, while oscillatory shear stress or low shear stress is one of the reasons for atherosclerosis. The purpose of our research is to screen the genes regulated by shear stress through machine learning and analyze their biological functions to help the diagnosis and treatment of atherosclerosis. Methods: GSE20739, GSE16706 and GSE1518 datasets were retrieved from the gene expression omnibus (GEO) database. We first used the limma package to screen the differential gene, then protein-protein interaction (PPI) network was used to delete some genes that do not interact with each other. Three machine learning algorithms, Least Absolute Shrinkage and Selection Operator (LASSO), Random Forest, Support Vector Machine (SVM) were utilized for screening signature genes. Receiver operating characteristic (ROC) curves estimating the diagnostic performance of signature genes. Results: Four signature genes were screened (ADM, RAMP2, GNA14, FABP5), which regulated by shear stress. The main biological functions of characteristic genes are related to AS. There was a strong correlation between ADM and GNA14 and oscillatory shear stress. There was gene co-expression between ADM and FABP5. Conclusion: Our study identified four genes that was regulated by shear stress and affect the occurrence and development of atherosclerosis, which may be helpful for the prevention and treatment of atherosclerosis in the future.

1. 引言

动脉粥样硬化(AS)是中风的主要原因，特别是在颈动脉分叉处[1]。AS是由多种因素引起的，其中血管剪切应力在AS的发生和发展中起着重要作用。作用于血管壁的剪切应力单位通常以dyn/cm²表示，它可以通过激活血管内皮细胞(EC)上的机械受体和调节内皮细胞基因表达来调节内皮结构和功能[2]。当内皮细胞处于层流中时，通常会受到高层流剪切应力(LS)或生理剪切应力(>15 dyn/cm²)的作用，此时血管不易形成动脉粥样硬化斑块。然而，在湍流中的血管内壁(在弯曲或分叉处)会受到振荡剪切应力(OS)或低剪切应力(<4 dyn/cm²)的作用，血管容易形成动脉粥样硬化[3] [4] (图1)。

一些研究发现，在低剪切应力下，动脉内皮细胞中抗动脉粥样硬化基因的表达受到抑制，而促动脉粥样硬化基因的表达显著增加[2] [5]，从而加速了AS的形成。目前，在体外实验中已经鉴定出一些相关的内皮基因。然而，目前通过生物信息学分析筛选与低剪切应力相关的基因的研究较少，尤其是在多数据集上的研究。

在本研究中，我们选择数据集中暴露于低或振荡剪切应力的样本和暴露于生理或高剪切应力的样本作为数据处理对象。我们通过多种方式减少差异基因(DEGs)的数量，最终筛选出一些可疑的特征基因。基于本研究，我们可以鉴定出一些受剪切应力调控并能影响AS的基因，同时也为心脑血管疾病的治疗提供新的治疗靶点。

Figure 1. Characteristic of carotid blood flow. The blood flow near the center of the vessel is laminar, the vessel wall is subjected to laminar shear stress or high shear stress (>15 dyn/cm²), while the lateral part of the vessel is turbulent, and the vessel wall is subjected to oscillatory shear stress or low shear stress (<4 dyn/cm²), which makes the lateral wall form atherosclerotic plaques easily

图1. 颈动脉血流的特征。靠近血管中心的血流是层流，血管壁受到层流剪切应力或高剪切应力(>15 dyn/cm²)的作用。而血管侧部的血流是湍流，血管壁受到振荡剪切应力或低剪切应力(<4 dyn/cm²)，这使得侧壁易形成动脉粥样硬化斑块

2. 材料与方法

2.1. 剪切应力数据集和预处理

从基因表达综合数据库(GEO数据库，http://www.ncbi.nlm.nih.gov/geo/)，下载了三个独立的剪切应力数据集，包括GSE20739、GSE16706和GSE1518 [6]。

在GSE20739数据集中，我们保留了所有样本，由于GSE16706和GSE1518中的一些样本不符合纳入标准，因此将其排除，包括GSM418610、GSM418611、GSM418614、GSM418615、GSM418618、GSM418619、GSM26287、GSM26288、GSM26289和GSM26290。最终包含了11个OS样本和11个LS样本(表1)。我们使用SVA包去除了数据集的批次效应。

Table 1. Summary of the shear stress datasets involved in our study

表1. 研究中涉及的切应力数据集的总结

GSE accession	Source types	Platform	Analysis type	Total number	OS	LS
GSE20739	HUVECs	GPL6947	Array	12	6	6
GSE16706	HUVECs	GPL6480	Array	6	3	3
GSE1518	HUVECs	GPL96	Array	4	2	2

OS：振荡剪切应力；LS：层流剪切应力；HUVECs：人脐静脉内皮细胞。

2.2. 差异基因的鉴定

我们使用Limma包筛选振荡剪切应力和层流剪切应力之间的差异基因(DEGs) [7]。DEGs的筛选标准如下：P值 < 0.05和绝对log2值的倍数变化(FC) > 1。火山图由R软件(版本4.3.1)生成。删除未上调和下调的基因后，通过tidyverse包筛选出25个上调基因和25个下调基因，并使用pheatmap包生成热图。通过tidyr、ggplot2和ggpubr包生成了50个基因的分类箱图。

2.3. 功能和通路富集分析

对50个DEG进行基因本体(GO)分析[8]，包括生物过程(BP)、细胞成分(CC)和分子功能(MF)。此外，使用京都基因与基因组百科全书(KEGG)数据库分析DEG的信号通路[9]。在R软件中使用Clusterprofiler包进行GO和KEGG分析，当P < 0.05时被认为具有统计显著性。

2.4. PPI网络的建立和分析

STRING (http://string-db.org/)用于建立PPI交互网络[10]，Cytoscape用于显示蛋白质–蛋白质相互作用的结果[11]，Cytoscape中的cytoHubba包用于分析PPI网络，分析结果生成的Degree值用于优化网络的可视化结果。

2.5. 特征基因的鉴定

我们使用三种机器学习算法——最小绝对收缩和选择算子(LASSO) [12]、随机森林[13]和支持向量机(SVM) [14]来筛选中心基因。最终我们确定了三种机器学习算法的交集，并通过ROC曲线验证这些基因的预测能力，使用R包“pROC”计算曲线下面积(AUC)。

2.6. 特征基因的功能和通路富集分析

通过GO分析4个中心基因，包括生物过程(BP)、细胞成分(CC)和分子功能(MF)。此外，使用KEGG数据库分析中心基因的信号通路。在R软件中使用Clusterprofiler包进行GO和KEGG分析，当P < 0.05时被认为具有统计显著性。

2.7. 列线图的建立

使用rms包将中心基因纳入列线图模型。ROC曲线和校准曲线用于评估列线图的准确性。

2.8. 基因共表达

我们使用future.apply包分析中心基因的共表达，当P < 0.05时被认为具有统计显著性。

2.9. 统计分析

本研究中的所有统计分析均在R软件(版本4.3.1)中实现。Wilcoxon用于分析两组之间的差异。除非另有说明，否则^*P < 0.05，^**P < 0.01和^***P < 0.001被认为具有统计显著性。

3. 结果

3.1. 鉴定OS和LS之间的差异基因(DEGs)

首先，通过使用sva包消除了三个数据集之间的批次效应(图2(A))。其次，使用limma包对数据集进行分析，并设置logFC为1，从而筛选出42个下调基因和52个上调基因，并绘制了火山图(图2(B))。热图展示了上调幅度最大的25个基因和下调幅度最大的25个基因(图2(C))。分类箱图显示OS组和LS组之间的50个差异基因存在显著差异(图2(D))。

Figure 2. (A) Principal component analysis (PCA) results and box plot for the combined expression profile before Combat and after Combat. (B) Volcanic map of differential genes (DEGs) distribution, blue nodes represent down-regulated genes, red nodes represent up-regulated genes, gray nodes represent genes with no significant change. (C) Heat map of DEGs, horizontal axis represents each sample, vertical axis represents each gene, blue represents low gene expression. Red represents high gene expression. (D) Box plot shows the difference in the expression of each DEG in laminar shear stress (LS) and oscillatory shear stress (OS). ^*P < 0.05, ^**P < 0.01, ^***P < 0.001

图2. (A) 主成分分析(PCA)结果和批次效应处理前后合并表达谱的箱线图。(B) 差异基因(DEGs)分布的火山图，蓝色节点代表下调基因，红色节点代表上调基因，灰色节点代表无显著变化的基因。(C) 差异基因的热图，横轴代表每个样本，纵轴代表每个基因，蓝色代表低基因表达，红色代表高基因表达。(D) 箱线图显示层流剪切应力(LS)和振荡剪切应力(OS)下每个DEG表达的差异。^*P < 0.05，^**P < 0.01，^***P < 0.001

3.2. 功能和通路富集分析

GO分析包括三个类别(图3(A))，即生物过程(BP)、细胞成分(CC)和分子功能(MF)。BP分析显示，细胞核分裂、细胞器裂变和有丝分裂核分裂显著富集。在CC分析中，含胶原蛋白的细胞外基质、纺锤体和有丝分裂纺锤体中区占据了前三位。此外，受体配体活性、激素活性以及糖胺聚糖结合在MF中起重要作用。KEGG分析显示，前三位富集的通路主要是矿物质吸收、白细胞跨内皮迁移以及血管平滑肌收缩(图3(B))。

Figure 3. Function enrichment analysis. (A) Gene ontology (GO) functional analysis of the differential genes (DEGs) between laminar shear stress (LS) and oscillatory shear stress (OS). (B) Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis of the DEGs between OS and LS

图3. 功能富集分析。(A) 层流剪切应力(LS)和振荡剪切应力(OS)之间差异基因(DEGs)的基因本体(GO)功能分析。(B) OS和LS之间差异基因的京都基因与基因组百科全书(KEGG)通路分析

3.3. PPI网络构建

Figure 4. Protein-protein interaction (PPI) networks, according to Cytohubba analysis, color and figure size represent Degree value

图4. 蛋白质–蛋白质相互作用(PPI)网络，根据Cytohubba分析，颜色和图形大小代表度值

我们将50个差异基因(DEGs)上传到STRING在线数据库，并将包含多个基因相互作用的TSV文件下载到Cytoscape软件，筛选出37个DEGs。然后，我们通过Cytohubba分析网络，排除了度数较小且位于网络外部的DEGs。最终，剩下32个DEGs。我们构建了这32个DEGs的PPI网络(图4)。

3.4. 通过LASSO、随机森林和SVM筛选核心基因

通过三种机器学习算法，筛选出受OS调控的特征基因。LASSO分析筛选出7个特征基因(图5(A)，图5(B))，基于SVM算法筛选出32个特征基因(图5(C))，随机森林分析筛选出9个特征基因(图5(D))。三种算法的结果交集，最终确定了四个特征基因(图5(E))，包括ADM、RAMP2、GNA14和FABP5。然后，通过ROC曲线评估算法的准确性(图5(F))。

Figure 5. Selection of signature genes regulated by shear stress in the combined GSE20739, GSE16706 and GSE1518 datasets. (A) Five-time cross-verification for tuning parameter selection in the least absolute shrinkage and selection operator (LASSO) model. Each curve corresponds to a single gene. (B) LASSO coefficient profiling. The solid vertical lines represent the partial likelihood deviance standard error (SE). The dotted vertical line is drawn at the optimal lambda. (C) Random forest shows the number of selected genes and accuracy. (D) Support vector machine-recursive feature elimination (SVM-RFE) algorithm for feature selection. (E) Venn diagram showing the characteristic genes shared by LASSO, Random Forest, and SVM-RFE algorithms. (F) The receiver operating characteristic (ROC) curves estimating the diagnostic performance of signature genes

图5. 在合并的GSE20739、GSE16706和GSE1518数据集中筛选剪切应力调控的特征基因。(A) 最小绝对收缩和选择算子(LASSO)模型中调节参数选择的五次交叉验证。每条曲线对应一个基因。(B) LASSO系数分析图。实线垂直线表示部分似然偏差标准误(SE)。虚线垂直线表示最佳lambda。(C) 随机森林显示所选基因的数量和准确性。(D) 支持向量机–递归特征消除(SVM-RFE)算法的特征选择。(E) 维恩图显示LASSO、随机森林和SVM-RFE算法共享的特征基因。(F) 受试者工作特征(ROC)曲线评估特征基因的诊断性能

3.5. 特征基因的功能富集分析

BP分析显示，肾上腺髓质素受体信号通路、胰淀素受体信号通路和降钙素家族受体信号通路显著富集。在CC分析中，胰淀素受体复合物2、肾上腺髓质素受体复合物和胰淀素受体复合物占据了前三位。此外，肾上腺髓质素受体结合、肾上腺髓质素结合以及肾上腺髓质素受体活性在MF中起重要作用(图6(A))。KEGG分析显示，前三位富集的通路主要是血管平滑肌收缩、PPAR信号通路以及查加斯病(图6(B))。

Figure 6. Function enrichment analysis. (A) Gene ontology (GO) functional analysis of the signature genes. (B) Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis of the signature genes

图6. 功能富集分析。(A) 特征基因的基因本体(GO)功能分析。(B) 特征基因的京都基因与基因组百科全书(KEGG)通路分析

3.6. 列线图模型构建

Figure 7. (A) Nomogram shows the correlation between the signature genes and the oscillatory shear stress. (B), (C) Receiver operating characteristic (ROC) and calibration curves show that the prediction ability of the nomogram model is accurate

图7. (A) 列线图显示特征基因与振荡剪切应力之间的相关性。(B)，(C) 受试者工作特征(ROC)曲线和校准曲线显示列线图模型的预测能力准确

Figure 8. The scatterplot showed that ADM and FABP5 were co-expressed (R = 0.74, P = 0.013)

图8. 散点图显示ADM和FABP5的共表达(R = 0.74, P = 0.013)

我们基于四个特征基因构建了一个列线图(图7(A))，以估计特征基因与OS的相关性。列线图中的每个基因向上投影到一个点，两个或多个变量的分数之和被转换为与OS的相关性，其中总分越高，与OS的相关性越高。ROC和校准曲线的结果显示，列线图模型的预测能力准确(图7(B)，图7(C))。

3.7. 特征基因的共表达分析

我们分析了四个特征基因的共表达，发现这四个基因中，只有ADM和FABP5存在基因共表达，表明它们之间存在一定的相关性(图8)。

4. 讨论

颈内动脉狭窄或闭塞是中风的重要原因之一，而动脉粥样硬化(AS)是颈内动脉狭窄的根本原因。除了高血压、糖尿病等风险因素外，低剪切应力或振荡剪切应力在AS的发生和发展中也起着重要作用[15]。一些研究表明，低剪切应力甚至可以影响斑块的稳定性[16]。

在本研究中，我们首先消除了三个数据集之间的批次效应，然后对合并数据集进行降维筛选。首先，通过R软件筛选出50个差异表达基因(DEGs)。GO和KEGG的结果显示，50个DEGs主要参与与AS相关的细胞增殖和白细胞穿内皮迁移等过程。在将50个DEGs导入String和Cytoscape后，排除18个基因，剩余的32个基因分别通过SVM、LASSO和随机森林进行筛选，最终选择三种方法的交集，并对剩余的四个特征基因进行功能富集分析。发现它们主要参与PPAR信号通路、血管平滑肌收缩等与AS相关的过程。列线图显示四个特征基因与振荡剪切应力(OS)相关，ADM和FABP5具有共表达关系。

G蛋白偶联受体由α、β和γ亚基组成的异源三聚体，是哺乳动物细胞中分布最广泛的膜表面蛋白，介导多种细胞功能[17]。PECAM-1(血小板内皮细胞粘附分子-1，也称CD31)是一种定位于内皮细胞间粘附部位的细胞粘附分子，流体流动会显著诱导PECAM-1发生酪氨酸磷酸化[18]。有研究表明Gαq/11与PECAM-1可以形成机械敏感复合物[19]，且剪切应力可以诱导Gαq-PECAM-1复合物解离，其分子机制尚未探索[18]。GNA14是Gαq/11亚家族的成员。GαQ家族可以调节多种内皮细胞功能的跨膜和细胞内信号转导系统，包括VEGF、Akt1、ERK1/2、磷脂酶C等。Zou等人的研究表明，GNA14可以影响HUVECs的细胞迁移，GNA14的上调可以消除由FGF2/GFP和VEGFA/GFP刺激的细胞增殖[20]。当将振荡剪切应力应用于HUVEC时，GNA14的表达下调，无法消除由FGF2/GFP和VEGFA/GFP刺激的细胞增殖，导致内皮细胞过度增殖，同时细胞迁移也受到了抑制，最终导致内皮细胞排列紊乱。以往研究表明，生理或高剪切应力可以使内皮细胞的排列更加稳定，而振荡或低剪切应力则使内皮细胞的排列更加紊乱[21]。然后，内皮细胞排列紊乱可能更容易导致内皮细胞损伤和低密度脂蛋白的侵入，最终发展为AS。

RAMP是一种具有单跨膜结构域的膜蛋白，具有三种亚型，分别为RAMP1、RAMP2和RAMP3，其中RAMP2对血管的完整性至关重要。Shindo等人的研究表明，当在内皮细胞中敲除RAMP2基因后，内皮细胞将从基底膜上分离，促进血管平滑肌细胞(VSMC)的增殖，并导致炎症细胞的积累[22]。ADM/RAMP2系统在维持血管结构、调节血管生成和保护血管免受损伤方面起关键作用[23]。ADM的血管生成功能依赖于RAMP2。在我们的研究中，发现RAMP2对剪切应力的调节有响应，在低剪切应力下RAMP2的表达下调，并且由于ADM/RAMP2系统的血管保护功能以及ADM在体内的半衰期较短，RAMP2也许可以成为治疗AS的新靶点。

肾上腺髓质素(ADM)是一种由52个氨基酸组成的血管扩张肽激素[24]，主要由血管内皮细胞产生，其受体是降钙素受体样受体(CLR) [25] [26]。CLR对ADM的特异性由RAMP调节[23]。研究表明，ADM具有抗炎、抗凋亡和抗氧化作用[27]-[29]，但在高血压和充血性心力衰竭患者中，ADM水平升高，且升高程度与疾病严重程度成正比[30] [31]，这证明ADM在一定程度上参与了这些疾病。ADM可以上调VEGF和eNOS的表达[32] [33]。Shindo等人的研究表明，ADM(-/-)小鼠由于血管结构异常而在子宫内死亡，表明ADM对血管发育非常重要[22]。还有研究表明ADM与AS有关。在冠状循环中，ADM在动脉粥样硬化性心脏病患者中减少[34]。有趣的是，我们的研究发现，在将振荡剪切应力应用于内皮细胞后，ADM在内皮细胞中的表达显著上调，而振荡剪切应力是AS的原因之一，这与上述研究结果相反，具体原因尚不清楚，需进一步研究证实。

脂肪酸结合蛋白(FABP)可以非共价结合长链脂肪酸，是代谢活动的关键介质[35]。FABP4和FABP5在apoE(-/-)小鼠中的表达在脂肪细胞和巨噬细胞中起着重要作用，影响脂肪酸运输、炎症和AS的发生。研究表明，在apoE(-/-)小鼠中，FABP4和FABP5的表达[36]，尽管FABP5在脂肪组织中的表达仅为FABP4的百分之一，但FABP5对抗AS的作用比FABP4更大，FABP5水平是与巨噬细胞胆固醇外流相关的AS风险的潜在标志物[37]。Yogi等人发现，FABP5在早期和晚期AS的泡沫细胞中显著上调，而FABP4的上调则不显著[38]。Chenwei Yu等人发现，FABP4和FABP5都有促血管生成作用，其中FABP4局限于微血管内皮，而FABP5可以扩展到更大的血管内皮细胞。FABP5在良好条件下可以促进内皮细胞增殖和迁移，而在不利条件下，如振荡剪切应力下，则会促进内皮细胞凋亡[39]。在我们的研究中，发现内皮细胞中的FABP5可以响应剪切应力，在振荡剪切应力下上调，与其他特征基因相比，但本研究的列线图显示FABP5与振荡剪切应力之间的相关性较低。我们使用的所有数据都是体外实验数据。如果在体内血管中，振荡剪切应力首先影响内皮细胞，导致内皮细胞的FABP5上调，促进内皮细胞凋亡；同时，一些其他促炎因子上调巨噬细胞的FABP5，加剧AS的进程。

本研究的优势在于通过机器学习筛选出响应OS的标志基因。以往对剪切应力的分析仅针对动物内皮细胞，数据较少。在我们的研究中，选择HUVECs作为分析对象，数据集数量为三个，是关于剪切应力研究中数据集最多的研究，包括22个样本。我们分析了每个标志基因与OS之间的相关性以及标志基因在体内的不同生物功能。所有特征基因都与AS的发生和发展有关，但ADM的结果与以往研究相反，需要进一步研究验证。由于数据都是体外实验数据，样本仅为HUVECs，因此无法分析OS作用下的免疫细胞浸润。

5. 结论

在本研究中，通过机器学习筛选出了由振荡剪切应力(OS)调控的四个特征基因。其中，GNA14和ADM与OS有很强的相关性，ADM和FABP5有共表达关系。所有筛选出的基因均与动脉粥样硬化的发生和发展有关。本研究的结果可能对未来剪切应力和动脉粥样硬化机制的研究有所帮助。

生成式人工智能在科学写作中的声明

在本文撰写过程中未使用人工智能进行辅助。

NOTES

^*通讯作者。