基于GEO数据库筛选动脉粥样硬化发生的关键差异circRNA

李 晶,张 芸＊，席建军，程文俊，鸿嘎鲁

（内蒙古医科大学附属医院全科医学科，内蒙古 呼和浩特 010050）

1 前言

心血管疾病(cardiovascular disease，CVDs)是导致全球发病和死亡的重要原因，每年影响数百万人。动脉粥样硬化（atherosclerosis，AS）是一种涉及不同细胞类型、多种细胞因子和粘附分子的慢性炎症，不稳定的动脉粥样硬化斑块破裂、血管狭窄或血小板聚集和血栓形成引起的闭塞可导致急性心血管疾病[1，2]。氧化低密度脂蛋白（OxLDL）具有高度的致AS作用，动脉粥样硬化的发病机制与巨噬细胞摄取低密度脂蛋白后转化为泡沫细胞密切相关[3]。因此，寻找与AS发展相关的生物标志物可以为临床研究提供理论依据。

环状RNAs(circRNAs)是一类广泛存在于各种生物细胞中具有调控基因表达功能的长链内源性非编码RNA，具有结构稳定和组织特异表达等特征[4]，它们通常参与包括肿瘤在内的各种疾病的发生和发展。大量的研究数据表明，circRNAs在多种肿瘤组织和癌旁正常组织中的表达存在显著性的差异[5~9]。一些特异性表达的circRNAs可能成为肿瘤等多种疾病诊断和预后的新的生物标志物。Chen L等人发现circRNA_100290可能作为一种竞争性的内源性RNA（ceRNA），通过吸收miR-29b家族成员来调控CDK6在口腔鳞状细胞癌中的表达[10]。尽管circRNAs在各种疾病中发挥作用，但是与AS相关的circRNAs的表达和功能鲜有研究，进一步研究有望用于预防或治疗AS。

在过去的几十年中，微阵列技术和生物信息学分析被广泛应用于差异基因基因筛选。在我们的研究中，利用基因表达综合数据库（gene expression omnibus，GEO）下载并分析微阵列数据集GSE107522，以获得经OxLDL处理得巨噬细胞THP-1实验组和THP-1对照组之间差异表达的circRNA(DEGs)。选择hsa_circ_0003645作为研究基因，并通过CSCD数据库预测hsa_circ_0003645 open reading frame（ORF）、RNA binding protein(RBP)和microRNA response element(MRE)。随后，对预测的靶基因进行Gene Ontology(GO)和Kyoto Encyclopedia of Genes and Genomes(KEGG)的通路富集分析，接着构建PPI网络，筛选关键基因，以帮助我们了解AS发展的分子机制，分析与关键基因相互作用的小分子药物。本文探讨了circRNAs可能成为AS的潜在标志物，为预防和治疗该疾病提供了理论基础。

2 材料与方法

2.1 GEO数据下载和预处理

在GEO数 据 库（https：//www.ncbi.nlm.nih.gov/geo/）中下载GSE107522数据集（3个巨噬细胞对照样本，3个氧化型低密度脂蛋白(OxLDL)处理的巨噬细胞样本）和GPL19978（Agilent-069978 Arraystar Human CircRNA microarray V1）平台注释文件，通过编写perl脚本用平台注释文件对GSE107522表达矩阵的探针进行注释获得circRNA表达矩阵，接着使用AgilentcircRNA注释文件将circRNA的ID环化为circBase格式的ID。

2.2 差异分析

R软 件（R version 3.6.1：https：//www.r-project.org/）加载limm包，以logFC>1或者logFC<1，P<0.05为筛选条件筛选差异表达的circRNA（DECs），接着利用pheatmap包对筛选出的DECs进行热图绘制。在circBase（http：//www.circbase.org/）数据库中检索DECs相关信息。接着在exorbase（http：//www.exorbase.org/）数据库中检索DECs在外泌体中的表达水平。

2.3 CircRNAs的ORF、RBP和MER预测

通过CSCD[11]数据库（Cancer-Specific CricRNA Database：http：//gb.whu.edu.cn/CSCD/）预测circRNAs的open reading frame（ORF）、RNA binding protein(RBP)和microRNA response element(MRE)。接着通过ORFfinder[7]数据库（https：//www.ncbi.nlm.nih.gov/orffinder/）和circbank数据库（http：//www.circbank.cn/index.html）预测circRNA的OFR以及编码的蛋白，并在UniProtKB数据库对编码的蛋白进行blast分析。最后我们在Circular RNA Interactome数据库（https：//circinteractome.nia.nih.gov/）预测circRNA的RBP，再 通 过venny2.1（https：//bioinfogp.cnb.csic.es/tools/venny/）在线网站与CSCD预测的RBP取交集。在starbase数据库（http：//starbase.sysu.edu.cn/）的通路模块下对两个数据库共有的RBP进行KEGG富集分析。

2.4 CircRNA MER的靶基因预测及GO和KEGG pathway分析

下载miRDB、miRTarbase和TargetScan三个数据的miRNA靶基因预测文件，对通过perl脚本得到circRNAs MER在三个数据共同预测到的靶基因进行GO和KEGG的富集分析。通过org.Hs.eg.db包将差异基因的Official Symbol转化为基因ID，cluster-Profiler包进行GO和KEGG富集分析，GO plot包[12]进行聚类分析，且p值和q值都小于0.05才具有显著差异。

2.5 PPI网络构建及关键基因筛选

使用STRING数据库（https：//string-db.org/）构建circRNA MER靶基因的PPI网络，并且综合数值大于0.9的相互作用才具有统计学意义的。Cyto-Hubba插件按degree方法筛选前20个关键基因。颜色越红色代表基因degree越大。

2.6 关键基因的KEGG和GO富集分析

利用DAVID在线数据库（https：//david.ncifcrf.gov/summary.jsp）进行生物学分析，P值小于0.05被认为具有统计学意义。通过Enrichr在线数据库(http：//amp.pharm.mssm.edu/Enrichr/)中的DrugMatrix分析预测可以与关键基因相互作用的小分子化合物。

3 结果

3.1 差异表达的circRNAs筛选

我们对GSE107522进行标准化处理后，得到了9个差异表达的DECs（见图1A）和DECs绘制的热图（见图1B）。与THP-1对照组相比，经oxLDL处理得巨噬细胞THP-1中有8个circRNAs上调，1个circRNAs下调（见表1）。

为了从9个差异表达的circRNA筛选到目标circRNA进行后续分析，我们在circBase数据中检索了这9个circRNA的相关信息。我们发现hsa_circ_0003645和hsa_circ_0005699的基因名称都是C16orf62且在许多样本中都有表达（见表2）。有趣的是通过核酸序列比较发现hsa_circ_0003645的核酸序列大于hsa_circ_0005699，因此我们选择hsa_circ_0003645进行后续分析。

3.2 CircRNA在人血外泌体中的表达水平

图1 DECs的鉴定及circRNAs的表达Fig.1 Identification of DECS and expression of circRNAs

表1 circRNAs在oxLDL巨噬细胞吞噬中的差异表达Tab.1 Differential expression of circRNAs in phagocytosis of oxLDL macrophages

我们通过在exorbase（http：//www.exorbase.org/）数据库中检索hsa_circ_0003645在人血外泌体中的表达水平，发现该circRNA在外泌体中有着比较高的表达水平，表达水平排在前10%~20%（expression ranK：10%~20%），且hsa_circ_0003645在冠心病和结肠癌患者血液外泌体中表达均下调（见图1C）。这些结果表明，hsa_circ_0003645可能通过外泌体运输到胞外调控相关靶细胞，预示着巨噬细胞可能通过分泌外泌circRNA调控AS进展。

3.3 CircRNA的ORF预测

通过CSCD在线数据库预测hsa_circ_0003645的ORF、RBP和MER（见图2A）。通过ORFfinder数据库预测hsa_circ_0003645的ORF（见表3），结果显示hsa_circ_0003645的ORF2可能编码包含108个氨基酸的蛋白。对hsa_circ_0003645的ORF2相关蛋白进行同源比对，发现与转录该circRNA的VPS35L基因（C16orf62）编码的蛋白同源性最高（见图2B）。在circbank数据库进一步验证hsa_circ_0003645的蛋白编码能力，结果进一步表明该circRNA具有较高的蛋白质编码能力（见图2C）。

表2 差异表达的circRNA在circbase数据库中检索信息Tab.2 Information retrieval of differentially expressed circRNA in circbase database

图2 预测hsa_circ_0003645的蛋白编码能力Fig.2 Forecasthsa_circ_Protein coding ability of 0003645

表3 has_circ_0003645的ORF预测Tab.3 has_circ_ORF prediction of 0003645

3.4 PBP预测与KEGG富集分析

我们通过CSCD数据库和circRNAs Interactome数据库共同对hsa_circ_0003645进行RBP预测分析，结果显示hsa_circ_0003645含有eIF4AIII和TDP-43的结合位点（见表4）。为了解eIF4AIII和TDP-43的功能，我们通过starbase数据库对eIF4AIII和TDP-43两个RBP进行KEGG富集分析，发现两者都主要富集在Spliceosome、Cell_Cycle和Wnt_Signaling_Pathway等相关通路上（见表5、6）。

表4 circRNAs的RBP分析Tab.4 RBP analysis of circRNAs

表5 TDP-43在KEGG途径中的富集分析Tab.5 Enrichment analysis of TDP-43 in KEGG pathway

表6 KEGG途径中eIF4AIII靶标的富集分析Tab.6 Enrichment analysis of eIF4AIII target in KEGG pathway

3.5 CircRNA MER靶基因预测及GO和KEGG pathway分析

同时，通过CSCD数据库MER分析发现hsa_circ_0003645与56个miRNA可能发生相互作用。通过perl脚本预测hsa_circ_0003645的MER在miRDB、miRTarbase和TargetScan三个数据中的共同靶基因，发现hsa_circ_0003645的MER含有578个靶基因，接着对这些靶基因进行GO和KEGG富集分析。Go分析结果显示（见图3A）：BP主要富集在myeloid cell differentiation、response to steroid hormone、cellular response to steroid hormone stimulus。CC主要富集在cytoplasmic region、transport vesicle、secretory granule membrane。MF主要富集在DNA-binding transcription activator activity，RNA polymerase II-specific、hormone receptor binding、nuclear hormone receptor binding。并对GO富集结果进行聚类分析。KEGG主要富集在Endocytosis、Focal adhesion、JAK-STAT signaling pathway、FoxO signaling pathway。并对KEGG富集结果进行聚类分析。

图3 GO和KEGG的富集分析Fig.3 Enrichment analysis of GOand KEGG

3.6 PPI网络构建及关键基因筛选和分析

为确定目标基因之间的关键基因，hsa_circ_0003645的靶基因通过STRING在线数据库构建PPI网络（见图4A），然后通过Cytoscape软件的NetworkAnalyzer plugin进行分析靶基因（见图4B），最后使用cytoHubba plugin的degree方法分析关键基因，排在前20的基因被鉴定为关键基因（见图4C）。接着为了解这些关键基因的功能，我们通过DVID数据库对hsa_circ_0003645的20个关键基因进行GO和KEGG富集分析。结果发现hsa_circ_0003645的20个关键基因的主要富集在protein ubiquitination、regulation of cell motility、cell periphery、FoxO signaling pathway、Focal adhesion、Wnt signaling pathway（见表7）。为了寻找与20个关键基因相互作用的小分子药物，我们通过Enrichr数据库的DrugMatrix分析，发现hsa_circ_0003645与新霉素、青霉素、氯雷他定三种小分子化合物相互作用（见表8）。这些结果表明这些小分子药物可能具有治疗动脉粥样硬化的潜在价值。

图4 PPI网络和筛选关键基因Fig.4 PPInetwork and screening of key genes

表7 关键基因的GO和KEGG通路富集分析Tab.7 Enrichment analysis of key genes through GO and KEGG pathways

4 讨论

心血管疾病(CVD)是当代社会普遍存在的问题。AS是导致人群心血管疾病(CVD)高死亡率的主要原因[13]。Liang Zong等人发现circRNA_102231在肺癌中表达上调，circRNA_102231沉默显著降低了细胞迁移和侵袭的能力，抑制了细胞增殖[14]。然而circRNAs在AS中的作用还有待进一步探讨。

CircRNAs是一类在真核转录组中高度表达的非编码RNA，可以作为miRNA海绵，从而降低其靶向mRNA的能力，还参与组织和器官的发育，并在各种疾病发生过程中发挥作用[15]。Yan Zhang等人研究发现circRNA KIAA1586作为ceRNA吸收了三种miRNA(hsa-miR-29b、hsa-miR-101、hsa-miR-15a)，其失调可能导致阿尔茨海默症相关生物学功能的异常[16]。Shen L等人研究发现circRNA 0044073在AS中表达上调，通过靶向miR107和激活JAK/STAT信号通路促进细胞增殖和侵袭，可能为AS的治疗策略提供新靶点[17]。这些研究结果有力地表明了circRNAs在多种疾病中起着重要作用。

表8 关键基因与16种小分子化合物相互作用分析Tab.8 Analysis of interaction between key genes and 16 small molecular compounds

本文研究通过微阵列数据集GSE107522获得经OxLDL处理的巨噬细胞THP-1实验组和THP-1对照组的DECs。我们发现共有9个差异表达的circRNAs，其中有8个表达上调，1个表达下调。结果显示，单核细胞衍生的巨噬细胞对oxLDL的吞噬作用以及随后巨噬细胞向泡沫细胞的分化是AS形成的关键。Chen L等人研究结果证实了circRNA_100290可作为治疗人类口腔鳞癌的靶点[10]。因此，这些DECs有可能为AS提供潜在的生物标志物及治疗靶点。我们根据circBase数据库检索9个差异表达的circRNA筛选出靶标hsa_circ_0003645后，通过exorbase数据库分析发现，hsa_circ_0003645在外泌体中的表达水平较高，且在CHD和CRC患者血液外泌体中表达均下调。而有研究表明由砷转化的L-02细胞衍生出的外泌体circRNA_100284可通过充当microRNA-217的海绵，加速细胞周期，促进肝细胞增殖[18]。因此，我们推测巨噬细胞可能通过分泌外泌circRNA调控AS进展。同时，我们还发现hsa_circ_0003645有eIF4AIII和TDP-43的结合位点，可能与56个miRNA发生相互作用。circRNAs可能通过miRNAs调控靶基因发挥作用。He JH等人的研究结果表明circRNA-ACAP2可作为miR-21-5p的海绵，通过调控Tiam1的表达，从而影响结肠癌SW480细胞的增殖、迁移和侵袭[19]。此外，通过starbase数据库对eIF4AIII和TDP-43两个RBP进行KEGG富集分析发现eIF4AIII和TDP-43两者均主要在Spliceosome、Cell_Cycle和Wnt_Signaling_Pathway等相关通路中富集。由此，我们推测hsa_circ_0003645可能与eIF4AIII和TDP-43相互作用，从而调控动AS相关功能。我们通过perl脚本预测hsa_circ_0003645的MER，发现hsa_circ_0003645的MER含有578个靶基因。有研究表明lncRNA和circRNA可作为ceRNA参与ceRNA网络的构建，可竞争结合相同miRNA，与mRNA之间相互作用，影响靶基因的翻译或稳定性，从而实现转录水平的基因调节，进而在肿瘤发病机制中发挥作用。有研究报道称circRNA 0072088可以通过作为miR-377-5p的竞争性ceRNA上调NOVA_2发挥作用，从而促进非小细胞肺癌的增殖和迁移[20]。这与我们的研究类似，因此，我们判断hsa_circ_0003645可能通过ceRNA的作用机制从而影响动脉粥样硬化的发展进程。随后，我们通过对hsa_circ_0003645的靶基因构建PPI网络并分析鉴定出20个关键基因，并对这20个关键基因进行GO和KEGG富集分析。分析结果发现这些关键基因主要富集在protein ubiquitination、regulation of cell motility、cell periphery、FoxO signaling pathway、Focal adhesion、Wnt signaling pathway。此外，我们通过Enrichr数据库的DrugMatrix分析还发现hsa_circ_0003645与新霉素、青霉素、氯雷他定三种小分子药物相互作用，这些结果表明这些小分子药物可能具有治疗动脉粥样硬化的潜在价值。

总之，本文研究鉴定了hsa_circ_0003645及其相关的RBP和MER的分析。研究发现了20个关键基因以及与hsa_circ_0003645相互作用的3种小分子化合物。该项研究可能为AS提供新的生物标记物以及具有潜在价值的小分子药物。随着基因测序技术的升级，越来越多的基因组数据库被建立，有许多在线数据库可以用于circRNAs的研究，目前仅有少数circRNAs被证明与AS有关。然而，hsa_circ_0003645在AS形成过程中可能的分子机制仍有待进一步探讨。此外，这些关键基因的生物学功能以及这些分子化合物在AS中的可能作用还需要进一步研究。