王 娟， 常 瑋

(南陽師范學(xué)院 生命科學(xué)與農(nóng)業(yè)工程學(xué)院，河南 南陽 473061)

單核苷酸多態(tài)性(single nucleotide polymorphism，SNP)是指在基因組上單個(gè)核苷酸的變異，包括轉(zhuǎn)換、顛換、缺失和插入.據(jù)報(bào)道，在人類基因組中，每1000～2000個(gè)堿基就有一個(gè)SNP，在一些高頻的區(qū)域，甚至能達(dá)到每300個(gè)堿基一個(gè)SNP[1-2].由于SNP在基因組中的高密度，使其成為一種理想的分子標(biāo)記，并已被用于構(gòu)建人類遺傳分析圖譜和人類復(fù)雜疾病性狀的研究[3-5].SNP在植物中也廣泛存在，例如在水稻基因組中的發(fā)生率接近1/268[6].根據(jù)SNP在基因中的位置，可將其分為基因編碼區(qū)SNP(coding SNP，cSNP)、基因內(nèi)含子區(qū)SNP和基因調(diào)控區(qū)SNPs(regulatory SNPs，rSNP).其中cSNP根據(jù)是否改變編碼的氨基酸類型又可分為同義cSNP(synonymous SNP，sSNP)和非同義cSNP(non-synonymous SNP，nsSNP).其中nsSNP由于能引起蛋白質(zhì)編碼序列的改變，因而以其作為研究基因功能的一種功能標(biāo)記.這使得nsSNP成為遺傳研究中應(yīng)用前景廣闊的理想分子標(biāo)記.

目前，開發(fā)SNP標(biāo)記(SNPs)的主要方法有：EST序列挖掘法[7]、基因芯片檢測(cè)法[8]、基因組比較法[9]、擴(kuò)增子重測(cè)序技術(shù)[10]以及基于下一代測(cè)序技術(shù)(Next-Generation Sequencing，NGS)的SNPs挖掘法[11].隨著測(cè)序技術(shù)的提高與成本的下降，NGS技術(shù)已經(jīng)成為SNPs挖掘的主要方法.例如，基于NGS技術(shù)的RNA-seq研究，由于可以直接獲得基因轉(zhuǎn)錄本的序列，也為挖掘nsSNP標(biāo)記(nsSNPs)提供了基礎(chǔ)[12].盡管基于NGS的RNA-seq方法可以獲得存在于細(xì)胞中所有的mRNA分子，但由于RNA剪接的復(fù)雜性會(huì)導(dǎo)致假陽性結(jié)果的出現(xiàn)[13]，因此如何獲得更為可靠的nsSNPs，已成為功能基因分子標(biāo)記開發(fā)、分子輔助育種的關(guān)鍵.隨著測(cè)序數(shù)據(jù)的不斷增長(zhǎng)和計(jì)算機(jī)技術(shù)的發(fā)展，目前已經(jīng)出現(xiàn)了許多降低SNP假陽性率的算法與程序，如基于決策樹的機(jī)器學(xué)習(xí)法[14]、ISMU(Integrated SNP Mining and Utilization)和GBS-SNP-CROP等程序包[15-16]以及具有交互式圖形用戶界面的軟件PLANET-SNP[17].這些算法及程序充分利用了SNP位點(diǎn)自身的信息進(jìn)行判斷，并能夠顯著降低SNP假陽性率，但在計(jì)算過程中沒有分析SNP位點(diǎn)側(cè)翼序列信息.趙輝等研究發(fā)現(xiàn)，基因組的堿基組分在不同程度上與堿基替換突變的發(fā)生相關(guān)[18].該研究同時(shí)表明水稻A&T堿基影響突變的范圍局限在突變位點(diǎn)上下游2個(gè)堿基內(nèi)，而擬南芥A&T堿基的影響范圍不超過上下游4個(gè)堿基.因此對(duì)SNP位點(diǎn)側(cè)翼序列成分進(jìn)行分析統(tǒng)計(jì)，可以進(jìn)一步降低SNP假陽性率[18].

大豆是重要的糧食作物，相比于玉米等作物，近年來產(chǎn)量提升緩慢，分子標(biāo)記輔助技術(shù)將是加快育種進(jìn)程的重要手段，而開發(fā)大豆nsSNPs是實(shí)現(xiàn)大豆分子育種的基礎(chǔ).本研究擬采用生物信息的方法對(duì)大豆ESTs數(shù)據(jù)中所包含的候選SNPs進(jìn)行挖掘，結(jié)合數(shù)據(jù)庫中公布的大豆SNPs，依據(jù)其編碼蛋白的情況挖掘其中的nsSNPs，并對(duì)這些nsSNPs分布規(guī)律及其附近的序列特征進(jìn)行統(tǒng)計(jì)分析，以期為降低nsSNPs假陽性率提供理論依據(jù).

1 材料與方法

1.1 材料

1.1.1 序列數(shù)據(jù)

從GenBank/dbEST(http：//www.ncbi.nlm.nih.gov/dbEST/)下載共1 453 823條大豆EST序列,大豆基因組序列及大豆cDNA序列來自phytozome數(shù)據(jù)庫(http://www.phytozome.net).

1.1.2 軟件工具

用于序列拼接的程序：CAP3[19]；用于本地化的序列比對(duì)程序：BLAST[20]；用于在Windows平臺(tái)上執(zhí)行Perl程序的工具軟件：ActivePerl；在Windows平臺(tái)上運(yùn)行的unix模擬環(huán)境：Cygwin；用于SNPs是評(píng)估的程序：MAVIANT[21]；用于非同義SNPs挖掘的工具：UltraCompare；用于序列分析的Perl程序工具包：BioPerl[22].

1.2 方法

1.2.1 SNPs挖掘及篩選

將全部1 453 823條EST序列用BLAST進(jìn)行聚類，E值小于10-10；將每一類的序列用CAP3進(jìn)行拼接，每一個(gè)重疊群必須有大于95%的相似度；將得到的拼接結(jié)果與相應(yīng)的質(zhì)量文件輸入MAVIANT，去除低質(zhì)量的結(jié)果，利用MAVIANT的SNP評(píng)估功能進(jìn)行候選SNPs的篩選.根據(jù)報(bào)道的主要作物SNP的出現(xiàn)頻率設(shè)置窗口大小為150bp，進(jìn)一步對(duì)所得的SNPs進(jìn)行篩選.

1.2.2 nsSNPs挖掘

將得到的包含候選SNPs的一致性序列與大豆cDNA序列進(jìn)行blastn分析，從而將SNPs定位在大豆cDNA序列上，記錄每個(gè)位點(diǎn)的位置；將包含cSNPs的編碼序列，按照野生型與突變型用BioPerl中的make_mrna_protein.pl腳本進(jìn)行翻譯，并將生成的兩蛋白序列文件用UltraCompare比較，找出其中的nsSNPs.

1.2.3 nsSNPs相鄰堿基組分分析

根據(jù)趙輝等對(duì)水稻和擬南芥基因組相鄰堿基組分與產(chǎn)生SNP的轉(zhuǎn)換或顛換的研究方法[18]，對(duì)大豆SNPs相鄰堿基組分進(jìn)行統(tǒng)計(jì)分析：首先統(tǒng)計(jì)基因組中對(duì)應(yīng)兩種堿基替換類型(轉(zhuǎn)換和顛換)SNP位點(diǎn)上下游各12個(gè)位置共24個(gè)位點(diǎn)組成的區(qū)域(稱為背景區(qū)域)堿基A&T的使用頻率，記為f0=(f01,f02)，其中f01對(duì)應(yīng)于堿基A的使用頻率，f02對(duì)應(yīng)于堿基T的使用頻率.然后統(tǒng)計(jì)SNP位點(diǎn)上下游各12個(gè)位置上A&T堿基的使用頻率，記為fi=(fi1,fi2)，i=-12,-11,…,-2,-1,1,2,…,11,12,其中fi1表示第i位置上堿基A的使用頻率，fi2表示第i位置上堿基T的使用頻率，SNP位點(diǎn)的位置記為0位，將5’端為負(fù)號(hào)位，3’端記為正號(hào)位.最后計(jì)算fd=fi-f0，并根據(jù)fd的值的大小，估計(jì)該位置A&T堿基對(duì)SNP發(fā)生概率的影響程度和影響范圍.

根據(jù)SNP位點(diǎn)上下游緊鄰的兩個(gè)位點(diǎn)上的A&T個(gè)數(shù)將其分為三類，即沒有A或T，只有一個(gè)A或T(另一個(gè)是G或C)和兩個(gè)都是A或T，分別用AT0，AT1和AT2記之.首先，統(tǒng)計(jì)這三類內(nèi)6種類型SNP的個(gè)數(shù)，然后再統(tǒng)計(jì)每一類(AT0，AT1和AT2)轉(zhuǎn)換(A/G和C/T型)和顛換(A/C，G/T，A/T和C/G型)的個(gè)數(shù)，從而得出每類的轉(zhuǎn)換和顛換的比值(Ts/Tv).

2 結(jié)果與分析

2.1 大豆EST序列中SNPs的挖掘

全部1 453 823條EST序列共拼接得到50 867個(gè)contigs和183 847個(gè)single ESTs，利用CAP3產(chǎn)生的50 867個(gè)contigs以及與之相對(duì)應(yīng)的質(zhì)量文件，經(jīng)篩選共得到16 093個(gè)候選SNPs.如表1所示：其中：[A/T]、[G/T]、[A/C]、[G/C]、[A/G]和[C/T]突變類型分別檢測(cè)到1 903、1 516、1 396、854、4 875和5 576次.其中發(fā)生轉(zhuǎn)換的次數(shù)為：[A/G]+[C/T]= 4 875+5 576=10 451次；發(fā)生顛換的次數(shù)為：[A/T]+[G/T]+[A/C]+[G/C]=1 903+1 516+1 369+854=5 642次.從結(jié)果中可以看出轉(zhuǎn)換發(fā)生的次數(shù)明顯高于顛換，在轉(zhuǎn)換中[C/T]突變類型高于[A/G]；而在顛換中，[A/T]突變類型最多，[G/C] 突變類型最少.

進(jìn)一步分析每種突變類型中兩種突變事件,如 [A/T]突變類型中的兩種突變事件：[A→T]與[T→A]的發(fā)生頻率，結(jié)果見表1.從表1中我們可以看出，在[A/G]轉(zhuǎn)換中，[G→A]發(fā)生的次數(shù)明顯高于[A→G]發(fā)生的次數(shù)；在[C/T]轉(zhuǎn)換中，[C→T]發(fā)生的次數(shù)明顯高于[T→C]發(fā)生的次數(shù).而其他四種顛換類型的兩種突變事件的發(fā)生次數(shù)之間沒有明顯的差別.

表1 大豆SNP類型統(tǒng)計(jì)表

2.2 編碼序列中nsSNPs的挖掘

將包含候選SNPs位點(diǎn)的16 093條一致性序列與大豆cDNA序列進(jìn)行blastn分析，共有14 297條能夠匹配相應(yīng)的cDNA序列.將這14 297條與植物蛋白數(shù)據(jù)庫比對(duì)，共獲得13 888個(gè)匹配結(jié)果，保留Score≥150、Identities≥65%的結(jié)果共7 085個(gè).搜尋這些cDNA序列上的ORF的起始位置，根據(jù)SNPs位點(diǎn)的位置，共獲得5’UTR-SNPs：63個(gè)、3’UTR-SNPs：1 508個(gè)以及cSNPs：5 514個(gè).根據(jù)UltraCompare的比較結(jié)果，其中sSNPs：3 748個(gè)；nsSNPs：1 766個(gè)(圖1).位于編碼區(qū)的SNPs占總數(shù)的77%，在這些SNPs中，sSNPs所占的比例約為nsSNPs的2倍，而在剩下的22%位于非編碼區(qū)內(nèi)的SNPs，約95%的位點(diǎn)位于3’UTR.分別統(tǒng)計(jì)分布于不同區(qū)域的SNPs中轉(zhuǎn)換Ts與顛換Tv的比值，5’UTR-SNPs、3’UTR-SNPs、sSNPs與nsSNPs的Ts/Tv分別為2、1.5、2.25、1.7，四者的關(guān)系為：sSNPs>5’UTR-SNPs>nsSNPs>3’UTR-SNPs.

圖1 大豆編碼序列SNP類型統(tǒng)計(jì)表

2.3 相鄰堿基組分與產(chǎn)生nsSNPs的轉(zhuǎn)換或顛換的分析

分析結(jié)果表明顯示，無論總SNPs還是不同分布的SNPs，均表現(xiàn)為Ts/Tv大于1，進(jìn)一步的分析結(jié)果還表明：大豆nsSNPs鄰近堿基中，A&T堿基在遠(yuǎn)離突變位點(diǎn)的其他位置的堿基組成無差異，而在與突變位點(diǎn)緊密相鄰的兩個(gè)堿基則表現(xiàn)出了一定的規(guī)律，即無論是轉(zhuǎn)換(圖2A)還是顛換(圖2B)，突變點(diǎn)5’端的A&T 堿基使用頻率偏差接近0；突變點(diǎn)3’端的A堿基使用頻率偏差約為-0.05；T堿基使用頻率偏差約為0.02，表明A&T 堿基的影響范圍局限在下游第一個(gè)堿基.大豆nsSNPs的AT0，AT1和AT2三類中 Ts/Tv 的比值依次為1.682、1.685、1.77.結(jié)果表明：只有在兩側(cè)堿基都為A或者T時(shí)，Ts/Tv 的比值才會(huì)升高.由于突變位點(diǎn)上游第一個(gè)堿基無明顯的頻率偏差，故推測(cè)下游第一個(gè)堿基的頻率偏差會(huì)導(dǎo)致Ts/Tv 的比值的升高，同時(shí)Ts中T的出現(xiàn)頻率(-0.052)小于Tv中T的出現(xiàn)頻率(-0.048)，說明在突變位點(diǎn)下游第一個(gè)堿基T頻率偏差與Ts/Tv的比值成反比.

3 討論

位點(diǎn)的可靠性是SNP挖掘的難點(diǎn)，一般來說可以從以下四方面進(jìn)行評(píng)價(jià)：第一，具有同一突變類型的片段越多，結(jié)果就越可靠，這也是本研究采用EST數(shù)據(jù)進(jìn)行SNP挖掘的主要原因；第二，由于測(cè)序的原因，相較于測(cè)序片段的兩端，讀序中段假陽性率更低；第三，由于連鎖不平衡現(xiàn)象的存在，相鄰兩個(gè)SNP位點(diǎn)應(yīng)緊密連鎖；第四，由于編碼區(qū)受選擇壓力更大，SNP的比例低于全基因組的水平[23].

圖2 nsSNPs相鄰堿基組分分布圖

本研究獲得的SNPs位點(diǎn)中，轉(zhuǎn)換變異10 451個(gè)，約為顛換變異的2倍.這是由于胞嘧啶易經(jīng)甲基化和脫氨作用轉(zhuǎn)換為胸腺嘧啶，致使置換發(fā)生率高于顛換[23].本研究發(fā)現(xiàn)[G→A]發(fā)生的次數(shù)明顯高于[A→G]發(fā)生的次數(shù)；[C→T]發(fā)生的次數(shù)明顯高于[T→C]發(fā)生的次數(shù)，與上述結(jié)論一致.此外，由于5’UTR中包含一些調(diào)控元件，因此相比較于3’UTR受到的選擇壓力更大一些，所以位于3’UTR的SNPs數(shù)目遠(yuǎn)遠(yuǎn)大于5’UTR；而由于編碼區(qū)受選擇壓影響，編碼區(qū)內(nèi)不改變氨基酸序列的候選SNP比改變氨基酸序列的候選SNP具有更高的存在概率，因此位于編碼區(qū)的SNPs中nsSNPs數(shù)目少于sSNPs的數(shù)目[24]，本研究也得到了相同的結(jié)果.

DNA復(fù)制酶的嚴(yán)謹(jǐn)性、各種化學(xué)物理誘變因素等都會(huì)影響基因組的堿基替換變異.我們的研究發(fā)現(xiàn)，大豆nsSNPs位點(diǎn)影響范圍局限在下游第一個(gè)堿基，且T的含量越低，越有利于轉(zhuǎn)換的發(fā)生.這一結(jié)果與趙輝等的研究一致[18].Radman和Wagner的研究對(duì)此現(xiàn)象的解釋為：當(dāng)某個(gè)位點(diǎn)周圍只有A或T堿基時(shí)，形成CG二聯(lián)碼的概率為0，因而不易發(fā)生轉(zhuǎn)換突變[25].

4 結(jié)論

在大豆SNP位點(diǎn)側(cè)翼序列中，A&T堿基的影響范圍局限在下游第1個(gè)堿基，且無論是轉(zhuǎn)換還是顛換，突變點(diǎn)5’端的A&T堿基使用頻率偏差均約為0；3’端的A堿基使用頻率偏差為-0.05；T堿基使用頻率偏差約為0.02.