基于生物信息学分析的肝癌circRNA-miRNA-mRNA调控轴的鉴定

摘要

近年来，环状RNA (circRNAs)被发现通过多种机制，特别是内源性竞争性RNA (ceRNA)机制，在肝细胞癌(HCC)中发挥重要的调节作用。因此，探索肝癌中的circrna具有重要意义。在本研究中，我们利用Cytoscape构建了ceRNA和存活网络。我们还利用R、Perl软件以及Gene Ontology (GO)、Kyoto Encyclopedia of Genes and Genomes (KEGG)等多个在线数据库和平台，对这些基因进行了总体存活、免疫细胞浸润、免疫检查点、通路活性、抗癌药物敏感性分析。最后进行受试者操作特征曲线(receiver operator characteristic curve, ROC)分析，确定基因的诊断价值。KEGG分析显示T细胞受体信号通路是主要富集途径。共筛选出29个与生存和预后相关的基因。结果提示ZNF544、WDR76、ACTG1、RASSF3、E2F3、ASRGL1和POGK与多水平免疫细胞浸润有关。此外，免疫检查点分析筛选出ACTG1、E2F3、RASSF3和WDR76。研究还发现，WDR76、E2F3、ASRGL1和POGK主要激活细胞周期和DNA损伤反应(DDR)途径。结果提示，对曲美替尼、瑞法美替尼(RDEA119)和塞鲁美替尼的敏感性与WDR76的表达相关。ROC分析显示，调控轴上所有基因的曲线下面积(area under the curve, AUC)均大于0.7。确定的hsa_circ_0000417/hsa_circ_0002688/hsa_circ_0001387—hsa-miR-199a-5p—WDR76调控轴可能为HCC的进展、临床诊断和治疗提供新的见解。

介绍

HCC是全球癌症相关死亡的第三大原因，发病率和死亡率高，特别是在亚洲1,2,3。它是隐匿的，容易转移，由于缺乏早期症状，经常在晚期发现。虽然目前肝癌的治疗方法，如手术切除、肝移植、靶向治疗、免疫治疗、局部治疗等技术已经相对成熟，但这些治疗方法对于晚期肝癌患者往往无效。这些治疗方法也无法延长这些患者的生存时间。因此，寻找新的、更有效的治疗靶点势在必行。非编码rna (ncRNAs)时代的到来为这一领域书写了新的篇章，circRNAs就是ncRNAs的新兴成员之一5,6。

环状rna最初是在一种植物病毒中发现的。然而，由于当时检测技术的限制，circrna主要被认为是传统线性剪接的副产品，其潜力未被发掘。多亏了测序技术，它们的潜力被发现了。环状RNA广泛存在于真核细胞中，具有没有5 '和3 '端的闭环结构，与线性RNA相比，具有对RNA酶的抗性和稳定性。环状rna由反向剪接前体mrna产生，由外显子、内含子或两者组成9,10。环状rna有几种形式，包括外显子环状rna (ecrna)、内含子环状rna (ciRNAs)和外显子内含子环状rna (eiciRNAs)。CiRNAs和eiciRNAs主要调控细胞核内的转录和剪接。此外，ecRNAs被转运出去执行特定的生物学功能，包括作为miRNA和蛋白质海绵或将核苷酸序列翻译成肽和蛋白质11,12,13。circrna也与多种疾病相关，包括心血管疾病、免疫系统疾病和肿瘤14,15,16。

已有证据表明，circRNAs和ceRNA机制参与了包括肝癌起始、肝癌干细胞的干性维持和肝癌进展在内的生物学过程。研究结果建议使用circRNA-miRNA-mRNA调控网络作为治疗肝细胞癌的生物标志物和治疗靶点17,18。近年来，生物信息学分析在科学研究领域越来越受欢迎。因此，本研究采用综合生物信息学分析方法。根据ceRNA机制，通过分析差异表达基因，构建了ceRNA网络。我们还构建了生存网络，分析了生存网络中上调mrna的免疫细胞浸润水平、免疫检查点和通路活性。本研究结果将为HCC的发病机制提供新的视角，提高HCC的诊断和治疗水平。图1给出了本研究的流程图。

图1

本研究流程图。

结果

下载数据及差异基因表达分析

根据确定的筛选条件，我们从Gene Expression Omnibus (GEO)数据库中下载了3个微阵列数据(GSE155949、GSE108724和GSE101728)和1个肝癌RNA-seq数据(TCGA)数据库中下载。数据集的细节如表1所示。差异基因表达分析揭示了36个circRNAs (DEcircRNAs)(图2A)、26个miRNAs (DEmiRNAs)(图2B)和2215个mrna (demmrnas)(图2C)的差异表达。这36个circrna被用于构建ceRNA网络。

表1 GEO和TCGA数据库的数据集信息。

图2

差异表达基因的鉴定。(A) DEcircRNAs热图，(B) DEmiRNAs热图，(C) demrna热图。

circrna和mirna的靶标预测

在线数据库共预测了2495个circRNAs靶基因，以及与DEmiRNAs交叉后的6个候选mirna(图3A)，分别是hsa-miR-375、hsa-miR-3188、hsa-miR-4270、hsa-miR-221-3p、hsa-miR-532-3p和hsa-miR-199a-5p。此外，数据库预测了6个候选mirna的2470个靶基因，与demrna重叠后鉴定出262个候选mrna(图3B)。

图3

miRNA (A)和mRNA (B)的重叠基因。

构建ceRNA网络

按照ceRNA机制，Perl筛选候选基因，获得构建ceRNA网络所需的输入文件，得到6个circrna、3个mirna和94个mrna。将输入的文件导入Cytoscape软件，构建包含103个节点和104条边的ceRNA网络(图4)。网格中节点的颜色表示基因表达，红色和蓝色分别表示高表达和低表达。不同的节点形状表示不同的基因类型:矩形代表circRNAs，三角形代表miRNAs，圆圈代表mrna。不同颜色的线表示节点之间的连接:紫色粗线表示circRNAs与miRNAs的关系，灰色细线表示miRNAs与mrna的关系。

图4

ceRNA网络。

GO和KEGG富集分析

对ceRNA网络中mrna的GO富集分析得出了28个GO条目。分析揭示了主要富集的生物过程，包括调节蛋白丝氨酸/苏氨酸激酶活性，细胞-基质粘附等(图5A,B)。KEGG富集分析表明，T细胞受体信号通路是主要富集途径(图5C,D)。

图5

ceRNA网络mrna的功能富集分析。(A,B) GO富集分析。(C,D) KEGG富集分析。

生存分析与生存网络构建

对ceRNA网络中的mrna进行生存分析，发现29个生存相关基因(图6)。使用Perl对这29个生存相关基因的上游基因进行匹配，确定它们之间的相互关系，获得构建生存网络所需的输入文件。使用Cytoscape创建了包含61个调控轴的存活网络(图7A)。存活网络中circRNAs高表达基因如表2所示。从存活网络中筛选出8个高表达mrna，分别为ZNF544、WDR76、ACTG1、RASSF3、E2F3、ASRGL1、SUCO和POGK进行后续分析。

图6

ceRNA网络中mrna的总体生存分析。

图7

生存网络构建及通路活性分析。(A)生存网络的构建。(B)存活网络中上调mrna的通路活性分析。

表2存活网络中高表达circrna的相关信息。

通路活性分析

通路活性分析表明，WDR76、E2F3、ASRGL1和POGK激活了细胞周期和DDR通路(图7B)。WDR76对细胞周期通路的激活作用最为显著。该基因还参与细胞凋亡和上皮-间质转化(EMT)途径的激活，并抑制RAS/MAPK途径。

免疫细胞浸润水平和免疫检查点分析

选取存活预后网络中高表达的8个mrna进行免疫细胞浸润分析。结果表明，ACTG1(图8A)、ASRGL1(图8B)、E2F3(图8C)、POGK(图8D)、RASSF3(图8E)、WDR76(图8F)、ZNF544(图8G)与肝脏肝细胞癌(LIHC)微环境中B细胞、CD4+ T细胞、CD8+ T细胞、中性粒细胞、巨噬细胞和树突状细胞的浸润水平存在显著相关性。但是，不能从数据库中检索SUCO。此外，免疫检查点分析显示，在LIHC中，ACTG1、E2F3、RASSF3和WDR76与程序性细胞死亡1 (PDCD1，图9A)、分化簇274 (CD274，图9B)和细胞毒性T淋巴细胞相关抗原4 (CTLA4，图9C)呈正相关。

图8

存活网络中上调mrna的免疫细胞浸润水平分析。(a) actg1， (b) asrgl1， (c) e2f3， (d) pogk， (e) rassf3， (f) wdr76， (g) znf544。

图9

免疫检查点，药物敏感性，ROC分析。ACTG1、E2F3、RASSF3和WDR76分别与PDCD1 (A)、CD274 (B)和CTLA4 (C)的免疫检查点分析。(D)药物敏感性分析。对hsa_circ_0000417 (E)、hsa_circ_0001387 (F)、hsa_circ_0002688 (G)、has-miR-199a-5p (H)和WDR76 (I)进行ROC分析。

药物敏感性及ROC分析

药敏分析显示，WDR76与曲美替尼、RDEA119、塞鲁美替尼的敏感性呈正相关(图9D)。研究结果表明，WDR76可能是存活网络中高表达基因中的枢纽基因。因此，我们从ceRNA网络中确定了hsa_circ_0000417/hsa_circ_0002688/hsa_circ_0001387-hsa_miR-199a-5p-WDR76调控轴。对鉴定出的调控轴上的基因进行ROC分析，以确定其对HCC诊断的价值。auc分别为0.738 (hsa_circ_0000417，图9E)、0.876 (hsa_circ_0001387，图9F)、0.785 (hsa_circ_0002688，图9G)、0.980 (hsa_miR-199a-5p，图9H)和0.844 (WDR76，图9I)。结果表明，这些基因可能用于HCC的诊断标记。并且，我们获得了circRNAs与miRNA之间的潜在结合位点(图10A)， hsa_circ_0000417(图10B)， hsa_circ_0002688(图10C)和hsa_circ_0001387(图10D)的结构。图10E显示了调节轴的模型。因此，本研究揭示了调节轴与LIHC患者的生存、免疫细胞浸润水平、免疫检查点、通路活性、药物敏感性和HCC诊断之间的联系。同时确定了其他可能参与调节肝细胞癌发生发展并对HCC的诊断和治疗具有参考价值的因素。

图10

circRNA-miRNA的调控轴模型及结合位点。(A) hsa_circ_0000417/0001387/0002688与has-miR-199a-5p之间的结合位点。结构为hsa_circ_0000417(B)/0002688 (C)/0001387 (D)， ORF打开读框。(E)调节轴模型。向上的红色箭头表示相应rna的表达上调，向下的箭头表示相反。

讨论

近年来，对circRNAs在肿瘤细胞中的各种机制，特别是ceRNA机制的研究不断扩大。ceRNA机制源于一种假设，即ceRNA可以竞争与靶向mirna结合，减少mirna的数量并增加下游mRNA靶标，从而调节疾病。有大量证据支持ceRNA机制在包括hcc21,22,23在内的多种肿瘤的发生和发展中发挥作用。本研究基于ceRNA机制和生物信息学分析，构建了ceRNA网络。多基因功能分析和鉴定数据库最终确定了HCC的关键调控轴。我们的研究结果将为HCC的发病机制、诊断和治疗提供有价值的参考。

已经确定了调节肿瘤发展的各种生物过程和信号通路。细胞基质粘附在包括浸润性乳腺癌在内的不同肿瘤细胞的侵袭性调节中起着至关重要的作用24,25。此外，T细胞受体信号通路是T细胞部分激活介导免疫应答的细胞内通路。研究发现分化146簇(CD146)联合淋巴细胞激酶(LCK)可促进小鼠T细胞受体信号的启动和抗肿瘤免疫应答26。几乎所有肿瘤的细胞周期通路都是异常的，这可能是肿瘤细胞无限复制的关键原因之一。研究还表明，细胞周期通路的遗传变异可能影响III期和iv期非小细胞肺癌(NSCLC)患者的生存27。细胞周期途径也可能调节前列腺癌的发生，这表明它们有可能作为乳腺癌治疗的靶点。DNA损伤反应是由于各种原因引起的DNA损伤后细胞进行修复的过程。受损的修复过程可能导致突变传播基因组的不稳定，并有利于肿瘤细胞的进展。一项研究发现apurinic apyrimidinic endodeoxyribonuclease 2 (APE2)可以调节DDR通路，维持胰腺癌基因组的完整性30。在胃癌中，DDR可以预测肿瘤进展和临床预后31。这些结果表明，我们研究中发现的ceRNA网络可能促进HCC的发生和HCC细胞的侵袭进展，从而恶化HCC患者的临床预后。生存分析的结果也支持这种临床预后价值。

mirna是一种短的ncrna，可以作为表观遗传调节剂促进或抑制恶性肿瘤(如HCC)的肿瘤发生和进展。在本研究中，我们通过生物信息学分析鉴定了hsa-miR-199a-5p。研究发现，靶向kelch样家族成员3 (KLHL3)的hsa-miR-199a-5p可能参与了NSCLC32的发生发展。发现长链ncRNA (lncRNA) AB209371的高表达有利于肝癌细胞的EMT，也可能是hsa-miR-199a-5p失活的直接原因。沉默AB209371联合过表达hsa-miR-199a-5p可抑制肝癌转移和肝癌细胞的EMT，提示hsa-miR-199a-5p可能作为抑制HCC转移的靶点33。这一发现与我们以前的研究结果一致。免疫细胞浸润是肿瘤微环境的重要组成部分。研究表明，B细胞、CD4+ T细胞、CD8+ T细胞、中性粒细胞、巨噬细胞、树突状细胞等浸润肿瘤的免疫细胞可影响肿瘤患者的免疫治疗效果和预后34,35。Gao等人发现SORT1基因与免疫细胞浸润之间存在关联，提示其有可能作为HCC预测的新生物标志物。干扰SORT1可以显著抑制HCC细胞的生长，提示开发抗癌策略的潜在新靶点，特别是HCC36。一般来说，免疫检查点作为免疫系统的调节器，维持和调节自身免疫和外周组织的免疫过程。然而，肿瘤细胞利用这种机制来逃避免疫37。研究发现多功能适配体(P1/C4-bi-apt)可以阻断CTLA4/B7和PDCD1/CD274信号通路，增强对HCC的免疫应答。该适体在HCC免疫治疗中的潜在应用此前已被探索过38。免疫检查点CD274上调可增加细胞分裂周期相关蛋白2 (CDCA2)的表达，预测HCC39不良预后。我们发现ceRNA网络中的ACTG1、E2F3、RASSF3和WDR76与免疫检查点密切相关，提示ceRNA网络调节HCC细胞的免疫逃逸。同时，药物敏感性分析显示WDR76与曲美替尼、RDEA119、塞鲁美替尼等抗癌药物的敏感性存在相关性。先前的研究表明，利法美替尼联合索拉非尼对失去手术机会的ras突变型HCC患者更有益40。然而，塞鲁美替尼联合索拉非尼治疗HCC的价值与其他单药治疗并无显著差异。已发现曲美替尼在治疗卵巢癌和高、低级别胶质瘤方面有显著效果42,43。然而，这种对HCC的影响尚未见报道。对调节轴上所有基因的ROC分析显示hsa_circ_0000417、hsa_circ_0002688、hsa_circ_0001387、hsa_miR-199a-5p和WDR76对HCC的潜在诊断价值。然而，我们的研究也存在局限性。本研究使用公共数据库、在线平台和其他生物信息学分析来确定特定临床诊断和治疗的价值。本研究缺乏细胞和动物实验来支持研究结果，应纳入未来的研究。

总之，通过分析差异表达的circrna、mirna和mrna，我们成功构建了ceRNA和存活网络。结果显示hsa_circ_0000417/hsa_circ_0002688/hsa_circ_0001387-hsa-miR-199a-5p-WDR76调节轴与HCC患者的生存预后、肿瘤浸润免疫细胞、免疫逃逸、通路活性和药物敏感性相关。研究结果表明，这种调节轴可以作为HCC临床诊断和治疗的新的潜在生物标志物和靶点。这些结果也为HCC的发生和发展机制提供了见解。

材料与方法

数据下载

在GEO数据库(http://www.ncbi.nlm.nih.gov/geo/)中以“肝癌”为关键词搜索相应的circRNA、miRNA和mRNA微阵列数据集。筛选和选择数据集的标准包括包含人类肝癌及其邻近正常肝组织的序列样本。标准组和肿瘤组的样本量均大于6例。据此筛选GSE155949 (circRNAs表达谱芯片)、GSE108724 (miRNAs表达谱芯片)、GSE101728 (mrna表达谱芯片)。筛选探针矩阵和平台文件进行进一步分析。从TCGA数据库(https://portal.gdc.Cancer.gov/)下载HCC表达和临床数据，用于后续的生存分析。

差异表达基因的分析

使用Perl (version 5.32.1001, https://strawberryperl.com/)对下载的芯片数据进行整理，得到相应的基因表达矩阵文件。根据探针矩阵文件划分肿瘤组和正常组。采用R软件(版本4.1.2,https://www.r-project.org/)和Limma软件包对数据进行分析，以调整后的P < 0.05为过滤阈值。最后，测定HCC组和正常组的DEcircRNAs、DEmiRNAs和demmrnas。这些decircrna被用作构建ceRNA网络的候选环状rna。对demirna和demmrna进行进一步分析。

从数据库中预测目标基因

decircrna的靶标使用三个在线数据库进行预测:CircBanK44 (http://www.circbank.cn/)、CircRNA Interactome45 (https://circinteractome.nia.nih.gov/)和Cancer-Specific CircRNA Database46 (http://gb.whu.edu.cn/CSCD/)。用于构建ceRNA网络的候选mirna是通过将三个数据库预测的靶标与demirna交叉获得的。使用miRDB47、48 (http://mirdb.org/)和TargetScanHuman49 (https://www.targetscan.org/vert_72/)数据库进行候选mirna的靶标预测。将预测的靶标与demrna重叠，以获得构建ceRNA网络的候选mrna。

ceRNA网络建设

使用Perl匹配识别的候选基因之间的关系，准备构建ceRNA网络所需的输入文件。候选基因之间的关系应符合ceRNA机制:(1)circRNA表达越高，靶mirna表达越低，下游mrna表达越高;(2) circRNA表达越低，靶mirna表达越高，下游mrna表达越低。使用Cytoscape软件50(版本3.8.0,https://cytoscape.org/)构建ceRNA网络。

GO和KEGG富集分析

我们使用R软件和' org. hs . egg .db '包将ceRNA网络中mrna的基因符号转换为功能富集分析所需的ID。途径富集分析基于Kanehisa实验室开发的KEGG数据库51,52,53。使用“clusterProfiler、ggplot2、org. hh . egg .db、enrichment plot、circlize、RColorBrewer、dplyr、ggpubr和complexHeatmap”软件包进行GO和KEGG富集分析和数据可视化。以P < 0.05为筛选阈值。

生存分析

从TCGA数据库中检索并下载肝细胞癌的转录组和临床数据。用R和Perl软件对数据进行处理，得到基因表达和生存时间文件。使用R软件和“Survival, survminer”软件包对ceRNA网络中的mrna进行生存分析，以筛选与HCC患者生存预后相关的mrna。P < 0.05提示mRNA与HCC患者生存预后有相关性。

构建生存网络

利用从上述步骤中获得的生存相关mrna，获得了它们的上游相关mirna和circRNA的信息。相应准备了网络关系、节点等相关输入文件。细胞景观构建存活网络。使用circprimer2.0软件检索存活网络中circrna的结构、位置和其他信息54。

通路活性分析

基因集癌症分析(GSCALite)55是基于网络和多数据库的癌症基因组分析平台(http://bioinfo.life.hust.edu.cn/web/GSCALite/)。使用该平台分析了生存网络通路中的8个高表达mrna。根据需要输入生存相关基因，选择TCGA表达数据，癌症类型选择LIHC。通过通路活性分析模块进行分析，得出最终结果。

免疫细胞浸润水平和免疫检查点分析

肿瘤免疫估计资源(TIMER)56是一个公共网络平台(https://cistrome.shinyapps.io/timer/)，用于分析不同类型癌症的免疫细胞浸润水平。选择存活网络中高表达的8个mrna，分析肝细胞癌免疫细胞浸润水平。这些基因被输入到搜索框中，LIHC被选为癌症类型。选择免疫浸润相关细胞(B细胞、CD4+ T细胞、CD8+ T细胞、中性粒细胞、巨噬细胞和树突状细胞)进行相关性分析。基因表达谱交互分析(GEPIA)57是一个多功能在线平台(http://gepia.cancer-pku.cn/about.html)，用于分析肿瘤和正常细胞之间基因的差异表达。该平台还可用于与所选癌症类型相关的生存分析、相似基因检测、相关分析、降维分析等。本研究利用GEPIA的相关分析模块，确定了8种mrna的表达与肿瘤免疫逃逸相关的3个经典免疫检查点PDCD1、CD274和CTLA4的相关性。癌症类型选择LIHC，相关系数采用Pearson。

药物敏感性分析

GSCALite分析平台整合了癌症药物敏感性基因组学(Genomics of drug sensitivity in cancer, GDSC)和癌症治疗反应门户网站(cancer Therapeutics Response Portal, CTRP)中癌细胞系的药物敏感性和基因表达谱数据。我们对8种与生存相关的高表达mrna进行了药物敏感性分析。肿瘤类型选择LIHC，药敏模块选择GDSC，肿瘤药敏数据集选择GDSC，对半数最大抑制浓度(IC50)进行Spearman相关分析。调整参数后，结果生成。

ROC分析

采用R软件和“pROC”软件包进行ROC分析，明确circRNA-miRNA-mRNA调控轴上基因对HCC的诊断价值。HCC组织和对照组的表达数据使用上述分析中调整的三个微阵列数据进行分析，AUC≥0.7被认为具有HCC的诊断价值。所有统计分析均使用R软件或在线平台进行。

结合位点的确定和调控轴的模型构建

利用circprimer2.0软件查询调控轴上3种circrna的相关信息，得到它们的序列。从miRWalk (version 3, http://mirwalk.umm.uni-heidelberg.de/)下载hsa-miR-199a-5p序列，并通过序列比对找到circRNAs与miRNA之间的潜在结合位点。同时，我们使用Microsoft PowerPoint软件(版本2021,https://www.microsoft.com/zh-cn/microsoft-365/powerpoint)绘制了监管轴的模型图，以显示这些元素之间的控制关系。

道德声明

本研究使用的数据均来自公共数据库，没有进行任何与动物或人体标本相关的实验，无需伦理审批。