基于全转录组测序分析多形性胶质母细胞瘤增殖及侵袭的相关机制

王仕超，刘 斌，赵 星，任晓敏，潘园园，李 霞，徐淑英

（1.内蒙古自治区呼和浩特市第一医院医学遗传实验室，内蒙古 呼和浩特 010030；2.内蒙古自治区人民医院神经内科，内蒙古 呼和浩特 010030）

多形性胶质母细胞瘤（glioblastoma multiforme，GBM）是原发性中枢神经系统肿瘤中最常见的一类恶性肿瘤，其预后较差[1]。随着对神经胶质瘤相关机制研究的逐步深入，有学者通过靶向阻断神经胶质瘤增殖、侵袭或迁移的关键信号通路，达到了治疗神经胶质瘤的目的[2-3]，这是一种理论上可行，且具有广阔应用前景的新型治疗方法[4]。在细胞特定状态下进行全转录组的研究，可以从转录层面系统地揭示生物学背后的调控规律。目前，学者们对非编码RNA（non-coding RNA，ncRNA）的研究主要集中在微小RNA（microRNA，miRNA）、长链非编码RNA（long non-coding RNA，lncRNA）和环状RNA（circular RNA，circRNA）上。多项研究证实，高频突变转录本在细胞周期、能量代谢和信号传导等方面发挥重要作用，或可作为潜在的生物标志物应用于临床肿瘤的靶向治疗[5-6]。近年来，Hh（Hedgehog）通路异常激活与各种恶性肿瘤之间的关系越来越受到重视[7]。BOC蛋白作为Hh通路上游的跨膜蛋白，能够通过增加局部Hh配体浓度或延长与Hh配体结合时间来放大Hh信号[8]，提示BOC蛋白可以影响整个中枢神经系统的Hh通路反应。本研究拟在转录组水平揭示GBM增殖、侵袭的发生、发展机制，探寻具体有哪些关键基因参与了GBM的进展，在此基础上构建circRNA-miRNA -mRNA、lncRNA-miRNA-mRNA的关联分析，进一步验证BOC蛋白对GBM的生物学作用机制，为临床制定诊疗方案提供更可靠的理论依据。

1 材料和方法

1.1 研究对象

选取2018年11月—2019年2月呼和浩特市第一医院GBM患者4例，其中男2例、女2例，年龄43～65岁；另选取2019年12月呼和浩特市第一医院65岁男性GBM患者1例。所有患者均依据《中国中枢神经系统胶质瘤诊断与治疗指南（2015）》[9]确诊，病理学分期为Ⅳ级。本研究经呼和浩特市第一医院伦理委员会批准，所有患者或家属均知情同意并签署知情同意书。

1.2 样本收集及处理

收集手术中切除的癌组织和癌旁组织（距癌组织边缘1 cm）。以0.5 cm为标准直径对癌组织样本及癌旁组织样本进行块状切割，用0.9%Nacl溶液冲洗后，用灭菌冻存管封存，并做好标记，立即储存于液氮罐内。

1.3 文库构建及转录组测序

样本测序工作由天津诺禾致源基因测序公司完成，其中miRNA和lncRNA均独立建库，circRNA采用lncRNA建库，去除样本中的核糖体RNA（ribosomal RNA，rRNA）后，对线性RNA和circRNA进行测序。将2018年11月—2019年2月收集的4例GBM患者分别作为A组、B组、C组、D组。4例患者的癌组织样本分别编号为AT、BT、CT、DT；对应的癌旁组织编号为AC、BC、CC、DC。4组混测的癌组织编号为Q1，癌旁组织编号为Q2。4组混测方法与单样本测序方法一致。

1.4 主要测序流程及内容（以miRNA为例）

主要分析内容包括样本总RNA检测、分类注释、已知miRNA分析、ncRNA分析、重复序列比对、新miRNA预测、小分子RNA（small RNA，sRNA）分类注释统计、miRNA分析、miRNA碱基编辑分析、miRNA家族分析、miRNA表达及差异分析、miRNA靶基因预测、差异miRNA靶基因富集分析、差异lncRNA表达共定位、基因本体（Gene Ontology，GO）富集（包括分子功能、生物学过程和细胞组分）及京都基因与基因组数据库（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路富集分析。

1.5 组织样本中BOC mRNA的检测

采用实时荧光定量聚合酶链反应（polymerase chain reaction，PCR）检测5例GBM患者组织样本中BOCmRNA的相对表达量。PCR试剂盒购自上海英拜生物有限公司，检测仪器为QuantStudio 6 Flex荧光定量PCR仪（美国ABI公司）。引物由上海吉凯生物有限公司合成。引物序列见表1。反应体系：2×Master Mix 5.0 μL，PCR引物（10 μmol/L）各0.3 μL，加ddH2O至总体积为9 μL。反应条件：95℃ 10 min；95℃ 15 s，60℃60 s，45个循环。以GAPDH为内参，采用2-ΔΔCt法计算BOCmRNA的相对表达量。

表1 引物序列

2 结果

2.1 mRNA的差异分析

高通量测序分析共筛选出3 088个癌组织与癌旁组织差异表达的mRNA，其中表达上调1 375个，表达下调1 713个，见图1（a）。表达上调居前10位的基因分别为CUL2、HDGF、BOC、PTCH1、ACY1、BCL2L13、CRYAB、WNT5B、PIK3R3、ZNF658。表达下调居前10位的基因分别为PKA、Slmb、CK1、MMP8、LSM4、NFATC1、ZNF211、WDR7、MMP13、SHTN1。

4组间共有差异基因33个，分别为HIVEP2、NUB1、OLIG1、BOC、SOX8、LLGL1、PTK2B、SEMA6D、PRUNE2、SPOP、IL11RA、MAP4K4、SPOCK3、GPRASP1、POGZ、HIP1R、CALM3、ANK2、DIO2、OLIG2、DCX、F5、CLU、FLNB、KCNIP4、AK5、ZEB2、POU2F1、KCNJ10、RABGAP1L、BMPR1B、VIPR2、RHOT1。A组中癌组织与癌旁组织差异基因1 669个，其中表达上调768个、表达下调901个；B组差异基因3 571个，其中表达上调1 717个、表达下调1 854个；C组差异基因4 354个，其中表达上调1 715个、表达下调2 639个；D组差异基因5 708个，其中表达上调2 750个、表达下调2 958个。见图1（b）。

对差异表达的mRNA进行KEGG通路富集分析，选择富集最显著的20条通路信息，涉及的KEGG通路包括癌信号途径、丝裂原活化蛋白激酶（mitogen-activated protein kinase，MAPK）信号通路、cAMP信号通路、轴突导向及钙信号途径等通路。见图1（c）。这些通路与免疫反应、氨基酸代谢、细胞增殖、蛋白转运等生理生化反应有关。

对3 088个差异表达基因进行GO富集分析。结果显示，有2 206个基因涉及生物学过程，其中1 082个表达上调、1 124个表达下调；有2 311个基因涉及细胞组分，其中901个表达上调、1 410个表达下调；有2 222个基因涉及分子功能，其中866个表达上调、1 356个表达下调。见图1（d）。

图1 mRNA测序分析结果

2.2 lncRNA的差异分析

测序结果显示，差异表达的lncRNA共有106个，其中表达上调52个，包括lnc_000266、lnc_000276、lnc_000417、lnc_000510、lnc_000565、lnc_000662、lnc_000853、lnc_001155、lnc_001315、lnc_001370、lnc_001788、lnc_001827等；表达下调54个，包括ENST00000572151、lnc_000009、lnc_000021、lnc_000036、lnc_000115、lnc_000410、lnc_000467等。见图2（a）。

4组间共有的差异lncRNA数目为0，见图2（b）。A组和B组筛选出的癌组织与癌旁组织差异表达的lncRNA均为151个，C组和D组分别筛选出221和290个。

KEGG通路富集结果显示，差异lncRNA共表达mRNA主要涉及的细胞通路有血小板活化及基础转录因子、环磷酸鸟苷蛋白激酶信号以及神经活性配体-受体互作等相关途径。见图2（c）。

对1 098个差异表达的lncRNA靶基因进行GO富集分析，其中685个基因表达上调、413个基因表达下调，见图2（d）。

图2 lncRNA测序分析结果

2.3 circRNA的差异分析

测序得到差异表达的circRNA共82个，其中表达上调32个，包括hsa_circ_0047720、hsa_circ_0050334、hsa_circ_0067774、hsa_circ_0068375、hsa_circ_0069397、hsa_circ_0069982、hsa_circ_0078241等；表达下调50个，包括hsa_circ_0004027、hsa_circ_0004101、hsa_circ_0004550、hsa_circ_0006323、hsa_circ_0006528、hsa_circ_0010459、hsa_circ_0010466、hsa_circ_0017673、hsa_circ_0027641等。见图3（a）。

4组间共有的差异circRNA为245个，其中已知的circRNA 156个，包括hsa_circ_0000814、hsa_circ_0000894、hsa_circ_0000970、hsa_circ_0001277、hsa_circ_0001289、hsa_circ_0001475等；新定义circRNA 89个，包括novel_circ_0000585、novel_circ_0000886、novel_circ_0001046、novel_circ_0002180、novel_circ_0002215、novel_circ_0002371等。见图3（b）。

KEGG通路富集结果显示，差异表达的circRNA涉及的细胞通路有新陈代谢、MAPK信号通路、细胞内嗜作用等相关途径。共有的差异circRNA的GO富集分类统计图见图3（d），细胞组分分析以胞液与细胞间组分最为显著。

图3 circRNA测序分析结果

2.4 miRNA的差异分析

差异表达的miRNA共15个，其中表达上调8个，分别为hsa-miR-3157-3p、hsa-miR-6761-5p、hsa-miR-3124-5p、hsa-miR-5009-5p、hsamiR-216b-3p、hsa-miR-296-3p、novel_912、novel_328；表达下调7个，分别为hsa-miR-135b-5p、hsa-miR-377-3p、hsa-miR-136-5p、hsa-miR-381-3p、hsa-miR-153-5p、hsa-miR-299-3p、hsamiR-4473。见图4（a）。

绘制多组比较的差异miRNA韦恩图，选择2个组、3个组或4个组比较得到的差异miRNA数进行统计。结果显示，4组间交集的miRNA数为0。见图4（b）。

差异表达miRNA的KEGG富集通路散点图见图4（c）。

图4 miRNA测序分析结果

4组癌组织样本（AT、BT、CT、DT）之间的r2值接近0.95，癌旁组织（AC、BC、CC、DC）之间的r2值接近0.93，表明不同样本之间表达模式的相似度较高。见表1。

表1 4组不同样本之间的相关性

2.5 mRNA-miRNA-lncRNA/circRNA的关联分析

在竞争性内源RNA（competitive endogenous RNA，ceRNA）理论的基础上，根据结合点位是否相同来寻找与miRNA同类的lncRNA与基因的配对等，完成ceRNA调控网络的构建。见图5、图6。

图5 mRNA-miRNA-circRNA的关联分析

图6 mRNA-miRNA-lncRNA的关联分析

利用miRBase数据库构建33个差异表达基因的miRNA-lncRNA-mRNA 的关联分析，见图7。其中BOC基因位于调控网络的关键节点，以BOC为核心进行拓扑分析，结果显示，BOC与has-miR-4763-5p、has-miR-6848-5p、has-miR-6860 miRNA具有潜在的靶向结合位点，这3类miRNA下游存在POU2F1、TP53等与肿瘤发生相关的基因转录本，见图8。

图7 mRNA-miRNA-lncRNA/circRNA的关联分析

图8 BOC与ncRNA、基因的关联分析

采用STRING数据库中的Confidence Score可靠指数对每个预测的相互关联性进行分析。结果显示，SHH、CDON、PTCH1等蛋白与BOC关联，主要定位在Hh通路上。见图9。

图9 BOC-蛋白互作分析图

2.6 癌组织与癌旁组织中BOC mRNA相对表达量的比较

5例GBM患者中有4例患者癌组织BOCmRNA水平明显高于癌旁组织（P＜0.05），1例患者癌组织BOCmRNA水平低于癌旁组织（P＜0.05）。见图10。

图10 5例GBM患者癌组织与癌旁组织BOC mRNA相对表达量比较

3 讨论

ncRNA在各种生物过程中起着重要的调控作用。本研究基于全转录组测序，同时分析5例GBM患者癌组织与癌旁组织中差异表达的mRNA、lncRNA、circRNA和miRNA，并且通过两两关联分析、三元关联分析、多元关联分析，深入挖掘相关转录调控机制。

本研究筛选出的差异表达的mRNA中已有多个被证实与胶质瘤的发生、发展密切相关，如Rab23[10]、RAS[11]、BRAF[12]、VEGF[13]、PTCH1[14]、MAPK[15-16]等基因。有研究结果显示，Rab23 mRNA及蛋白在胶质瘤组织中高表达，其表达水平与病理分级相关[10-11]。本研究结果显示，癌组织RAS家族中与增殖相关的蛋白（ERAS、GEM、HRAS、KRAS、MRAS）的mRNA相对表达量显著高于癌旁组织，与文献报道[13]一致。

BRAF基因编码的蛋白属于raf/mil家族的丝氨酸/苏氨酸蛋白激酶，参与调控MAPK/细胞外调节蛋白激酶信号通路，通过细胞分裂、分化、分泌等功能促进肿瘤形成[12]。MAPK是在真核生物中非常保守的丝氨酸/苏氨酸蛋白激酶，参与了细胞增殖、分化、运动及凋亡等生物过程。本研究结果显示，胶质瘤患者癌组织BRAFmRNA、MAPKmRNA呈高表达，提示BRAF、MAPK等癌基因参与了GBM的恶性增殖、侵袭等生物学过程。

本研究筛选得到的差异表达mRNA的转录本富集通路包括Shh（Sonic hedgehog）信号途径、胶质瘤信号途径、胰腺癌信号途径、轴突导向、TNF信号途径等通路，这些通路与免疫反应、氨基酸代谢、细胞增殖、蛋白转运等生理生化反应相关。BOC作为跨膜蛋白定位在Shh信号途径上游。本研究结果显示，癌组织BOCmRNA水平高于癌旁组织（P＜0.05），与文献报道[17]一致。SHERGALIS等[17]利用TCGA数据库发现胶质瘤组织BOC高于正常脑组织（P＜0.05），高表达BOC的胶质瘤患者生存期明显缩短。提示BOC激活了Hh信号途径，可促进GBM的进展，因此BOC有成为药物靶点的潜力。

GBM的发生、发展与mRNA、lncRNA、miRNA、circRNA之间复杂的调控关系密切相关。本研究构建GBM mRNA-miRNA-lncRNA/circRNA调控网络后分析了网络拓扑特性，发现BOC、miR-4763-5P、miR-6848-5P、miR-6880是关键节点。进一步对这些关键节点及其相互作用的RNA进行功能富集分析和生存分析，结果显示，与miR-6880具有潜在结合位点的lncDUBR、lnc003875已被证实与多种肿瘤发生相关[18]，提示BOC可能通过与miRNA结合，调控下游基因来影响肿瘤的进展。lnc003875/miR-363/EGR1在癌相关成纤维细胞调控网络中的作用是控制人胎盘部位滋养层肿瘤的血管生成。另外，lnc003745、lnc003874、lnc001254、lnc001763、lnc002545、lnc002436、lnc003671、lnc3371、lnc00798、DCTN1-AS1、RP11-1055B83、lnc00653、c20orf166-AS1共计13个lncRNA也出现在BOC相关的互作网络中，其中C20orf166-AS1与抗原处理、表达和细胞黏附分子通路相关。目前，C20orf166-AS1已被证实与MAPK信号通路相关，可导致前列腺癌、胃癌的发生[19]。lnc001763、lnc002436与病灶粘连、细胞外基质受体相互作用、MAPK信号通路相关，共涉及到12个基因。

关于circRNA在恶性肿瘤发生、转移、侵袭性生长和肿瘤血管新生等过程中的调控作用是目前肿瘤研究的热点。本研究结果显示，circRNA0000543、circRNA0003666、circRNA0047270、circRNA0010465等8个circRNA参与了BOC的互相作用。既往研究发现，circSMARCA5和circCFH与胶质瘤的特异性表达模式有关，这些circRNA或可作为肿瘤的生物标志物[20]。此外，一些circRNA已被发现在GBM中发挥致癌作用，如circNFIX、circNT5E；而另一些circRNA则具有抑癌作用，如circFBXW7、circSHPRH[21]。本研究还发现了多个差异表达且尚未报道的circRNA，其对GBM的调控机制还有待进一步验证。

本研究全转录组多元关联分析结果显示，BOC位于调控网络关键节点位置，其定位的Hh信号途径在被异常激活的情况下会促进多种恶性肿瘤的进展，提示BOC可能受相关miRNA、circRNA、lncRNA的调控，影响下游靶基因的表达，参与了肿瘤增殖与侵袭的过程。本研究还对5例GBM患者癌组织和癌旁组织中BOCmRNA的相对表达量进行了检测，结果显示，有4例患者癌组织BOCmRNA水平明显高于癌旁组织（P＜0.05），1例患者癌组织BOCmRNA水平低于癌旁组织（P＜0.05）。后续研究将重点针对BOC是否在GBM的致病机制中发挥癌基因作用，进一步在细胞层面进行验证。

综上所述，本研究利用二代测序技术对4例GBM患者的肿瘤组织及癌旁组织进行全转录组测序分析，筛选后得到33个共有差异基因，深入分析后挖掘到1个可能与GBM增殖、侵袭、迁移密切相关且定位于Hh通路的关键基因BOC。进一步对mRNA-miRNA-lncRNA/circRNA多元关联网络拓扑分析发现，BOC位于调控网络的关键节点位置，且其定位的Hh信号途径在异常激活的情况下会促进多种肿瘤的恶性进程，提示BOC受到相关miRNA、circRNA、lncRNA的调控，影响下游靶基因的表达，从而导致GBM的发生。