基于全轉錄組測序分析多形性膠質母細胞瘤增殖及侵襲的相關機制

王仕超，劉 斌，趙 星，任曉敏，潘園園，李 霞，徐淑英

（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的發生。