mRNA序列與相應(yīng)內(nèi)含子序列匹配的普適性分析

張 強，趙小慶，薄素玲，曹艷娟 ，蘇文霞

(1．內(nèi)蒙古農(nóng)業(yè)大學(xué) 理學(xué)院，呼和浩特 010018; 2．內(nèi)蒙古醫(yī)科大學(xué) 計算機信息學(xué)院， 呼和浩特 010110；3．內(nèi)蒙古自治區(qū)農(nóng)牧業(yè)科學(xué)院，呼和浩特 010031)

近幾年，諸多研究都開始重視內(nèi)含子對基因表達的影響[1]。大量研究表明，內(nèi)含子是一類具有生物學(xué)功能的序列。許多基因表達調(diào)控元件，例如基因轉(zhuǎn)錄和mRNA加工(尤其是可變剪接)，都屬于內(nèi)含子序列的一部分。各種非編碼RNA，例如microRNA和snoRNA，也屬于內(nèi)含子[2]。對內(nèi)含子序列而言，它的丟失和獲得對非編碼RNA變異和基因重組有影響，這是關(guān)系到真核基因進化的主要方面[3-7]。許多疾病的發(fā)生是由于內(nèi)含子的突變造成的[8-9]。內(nèi)含子兩端和中間序列的突變，都是因為激活隱性切割位點進而導(dǎo)致的疾病。雖然部分基因的表達并不需要內(nèi)含子序列的參與甚至不存在內(nèi)含子，然而內(nèi)含子在許多情況下能夠最大限度地增強轉(zhuǎn)基因生物的基因表達[10-11]。內(nèi)含子序列已經(jīng)成為對轉(zhuǎn)基因生物中外源基因表達進行改善不可或缺的一部分[12-15]。以上所有情況間接表明內(nèi)含子序列的存在與否對基因表達有顯著差別。同時，在pre-mRNA上，所有剪接過程的發(fā)生，都可以看作內(nèi)含子和外顯子序列相互作用后，所產(chǎn)生的必然結(jié)果。序列匹配是體現(xiàn)上述相互作用最基本的形式[16-19]。那么必定只能由更加復(fù)雜的序列結(jié)構(gòu)才能組成內(nèi)含子序列并執(zhí)行多種生物功能。比如，內(nèi)含子通過與mRNA的序列匹配，調(diào)節(jié)結(jié)合蛋白因子與mRNA的相互作用，才能進行調(diào)節(jié)[20-22]。在mRNA序列上存在功能區(qū)域，例如翻譯起始與終止位點，還存在與內(nèi)含子匹配的特殊形式，這些形式對基因的表達調(diào)控具有至關(guān)重要的作用[23]。因此，尋找mRNA序列與相應(yīng)內(nèi)含子的最佳匹配區(qū)域是解決問題的思路。

基于此思路，以13個生物基因的編碼序列作為研究樣本，目的在于找到mRNA序列上，能夠與內(nèi)含子序列存在相互作用關(guān)系的最佳匹配片段，進而分析這些片段的序列特征，討論在mRNA序列上匹配頻率的分布特點，研究該分布在不同物種中的進化規(guī)律。

1 數(shù)據(jù)與方法

1.1 基因序列數(shù)據(jù)

13個物種1號染色體編碼基因序列取自The Exon-Intron Database (EID)數(shù)據(jù)庫。13個物種中包括4個植物，2個無脊椎動物和7個有脊椎動物。在選取樣本(見表1)中去掉了含有重復(fù)元件和ncRNA等已知非剪接功能的基因。

1.2 比對方法

mRNA序列與相應(yīng)內(nèi)含子序列存在堿基互補匹配片段，該片段能夠反映它們之間的相互特點。因此它們之間的相互作用可由最佳匹配片段表征出來。首先把內(nèi)含子進行轉(zhuǎn)化，得到它的互補序列，其次通過軟件Smith-Waterman，對其進行局部相似性比對，再次比對互補序列與其所對應(yīng)的mRNA序列，進而找出最佳相似片段，最后通過轉(zhuǎn)化得到兩者最佳匹配片段。

1.3 最佳匹配頻數(shù)分布

定義1：匹配打分函數(shù)

對給定序列的每個堿基位點賦予一個分值，在最佳匹配區(qū)域內(nèi)，則堿基位點賦值為1；如果不在最佳匹配區(qū)域內(nèi)，則賦堿基位點值為0。定義匹配打分函數(shù)如下：

(1)

公式中，i代表第i條序列，j是第i條序列中，第j個堿基位點(j=1，2，…，Li)；那么Li是第i條序列的長度，則Nis和Nie分別是第i條序列上最佳匹配片段的起始和終止堿基位點。這樣一條序列就轉(zhuǎn)換成由0和1組成的數(shù)字串，1代表最佳匹配片段位置。

表1 13種真核生物基因Table 1 Protein coding genes of 13 eukaryotes

定義2：匹配頻率函數(shù)

在所分析的m條序列上，統(tǒng)計第j堿基位點上出現(xiàn)1的次數(shù)，除以m就得到該位點的匹配頻率。匹配頻率函數(shù)F(j)定義如下：

(2)

函數(shù)F是匹配頻率，它代表的是比對序列的第j個堿基位點，與被比對序列兩者的匹配強度，或者是它們相互作用的強度。由此可以推出，根據(jù)F值的大小，可以對此位點參與最佳匹配的概率進行評價。如果F(j)等于1，那么就能夠說明，所有m條序列的第j個堿基位點，應(yīng)該全都位于最佳匹配區(qū)域內(nèi)。由此得到匹配頻率F值在序列上的分布。

定義3：序列長度標(biāo)準(zhǔn)化

實際序列可以是編碼序列，或者是內(nèi)含子序列，針對它們長度各異這一情況，序列的最佳匹配片段要進行相對位置分布比較，為使其比較起來更方便，則需要把序列長度標(biāo)準(zhǔn)化，得到長度為L的目標(biāo)序列。長度標(biāo)準(zhǔn)化采用如下方法：

(3)

分析以上公式，將其中第i條比對序列長度設(shè)置為Li；那么位于第i條比對序列上，第j個堿基位點是Nij；在長度標(biāo)準(zhǔn)化之后，與第i條比對序列相對應(yīng)，第j個相對位點是Nij。高斯取整函數(shù)用方括號代表，要求是，取一個實數(shù)的整數(shù)部分。如此，就在m條長度各不相同的序列經(jīng)過轉(zhuǎn)化之后，得到L長度的目標(biāo)序列。實際上，序列長度標(biāo)準(zhǔn)化的過程就是，按照一定的比例將長度各異的多條序列縮放成等長的序列，便于更好地對最佳匹配區(qū)域的相對位置分布進行分析。

2 結(jié) 果

分析最佳匹配片段的序列構(gòu)成，有助于深入研究內(nèi)含子與相應(yīng)mRNA序列之間相互作用。通過對mRNA序列和相應(yīng)內(nèi)含子序列之間的比對，得到各物種mRNA序列上全部最佳匹配片段集合。要對所有最佳匹配片斷進行序列特征分析，那么就要用到最佳匹配片段中的3個特征參數(shù)，它們分別為最佳匹配片段長度、最佳配對率及其中的GC含量。

2.1 最佳匹配長度分布

通過mRNA序列與相應(yīng)內(nèi)含子序列之間的比對分析，得到了13個物種mRNA序列上的最佳匹配長度分布(見圖1)。考慮這13類生物的最佳匹配片段長度分布形式是一樣的，與α>1的伽馬分布相接近。長度大多在10～80 bp之內(nèi)分布。無脊椎動物、植物，它們的長度平均為20 bp左右；而有脊椎動物，長度大約平均為30 bp。我們發(fā)現(xiàn)，對于高等真核生物，其最佳匹配片段的長度，平均起來要長于低等真核生物。 這就是說生物進化的同時，最佳匹配片段長度，亦隨之增加。那么也就是說，物種愈復(fù)雜，內(nèi)含子與其相應(yīng)mRNA序列二者之間的相互作用模式，也會愈復(fù)雜。

2.2 最佳匹配片段的配對率

對配對率給出如下定義：把L設(shè)定為最佳匹配片段長度，如果k個堿基可以進行完全配對，也就是存在C-G和U-A，那么，k/L即是該最佳匹配片段的配對率。最佳匹配片段的配對率分布圖(見圖2)。不同的物種的配對率分布基本無差別，大多分布在60%～85%，甚至極少部分最佳匹配片段的配對率能達到100%，這些片段的長度都很短，低于15 bp。配對率分布曲線上有一些小的峰值出現(xiàn)，比較明顯的峰值出現(xiàn)在配對率為0.75, 0.8, 0.85 和0.95處，這些峰值分布對各個物種都是一樣的，說明了最佳匹配區(qū)域序列匹配是有規(guī)律的。最佳匹配片段配對率分布在各個物種中具有普適性，反映了配對率的保守性。這表明不同進化水平的物種其mRNA與相應(yīng)內(nèi)含子的相互作用或序列匹配方式遵循同一個匹配機制。

圖1 最佳匹配片段長度分布Fig.1 Length distributions of optimal matched segments

2.3 最佳匹配片段GC含量分布

對mRNA的序列特征分析時，發(fā)現(xiàn)，GC含量是非常重要的參數(shù)。13類生物基因的mRNA序列上的最佳匹配片段的GC含量分布情況(見圖3)。四種植物最佳匹配片段的GC含量分布在一個比較窄的范圍(0.2～0.6)之內(nèi)，最概然分布是0.35，兩個脊椎動物雞和斑馬魚的分布范圍比植物要寬，在0.1～0.7，最概然GC含量是0.45。線蟲的分布與植物相似，果蠅的分布與脊椎動物的分布相似。哺乳動物最佳匹配區(qū)域的GC含量分布各不相同，大鼠和小鼠的最概然GC含量更大一些(約為0.55)，牛的最概然GC含量很小為0.3。但它們的分布范圍很寬，主要在0.1～0.8。人類和狗的分布與其它物種的分布不同，不僅GC含量分布更寬，在0.1～0.9，而且最佳匹配區(qū)域的GC含量分布是雙峰分布。第一個峰的最概然GC含量值是0.3，第二個峰的最概然GC含量值是0.7。

圖2 最佳匹配片段配對率分布Fig.2 Matching rate distributions of optimal matched segments

圖3 最佳匹配片段的GC含量分布Fig.3 GC content distribution of optimal matched segments

可見最佳匹配片段的最概然GC含量隨著物種進化在逐漸增加，這與基因組序列GC含量隨物種進化而增加的現(xiàn)象是一致的。值得注意的是最佳匹配片段的GC含量分布范圍很廣，在0.1～0.9廣泛分布。由此說明最佳匹配片段的發(fā)生不受GC含量水平的影響，有的GC含量高，有的GC含量低。

對于大多數(shù)低等生物，最佳匹配片段中GC含量均偏低，接近內(nèi)含子的GC含量。這表明基因中mRNA序列上，最佳匹配多發(fā)生在低GC區(qū)域。而且內(nèi)含子與mRNA之間存在一種弱相互作用，從而形成弱雙鏈結(jié)構(gòu)。然而對高等生物而言，最佳匹配片段中GC含量多數(shù)偏高，與其外顯子中GC含量接近。在人類和狗中，最佳匹配片段中GC含量分布特征與線蟲和小鼠的分布一致。盡管人類與小鼠基因在進化上很接近，但最佳匹配片段中GC含量分布卻有差別。這說明內(nèi)含子與mRNA序列之間的最佳匹配序列組成在一定程度上反映了表觀遺傳的進化差異。

2.4 最佳匹配片段在mRNA序列上分布特征

針對13個物種基因序列，把它們的內(nèi)含子和對應(yīng)的成熟mRNA放在一起，進行局域比對。在成熟mRNA上，出現(xiàn)最佳匹配區(qū)域位置。因成熟mRNA序列長度不一致，為研究mRNA序列上各個相對位置的匹配強度，所以把成熟mRNA序列進行長度標(biāo)準(zhǔn)化，長度標(biāo)準(zhǔn)化到100 bp。發(fā)現(xiàn)成熟mRNA上出現(xiàn)相對匹配頻率的分布情況(見圖4)。結(jié)果如下：(1)所有物種的匹配強度分布趨勢一致，在序列UTR區(qū)(兩端)的匹配頻率明顯高于編碼序列(中部)，特別是在3’UTR區(qū)出現(xiàn)極大值分布。(2)無脊椎動物在3’UTR區(qū)的峰值最大，比如果蠅和線蟲的極大值大約是它們極小值的7倍(見圖4b)，植物的匹配頻率的極大值略高于有脊椎動物的峰值。(3)在mRNA序列的5’UTR區(qū)，有3個有脊椎動物(人類，狗和大鼠)的匹配F值高于它們的CDS區(qū)(見圖4a)，而其它有脊椎物種則低于CDS區(qū)。單子葉植物(水稻)的F值低于雙子葉植物的F值，也低于CDS區(qū)的F值(見圖4c)。(4)所有物種在CDS區(qū)匹配強度大小較一致，并且F值都比較低。從整體上來說，植物物種在CDS區(qū)的F值相對最高，無脊椎動物的F值相對最低，而有脊椎動物F值介于它們之間。

圖4 mRNA序列和相應(yīng)的內(nèi)含子序列局域比對的相對頻率分布Fig.4 Distributions of matching frequency between mRNA sequences and cooresponding intron sequences

圖3中，為找到最佳匹配片段的GC含量分布情況，把最佳匹配片段按照GC含量的大小分為2組：GC含量高于0.5的高GC片段和GC含量小于0.3的低GC片段。得到它們在mRNA序列的上分布(見圖5、圖6)。分析發(fā)現(xiàn)，在mRNA序列的UTR區(qū)，低GC片段在其匹配強度更高，而在編碼區(qū)對應(yīng)的F值則更低。無脊椎動物和植物在3’UTR區(qū)有顯著的峰值(見圖5b,5c)，對有脊椎動物來說，3’UTR區(qū)出現(xiàn)了多個極值的匹配區(qū)域(見圖5a)。4個植物物種和果蠅仍在mRNA序列的5’UTR區(qū)出現(xiàn)極大值分布，且約為極小值的3倍(見圖5c)，但是其他物種在該區(qū)域的匹配頻率極低。

對于高GC片段在mRNA序列上的匹配強度分布，我們發(fā)現(xiàn)，無脊椎動物在mRNA序列上的各個位點的F值大小相近(見圖6b)，然而植物和有脊椎動物在5’UTR匹出現(xiàn)極大值，其中有脊椎動物出現(xiàn)顯著的極大值約為極小值的5.5倍(見圖6a,6c)。另外在CDS區(qū)和3’UTR區(qū)的匹配頻率一致，且F值均較低。

mRNA序列上存在許多與內(nèi)含子序列的匹配區(qū)域，低GC片段偏好與UTR區(qū)作用，特別是在mRNA序列的3’UTR區(qū)。高GC片段偏好與mRNA序列的5’UTR區(qū)匹配。高GC片段在mRNA序列的3’UTR區(qū)和編碼區(qū)分布沒有顯著的差別，說明在內(nèi)含子上還存在一些高GC含量區(qū)域，它們與整個mRNA序列都有作用。

圖5 GC含量在0-0.3內(nèi)mRNA序列匹配頻率Fig.5 Distribution of F value with GC content between 0.0 and 0.3

圖6 GC含量大于0.50的mRNA序列匹配頻率Fig.6 Distribution of F value with GC content greater than 0.5

3 結(jié)果與討論

通過序列比對獲得了13個物種基因中內(nèi)含子與mRNA二者的最佳匹配片段。通過這些最佳匹配片段序列特征及匹配頻率的分布規(guī)律，我們發(fā)現(xiàn)，siRNA和miRNA的結(jié)合特征與得到的最佳匹配片段的平均長度和配對率分布一致；mRNA上的UTR區(qū)偏好與內(nèi)含子相互作用，而CDS區(qū)域與內(nèi)含子的匹配程度較低。結(jié)論表明內(nèi)含子與成熟mRNA序列存在相互作用。

目前，人們對內(nèi)含子與mRNA序列二者的最佳匹配片段還缺乏一定的了解。小干擾RNA(siRNA)、微小RNA(miRNA)以及piRNA(Piwi-interactingRNA)這三個非編碼序列承擔(dān)了RNA干涉和RNA抑制的功能。一方面，Dicer可以把長度在21～25 bp之間的siRNA加工形成雙鏈RNA，然后雙鏈RNA與靶mRNA進行嚴格地互補再指導(dǎo)mRNA沉默[24]。另一方面，要形成單鏈RNA，Dicer還會在18～25 bp之間，選出miRNA進行加工，單鏈RNA與目標(biāo)mRNA這二者之間，在不同程度地進行完成互補，之后對靶mRNA轉(zhuǎn)錄和表達進行干涉和抑制[25]。通過統(tǒng)計數(shù)據(jù)，可以看到miRNA與靶mRNA的匹配率在65%～95%的范圍內(nèi)，那么我們認為miRNA在調(diào)控發(fā)育過程中發(fā)揮著關(guān)鍵作用。piRNA是一類長度為26～31 bp單鏈的小RNA，大部分集中在29～30 bp之間，而且只有通過小RNA與屬于PIWI蛋白的家族成員進行結(jié)合，piRNA才能發(fā)揮調(diào)控作用。piRNA的發(fā)現(xiàn)對非編碼小分子RNA的研究開拓出新的領(lǐng)域，因此Science把該研究稱為2006年十大科技進展之一[26]。

把內(nèi)含子與mRNA之間相互作用的最佳匹配片段的長度以及配對率的數(shù)據(jù)區(qū)間范圍與miRNA和siRNA進行比較，結(jié)果顯示雙方的長度和配對率范圍的特性，竟然相當(dāng)?shù)囊恢隆?nèi)含子在經(jīng)過剪切之后，和對應(yīng)的mRNA之間，相互作用的程度高，無論是mRNA與內(nèi)含子作用的最佳匹配片段還是siRNA、miRNA和piRNA，從生物選擇的功能片段長度到配對率，它們所遵循的生物學(xué)機制是一樣的。由此推論，在生物基因組序列上，還存在大量的各種形式的類似miRNA的功能片段，它們在基因表達調(diào)控和實現(xiàn)表觀遺傳多樣性方面起著決定性的作用。在不同程度上，通過和mRNA序列之間的互補，一些存在于內(nèi)含子序列上的區(qū)域，才能夠?qū)蛘{(diào)控和表達產(chǎn)生深刻影響。對配對率的分布進行研究之后，我們發(fā)現(xiàn)，絕大部分的片段所進行的配對并不嚴格，然而有極少的片段，嚴格進行配對，但這些片段的長度很短，遠小于20 bp。按照RNA干涉理論，完全匹配對基因表達是致死的，因此在mRNA與內(nèi)含子之間存在的相互作用片段，顯然避開了這類致死的匹配片段。從這個角度來看，mRNA-內(nèi)含子的相互作用理論是合理的。分析認為，對于真核基因這一基因種群，他們自己就具備基因調(diào)控所需的原件。因為內(nèi)含子是非編碼序列，所以它可以完成剪接與可變剪接這兩個任務(wù)，而且能夠?qū)嵺`基因表達調(diào)控這一重要功能。內(nèi)含子序列作為具有功能的一類RNA-RNA相互作用的非編碼序列集合，研究人員需要對此重要特征予以關(guān)注。

發(fā)現(xiàn)內(nèi)含子與mRNA序列的匹配片段在mRNA序列上有分布，其匹配頻率，發(fā)生在兩端非編碼序列區(qū)域的值較大，在中間編碼序列區(qū)域，卻有較低值。同時，內(nèi)含子又出現(xiàn)偏好和mRNA序列上3’UTR區(qū)相互作用。在mRNA序列的3’UTR區(qū)，GC值比較小的片段，擁有相對高的匹配強度，匹配強度卻在編碼區(qū)較低。那么就說明內(nèi)含子與相應(yīng)mRNA之間存在的相互作用主要是弱鍵，也就是所謂的AU匹配，同時還包括GC值大的匹配。但是，大家有可能提出疑問的是，基因序列進化時，UTR類似于內(nèi)含子的GC值，很清楚的顯示出不同于編碼序列的GC值，這一現(xiàn)象的原因是什么呢？想要出現(xiàn)內(nèi)含子和UTR序列發(fā)生較強的相互作用，然后通過該作用對基因進行調(diào)控，從而使以上兩類序列的進化趨于一致。在UTR和編碼區(qū)上GC含量較高的最佳匹配片段分布幾乎不存在差別，這也表示一些高GC區(qū)域存在于內(nèi)含子序列上。可以看出，要想進一步探究內(nèi)含子對基因的表達與調(diào)控，還必須繼續(xù)挖掘mRNA序列上匹配頻率分布的內(nèi)涵。

綜上所述，所有的研究結(jié)論，有力地證明了，內(nèi)含子和mRNA之間發(fā)生的相互作用，是真實存在的。在維持基因組正常運轉(zhuǎn)的過程中，內(nèi)含子可能起到了比較關(guān)鍵的調(diào)控作用[27]。然而尚需進一步對真核生物內(nèi)含子與相應(yīng)的mRNA序列的相互作用進行深入探討, 本文雖然揭示了一些有意思的論點，然而一些結(jié)論尚需通過實驗進一步驗證。