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

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

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

可變剪切是一種關(guān)鍵的轉(zhuǎn)錄后基因表達調(diào)控方式，有助于擴大蛋白質(zhì)種類復(fù)雜性和調(diào)控mRNA代謝［6-7］，研究報道m(xù)RNA的可變剪切在多種腫瘤中能夠改變蛋白質(zhì)組的多樣性，異常的可變剪切參與腫瘤細胞增殖和轉(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ā)展促進了關(guān)于腫瘤基因變異的相關(guān)研究，該技術(shù)為研究者在基因表達差異、可變剪切和新轉(zhuǎn)錄本鑒定等研究提供極大的便利［12-13］。盡管高通量測序以前所未有的準(zhǔn)確度和深度提供了巨大的測序能力，然而，傳統(tǒng)的二代測序技術(shù)由于短讀長的特點，其測序數(shù)據(jù)需要通過生物信息學(xué)算法拼接來分析可變剪接，因此無法完全準(zhǔn)確地反映全長剪接轉(zhuǎn)錄本的序列和定量表達情況［14-15］，這一缺點將限制其在腫瘤中對可變剪接的大規(guī)模研究和轉(zhuǎn)化應(yīng)用。最新的三代測序技術(shù)克服了二代測序技術(shù)的不足，比如納米孔測序技術(shù)，具有長讀長的優(yōu)勢，平均測序長度能夠達到10 kb［16］。長讀長測序技術(shù)能夠更加全面準(zhǔn)確地識別新轉(zhuǎn)錄本、可變剪切和基因融合等結(jié)構(gòu)變異，并進行轉(zhuǎn)錄本定量分析［17-18］。

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

1 材料與方法

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

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

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

1.3.1 轉(zhuǎn)錄本表達量分析轉(zhuǎn)錄組測序可以模擬成一個隨機抽樣的過程，為了讓片段數(shù)目能真實地反映轉(zhuǎn)錄本表達水平，需要對樣品中Mapped Reads的數(shù)量進行歸一化。采用CPM（counts per million）［19］作為衡量轉(zhuǎn)錄本或基因表達水平的指標(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）進行兩株細胞間的差異表達分析。使用Benjamini和Hochberg軟件計算PValue和FDR（False discovery rate），將Pvalue＜0.01和差異倍數(shù)≥1.5的基因認定為差異表達基因。

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

1.3.4 融合基因分析融合基因是指將兩個或多個基因的編碼區(qū)首尾相連，置于同一套調(diào)控序列（包括啟動子、增強子、核糖體結(jié)合序列及終止子等）控制之下，構(gòu)成的嵌合基因。融合基因的表達產(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）兩位點間距離必須達到10 kb以上。

2 結(jié) 果

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

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

2.2 差異基因與轉(zhuǎn)錄本表達分析本研究中兩株細胞共注釋到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個差異表達基因，包括9個新注釋的基因。這些差異表達基因中，有116個基因在MHCC97H中高表達，剩余177個基因則在MHCC97H中低表達（圖1A）。同時在轉(zhuǎn)錄本水平比較分析發(fā)現(xiàn)兩株細胞共有74條差異轉(zhuǎn)錄本，其中10條為新發(fā)現(xiàn)的轉(zhuǎn)錄本。這些差異表達的轉(zhuǎn)錄本中有22條在MHCC97H中高表達，剩余52條則在MHCC97H中低表達（圖1B），其中ITM2A的兩個轉(zhuǎn)錄本僅在MHCC97H細胞中表達，而IL24的四個轉(zhuǎn)錄本僅在MHCC97L細胞中表達（圖2）。

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

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

差異表達轉(zhuǎn)錄本的KEGG通路注釋發(fā)現(xiàn)，大部分差異表達轉(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細胞中差異表達轉(zhuǎn)錄本的KEGG數(shù)據(jù)庫注釋結(jié)果

2.3 可變剪切轉(zhuǎn)錄本分析鑒定在MHCC97H和MHCC97L細胞中共鑒定到1 008次可變剪切事件，其中MHCC97H細胞中鑒定到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細胞中鑒定到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）。兩種細胞中各種可變剪切事件的占比基本一致（圖4），可變剪切事件的具體信息詳見表2。

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

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

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

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

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

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