張 強(qiáng)，趙小慶，薄素玲，曹艷娟 ，蘇文霞

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

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

基于此思路，以13個生物基因的編碼序列作為研究樣本，目的在于找到mRNA序列上，能夠與內(nèi)含子序列存在相互作用關(guān)系的最佳匹配片段，進(jìn)而分析這些片段的序列特征，討論在mRNA序列上匹配頻率的分布特點(diǎn)，研究該分布在不同物種中的進(jìn)化規(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)含子序列存在堿基互補(bǔ)匹配片段，該片段能夠反映它們之間的相互特點(diǎn)。因此它們之間的相互作用可由最佳匹配片段表征出來。首先把內(nèi)含子進(jìn)行轉(zhuǎn)化，得到它的互補(bǔ)序列，其次通過軟件Smith-Waterman，對其進(jìn)行局部相似性比對，再次比對互補(bǔ)序列與其所對應(yīng)的mRNA序列，進(jìn)而找出最佳相似片段，最后通過轉(zhuǎn)化得到兩者最佳匹配片段。

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

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

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

(1)

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

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

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

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

(2)

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

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

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

(3)

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

2 結(jié) 果

分析最佳匹配片段的序列構(gòu)成，有助于深入研究內(nèi)含子與相應(yīng)mRNA序列之間相互作用。通過對mRNA序列和相應(yīng)內(nèi)含子序列之間的比對，得到各物種mRNA序列上全部最佳匹配片段集合。要對所有最佳匹配片斷進(jìn)行序列特征分析，那么就要用到最佳匹配片段中的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)，對于高等真核生物，其最佳匹配片段的長度，平均起來要長于低等真核生物。 這就是說生物進(jìn)化的同時，最佳匹配片段長度，亦隨之增加。那么也就是說，物種愈復(fù)雜，內(nèi)含子與其相應(yīng)mRNA序列二者之間的相互作用模式，也會愈復(fù)雜。

2.2 最佳匹配片段的配對率

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

圖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含量隨著物種進(jìn)化在逐漸增加，這與基因組序列GC含量隨物種進(jìn)化而增加的現(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含量分布特征與線蟲和小鼠的分布一致。盡管人類與小鼠基因在進(jìn)化上很接近，但最佳匹配片段中GC含量分布卻有差別。這說明內(nèi)含子與mRNA序列之間的最佳匹配序列組成在一定程度上反映了表觀遺傳的進(jìn)化差異。

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

針對13個物種基因序列，把它們的內(nèi)含子和對應(yīng)的成熟mRNA放在一起，進(jìn)行局域比對。在成熟mRNA上，出現(xiàn)最佳匹配區(qū)域位置。因成熟mRNA序列長度不一致，為研究mRNA序列上各個相對位置的匹配強(qiáng)度，所以把成熟mRNA序列進(jìn)行長度標(biāo)準(zhǔn)化，長度標(biāo)準(zhǔn)化到100 bp。發(fā)現(xiàn)成熟mRNA上出現(xiàn)相對匹配頻率的分布情況(見圖4)。結(jié)果如下：(1)所有物種的匹配強(qiáng)度分布趨勢一致，在序列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ū)匹配強(qiáng)度大小較一致，并且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片段在其匹配強(qiáng)度更高，而在編碼區(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序列上的匹配強(qiáng)度分布，我們發(fā)現(xiàn)，無脊椎動物在mRNA序列上的各個位點(diǎn)的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進(jìn)行嚴(yán)格地互補(bǔ)再指導(dǎo)mRNA沉默[24]。另一方面，要形成單鏈RNA，Dicer還會在18～25 bp之間，選出miRNA進(jìn)行加工，單鏈RNA與目標(biāo)mRNA這二者之間，在不同程度地進(jìn)行完成互補(bǔ)，之后對靶mRNA轉(zhuǎn)錄和表達(dá)進(jìn)行干涉和抑制[25]。通過統(tǒng)計數(shù)據(jù)，可以看到miRNA與靶mRNA的匹配率在65%～95%的范圍內(nèi)，那么我們認(rèn)為miRNA在調(diào)控發(fā)育過程中發(fā)揮著關(guān)鍵作用。piRNA是一類長度為26～31 bp單鏈的小RNA，大部分集中在29～30 bp之間，而且只有通過小RNA與屬于PIWI蛋白的家族成員進(jìn)行結(jié)合，piRNA才能發(fā)揮調(diào)控作用。piRNA的發(fā)現(xiàn)對非編碼小分子RNA的研究開拓出新的領(lǐng)域，因此Science把該研究稱為2006年十大科技進(jìn)展之一[26]。

把內(nèi)含子與mRNA之間相互作用的最佳匹配片段的長度以及配對率的數(shù)據(jù)區(qū)間范圍與miRNA和siRNA進(jìn)行比較，結(jié)果顯示雙方的長度和配對率范圍的特性，竟然相當(dāng)?shù)囊恢?。?nèi)含子在經(jīng)過剪切之后，和對應(yīng)的mRNA之間，相互作用的程度高，無論是mRNA與內(nèi)含子作用的最佳匹配片段還是siRNA、miRNA和piRNA，從生物選擇的功能片段長度到配對率，它們所遵循的生物學(xué)機(jī)制是一樣的。由此推論，在生物基因組序列上，還存在大量的各種形式的類似miRNA的功能片段，它們在基因表達(dá)調(diào)控和實現(xiàn)表觀遺傳多樣性方面起著決定性的作用。在不同程度上，通過和mRNA序列之間的互補(bǔ)，一些存在于內(nèi)含子序列上的區(qū)域，才能夠?qū)蛘{(diào)控和表達(dá)產(chǎn)生深刻影響。對配對率的分布進(jìn)行研究之后，我們發(fā)現(xiàn)，絕大部分的片段所進(jìn)行的配對并不嚴(yán)格，然而有極少的片段，嚴(yán)格進(jìn)行配對，但這些片段的長度很短，遠(yuǎn)小于20 bp。按照RNA干涉理論，完全匹配對基因表達(dá)是致死的，因此在mRNA與內(nèi)含子之間存在的相互作用片段，顯然避開了這類致死的匹配片段。從這個角度來看，mRNA-內(nèi)含子的相互作用理論是合理的。分析認(rèn)為，對于真核基因這一基因種群，他們自己就具備基因調(diào)控所需的原件。因為內(nèi)含子是非編碼序列，所以它可以完成剪接與可變剪接這兩個任務(wù)，而且能夠?qū)嵺`基因表達(dá)調(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值比較小的片段，擁有相對高的匹配強(qiáng)度，匹配強(qiáng)度卻在編碼區(qū)較低。那么就說明內(nèi)含子與相應(yīng)mRNA之間存在的相互作用主要是弱鍵，也就是所謂的AU匹配，同時還包括GC值大的匹配。但是，大家有可能提出疑問的是，基因序列進(jìn)化時，UTR類似于內(nèi)含子的GC值，很清楚的顯示出不同于編碼序列的GC值，這一現(xiàn)象的原因是什么呢？想要出現(xiàn)內(nèi)含子和UTR序列發(fā)生較強(qiáng)的相互作用，然后通過該作用對基因進(jìn)行調(diào)控，從而使以上兩類序列的進(jìn)化趨于一致。在UTR和編碼區(qū)上GC含量較高的最佳匹配片段分布幾乎不存在差別，這也表示一些高GC區(qū)域存在于內(nèi)含子序列上?？梢钥闯?，要想進(jìn)一步探究內(nèi)含子對基因的表達(dá)與調(diào)控，還必須繼續(xù)挖掘mRNA序列上匹配頻率分布的內(nèi)涵。

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