劉秋芳,李 華

(上海交通大學(xué)生物醫(yī)學(xué)工程學(xué)院,上海200240)

真核生物基因能夠產(chǎn)生大量不同的mRNA產(chǎn)物，例如，在人類已注釋的基因中，每個基因平均能夠產(chǎn)生5個已知的轉(zhuǎn)錄本異構(gòu)體，有些基因甚至能夠產(chǎn)生多達80個轉(zhuǎn)錄本異構(gòu)體(ensembl release 75)[1]。這些轉(zhuǎn)錄本異構(gòu)體是由選擇性轉(zhuǎn)錄起始、可變剪切以及多聚腺苷酸化等轉(zhuǎn)錄調(diào)控過程的共同作用產(chǎn)生的[2-4]。在動物、植物以及其他真核生物中，通過剪接的方式，將同一個基因的外顯子按照多種組合方式拼接在一起形成許多不同的mRNAs[1]。mRNA的表達水平和對應(yīng)蛋白的水平并不一致,目前，有很多科研工作人員致力于研究造成這種差異的原因[5-7]。研究人員在人類的研究中發(fā)現(xiàn)不同的可變剪切形式對蛋白的翻譯水平有顯著的影響，他們發(fā)現(xiàn)結(jié)合到多核糖體上的轉(zhuǎn)錄本異構(gòu)體的3′UTR序列更短[8]?？勺兗羟袑虮磉_的影響目前比較公認的是通過2種方式進行的[8]：一是產(chǎn)生不同的可變剪切轉(zhuǎn)錄本異構(gòu)體，從而增加蛋白的多樣性；二是產(chǎn)生包含提前終止密碼子的轉(zhuǎn)錄本，啟動NMD(nonsense-mediated decay)通路降解該轉(zhuǎn)錄本，從而降低基因的表達水平。核糖體是蛋白質(zhì)的加工場所，能夠?qū)RNA序列通過特定的編碼方式翻譯成蛋白質(zhì)。以前的研究主要著眼于基因水平上的mRNA豐度與蛋白質(zhì)水平的關(guān)系，而忽視了轉(zhuǎn)錄本異構(gòu)體本身與核糖體的結(jié)合在翻譯中所起的作用。利用傳統(tǒng)的核糖體圖譜法可以收集并提取與核糖體結(jié)合的全長轉(zhuǎn)錄本，然后對這些與核糖體結(jié)合的轉(zhuǎn)錄本進行測序，可以得到翻譯組數(shù)據(jù)[9]。通過對細胞質(zhì)和核糖體結(jié)合RNA分別進行測序，可以得到這兩種不同細胞組分的RNA。

本文中，筆者收集果蠅早期胚胎0～4 h的轉(zhuǎn)錄組(細胞質(zhì)RNA)和翻譯組(核糖體結(jié)合RNA)數(shù)據(jù)，并通過系統(tǒng)的比較分析來揭示果蠅不同細胞組分轉(zhuǎn)錄本多樣性。

1 材料與方法

1.1 數(shù)據(jù)

本文中使用的果蠅早期胚胎0～4 h的轉(zhuǎn)錄組和翻譯組數(shù)據(jù)來源于Li等[10]公布的數(shù)據(jù)。

1.2 方法

1.2.1 高通量數(shù)據(jù)分析

對于來自果蠅不同細胞組分的高通量測序數(shù)據(jù)，首先，使用TopHat 2(v2.0.9)分別將轉(zhuǎn)錄組和翻譯組數(shù)據(jù)比對到果蠅的參考基因組上去[11]。接下來，使用Cufflinks(v2.2.1)的Cuffnorm模塊，利用唯一比對的數(shù)據(jù)計算果蠅已知轉(zhuǎn)錄本的豐度，Cuffnorm計算轉(zhuǎn)錄本的表達水平時會將比對上參考基因組的reads數(shù)轉(zhuǎn)化為FPKM值(fragments per kilobase of exon model per million mapped reads)[12]。最后，為了驗證數(shù)據(jù)的可靠性，對來自不同細胞組分的數(shù)據(jù)進行相關(guān)性分析。

1.2.2 主要表達轉(zhuǎn)錄本分析

首先，為了得到不同細胞組分中主要表達的轉(zhuǎn)錄本，筆者分別對轉(zhuǎn)錄組和翻譯組數(shù)據(jù)進行如下處理：一是刪除低表達的基因，在某一細胞組分中，如果一個基因的任意一個轉(zhuǎn)錄本的表達值大于1,則在該細胞組分中保留該基因[13]，否則刪除該基因；二是計算不同細胞組分中主要表達的轉(zhuǎn)錄本；三是比較不同細胞組分中主要表達的轉(zhuǎn)錄本，找出主要表達轉(zhuǎn)錄本發(fā)生改變的基因。其次，對于不同細胞組分中主要表達轉(zhuǎn)錄本發(fā)生改變的基因，根據(jù)這些基因在不同細胞組分中主要表達轉(zhuǎn)錄本的非翻譯區(qū)(UTR)和編碼序列(CDS)是否發(fā)生變化進行分類。最后，將不同細胞組分中主要表達轉(zhuǎn)錄本CDS發(fā)生變化的基因用DAVID(v6.8)進行功能富集分析[14]。

2 結(jié)果與討論

2.1 數(shù)據(jù)比對結(jié)果

為了得到不同細胞組分的轉(zhuǎn)錄本豐度，首先，使用TopHat 2(--library-type=fr-firststrand -G，其他為默認參數(shù))分別將轉(zhuǎn)錄組和翻譯組數(shù)據(jù)比對到果蠅的參考基因組上去。由于測序數(shù)據(jù)會比對到基因組的多個位置，這些比對結(jié)果可能是錯誤的比對，會對后續(xù)的分析造成一定的影響，因此，為了提高后續(xù)轉(zhuǎn)錄本豐度計算的準確性，去除掉多比對的數(shù)據(jù)，后續(xù)分析只使用唯一比對(unique mapping)，數(shù)據(jù)結(jié)果見表1。由表1可以看到細胞質(zhì)和核糖體結(jié)合RNA數(shù)據(jù)比對上的比率分別為94.49%和93.55%，唯一比對數(shù)據(jù)的比率分別為92.05%和90.55%，這說明所選擇的數(shù)據(jù)具有較高的質(zhì)量，能夠很好地比對到參考基因組上去，可以用于后續(xù)的分析。

表1 比對前和比對后的數(shù)據(jù)

注：Cyto表示細胞質(zhì)RNA數(shù)據(jù)(即轉(zhuǎn)錄組數(shù)據(jù))，Poly表示核糖體結(jié)合RNA數(shù)據(jù)(即翻譯組數(shù)據(jù))。

2.2 相關(guān)性分析結(jié)果

為了驗證數(shù)據(jù)的可靠性，對早期胚胎細胞轉(zhuǎn)錄組和翻譯組數(shù)據(jù)進行相關(guān)性分析。在相關(guān)性分析中，要求轉(zhuǎn)錄本的表達水平在2個細胞組分中的FPKM都高于1，以減少測序噪音對相關(guān)性的影響[13]。利用R語言中的cor.test函數(shù)計算轉(zhuǎn)錄組和翻譯組數(shù)據(jù)的pearson相關(guān)系數(shù)，結(jié)果見圖1。

由圖1可知：早期胚胎細胞的轉(zhuǎn)錄組和翻譯組的相關(guān)系數(shù)為0.85，這說明來自不同細胞組分的測序數(shù)據(jù)具有很好的相關(guān)性，符合預(yù)期，因此一定程度上支持了數(shù)據(jù)的可靠性。

圖1 果蠅0～4 h早期胚胎轉(zhuǎn)錄組和翻譯組數(shù)據(jù)相關(guān)性Fig.1 Correlation between transcriptome data and translatome data of 0～4 h embryo cell in Drosophila melanogaster

2.3 不同細胞組分主要表達轉(zhuǎn)錄本分析結(jié)果

真核生物的基因表達水平與蛋白水平并不具有完全的一致性，有些基因通過選擇性剪切能夠轉(zhuǎn)錄得到多個轉(zhuǎn)錄本，每個轉(zhuǎn)錄本可能對應(yīng)著不同的轉(zhuǎn)錄水平，不同的轉(zhuǎn)錄本可能具有不同的翻譯效率。Floor等[1]研究發(fā)現(xiàn)，即使基因水平的表達豐度不變，在轉(zhuǎn)錄本水平上發(fā)生小小的改變也會對蛋白的水平產(chǎn)生巨大的影響。核糖體是轉(zhuǎn)錄本翻譯成蛋白質(zhì)的場所，因此，筆者想要研究是否是細胞質(zhì)和核糖體結(jié)合RNA中轉(zhuǎn)錄本組成差異造成基因與蛋白表達水平的不一致性。

果蠅已知基因數(shù)是17 746個，已知轉(zhuǎn)錄本數(shù)是35 113個，平均每個基因含有2個轉(zhuǎn)錄本。為了得到不同細胞組分中主要表達的轉(zhuǎn)錄本，首先分別對果蠅轉(zhuǎn)錄組和翻譯組數(shù)據(jù)的基因進行篩選，結(jié)果見圖2。

圖2 果蠅早期胚胎不同細胞組分表達的基因重合數(shù)Fig.2 Overlap of expressed gene number of different cellular compartments in Drosophila melanogaster

由圖2可知：在轉(zhuǎn)錄組和翻譯組數(shù)據(jù)中剩余的基因數(shù)分別是7 525個(占所有已知基因的42.4%)和7 525個(占所有已知基因的42.4%)，不同細胞組分的基因重合數(shù)是7 207個(其中表達2個及2個以上轉(zhuǎn)錄本的基因數(shù)為3 650個，占重合基因總數(shù)的50.65%)，表明不同細胞組分中表達的基因具有較高的一致性。然后，通過計算和比較分析，在表達2個及2個以上轉(zhuǎn)錄本的重合基因中找出了主要表達的轉(zhuǎn)錄本發(fā)生變化的基因數(shù)766個。接下來，對這些基因按照不同細胞組分中主要表達的轉(zhuǎn)錄本UTR和CDS是否發(fā)生變化進行分類，UTR和CDS的信息來源于Flybase數(shù)據(jù)庫發(fā)布的r6.13版本的果蠅已知注釋信息(ftp://ftp.flybase.net/genomes/Drosophila_melanogaster/dmel_r6.13_FB2016_05/gtf/)[15]。

經(jīng)分類計算，發(fā)現(xiàn)400個基因的UTR發(fā)生變化，366個基因的CDS發(fā)生變化。最后對CDS發(fā)生變化的基因用DAVID進行GO功能富集分析和KEGG通路分析，結(jié)果見表2～3。由表2可知，這些基因與蛋白結(jié)合(protein binding)、GTP酶激活活性(GTPase activator activity)等分子功能相關(guān)。由表3可知，這些基因參與到果蠅早期胚胎背腹軸的形成過程。

以上結(jié)果表明，基因在細胞質(zhì)和核糖體結(jié)合RNA中主要表達的轉(zhuǎn)錄本發(fā)生了改變。核糖體是蛋白質(zhì)的生產(chǎn)場所，這可以在一定程度上解釋蛋白質(zhì)表達水平和轉(zhuǎn)錄本表達水平的不一致性。

2.4 討論

筆者收集了果蠅早期胚胎細胞的轉(zhuǎn)錄組和翻譯組數(shù)據(jù)，并通過系統(tǒng)的分析展示了不同細胞組分中RNA的組成差異。具體而言，首先利用TopHat 2軟件將轉(zhuǎn)錄組和翻譯組數(shù)據(jù)比對到果蠅的參考基因組上，得到不同細胞組分中的轉(zhuǎn)錄本表達豐度。然后，對于表達多個轉(zhuǎn)錄本的基因，在細胞質(zhì)和核糖體結(jié)合RNA數(shù)據(jù)中分別找到其主要表達的轉(zhuǎn)錄本。最后，對于在2個樣本中都表達的基因，計算了這些基因主要表達轉(zhuǎn)錄本的組成異同。最終，找到了766個主要表達轉(zhuǎn)錄本發(fā)生差異的基因，這表明在果蠅中大量基因的不同轉(zhuǎn)錄本的翻譯效率差異可能很大。通過對基因序列特征的分析，將主要表達轉(zhuǎn)錄本發(fā)生變化的基因分為兩類：UTR變化類(400個基因)和CDS變化類(366個基因)。轉(zhuǎn)錄本CDS的改變意味著編碼得到的蛋白也發(fā)生了改變，這可以一定程度上解釋蛋白質(zhì)表達水平和轉(zhuǎn)錄本表達水平的不一致性。

由于果蠅的基因與人類的基因存在比較大的同源性，對果蠅胚胎發(fā)育的研究成果可以一定程度地應(yīng)用在人類身上，為了解人類的胚胎發(fā)育過程提供一定的依據(jù)，因此，對果蠅不同細胞組分轉(zhuǎn)錄本多樣性的研究是很有價值的。

表2 果蠅早期胚胎細胞中CDS發(fā)生變化的基因功能富集分析結(jié)果

表3 果蠅早期胚胎細胞中CDS發(fā)生變化的基因KEGG通路分析結(jié)果

3 結(jié)論

筆者對果蠅不同細胞組分轉(zhuǎn)錄本組成多樣性進行了初步的研究，對CDS發(fā)生變化的基因進行功能富集分析和KEGG通路分析，發(fā)現(xiàn)這些基因顯著富集到某些分子功能上并在特殊的通路中發(fā)揮功能。另外，UTR是序列上的非編碼區(qū)域，對轉(zhuǎn)錄本的翻譯起著重要的調(diào)控作用。因此，將來有必要對UTR上的調(diào)控元件做進一步的研究。

(責任編輯 荀志金)