乳腺癌中基于開放染色質的環狀RNA互作網絡構建及功能分析

張蘊顯 王雅梅 周 萍*

(1.首都醫科大學生物醫學工程學院生物醫學信息學系，北京 100069；2.首都醫科大學基礎醫學院生物化學與分子生物學系，北京 100069)

女性乳腺由皮膚、纖維組織、乳腺腺體和脂肪組成。乳腺癌是一種發生在乳腺上皮組織中的惡性腫瘤。目前，乳腺癌已成為威脅女性身心健康的常見腫瘤。據估計，2017年美國有255 180例新發病例和41 070例乳腺癌死亡病例[1]。乳腺癌常見的3種亞型是雌激素受體陽性(estrogen receptor positive,ER+)、人類表皮生長因子受體2陽性(HER2+)和三陰型。相關研究[2-5]表明肥胖、雌激素和孕激素的使用、高齡初產、酗酒等因素能夠增加乳腺癌的發生風險。除了以上幾種因素外，基因突變、表觀遺傳機制在乳腺癌的發生、發展、轉移中也起到重要作用。過去幾年中，很多研究[6-8]表明非編碼RNA(non-coding RNA，ncRNA)參與到致癌過程中，例如微小RNA(micro RNA，miRNA)，長非編碼RNA(long non-coding RNA，lncRNA)等[6]，可以作為乳腺癌識別、預后的生物標志物[7]。作為新的內源性的ncRNA，環狀RNA(circular RNA，circRNA)對腫瘤發生發展過程的調控作用已經成為研究熱點之一[8]。

circRNA是一種新的內源性ncRNA，其長度為幾百到幾千個核苷酸[9]。與線性RNA轉錄產物不同的是，circRNA通過RNA的3′端和5′端的共價鍵形成，具有穩定的結構[10-11]，但circRNA的生物發生過程尚不清楚。盡管如此，越來越多的研究[6]證明circRNA參與到了一系列的病理發生過程中。circRNA參與疾病發生過程的機制有以下幾種途徑：①circRNA作為競爭性內源RNA(competing endogenous RNAs,ceRNA)和“海綿吸附miRNA”(miRNA sponges)參與致癌的調節過程。miRNA通過與mRNA的3′端結合負向調節mRNA，有關研究[12]顯示，其他具有miRNA靶點的RNA也同樣可以與mRNA競爭miRNA。研究[13]顯示含有miRNA的反應元件(miRNA reponse element，MER)的circRNA可以與miRNA相互作用并充當“海綿吸附miRNA”，以此削弱miRNA對靶基因的抑制作用。②circRNA可以特異性結合蛋白質，直接或通過RNA間接阻斷蛋白質的功能[14-15]。③調節基因轉錄，雖然大部分circRNA作為ceRNA和“海綿吸附miRNA”調節miRNA，但相關研究[16]表明，一些circRNA可以順式或反式的調節基因的轉錄。④普遍認為circRNA是一種特殊的內源性的ncRNA，不能直接翻譯為蛋白質。但最近的研究[17-18]顯示，一些動物體內circRNA(如小鼠、果蠅)可以直接翻譯蛋白質并調控其生理過程。circRNA在乳腺癌患者體內存在差異表達，可能是乳腺癌細胞增生的生物標志物且與乳腺癌亞型相關[19-20]。

轉座酶可接近性染色質測序技術(assay for transposase accessible chromatin with high-throughput sequencing, ATAC-seq)，可以只利用少量細胞快速得到調控的多維信息。ATAC-seq可以同時獲得“開放”染色質的位置、轉錄因子的結合位點、核小體的調控區域和染色質狀態等信息，在表觀遺傳機制研究領域有廣闊的應用前景[21-22]。脫氧核苷酸酶Ⅰ超敏感位點測序技術(DNase Ⅰ hypersensitive sites sequencing，DNase-seq)是利用脫氧核苷酸酶Ⅰ(DNase Ⅰ) 超敏感位點捕獲技術，通過DNase Ⅰ對裸露在蛋白質包裹之外的DNA進行酶切，捕獲基因活化區域DNA片段，然后進行高通量測序[23]。

本文利用ATAC-seq和DNase-seq數據，全基因范圍分析乳腺癌MCF-7、T-47D和MDA-MB-231細胞系中circRNA的轉錄調控，進而構建了乳腺癌中活化circRNA基因的功能互作網絡，從而挖掘出乳腺癌中的主要功能組circRNA和關鍵活化circRNA基因。

1 材料與方法

1.1 數據

本研究使用到乳腺癌的MCF-7、T-47D、MDA-MB-231細胞系的ATAC-seq和DNase-seq測序數據和人類circRNA數據庫。①MCF-7、MDA-MB-231細胞系的ATAC-seq實驗數據下載于基因表達數據庫(Gene Expression Omnibus, GEO)(網址：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE97583)的GSE97583[24]實驗數據。②MCF-7、T-47D細胞系的DNase-seq實驗數據下載于UCSC基因瀏覽器中Table Browser的Duke DNase-seq HS實驗數據(網址：https://genome.ucsc.edu/cgi-bin/hgTables)，和GEO數據庫的GSE108167(網址：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE108167)實驗數據。③circRNA的高通測序數據下載于circBase中所有已發現的人類circRNA[25-27]實驗數據(網址：http://www.circbase.org/)。

1.2 軟件

本文所用軟件有：Windows 10，Access 2010數據庫，Cytoscape 3.6.1(網址：https://cytoscape.org/)版本的圖形顯示軟件[28]，Cytoscape3.6.1插件ClueGo V2.5.2[29]，KEGG pathway數據庫(網址：https://www.kegg.jp/kegg/mapper.html)。

1.3 方法

本文工作主要環節為特異circRNA篩選、網絡構建和功能分析，方法流程圖如圖1所示。

圖1 分析方法流程圖Fig.1 The flow chart of analyzing method

1.3.1 于Access數據庫的數據挖掘

首先利用Access 2010構建“乳腺癌開放染色質高通測序數據庫”，該數據庫基礎表由6個高通測序樣本、人類基因hg19(網址：http://genome.ucsc.edu/)和circRNA(網址：http://www.circbase.org/)構成，基因版本配準采用在線轉換工具：http://genome.ucsc.edu/cgi-bin/hgLiftOver。然后利用SQL語句對數據進行篩選挖掘。數據篩選實例詳見圖2(T-47D的DNase-seq實驗數據數據和circRNA數據為例)。

圖2 數據篩選實例Fig.2 Examples of data mining

以此初步篩選出與乳腺癌相關的circRNA，以備進一步分析。

1.3.2 活化circRNA基因KEGG pathway富集篩選

將篩選出的活化circRNA對應的基因放入KEGG數據庫中，利用KEGG Mapper對這些特異基因進行功能分析，并進一步篩選出與腫瘤(特別是與乳腺癌)相關通路的活化circRNA基因。在此結果中最終篩選出3個細胞系共有活化circRNA對應的基因，依據Gene Ontology(GO)的分子生物過程構建乳腺癌活化circRNA基因互作網絡。

1.3.3 基于GO生物過程的乳腺癌活化circRNA基因互作網絡構建

在Cytoscape 3.6.1平臺中，使用ClueGo 2.5.2將經過數據挖掘篩選出的活化circRNA基因的集合和KEGG pathway富集篩選出的3個細胞系共有的活化circRNA基因的集合依據GO 生物過程構建乳腺癌活化circRNA基因互作網絡，GO版本為Biological Process-EBI-UniProt-GOA(日期：2018-04-09)，網絡構建方法及參數選擇：依據GO 術語進行分組，網絡特異性為中等，P<0.05，為了細化模組功能，GO 術語層級設置為10～20級，模組最少基因數為3，占所在功能或通路基因的百分比為4%以上，網絡連接度的Kappa 評分為0.5，富集分析為雙尾超幾何分布檢驗，P值校正方法為“Bonferroni step down”。

2 結果

2.1 MCF-7高敏感位點circRNA構成比

本文分析MCF-7細胞捕獲高敏感位點circRNA為7 458個，含有utr區的基因是10 163個，捕獲不含3utr和5utr的基因的circRNA是3 501個。由此可見在乳腺癌中circRNA的表達調控是非?；钴S的。根據circRNA綜合數據庫CIRCpedia v2(網址：http://www.picb.ac.cn/rnomics/circpedia/)的研究，線性RNA中也包含一定比例的circRNA[30]。為此本研究以DNase-seq和ATAC-seq所有捕獲的circRNA為研究對象，即不考慮該基因是否包含utr區。

2.2 3種乳腺癌亞型對比

T-47D捕獲circRNA最高，其次是MCF-7，最少的是MDA-MB-231。就兩種實驗而言，DNase-seq捕獲率遠高于ATAC-seq。將3種細胞捕獲的活化circRNA基因分別構建了互作網絡，從結果看它們拓撲結構基本相似，只是每個功能模組基因構成比例略有差異。

綜上所述，萊考夫非常強調概念隱喻在數學認知中的作用，認為概念隱喻是理解復雜數學思想的核心認知機制。概念隱喻是算術所需要的一種重要的認知能力。人們正是通過隱喻對數學中的概念進行認知和理解，概念隱喻使人們超越了極少的先天算術和簡單的數數能力，延伸了人們的認知能力，獲得了進一步的算術能力。物體集合隱喻、對象建構隱喻、量尺隱喻和沿路線運動隱喻是四種基本的基礎隱喻，也是人們擴展算術的重要的隱喻能力。這四種基本的基礎隱喻都是以隱喻描述了人們的日常經驗與數字之間的映射關系。

本研究結果顯示，3種細胞捕獲的活化circRNA基因主要富集在細胞代謝的生物過程中，其中DNase-seq實驗捕獲的MCF-7和T-47D細胞的活化circRNA基因富集在有機氮化合物代謝過程(分別占11.58%、10.67%)；ATAC-seq實驗捕獲的MDA-MB-231和MCF-7細胞的活化circRNA基因分別主要富集在細胞大分子代謝過程(15.43%)和細胞蛋白質代謝過程(18.35%)。

除了細胞代謝的生物過程，本研究捕獲的活化circRNA基因也顯著富集在其他與腫瘤相關的生物過程中，包括：DNase-seq實驗捕獲的MCF-7細胞的活化circRNA基因顯著富集在細胞正調節(7.61%)、細胞內轉導(3.65%)、細胞組分生物合成調控(3.57%)、生物過程負調控(3.01%)、細胞有絲分裂周期(2.3%)；DNase-seq實驗捕獲的T-47D細胞的活化circRNA基因顯著富集在細胞正調節(7.26%)、超分子纖維組織(3.48%)、細胞蛋白質定位(3.63%)、生物過程負調控(2.74%)、細胞有絲分裂周期(2.3%)；ATAC-seq實驗捕獲的MDA-MB-231細胞的活化circRNA基因顯著富集在細胞組分組織正調節(10.8%)、細胞分解代謝過程(8.49%)、細胞內轉導(4.94%)、細胞骨架組織(4.01%)、蛋白質定位(3.24%)、細胞有絲分裂周期(3.09%)；ATAC-seq實驗捕獲的MCF-7細胞活化circRNA基因顯著富集在細胞過程正調節(9.8%)、細胞器蛋白質定位(9.24%)、細胞蛋白質定位(5.6%)、細胞組分裝配(4.06%)、蛋白質定位(2.94%)、有絲分裂細胞周期過程(2.8%)。為此本研究用3種細胞系共有的circRNA進行進一步研究。

2.3 KEGG功能分析

對3種細胞捕獲的活化circRNA基因進行KEGG功能分析，由結果可見比對到的關鍵信號通路皆與腫瘤有關。本研究中選取了對腫瘤起關鍵作用的18個信號通路(表1)。根據這18個信號通路最終篩選出了326個共有活化circRNA基因，并據此進行共有基因的表達調控互作網絡構建及GO分子功能富集分析。

表1 活化circRNA基因所屬關鍵信號通路Tab.1 Key signaling pathways for activated circRNA gene

2.4 互作網絡構建及分子功能分析

為了全面探索3個細胞系共有活化circRNA基因的生物學功能，本研究進行分子功能富集分析及生物過程的表達調控互作網絡構建。分子功能富集分析結果如圖3所示，由結果可見,這些活化circRNA基因顯著富集在結合蛋白質和刺激活化功能區。

圖3 活化circRNA基因分子功能Fig.3 Molecular function of activated circRNA gene

利用Cytoscape 3.6.1的ClueGo插件構建的生物過程表達調控互作網絡結果如圖4所示。

圖4 活化circRNA基因功能互作網絡Fig.4 Functional interaction network of activated circRNA genes

根據構建的生物過程表達調控互作網絡，其功能模組主要分布于細胞代謝、重要信號通路、生物過程的調控、生物過程反應等生物過程中，結果如圖5所示。根據主要的功能模組最終本研究識別出了PTK2B、PDCD6IP、ABL1、EGFR、RHOA、MTOR、NRP1、ATR、CTNNB1、ILK、NF1、PRKCA、HNRNPK、MEF2C、PMAIP1、LYN、NR4A3、SRF、AREG、PML、PTK2、ROCK2、TGFBR2、VEGFA、DLG1、HRAS、ITGA2、CREB1、STAT3、RUNX2、TIAM1等31個hub基因。這些基因所對應的捕獲到的活化circRNA如圖6所示。

圖5 活化circRNA基因重要功能組Fig.5 Important functional groups of activated circRNA genescircRNA:circular RNA.

圖6 hub基因的活化circRNA列表Fig.6 Activated circRNA list of hub genescircRNA:circular RNA.

從圖6中看出，PTK2基因對應的活化circRNA最多，PTK2在ErbB信號通路中，對細胞粘連遷移起直接作用，由于PTK2活化circRNA多，因此circRNA吸附的miRNA的數量和種類也會隨之增多，從而降低了miRNA對PTK2基因的mRNA調控。其次是MOTR基因對應的活化circRNA，MOTR基因在乳腺癌通路中直接對細胞的增生和存活起作用。由此可見上述circRNA在乳腺癌中起重要作用。

由于circRNA可以直接影響其對應基因的蛋白質表達水平， circRNA活化程度越高其對應基因的蛋白質表達水平也會增高，圖5中顯示了依據活化circRNA構建的基因互作網絡功能組，更顯示出了circRNA在乳腺癌中的作用。

3 討論

隨著circRNA領域的發展，已經建立了許多關于circRNA的數據庫以促進circRNA的分析。circBase，CIRC pedia v2和Deepbase 2.0數據庫包含許多關于不同物種的circRNA，并提供詳細信息。CSCD和CircInteractome數據庫可用作預測miRNA反應元件和RNA結合蛋白質的工具。CirclncRNAnet數據庫提供了一種分析測序結果的簡便方法。 ExoRBase數據庫提供了人類血液外泌體中存在的58 330個circRNA。 circRNA Disease數據庫記錄了多種疾病中經過實驗驗證的circRNA[31]。Starbase v2.0數據庫是一個提供miRNA-ncRNA互作網絡和生物信息學應用的數據庫，該數據庫的開發者利用108個CLIP-seq(PAR-CLIP，HITS-CLIP，iCLIP，CLASH)數據集鑒定出了9 000多個miRNA-circRNA互作關系[32]。最近有研究人員[33]使用Starbase v2.0數據庫鑒定CLIP-seq數據中miRNA-circRNA的互作關系，并開發了circSeeker，circAnno和clipSearch生物信息學分析軟件用于注釋和識別circRNA及其與miRNA的相互作用。目前對circRNA生物信息學分析大多基于表達譜數據篩選在腫瘤細胞內差異表達的circRNA，并利用上述公開數據庫和軟件分析預測circRNA-miRNA的互作網絡，例如有研究人員[34]利用表達譜數據篩選出了肝癌細胞中差異表達的circRNA，并利用CSCD和CircInteractome數據庫預測差異表達的circRNA的潛在靶miRNA，構建circRNA-miRNA-mRNA互作網絡，揭示circRNA對肝癌的調控作用。本研究立足于circRNA的表達調控，借助Access數據庫使用SQL語句對MCF-7、T-47D、MDA-MB-231細胞系的ATAC-seq和DNase-seq實驗數據進行分析挖掘，鑒定出在乳腺癌細胞中可能出現circRNA的活化染色質位點，篩選出3種乳腺癌共有的活化circRNA基因，并從功能角度構建活化circRNA互作網絡。由于目前對circRNA的注釋信息尚不完善，而對基因的注釋比較完整，因此筆者利用ClueGo插件對活化circRNA基因進行生物學分析，并構建活化circRNA基因的生物過程調控互作網絡，以此預測篩選出的活化circRNA對乳腺癌的調控作用。

目前相關研究[35]表明，癌細胞的一個新興特點是改變代謝。腫瘤的發生發展依賴細胞代謝重編程，細胞代謝重編程也是致癌性病變的直接或間接結果[36]。本研究的結果顯示，乳腺癌細胞染色質中處于活化狀態的circRNA對應的基因，其參與的生物過程主要是細胞代謝，正、負調控，細胞周期，細胞凋亡、增生、分化，重要信號通路，細胞反應等重要的細胞生物過程，因此表明circRNA在乳腺癌的發生發展過程中起到重要作用。尤其是發現的hub基因的circRNA起到至關重要的作用。

本研究發現盡管一個基因編碼區可能包含多個circRNA，但并非捕獲到基因編碼區上全部circRNA。例如PTK2B基因在其編碼區內有8個circRNA：hsa_circ_0083760、hsa_circ_0083759、hsa_circ_0083765、hsa_circ_0083764、hsa_circ_0083761、hsa_circ_0083762、hsa_circ_0083758和hsa_circ_0083763，但在本研究中僅捕獲到hsa_circ_0083759和hsa_circ_0083762，根據miRTarBase(miRNA靶基因預測數據庫)數據庫hsa_circ_0083759和hsa_circ_0083762可能作為hsa-miR-23b-3p、hsa-miR-517c-3p和hsa-miR-517a-3p的競爭性內源RNA和“海綿吸附miRNA”，這說明在乳腺癌中circRNA的啟動是有針對性的，很可能是針對某些特定的miRNA作為競爭性內源RNA和“海綿吸附miRNA”。

在本研究中連接度最高的circRNA對應的基因是PTK2B，該基因編碼細胞質蛋白質酪氨酸激酶。最新研究[37]表明，PYK2B水平可以預測結腸腺癌手術切除后的預后。各種治療模型中circRNA的表達譜尚未得到廣泛研究，但已經報道了在放射抗性癌細胞中cirRNA的差異表達[38]。雖然目前乳腺癌還沒有廣泛開展PYK2B基因的circRNA研究，但是本文研究表明，PYK2B基因的circRNA很可能對乳腺癌的診斷治療起重要作用。最近的研究[39]已經證明外泌體介導腫瘤微環境中不同類型細胞之間的相互作用，它們調節腫瘤生長、轉移、耐藥性(通過轉運腫瘤相關mRNA，miRNA和蛋白質)、血管生成、免疫逃逸和其他過程。PDCD6IP基因對外泌體的發生有著重要的調控作用，能直接介導腫瘤微環境中對腫瘤的重要調控作用。ABL1基因是一種蛋白質酪氨酸激酶和有效的致癌基因，相關研究[40]證實，它是miR-203(一種腫瘤抑制性miRNA)的真正靶標。ABL1基因通過與miR-203結合導致miR-203腫瘤抑制能力減弱或者喪失，使正常抑制的致癌基因不受表達，從而達到調控腫瘤發生發展和轉移。雖然ABL1基因的circRNA對乳腺癌的作用目前尚未有相關研究，但ABL1基因的circRNA極有可能作為ceRNA和“海綿吸附miRNA”參與到乳腺癌致癌的調節過程。最近的相關研究[40]證明circHIPK3可以使內源性miR-7結合以隔離和抑制miR-7活性，從而導致表皮生長因子受體(epidermal growth factor receptor，EGFR)表達增加達到調控直腸癌生長和轉移的作用[41]。雖然尚未有EGFR基因的circRNA對乳腺癌的作用的相關研究，但本研究證明EGFR基因的circRNA很可能對乳腺癌的診斷及預后有重要作用。RhoA表達水平在多種腫瘤組織中增加并且與腫瘤惡性相關，表明RhoA在腫瘤發生和轉移中的重要作用。相關研究[42]證實miR-200b可介導RhoA基因以形成內源競爭并在包括肝細胞癌在內的多種腫瘤的發生中具有重要的調節作用。研究[42]顯示，通過多種信號轉導途徑，RhoA參與腫瘤的惡性轉化、侵襲、轉移、血管生成等。RhoA基因的circRNA極有可能作為ceRNA和“海綿吸附miRNA”調控RhoA基因的表達并參與到乳腺癌致癌的調節過程中，并作為用于治療癌癥轉移的潛在的治療標靶。

本文從生物信息學角度對高通測序數據進行深度挖掘，發現的hubcircRNA基因及其對應的circRNA對乳腺癌基礎臨床研究具有一定的參考價值。

本研究通過3種乳腺癌細胞的DNase-seq數據和ATAC-seq數據對3種乳腺癌細胞內circRNA捕獲率的分析，發現DNase-seq實驗的捕獲率及精度要高于ATAC-seq實驗。通過對捕獲的活化circRNA基因的生物信息學分析發現，乳腺癌細胞內的活化circRNA基因大多參與到細胞代謝，正、負調控，細胞的凋亡、增生、分化，以及目前已證實的跟腫瘤相關的一些重要信號通路中，因此表明活化circRNA在乳腺癌的發生發展過程中起到重要作用。本研究顯示基因編碼區包含多個circRNA，但在本研究中僅捕獲基因編碼區中對應的部分circRNA，表明這些基因對應的circRNA可能針對某些特定的miRNA作為ceRNA和“海綿吸附miRNA”參與到乳腺癌致癌的調節過程。通過對本研究識別出的hubcircRNA基因進行生物信息學分析，發現該基因參與到了腫瘤的惡性轉化、侵襲、轉移等過程中，這表明hubcircRNA在此過程中起重要作用，雖然目前尚未廣泛研究這些hubcircRNA對乳腺癌的作用，但本研究的結果表明這些hubcircRNA可能作為用于治療癌癥轉移的潛在靶點。