Profiling expression of coding genes, long noncoding RNA, and circular RNA in lung adenocarcinoma by ribosomal RNA-depleted RNA sequencing

用去除核糖體RNA（rRNA-depleted）的RNAseq技術研究肺腺癌mRNA, lncRNA, circRNA表達譜

期刊：FEBS Open Bio 影響因子：2.143

發(fā)表單位：江蘇省腫瘤醫(yī)院放射科

1 摘要

    非編碼RNA在多種生物學過程、疾病包括腫瘤中發(fā)揮著重要作用。肺腺癌lung adenocarcinoma（簡稱LUAD），還沒有進行過circRNA的系統(tǒng)研究。本研究，我們使用核糖體RNA去除的RNAseq技術對配對肺腺癌樣品做了轉(zhuǎn)錄組表達譜研究包括編碼基因、lncRNA, circRNA。測序得到的reads首先比對到人基因組分析mRNA和lncRNA。未能比對的reads使用circRNA預測算法鑒定候選circRNA。一共鑒定到1282個差異表達的編碼基因。檢測到19023個lncRNA，其中244個差異表達。AFAP1-AS1, BLACAT1, LOC101928245, and FENDRR 差異變化最高。大于2個backspliced reads的circRNA共鑒定到9340個，包括3590個新的circRNA，circRNA長度中值是530nt。circRNA豐度較低，超過一半的circRNA的backspliced reads <10。使用瓊脂糖電泳和sanger法測序驗證4個候選的circRNA，發(fā)現(xiàn)確實是環(huán)狀的。本研究概述了肺腺癌mRNA，lncRNA，circRNA的表達譜。鑒定到9340個circRNA，提示circRNA在肺腺癌中廣泛表達。

    注：原始測序數(shù)據(jù)存儲在GEO數(shù)據(jù)庫，編號為GSE104854 。

2 材料方法

1 病人及組織樣品

    本研究獲得江蘇省腫瘤研究所倫理會的支持。配對的腺癌樣品及癌旁樣品來自于非小細胞肺癌（NSCLC）手術樣品，時間是2012-2014，來源于江蘇省腫瘤醫(yī)院胸外科。由病理診斷豐富的病理醫(yī)生判讀配對的癌與癌旁樣品。收集病人的臨床和病理資料，同時獲得病人的知情同意書。共收集9對肺腺癌樣品，包括癌與癌旁組織。RNA樣品混合測序共測得6個樣品（3個癌組織，3個癌旁組織）。臨床資料及RNA樣品信息如Table S1。RNAseq的數(shù)據(jù)也已上傳至GEO數(shù)據(jù)庫。

2 細胞培養(yǎng)及PCR

    用RPMI 1640 (GIBCO, Invitrogen, Shanghai, China) 培養(yǎng)基培養(yǎng)A549細胞，培養(yǎng)基加入10% FBS (GIBCO-BRL, Invitrogen, Carlsbad, CA, USA), 100U.mL-1 penicillin, and 100mg.mL-1 streptomycin (Key-GEN)，在空氣潮濕，溫度37 °C，5% CO2環(huán)境中培養(yǎng)。使用TRIzol (Invitrogen, Grand Island, NY, USA) 提取A549的細胞的totalRNA。使用PrimerScript RT Master Mix (Takara, Dalian, China)將1ug totalRNA在20ul體系中反轉(zhuǎn)錄。引物如Table S2。產(chǎn)物使用瓊脂糖凝膠電泳分離并進行Sanger測序。

3 樣品準備

    使用TRIzol (Invitrogen, Shanghai,China)提取肺腺癌樣品及癌旁樣品total RNA。Bioanalyzer 2200分析RNA的完整性，并保存在-80℃。使用GeneRead rRNA Depletion Kit (QIAGEN, Venlo, Germany)去除rRNA。

4 文庫構(gòu)建及RNA測序

    使用Ion Total RNA-Seq Kit v2.0 (Thermo Fisher Scientific, Inc., Waltham, MA, USA)準備單端cDNA文庫。cDNA文庫構(gòu)建根據(jù)試劑盒說明書進行。測序儀器使用 proton sequencers（此儀器建庫上機過程可查相關資料，并且此儀器目前很少使用，此處不多介紹）。

5 RNA測序及reads比對

    對原始下機數(shù)據(jù)reads去除接頭序列，去除“N”超過5%的序列，去除堿基質(zhì)量值<13含量大于20%的reads。MAPSPLICE 軟件 (v2.1.8) 將clean reads比對到人基因組(version: GRCh37)，為了獲得更多的可變剪切轉(zhuǎn)錄本，設置參數(shù)為(–s 22 –p 12 –ins 6 –del 6 –noncanonical)。

6 circRNA鑒定及定量

    使用acfs流程https://code.google.com/p/acfs/鑒定circRNA， 主要包括以下幾步：

    1）獲得未能比對到基因組的reads。使用BOWTIE2 version 2.2.5將reads比對到參考基因組[GRCH37.p13 NCBI]。

    2）circRNA鑒定。未能比對的reads，使用BWA mem (bwa mem –t 1 –k 16 –T 20) 再次比對。部分比對的reads（a）兩部分比對到同一個染色體上，距離不超過1MB；（b）同一條鏈；（c）比對方向相反，支持head-to-tail連接。將這類reads的連接點進一步用MaxEntScan33評估。根據(jù)受體與供體點的最高得分來確定準確連接點。支持head-to-tail連接點的reads數(shù)量≥2，連接點的得分值≥10。

    3）表達量分析：為了得到circRNA的表達量。我們將unmapped reads使用 BWA mem 參數(shù) (bwa mem –t 1 –k 16 –T 20) 重新比對到候選circRNAs。由于沒有直接的證據(jù)能夠確定circRNA的序列，我們通過中間存在的外顯子注釋獲得。序列從5’端延伸到3’端以形成環(huán)狀。能夠比配到連接點的reads（至少重疊6nt）進行計數(shù)。

7 差異表達分析

    使用DESeq算法過濾差異表達的基因，F(xiàn)C和FDR作為篩選標準，標準如下(a) FC > 2 or < 0.5; (b) FDR< 0.05。

8 GO富集分析

    對差異表達基因（mRNA）進行GO富集分析，以了解基因參與生物學過程。首先下載GO注釋信息，NCBI (http://www.ncbi.nlm.nih.gov/)、UniProt (http://www.uniprot.org/)、GO (http://www.geneontology.org/)，將這些基因的注釋信息進行整合，并作為背景信息。GO分析使用Fisher’s exact test計算p-value，使用FDR值對p-value進行校對。

9 通路分析

    使用KEGG數(shù)據(jù)庫進行通路分析以鑒定有意義的通路，同樣使用Fisher’s exact test計算p-value，使用FDR值對p-value進行校對。FDR<0.05的認為有意義。

10 GO tree

    GO tree是一種有向無環(huán)圖，條目之間彼此關聯(lián)，可以提供友好的數(shù)據(jù)整理和展示。我們挑選了有意義的GO條目（P-value<0.01）構(gòu)建GO tree來了解實驗影響的功能。

11 通路網(wǎng)絡

    KEGG數(shù)據(jù)庫包含代謝、膜轉(zhuǎn)運、信號轉(zhuǎn)導、細胞周期通路。我們選擇富集到的生物學過程的基因，使用CYTOSCAPE(version 3.2.1)進行通路網(wǎng)絡化展示。

12 共表達網(wǎng)絡

    我們使用基因共表達網(wǎng)絡，尋找基因之間的關系，基因共表達網(wǎng)絡構(gòu)建基于標準化的基因表達量。對于每兩個基因，我們計算Pearson相關性系數(shù)，選擇相關性高的成對基因構(gòu)建網(wǎng)絡。在網(wǎng)絡分析中，degree cetrality（度中心性，一個節(jié)點連接其他節(jié)點的數(shù)量為度）是最簡單也是最重要的一種分析方法。度中心性是某個點與相連接點的數(shù)量。為了研究網(wǎng)絡的權重，使用W-core理論來進行圖形拓撲分析，W-core網(wǎng)絡是一個子網(wǎng)絡，在這個網(wǎng)絡中節(jié)點至少與W個其他節(jié)點相連接。W-core網(wǎng)絡是一個protein-protein網(wǎng)絡，通常是緊密相關的蛋白。網(wǎng)絡分析的目的是為了確定核心調(diào)控因子，在一個網(wǎng)絡中，核心調(diào)控因子與許多其他基因相連接，并有最高的度。考慮到不同的網(wǎng)絡，核心調(diào)控因子要考慮在多個網(wǎng)絡中存在，并具有較高的差異。

13 LncRNA-protein關系

    使用NPInter(http://www.bioinfo.org/NPInter)建立lncRNA與蛋白質(zhì)之間的網(wǎng)絡，這個數(shù)據(jù)庫包含了經(jīng)過驗證的noncoding RNA(except tRNA and rRNA) 與生物分子(proteins, RNA and DNA)的關系。這些關系是根據(jù)同行評審的論文人工收集的，許多來源于NONCODE (http://www.noncode.org/),miRBase (http://www.mirbase.org/), and UniProt (http://www.uniprot.org/)。NPInter 根據(jù)ncRNA參與功能的過程對功能關系進行了分類。同時也提供高效的搜索功能，包括關系的搜索，相關論文及其他的信息。

3 結(jié)果

1 RNAseq概述

    數(shù)據(jù)分析流程如Fig. 1。肺腺癌組織及癌旁組織提取total RNA，然后去除 rRNA。首先將reads比對到基因組分析mRNA, lncRNA的表達量，然后將 unmapped reads進行circRNA預測以鑒定circRNA表達譜，然后分析lncRNA, circRNA, mRNA之間的關系。最后，對差異表達的基因，進行功能富集分析，以了解其功能。大多數(shù)reads比對到exon, intron, and TSS，5’UTR的最少（fig. 2A）。Chr18,chr14,chr1比對到的reads數(shù)最多，同時有些reads比對到線粒體（fig. 2B）。figure 2C概括了lncRNA, circRNA, mRNA差異表達的情況。

2 編碼基因表達譜

    共檢測到20384個mRNA, 1282個顯著差異表達(Fig. 3A)。GO分析顯示，差異表達基因主要涉及到分子功能如calciumion biding（鈣離子結(jié)合）, receptor activity（受體激活）, transmembrane signaling receptor activity（跨膜信號受體激活）, 及其他重要功能(Fig. S1A)。GO細胞組分分析主要包括，extracellular matrix （細胞外基質(zhì)）and plasma membrane（質(zhì)膜）(Fig. S1A)。GO生物學過程顯著富集到cell adhesion（細胞黏著）, angiogenesis（血管生成）, extracellular matrix organization（胞外機制構(gòu)成）, and cell surface receptor signaling pathways（細胞表面受體通路）(Fig. 3B)。上述過程均與癌轉(zhuǎn)移相關。對生物學過程進一步GO tree分析發(fā)現(xiàn)，細胞-細胞黏著，血管生成是核心過程(Fig. S1B)。KEGG顯著富集到cGMP-PKG signaling pathway, cell adhesion molecules (CAMs), neuroactive ligand–receptor interaction(Fig. 3C)。Calcium signaling pathway, cell adhesion molecules, PI3K-Akt signaling pathway, and focal adhesion是核心通路(Fig. 3D)。

3 肺腺癌lncRNA表達

    lncRNA共鑒定到19023個，顯著差異的有244個(Fig. 4A; Table S3)。比對到轉(zhuǎn)錄組的reads類型如Fig. 4B, 可以看出ncRNA reads大概占據(jù)所有比對reads的三分之一。AFAP1-AS1, BLACAT1, and LOC101928245是上調(diào)最顯著的lncRNA。FENDRR, LOC400406, and LINC00842是下調(diào)最顯著lncRNA。我們組及其他研究人員曾報道AFAP1-AS1、BLACAT1在肺癌中扮演重要角色。FENDRR在肺癌中的作用也曾報道過。LncRNA的異常表達，常與臨床特征，病人生存率密切相關。高表達的AFAP1-AS1、BLACAT1與肺癌較差預后密切相關。通過一些網(wǎng)上工具，我們發(fā)現(xiàn)POU6F2-AS1,LOC101929398, and LOC101928612與肺癌的生存率相關(Fig. 4C-E)。提示lncRNA可能在肺癌中扮演重要角色。

    我們對lncRNA-protein關系進行預測，結(jié)果提示lncRNA與蛋白質(zhì)ELAVL1(Fig. S2A)相關。構(gòu)建lncRNA-mRNA網(wǎng)絡進而預測lncRNA的生物學過程?？梢园l(fā)現(xiàn)肺癌的共表達網(wǎng)絡與癌旁共表達網(wǎng)絡是不同的。lncRNA BALCAT1, LOC101929398, and MYO16-AS1是關鍵的節(jié)點，而且擁有網(wǎng)絡最高的度?？紤]到BALCAT1 and LOC101929398與肺癌預后的關系，推測，網(wǎng)絡中度較高的lncRNA，可能在肺癌發(fā)生過程起重要的功能。

4 肺腺癌circRNA表達

    circRNA預測算法共鑒定到9340個circRNA(≥ 2 backspliced reads)；與數(shù)據(jù)庫報道的相比，5750個circRNA是已知的，3590個novel circRNA是本次研究鑒定到的。新鑒定到的circRNA，在Appendix S1。大部分的circRNA長度<1000nt（6579個circRNA的長度小于1000nt）,中值是530nt(Fig. 5A)，與以往的報道一致。circRNA表達豐度要小于mRNA和lncRNA。鑒定到112543 backspliced reads，其中超過一半的circRNA的backspliced reads <10(Fig. 5B)。一個來源基因可以產(chǎn)生多個circRNA(Fig. 5C), MACF1可以產(chǎn)生多達26個circRNA。另外，許多基因只產(chǎn)生一個circRNA, 約有1/4的circRNA是其來源基因僅有的環(huán)狀RNA產(chǎn)物(2307 of 9340 circRNA, Fig. 5C)。

    為了確定backspliced reads是否確實是環(huán)狀的，而非普通轉(zhuǎn)錄本剪切的產(chǎn)物，對四個差異表達交大的circRNA，設計divergent primers（發(fā)散性引物），擴增A549 cell cDNA(Fig. 6A)，使用Sanger測序確定其連接位點(Fig. 6B)。Fold change (FC) >2 或 <0.5, flase discovery rate (FDR) <0.05, 共發(fā)現(xiàn)56個差異表達的circRNA(Table S4)。circRNA與mRNA共表達分析以預測circRNA潛在生物學功能。功能相似或位于相同通路的基因通常具有相似的表達模式。一些重要的癌基因與circRNA共表達，如Wnt3A, CDK1,and BUB1，說明這些circRNA可能參與肺癌的發(fā)生過程(Fig. S3)。

4 討論

    雖然發(fā)現(xiàn)了幾十年，直到近些年，非編碼RNA才受到廣泛關注。肺癌的lncRNA表達譜也被廣泛報道，但肺癌中circRNA的表達譜卻研究的很少。本研究，我們使用RNAseq rRNA去除的方式深入分析了肺腺癌中的mRNA, lncRNA, circRNA表達。發(fā)現(xiàn)了244個差異表達的lncRNA，及3590個新的circRNA, 并有56個差異表達的circRNA。對于circRNA形成，人們提出過幾種機制。Bachmayr et al.發(fā)現(xiàn)，circRNA的表達豐度與增值呈負相關，推斷circRNA更容易在非增值細胞中發(fā)生積累。在研究上皮間質(zhì)轉(zhuǎn)化時，Conn發(fā)現(xiàn)許多高豐度circRNA，受到剪切因子的調(diào)控，提示circRNA的形成，會受到RNA結(jié)合蛋白的調(diào)控。對于外顯子來源circRNA, 互補序列，比如在兩端內(nèi)含子上Alu elements，上下游內(nèi)含子會與蛋白質(zhì)結(jié)合，促使circRNA環(huán)化。本研究我們發(fā)現(xiàn)大部分的circRNA來源于前體mRNA外顯子和內(nèi)含子的環(huán)化，這與之前報道的大部分circRNA來源于編碼序列的報道一致。對circRNA 及來源基因的數(shù)量分析，發(fā)現(xiàn)一個基因可以產(chǎn)生多個circRNA。我們觀察到MACF1基因可以產(chǎn)生26個不同的circRNAs，之前有報道PTK2基因能夠產(chǎn)生47個不同的circRNA。大部分circRNA的長度在200-800nt之間，與之前報道的circRNA中值的長度500nt一致。眾多證據(jù)表明，circRNA在多種生物學過程和疾病中發(fā)揮作用，在大部分功能中circRNA發(fā)揮microRNA海綿吸附作用。一個非常夸張的例子是，環(huán)狀RNA CDR1as在腦組織中具有miR-7超過60個保守結(jié)合位點。為了預測高豐度表達circRNA的功能，我們構(gòu)建circRNA-mRNA網(wǎng)絡。許多癌基因與circRNA共表達，比如Wnt3a, and CDK1。由此推斷circRNA或許參與這些生物學過程。對真核RNA的研究，通常會先進行polyA富集。這會失去一部分沒有polyA尾的lncRNA。本研究我們使用rRNA去除RNAseq建庫方式，來研究肺腺癌lncRNA表達譜，這有助于我們檢測到?jīng)]有polyA尾的lncRNA。據(jù)我們所知，這是第一次使用rRNA去除的RNAseq技術研究肺腺癌。AFAP1-AS1, BLACAT1, and FENDRR是肺腺癌差異最明顯的lncRNA, 這與我們之前的工作和別人的報道一致，說明RNAseq技術是比較可靠的。構(gòu)建lncRNA-mRNA共表達網(wǎng)絡來預測lncRNA的功能。在網(wǎng)絡中，BALCAT1 and LOC101929398擁有的度最高，位于網(wǎng)絡的中心。BALCAT1之前報道是肺癌的一個致癌基因。Kaplan–Meier曲線顯示高表達的LOC101929398肺癌差預后相關。這些說明，LOC101929398可能是一個原癌基因。通過共表達網(wǎng)絡，我們可能鑒定到在肺腺癌中發(fā)揮作用的lncRNA。

5 結(jié)論

    通過rRNA去除的RNAseq技術，我們分析了肺腺癌mRNA, lncRNA, circRNA表達譜。AFAP1-AS1，BLUADAT1，LOC101929398是上調(diào)最明顯的lncRNA。我們檢測到9340個circRNA，包括3570個新的circRNA，說明circRNA在肺腺癌中廣泛表達。

微信掃描二維碼關注公眾號獲取更多相關文獻信息