袁杰，李敏，黃詩圣，舒文杰，任超

1.軍事醫(yī)學研究院 輻射醫(yī)學研究所，北京 100850；2.廣州大學 生命科學學院精準基因編輯工程中心，廣東 廣州 510006；3.上?？萍即髮W 生命科學與技術學院，上海 201210

盡管高通量的短讀長測序極大地促進了轉錄組學研究[1]，然而基于組裝的短讀長測序不足以精確地進行生物信息學分析[2]。此外，先前研究表明真核生物的轉錄組非常復雜。包括可變剪切在內的前體mRNA的轉錄后加工極大地增加了轉錄組的多樣性[3-4]。由于高通量測序的長度限制，它無法完全解析真核生物的轉錄組，特別是與真核生物中普遍存在的新型可變剪切事件相關聯(lián)的復雜性[5-6]。近年出現(xiàn)的第三代測序技術，即長讀長測序，通過使用長讀長技術來實現(xiàn)堿基序列的實時讀取，縮短了測序時間。此外，可以通過直接獲得錯誤率小于1%的全長轉錄本來克服二代短讀長轉錄組測序技術（RNA-seq）的測序長度限制[7]。最近，長讀長測序被用來證明，即使是來自人體器官的高度表征過的轉錄組在基因和同源異構體層面也是不完整的[8-9]。相比之下，使用PacBio單分子實時（single-molecule real-time，SMRT）測序對玉米和草莓的轉錄組進行分析，發(fā)現(xiàn)了許多新穎的剪切同源異構體、長鏈非編碼RNA（long noncoding RNAs，LncRNA）、融合轉錄本，以及新穎的可變剪切事件[10-11]。

等位基因特異性表達（allele-specific expression，ASE）是指二倍體生物體中來自2個等位基因的轉錄本的相對表達水平。等位基因特異性表達可能是由于轉錄速率、mRNA穩(wěn)定性或其他影響轉錄本豐度的機制的不同所造成的[12]。在小鼠早期胚胎發(fā)育過程中，它會經(jīng)歷大規(guī)模重編程過程，以完成母源mRNA的降解和合子基因組激活（zygote genome activation，ZGA）[13]。這些重編程過程可以幫助調節(jié)胚胎基因組轉錄的激活，并為隨后的胚胎發(fā)育和分化奠定基礎[14]。

早期胚胎的起源效應分析可以鑒定出某些顯示出ASE的基因，從而大大增強了我們對早期胚胎發(fā)育過程中重編程的理解。先前的研究表明，包括等位基因特異的剪切和甲基化不對稱性在內的轉錄組和表觀組的起源效應極大地影響著早期胚胎的發(fā)育[15]。然而，使用二代短讀長測序技術從2個等位基因的相對豐度推斷轉錄本來源的方法具有較大的局限性。在最近的研究中，我們利用三代測序技術生成了包括新穎轉錄本在內的更加完整的轉錄組[16]。在本研究中，我們應用三代測序技術描繪了小鼠早期胚胎的剪切圖譜和ASE。我們還將三代與二代RNA-seq數(shù)據(jù)結合在一起，以獲得完整的剪切圖譜，并且利用更完整的剪切信息探究起源效應。我們的目標是更好地表征剪切圖譜和轉錄組的等位基因特異性的注釋信息，以增強對小鼠胚胎發(fā)育的了解。

1 材料和方法

1.1 材料

本研究用到的小鼠早期胚胎的二代測序數(shù)據(jù)、三代測序數(shù)據(jù)、全長轉錄本及三代長讀長測序數(shù)據(jù)識別的可變剪切事件和差異可變剪切事件均來自GSE138760。

1.2 RNA-seq數(shù)據(jù)預處理

首先用 TrimGalore（v0.6.1）修剪 RNA-seq數(shù)據(jù)，隨后用STAR（v2.5.0a）[17]將修剪后的數(shù)據(jù)比對到mm10參考基因組（參數(shù)：--two pass Mode Basic-outSAM typeBAM Unsorted-outSAM strand Field intron Motif）。用 Cufflinks（v2.2.1）[18]拼接各階段的短讀長轉錄本，并且過濾掉不包含正負鏈信息和每千個堿基的轉錄每百萬映射讀取的片段（fragments per kilobase million，F(xiàn)PKM）的值小于1的轉錄本。

1.3 等位基因特異性識別

為了將雜交小鼠（♂ DBA/2×♀C57/BL6J）的序列無偏映射到參考基因組，通過SNPsplit（v0.3.4），采用來自桑格小鼠基因組計劃數(shù)據(jù)庫（dbSNP142）的DBA/2特異的單核苷酸多態(tài)性（single nucleotide polymorphisms，SNPs）生成一套偽基因組。隨后，用GMAP[19]將高質量的三代長讀長轉錄本比對到N掩蓋（N-masked）的基因組，并使用STAR將二代RNA-seq序列比對到N-masked的基因組。根據(jù)在N堿基處的等位基因特異性錯配的SNP數(shù)量判定三代長讀長轉錄本的來源?；贜位點處的堿基，使用SNPsplit將唯一比對上（uniquely mapped）的二代RNA-seq序列分為DBA/2特異組、C57/BL6J特異組和不可分配組。

為了識別等位基因特異的二代短讀長轉錄本，首先過濾掉在任一生物重復樣本中C57特異序列（C57 reads）+DBA特異序列（DBA reads）小于5的轉錄本。等位基因比率計算公式為ReadsC57/（ReadsC57+ReadsDBA）。然后，應用卡方檢驗計算每個重復樣本中等位基因特異的序列偏性的P值，并且用fisher方法將P值合并。依據(jù)等位基因比率和P值來定義等位基因平衡和等位基因不平衡的轉錄本。│等位基因比率-0.5│＜0.16或P≥0.05的轉錄本被定義為等位基因平衡的轉錄本，│等位基因比率-0.5│≥0.16且P＜0.05的轉錄本被定義為等位基因不平衡的轉錄本。

采用2種不同的分析策略，比較三代測序數(shù)據(jù)與二代測序數(shù)據(jù)中的等位基因特異性。長讀長轉錄本的起源與等位基因特異的SNP一致，而短讀長轉錄本的等位基因特異性是由等位基因的比率和偏性決定的。

1.4 計算等位基因特異的基因與同源異構體之間表達模式的相關性

通過 StringTie（v1.3.3b）[20]，使用二代 RNA-seq數(shù)據(jù)計算長讀長基因和同源異構體的總表達量（FPKM）和等位基因特異的表達量（FPKM）。計算等位基因特異的表達量時，輸入的數(shù)據(jù)是等位基因特異的序列。用皮爾遜相關性檢驗計算基因和同源異構體的總表達量以及等位基因特異表達量之間的相關性，相關性系數(shù)（cor）＞0且P＜0.05的模式被定義為一致，其余模式被定義為不一致。

1.5 基因本體（gene ontology，GO）功能富集分析

使用 PANTHER（v14）[21]，對新發(fā)現(xiàn)的等位基因特異的轉錄本和相關性模式為不一致的基因進行GO功能富集分析。

1.6 識別可變剪切事件與差異可變剪切事件

使用SUPPA2（v2.2.1）[22]，在7個階段識別7種類型的可變剪切事件，包括外顯子跳躍（skipping exons，SE）、可 變 的 5′端（alternative 5′splice sites，A5）、可變的3′端（alternative 3′splice sites，A3）、內含子保留（retained introns，RI）、互斥外顯子（mutually exclusive exons，MX）、可變起始外顯子（alternative first exons，AF）、可變末端外顯子（alternative last exons，AL）。由注釋的轉錄本生成的事件被定義為注釋事件，其余事件被定義為新穎事件。差異可變剪切事件是從我們最近發(fā)表的研究中獲取的[16]。我們通過STAR使用二代數(shù)據(jù)計算比對到剪切結上的序列（reads）。與等位基因特異的轉錄本相關聯(lián)的可變剪切事件或差異可變剪切事件被定義為等位基因特異的可變剪切事件或差異可變剪切事件。

1.7 可變剪切事件的累積數(shù)量分析

為了評估三代測序數(shù)據(jù)識別可變剪切事件的潛力，比較了通過二代數(shù)據(jù)和三代與二代數(shù)據(jù)的組合來識別可變剪切事件的累積數(shù)量。首先隨機選取一個階段的二代數(shù)據(jù)識別的轉錄本，接著使用Cuffmerge（v2.2.1）[18]合并單個階段的轉錄本與前面所有階段的轉錄本，隨后用SUPPA2從合并的轉錄本中識別可變剪切事件，當7個階段都被合并時，完成一次循環(huán)。將上述過程重復100次，在每個點計算可變剪切事件數(shù)量的平均值和99%置信區(qū)間。當我們對三代數(shù)據(jù)與二代數(shù)據(jù)的組合進行累積數(shù)分析時，與上述策略唯一的區(qū)別在于分析之前每個階段的長讀長與短讀長轉錄本都已經(jīng)被合并。

2 結果

2.1 使用二代測序數(shù)據(jù)和三代測序數(shù)據(jù)識別等位基因特異的轉錄本

為了表征PacBio三代測序技術在研究早期胚胎轉錄組中等位基因特異的轉錄本和可變剪切事件的優(yōu)勢，我們從我們最近發(fā)表的研究中收集了小鼠早期胚胎7個階段［精子（sperm，SP）、卵母細胞（oocyte，Oo）、1細胞（1-cell，1C）、2細胞（2-cell，2C）、4細胞（4-cell，4C）、8細胞（8-cell，8C）和囊胚（blastocyst，BL）］的測序數(shù)據(jù)和注釋數(shù)據(jù)。這些數(shù)據(jù)包括二代測序數(shù)據(jù)、三代測序數(shù)據(jù)、來自7個階段的全長轉錄本，以及7個階段合并的全長轉錄本。這些數(shù)據(jù)被用來鑒別等位基因特異的轉錄本、可變剪切事件和轉錄組的起源效應（圖1A）。借助三代測序技術的優(yōu)勢，可以將等位基因特異的單核苷酸多態(tài)性SNP精確定位在單個轉錄本中，從而有助于準確識別轉錄本的來源。絕大多數(shù)（97%）的長讀長轉錄本中至少存在1個等位基因特異的SNP，接近一半的轉錄本（46%）至少包含3個等位基因特異的SNP。因此，我們把至少包含3個等位基因特異的SNP的轉錄本定義為等位基因特異的轉錄本。然后，我們基于每個階段獲得的長讀長轉錄本中存在的等位基因特異的SNP來確定長讀長轉錄本的起源。這些轉錄本被分為可以區(qū)分來源的轉錄本［包括C57特異（母源）、DBA特異（父源）以及雙等位基因的轉錄本］和無法區(qū)分來源的轉錄本（圖1B）。與先前研究結果一致[13]，從1細胞到囊胚階段，C57特異的轉錄本所占比例逐漸上升，DBA特異的轉錄本逐漸下降（圖1B）。在7個階段中，我們鑒別出734～1288個等位基因特異的轉錄本（圖1C）。在卵母細胞和精子期，我們觀察到大量的C57特異或DBA特異的轉錄本，然而隨著胚胎發(fā)育，雙等位基因的轉錄本比例逐漸上升（圖1C）。通過與GENCODE（vM20）注釋進行比較，每個階段平均鑒定出532個C57特異的轉錄本和397個DBA特異的轉錄本，參考我們先前的發(fā)現(xiàn)，每個階段平均存在650個注釋的轉錄本和234個新基因和同源異構體（圖1D）[16]。

接著我們比較了三代測序數(shù)據(jù)和二代測序數(shù)據(jù)鑒別的等位基因特異的轉錄本。盡管由于測序深度的原因，三代數(shù)據(jù)中識別的等位基因特異的轉錄本要少于二代數(shù)據(jù)，但不管是對于三代數(shù)據(jù)（圖1B、C）還是二代數(shù)據(jù)，等位基因特異的轉錄本的數(shù)量和比例都是從1細胞到囊胚階段逐漸減少。然而兩者之間的重疊率逐漸降低，在囊胚階段兩者識別等位基因特異的轉錄本僅6%是一致的（圖1E）。在7個階段中，僅被三代數(shù)據(jù)所識別的新發(fā)現(xiàn)的等位基因特異的轉錄本數(shù)量范圍為378～872（圖1F）。我們對這些新發(fā)現(xiàn)的等位基因特異的轉錄本做了GO分析，發(fā)現(xiàn)這些轉錄本與細胞代謝過程、感覺知覺和細胞周期過程存在關聯(lián)。

我們進一步表征了從三代數(shù)據(jù)中鑒別出的新發(fā)現(xiàn)的等位基因特異的轉錄本，將其分為3種類型，即沒有特異的二代數(shù)據(jù)支持（未被等位基因特異的短讀長序列所識別）、等位基因特異的偏性相反（與三代數(shù)據(jù)識別的結果相比，二代數(shù)據(jù)識別結果的親代定位是相反的）、雙等位基因（存在等位基因特異的SNP但等位基因的表達水平無差異）（圖1F）。例如，在囊胚階段識別的基因PB.2249的2個異構體在三代數(shù)據(jù)中被識別為品系特異，但在二代數(shù)據(jù)中被錯誤地鑒定為雙等位基因（圖1G）。這些結果證明了使用三代數(shù)據(jù)鑒別等位基因特異的轉錄本的優(yōu)勢。

圖1 三代測序數(shù)據(jù)與二代測序數(shù)據(jù)定義的等位基因特異的轉錄本

2.2 比較三代測序和二代測序數(shù)據(jù)中的可變剪切事件

為了表征三代數(shù)據(jù)識別可變剪切事件的優(yōu)勢，分別使用三代數(shù)據(jù)和二代數(shù)據(jù)在7個階段識別可變剪切事件?；贕ENCODE注釋，這些事件被分為注釋的事件和新穎的事件（圖2A、B）。與二代數(shù)據(jù)相比，三代數(shù)據(jù)能鑒別出更多罕見的可變剪切事件，例如AF（P=1.3e-18，卡方檢驗），從而證明了三代測序技術捕獲復雜可變剪切事件的能力。此外，從三代數(shù)據(jù)中識別的新穎事件的比例要顯著大于二代數(shù)據(jù)中識別的新穎事件的比例（圖2C）。這些結果證明，三代測序技術對于分析可變剪切事件極具價值。

接著，比較了分別從長讀長轉錄本和短讀長轉錄本中提取的剪切結（splicing junction，SJ）。大多數(shù)剪切結在兩者之中都存在，但仍有數(shù)千個剪切結（平均每個階段7058個）只能被三代數(shù)據(jù)所識別（圖2D）。在這些僅能被三代數(shù)據(jù)所識別的特有的剪切結中，平均每個階段有87%的剪切結能被至少5個短讀長序列所支持（圖2E）。這些結果表明借助三代測序技術，我們能識別出大量高準確度且特異的剪切結。

圖2 使用三代測序數(shù)據(jù)和二代測序數(shù)據(jù)識別可變剪切事件

為了評估二代數(shù)據(jù)與三代數(shù)據(jù)的組合識別可變剪切事件的能力，我們比較了二代數(shù)據(jù)和二代與三代的組合數(shù)據(jù)識別可變剪切事件的累積數(shù)量。在當前的測序深度之下，從組合數(shù)據(jù)中識別出的事件顯著地比僅從短讀數(shù)據(jù)中識別的事件的數(shù)量多（P＜1e-100，威爾遜配對秩和檢驗）（圖2F）。這個結果顯示出利用二代數(shù)據(jù)和三代數(shù)據(jù)的組合識別可變剪切事件的優(yōu)勢。

2.3 小鼠早期胚胎中與起源效應相關的同源異構體的表達

由于早期胚胎中親本轉錄本的動態(tài)變化，可變剪切事件與等位基因特異的轉錄本存在關聯(lián)。為了研究特定于起源的剪切異構體，我們對具有不少于2個轉錄本的基因以及關聯(lián)轉錄本的總表達和等位基因特異的表達進行了定量。通過計算基因與同源異構體表達量之間的相關性，發(fā)現(xiàn)大多數(shù)異構體（73%）的表達模式與對應的基因一致。相反，有27%的異構體與基因表達模式不一致（圖3A）?；诋悩嬻w與異構體等位基因特異性表達的相關性，我們將后者分為4類。我們觀察到大多數(shù)異構體的表達模式與它們在等位基因特異性水平的表達不一致（圖3B、C）。例如，Hsd17b6基因的C57特異性表達與一個新鑒定的Hsd17b6的異構體的表達一致，但是與其在整個基因層面的表達模式不一致；同樣，Trim43a基因的DBA特異性表達與異構體的表達和其在整個基因層面的表達模式都不一致；然而，在ZGA過程中，Mcph1基因的異構體的DBA特異性表達與C57特異性表達模式相反（圖3D、F）。我們對4組表達模式不一致的異構體進行了GO分析，這些基因參與了細胞代謝、基因表達及RNA加工等過程（圖3G）。因此，我們可以通過整合三代測序數(shù)據(jù)和二代測序數(shù)據(jù)來識別階段特異和等位基因特異的轉錄本。

圖3 等位基因特異的基因與轉錄本的表達模式的相關性

2.4 早期胚胎中基于三代測序數(shù)據(jù)的等位基因特異的可變剪切的分析

根據(jù)上述結果（圖1、2），我們鑒定了一些新發(fā)現(xiàn)的等位基因特異的轉錄本，這些轉錄本可能是由等位基因特異的可變剪切事件產(chǎn)生的。因此，我們在三代測序數(shù)據(jù)中分析了等位基因特異的剪切事件的發(fā)生。每個階段平均鑒定出230個等位基因特異的剪切事件（圖4A）。我們觀察到一些包含DBA特異與C57特異的mRNA的剪切事件。例如，在2細胞階段，Tor1aip1基因的一個外顯子跳躍事件產(chǎn)生了DBA特異的轉錄本TCONS_00001249和C57特異的轉錄本TCONS_00001250（無跳躍外顯子）（圖4B）。此外，我們結合三代和二代數(shù)據(jù)來檢測早期胚胎發(fā)育過程中的等位特異的差異可變剪切事件（圖4C）。等位基因特異的差異剪切事件在1細胞到2細胞階段以及4細胞到8細胞階段發(fā)生的頻率更高。仍以Tor1aip1為例，在1細胞階段，僅僅只有C57特異的轉錄本TCONS_00001250存在表達，然而DBA特異的轉錄本TCONS_00001249在1細胞到2細胞的轉變中被激活（圖4D、E）。隨著一個新發(fā)現(xiàn)的DBA特異的差異可變剪切事件的出現(xiàn)（圖4D、E），Tor1aip1基因中的DBA特異的轉錄本和C57特異的轉錄本的表達量都有所上升（圖4F～H）?？偟膩碚f，這些數(shù)據(jù)表明某些親本特異的轉錄本表現(xiàn)為親本特異的剪切異構體，這些等位基因特異的剪切事件和差異可變剪切事件的發(fā)現(xiàn)將改善我們對早期胚胎發(fā)育過程的理解。

圖4 基于三代測序數(shù)據(jù)的等位基因特異的可變剪切

3 討論

在本研究中，我們使用了最近發(fā)表的研究中的高分辨率的轉錄組信息來分析起源效應，并且比較了小鼠早期胚胎的三代測序數(shù)據(jù)和二代測序數(shù)據(jù)識別的剪切事件[16]。我們的目的是探索三代測序技術在研究轉錄組起源效應方面的優(yōu)勢，并且鑒定新穎的可變剪切事件。

由于包含2個長距離的SNP的片段很難被擴增以及無法使用Sanger在整個序列中對800 bp以上的片段進行測序等技術限制，目前暫時無法對三代測序技術新識別的等位基因特異的轉錄本進行準確的實驗驗證。我們建立了一套生物信息學的流程來討論三代測序技術在識別等位基因特異性方面的優(yōu)勢，并且借助這一優(yōu)勢來識別新的等位基因特異的轉錄本，在同源異構體的層面上研究了早期胚胎中等位基因特異的轉錄激活。結果顯示，與基于組裝的二代測序技術相比，三代測序技術能夠更準確地識別等位基因特異的轉錄本和剪切事件。盡管二代數(shù)據(jù)能識別大量等位基因特異的轉錄本，但是考慮到這種基于組裝的策略的可靠性與準確性的不足，我們僅關注于利用三代數(shù)據(jù)的準確性和優(yōu)勢來發(fā)現(xiàn)等位基因特異的轉錄本。隨著母源mRNA在ZGA過程中降解[23]，我們觀察到在1細胞到2細胞轉變過程中C57特異的轉錄本數(shù)量明顯下降。但是即便是到囊胚期，仍有數(shù)百個轉錄本保持著等位基因失衡的狀態(tài)，這表明等位基因特異的轉錄本是胚胎發(fā)育過程中的常規(guī)產(chǎn)物[24]。

三代測序數(shù)據(jù)具有極高的識別新穎可變剪切事件的潛力。我們對可變剪切事件的累積數(shù)量分析表明，將三代數(shù)據(jù)和二代數(shù)據(jù)合并可以快速增加可變剪切事件的數(shù)量。此外，三代數(shù)據(jù)有助于識別大量二代數(shù)據(jù)無法識別的新穎的剪切結。這突顯了三代數(shù)據(jù)在識別新的剪切事件方面的優(yōu)勢。

我們還證明了，很大一部分基因的表達模式與它對應的異構體或它在親本水平的表達都不一致。我們推斷某些DBA或C57偏性的轉錄本在特定階段被特異性激活，它們可能在功能上參與了胚胎發(fā)育過程。因此，迫切需要開發(fā)適當和有效的工具來闡明等位基因特異的轉錄本的功能[25-26]。最近的一項研究通過靶向非甲基化基因座實現(xiàn)了印記基因的等位基因特異性編輯[26]，這可以用于進一步研究等位基因特異性基因在發(fā)育過程中的功能。

當依據(jù)等位基因特異的SNP來分離新穎的剪切異構體時，我們觀察到數(shù)百個品系特異的剪切事件?？傮w而言，與DBA特異的轉錄本相比，識別出更多的C57特異的轉錄本，尤其是在1細胞和2細胞階段。這些發(fā)現(xiàn)與已有的母源基因編碼的轉錄因子激活了ZGA過程的推測一致[24]。在啟動子區(qū)，基因間區(qū)和遠端區(qū)域的等位基因特異的甲基化和母源組蛋白H3第27位賴氨酸三甲基化（H3 lysine 27 trimethylation，H3K27me3）修飾[27-29]，在等位基因特異的可變剪切的分子調節(jié)過程中成為不依賴于DNA甲基化的印記機理[30-32]。這一過程可能控制了等位基因特異的基因表達和可變剪切事件。因此，利用三代測序技術識別父母源特異的DNA甲基化過程或H3K27me3印記現(xiàn)象，有助于進一步解釋等位基因特異的轉錄調控。

總之，我們建議借助三代測序技術的優(yōu)勢準確識別等位基因特異的轉錄本與剪切事件。通過三代測序技術，我們可以獲得更多等位基因特異的轉錄本、新穎可變剪切事件以及更多的剪切結，并且我們還報道了早期胚胎過程中的等位基因特異的可變剪切和差異可變剪切。這些發(fā)現(xiàn)可以加深我們對于早期胚胎發(fā)育過程中轉錄組起源效應的理解。