全長轉(zhuǎn)錄組測序技術(shù)解析不同轉(zhuǎn)移性能肝癌細(xì)胞系的轉(zhuǎn)錄本表達(dá)譜與結(jié)構(gòu)變異

楊淑欣，彭民武，羅娟娟，徐景祥，萬紹貴

（贛南醫(yī)學(xué)院基礎(chǔ)醫(yī)學(xué)院分子病理中心，江西 贛州 341000）

肝癌是一類發(fā)病率高、死亡率高的惡性腫瘤，也是全球范圍內(nèi)最常見的癌癥之一，而肝細(xì)胞癌是原發(fā)性肝癌中最主要的亞型［1］。中國是全球肝癌發(fā)病率最高的國家［2］，同時肝癌在我國常見腫瘤中發(fā)病率位居第二［3］。盡管肝細(xì)胞癌的治療取得了巨大進(jìn)展，但由于術(shù)后復(fù)發(fā)和轉(zhuǎn)移率高，5年生存率僅為15%~30%［4］。常見的肝內(nèi)和肝外轉(zhuǎn)移是導(dǎo)致肝細(xì)胞癌患者臨床預(yù)后不良的主要原因［5］。因此，揭示肝癌轉(zhuǎn)移的分子機制對改善其臨床治療具有重要意義。

可變剪切是一種關(guān)鍵的轉(zhuǎn)錄后基因表達(dá)調(diào)控方式，有助于擴大蛋白質(zhì)種類復(fù)雜性和調(diào)控mRNA代謝［6-7］，研究報道m(xù)RNA的可變剪切在多種腫瘤中能夠改變蛋白質(zhì)組的多樣性，異常的可變剪切參與腫瘤細(xì)胞增殖和轉(zhuǎn)移的調(diào)控過程，在腫瘤發(fā)生、發(fā)展、藥物治療及耐藥等方面發(fā)揮重要作用［8-10］。ARHGEF2基因的可變剪切轉(zhuǎn)錄本V1和V3變體與肝癌的轉(zhuǎn)移和預(yù)后密切相關(guān)，提示可變剪切調(diào)控可能是肝癌轉(zhuǎn)移的重要分子機制［11］。

高通量測序技術(shù)的發(fā)展促進(jìn)了關(guān)于腫瘤基因變異的相關(guān)研究，該技術(shù)為研究者在基因表達(dá)差異、可變剪切和新轉(zhuǎn)錄本鑒定等研究提供極大的便利［12-13］。盡管高通量測序以前所未有的準(zhǔn)確度和深度提供了巨大的測序能力，然而，傳統(tǒng)的二代測序技術(shù)由于短讀長的特點，其測序數(shù)據(jù)需要通過生物信息學(xué)算法拼接來分析可變剪接，因此無法完全準(zhǔn)確地反映全長剪接轉(zhuǎn)錄本的序列和定量表達(dá)情況［14-15］，這一缺點將限制其在腫瘤中對可變剪接的大規(guī)模研究和轉(zhuǎn)化應(yīng)用。最新的三代測序技術(shù)克服了二代測序技術(shù)的不足，比如納米孔測序技術(shù)，具有長讀長的優(yōu)勢，平均測序長度能夠達(dá)到10 kb［16］。長讀長測序技術(shù)能夠更加全面準(zhǔn)確地識別新轉(zhuǎn)錄本、可變剪切和基因融合等結(jié)構(gòu)變異，并進(jìn)行轉(zhuǎn)錄本定量分析［17-18］。

本研究利用納米孔三代測序技術(shù)對兩株不同轉(zhuǎn)移潛能的肝癌細(xì)胞系MHCC97H和MHCC97L進(jìn)行全長轉(zhuǎn)錄組測序分析，旨在從轉(zhuǎn)錄本水平上系統(tǒng)地揭示肝癌轉(zhuǎn)移可能存在的分子機制。

1 材料與方法

1.1 細(xì)胞系本實驗使用的高轉(zhuǎn)移潛能肝癌細(xì)胞系（MHCC97H）及低轉(zhuǎn)移潛能肝癌細(xì)胞系（MHCC97 L）由贛南醫(yī)學(xué)院血管生物學(xué)研究平臺提供。

1.2 全長轉(zhuǎn)錄組文庫構(gòu)建及上機測序利用RNA simple總RNA提取試劑盒（天根生化科技北京有限公司）分別提取MHCC97H及MHCC97L細(xì)胞系總RNA，每株細(xì)胞做3個重復(fù)，共6個樣本。然后在Qsep-100 Advance毛細(xì)管電泳儀上使用RNA卡夾對RNA樣本進(jìn)行完整度質(zhì)控檢測。利用Maxima H Minus Reverse Transcriptase試劑盒（Thermo Fisher公司，美國）進(jìn)行反轉(zhuǎn)錄，并使用cDNA-PCR測序試劑盒（SQK-PCB109，Oxford Nanopore Technologies公司）制備cDNA測序文庫，具體步驟按照試劑盒說明書進(jìn)行，簡要過程如下：逆轉(zhuǎn)錄得到的cDNA加Switch Oligo，再合成互補鏈，然后經(jīng)DNA末端修復(fù)加A尾，AMPure XP磁珠純化并添加測序接頭。將上述cDNA測序文庫上樣到Nanopore測序芯片（FLO-MIN106D，R9.4.1）中，在Nanopore MinION平臺上通過MinKNOW2.2軟件操作進(jìn)行全長轉(zhuǎn)錄組測序。

1.3 數(shù)據(jù)過濾與生信分析Nanopore MinION測序下的原始數(shù)據(jù)（Raw data）格式為包含所有原始測序信號的FAST5格式，通過MinKNOW 2.2軟件包中的Guppy軟件對原始數(shù)據(jù)進(jìn)行電流信號到堿基序列信息的轉(zhuǎn)換，并將數(shù)據(jù)轉(zhuǎn)換為FASTQ格式。委托北京百邁客生物科技有限公司對上述FASTQ格式數(shù)據(jù)進(jìn)行過濾與生信分析，進(jìn)一步過濾短片段和低質(zhì)量的原始數(shù)據(jù)，以及去除接頭序列得到有效數(shù)據(jù)（Clean data），過濾條件為所測序列長度＞500 bp，測序質(zhì)量Q score＞7。基于Oxford Nanopore Technologies單分子實時測序技術(shù)的全長轉(zhuǎn)錄組測序無須打斷RNA片段，反轉(zhuǎn)錄得到全長cDNA。該平臺的超長讀取包含了單條完整轉(zhuǎn)錄本序列信息，后期分析無需組裝，所測即所得。根據(jù)所測轉(zhuǎn)錄本與參考基因組的進(jìn)行比對分析，將比對到參考基因組gtf文件未注釋區(qū)域的序列定義為新基因。真核生物的基因往往具有多個轉(zhuǎn)錄本，將由可變剪切等導(dǎo)致的結(jié)構(gòu)不同的轉(zhuǎn)錄本定義為新轉(zhuǎn)錄本。

1.3.1 轉(zhuǎn)錄本表達(dá)量分析轉(zhuǎn)錄組測序可以模擬成一個隨機抽樣的過程，為了讓片段數(shù)目能真實地反映轉(zhuǎn)錄本表達(dá)水平，需要對樣品中Mapped Reads的數(shù)量進(jìn)行歸一化。采用CPM（counts per million）［19］作為衡量轉(zhuǎn)錄本或基因表達(dá)水平的指標(biāo)，CPM計算公式如下：CPM=reads mapped to transcript/total reads aligned in sample×1 000 000（“reads mapped to transcript”表示比對到某一轉(zhuǎn)錄本上的reads數(shù)，“total reads aligned in sample”表示比對到參考轉(zhuǎn)錄組的片段總數(shù)）。使用edgeR R package（3.8.6）進(jìn)行兩株細(xì)胞間的差異表達(dá)分析。使用Benjamini和Hochberg軟件計算PValue和FDR（False discovery rate），將Pvalue＜0.01和差異倍數(shù)≥1.5的基因認(rèn)定為差異表達(dá)基因。

1.3.2 差異表達(dá)轉(zhuǎn)錄本KEGG注釋對差異表達(dá)基因的信號通路注釋分析有助于進(jìn)一步解讀基因的功能。KEGG（Kyoto Encyclopedia of Genes and Genomes，http：//www.genome.jp/kegg/）是系統(tǒng)分析基因功能、基因組信息的數(shù)據(jù)庫，它有助于研究者把基因及表達(dá)信息作為一個整體網(wǎng)絡(luò)進(jìn)行研究。本研究利用KOBAS［20］軟件進(jìn)行差異表達(dá)轉(zhuǎn)錄本的KEGG信號通路注釋。

1.3.3 可變剪切分析通過Astalavista軟件［21］獲取每個樣品存在的可變剪切類型，主要的基因可變剪切類型主要包括外顯子跳躍（Exon skipping，ES）、內(nèi)含子保留（Intron retained，IR）、可變5?剪切位點（Alternative 5?splice-site，A5S）、可 變3?剪切 位 點（Alternative 3?splice-site，A3S）、外 顯 子 互 斥（Mutually exclusive exons，MEE），從Astalavista軟件分析結(jié)果中，對轉(zhuǎn)錄本發(fā)生上述5種可變剪切事件情況進(jìn)行統(tǒng)計。

1.3.4 融合基因分析融合基因是指將兩個或多個基因的編碼區(qū)首尾相連，置于同一套調(diào)控序列（包括啟動子、增強子、核糖體結(jié)合序列及終止子等）控制之下，構(gòu)成的嵌合基因。融合基因的表達(dá)產(chǎn)物為融合蛋白。使用Tofu（版本：13.0.0；參數(shù)：default）比對及尋找融合轉(zhuǎn)錄本，檢測融合轉(zhuǎn)錄本的分析原理有：（1）比對到2個或多個位點；（2）每個位點必須比對上至少5%的轉(zhuǎn)錄本長度，最小比對長度為1 bp；（3）所有位點比對到的總長度必須占轉(zhuǎn)錄本總長度的95%以上；（4）兩位點間距離必須達(dá)到10 kb以上。

2 結(jié) 果

2.1 全長轉(zhuǎn)錄組數(shù)據(jù)質(zhì)控及分析通過Nanopore MinION平臺測序，在兩株細(xì)胞系（各3個重復(fù)樣本）的6個文庫中共得到5 801 060條有效reads，N50和平均讀長分別為985 bp和815 bp，最大讀長達(dá)到134 443 bp。過濾核糖體RNA后的reads數(shù)為4 284 539條，全長序列為3 796 820條，全長序列占總mRNA reads數(shù)的88.6%，詳細(xì)信息見表1。

表1 全長序列數(shù)據(jù)統(tǒng)計表

2.2 差異基因與轉(zhuǎn)錄本表達(dá)分析本研究中兩株細(xì)胞共注釋到9 807個基因，包括643個已有參考基因組中未注釋的新基因。獲得28 532條轉(zhuǎn)錄本，其中已知轉(zhuǎn)錄本為26 061條，及與已有轉(zhuǎn)錄本具有結(jié)構(gòu)差異的轉(zhuǎn)錄本2 471條。在MHCC97H vs.MHCC97L的比較分析中發(fā)現(xiàn)293個差異表達(dá)基因，包括9個新注釋的基因。這些差異表達(dá)基因中，有116個基因在MHCC97H中高表達(dá)，剩余177個基因則在MHCC97H中低表達(dá)（圖1A）。同時在轉(zhuǎn)錄本水平比較分析發(fā)現(xiàn)兩株細(xì)胞共有74條差異轉(zhuǎn)錄本，其中10條為新發(fā)現(xiàn)的轉(zhuǎn)錄本。這些差異表達(dá)的轉(zhuǎn)錄本中有22條在MHCC97H中高表達(dá)，剩余52條則在MHCC97H中低表達(dá)（圖1B），其中ITM2A的兩個轉(zhuǎn)錄本僅在MHCC97H細(xì)胞中表達(dá)，而IL24的四個轉(zhuǎn)錄本僅在MHCC97L細(xì)胞中表達(dá)（圖2）。

圖1 MHCC97H和MHCC97L細(xì)胞的（A）差異表達(dá)基因熱圖和（B）差異表達(dá)轉(zhuǎn)錄本熱圖

圖2 部分差異表達(dá)轉(zhuǎn)錄本在MHCC97H和MHCC97L細(xì)胞中的表達(dá)差異

差異表達(dá)轉(zhuǎn)錄本的KEGG通路注釋發(fā)現(xiàn)，大部分差異表達(dá)轉(zhuǎn)錄本被注釋在Human Diseases這個類別的腫瘤相關(guān)信號通路中（圖3），其中有7條轉(zhuǎn)錄本注釋到癌癥蛋白多糖（Proteoglycans in cancer）通路，分別有5條轉(zhuǎn)錄本注釋到結(jié)直腸癌（Colorectal can?cer）、癌癥信號通路（pathways in cancer）及膠質(zhì)瘤（Glioma）等通路。

圖3 MHCC97H和MHCC97L細(xì)胞中差異表達(dá)轉(zhuǎn)錄本的KEGG數(shù)據(jù)庫注釋結(jié)果

2.3 可變剪切轉(zhuǎn)錄本分析鑒定在MHCC97H和MHCC97L細(xì)胞中共鑒定到1 008次可變剪切事件，其中MHCC97H細(xì)胞中鑒定到619次可變剪切事件，包括83次3’端可變剪切（Alternative 3'splice site，A3S）、87次5’端可變剪切（Alternative 5'splice site，A5S）、403次外顯子跳躍（Exon skipping，ES）、28次內(nèi)含子保留（Intron retention，IR）及18次外顯子互斥（Mutually exclusive exon，MEE）；MHCC97L細(xì)胞中鑒定到870次可變剪切事件，包括120次3’端可變剪切（Alternative 3'splice site，A3S）、107次5’端可變剪切（Alternative 5'splice site，A5S）、572次外顯子跳躍（Exon skipping，ES）、42次內(nèi)含子保留（Intron retention，IR）及29次外顯子互斥（Mutually exclusive exon，MEE）。兩種細(xì)胞中各種可變剪切事件的占比基本一致（圖4），可變剪切事件的具體信息詳見表2。

表2 部分基因的可變剪切信息概要

圖4 MHCC97H和MHCC97L細(xì)胞可變剪切分析結(jié)果

2.4 融合基因分析在MHCC97H和MHCC97L細(xì)胞系中，共鑒定出5條融合轉(zhuǎn)錄本（表3），這5條融合轉(zhuǎn)錄本均為新發(fā)現(xiàn)的融合突變類型，之前未見文獻(xiàn)報道。其中包括肝癌中常見異?；駽TNNB1的融合轉(zhuǎn)錄本及兩個新轉(zhuǎn)錄本ONT.5884與ONT5894形成的融合基因。

表3 融合轉(zhuǎn)錄本統(tǒng)計結(jié)果

3 討論

納米孔測序技術(shù)因具有超長讀長的特點，使得其在全長轉(zhuǎn)錄組測序方面表現(xiàn)出獨特優(yōu)勢。本研究通過納米孔三代測序技術(shù)對具有不同轉(zhuǎn)移潛能的肝細(xì)胞癌細(xì)胞系MHCC97H和MHCC97L進(jìn)行全長轉(zhuǎn)錄組分析，共得到了5 801 060條Clean reads，N50和平均讀長分別為985 bp和815 bp，最大讀長達(dá)到134 443 bp，體現(xiàn)出納米孔測序技術(shù)在鑒定肝癌細(xì)胞全長轉(zhuǎn)錄本方面的優(yōu)勢。通過生物信息學(xué)分析共注釋到9 807個基因，包括643個新基因，獲得28 532條轉(zhuǎn)錄本，其中已知轉(zhuǎn)錄本26 061條，本研究中鑒定出的新轉(zhuǎn)錄本2 471條。在這兩種不同轉(zhuǎn)移潛能的肝癌細(xì)胞系中存在著74條差異轉(zhuǎn)錄本，包括10條新發(fā)現(xiàn)的轉(zhuǎn)錄本。在結(jié)構(gòu)變異分析中發(fā)現(xiàn)了1 008次可變剪切轉(zhuǎn)錄事件，其中占比最多的是外顯子跳躍，這個結(jié)果與CHEN等利用Pacibio三代測序技術(shù)在肝細(xì)胞癌樣本和MIHA細(xì)胞中的研究結(jié)果一致［11］。本研究結(jié)果為在轉(zhuǎn)錄本水平的基因表達(dá)調(diào)控參與肝細(xì)胞癌轉(zhuǎn)移分子機制提供了初步的研究思路與線索。

由于二代測序技術(shù)短讀長的限制，無法對轉(zhuǎn)錄本進(jìn)行結(jié)構(gòu)精確定量和差異表達(dá)分析。同一個基因轉(zhuǎn)錄而來的前體mRNA（pre-mRNA）通過可變剪切可形成不同的剪接異構(gòu)體，最終形成不同的蛋白質(zhì)而發(fā)揮不同的功能［22］。三代測序技術(shù)為深入研究轉(zhuǎn)錄本結(jié)構(gòu)提供了強大的工具［23］，基于二代短讀長測序數(shù)據(jù)只能進(jìn)行基因表達(dá)量的計算和差異表達(dá)分析，但基于三代長讀長測序數(shù)據(jù)不僅能夠同時進(jìn)行基因和轉(zhuǎn)錄本表達(dá)量的計算和差異表達(dá)分析，還能對基因和轉(zhuǎn)錄本的結(jié)構(gòu)進(jìn)行精確的可變剪切分析?？勺兗羟凶鳛橐环N轉(zhuǎn)錄后的調(diào)控機制，在腫瘤的發(fā)生發(fā)展中經(jīng)常出現(xiàn)異常［24］，比如促進(jìn)腫瘤進(jìn)展和轉(zhuǎn)移的特異性可變剪切體［25-26］。與正常肝臟相比，肝細(xì)胞癌腫瘤組織中存在高度差異性可變剪切，其中許多可變剪切差異與肝細(xì)胞癌患者的存活率密切相關(guān)［27-28］。CHEN等通過Pacbio三代長讀長測序技術(shù)對肝細(xì)胞癌患者樣本和MIHA細(xì)胞系進(jìn)行全長轉(zhuǎn)錄組分析鑒定出了肝癌細(xì)胞特有的isoform，研究結(jié)果體現(xiàn)了三代測序技術(shù)在鑒定可變剪切事件上的優(yōu)勢［11］。本研究通過納米孔測序，在全長轉(zhuǎn)錄組水平全面系統(tǒng)地刻畫了MHCC97H及MHCC97L細(xì)胞的可變剪切事件，為后續(xù)進(jìn)一步研究可變剪切轉(zhuǎn)錄本的功能及肝癌轉(zhuǎn)移分子機制提供了重要線索和數(shù)據(jù)信息。

本研究通過全長轉(zhuǎn)錄組測序更加精確地比較不同轉(zhuǎn)錄本在樣本間的差異。MHCC97H與MHCC97L細(xì)胞是從MHCC97人肝癌細(xì)胞株再次分離培養(yǎng)篩選得到的具有高低轉(zhuǎn)移潛能差異的兩株細(xì)胞［29-30］。由于這兩株細(xì)胞均來源于MHCC97，雖然具有不同的轉(zhuǎn)移潛能，但總體差異較小，這可能是本研究僅篩選到74個差異轉(zhuǎn)錄本的原因。在這些差異表達(dá)轉(zhuǎn)錄本中，ITM2A基因的兩個轉(zhuǎn)錄本僅在高轉(zhuǎn)移肝癌細(xì)胞株MHCC97H中表達(dá)，而IL24基因的四個轉(zhuǎn)錄本僅在低轉(zhuǎn)移肝癌細(xì)胞株MHCC97L中表達(dá)。有研究報道ITM2A表達(dá)與肝癌的預(yù)后密切相關(guān)，可能作為肝癌臨床預(yù)后的標(biāo)志物［31］，IL24基因能夠抑制肺癌細(xì)胞的遷移和侵襲［32］。

本研究利用納米孔三代測序技術(shù)對兩種不同轉(zhuǎn)移潛能的肝癌細(xì)胞系進(jìn)行全長轉(zhuǎn)錄組測序，在轉(zhuǎn)錄本水平比較了兩種細(xì)胞的表達(dá)差異，同時揭示了其可變剪切等基因結(jié)構(gòu)變異。研究結(jié)果進(jìn)一步體現(xiàn)出三代長讀長測序技術(shù)在鑒定可變剪切事件上的優(yōu)勢，可變剪切轉(zhuǎn)錄本可能作為肝細(xì)胞癌治療的新型潛在分子靶標(biāo)，同時為后續(xù)進(jìn)一步揭示肝癌轉(zhuǎn)移的分子機制提供了新的思路和線索。