王 璐

(西安建筑科技大學(xué)理學(xué)院數(shù)學(xué)系 陜西西安 710055)

隨著人類基因組計(jì)劃的完成，生物基因數(shù)據(jù)呈指數(shù)形式增長(zhǎng)，找出蛋白質(zhì)編碼基因，是進(jìn)行基因組分析的基礎(chǔ)，在生物信息處理中占有非常重要的地位[1].通常的基因識(shí)別方法大致可以分為如下3類：序列相似性方法、從頭預(yù)測(cè)方法、序列相似性和從頭預(yù)測(cè)方法相結(jié)合的第三類方法[2].由于物種的多樣性，生物基因數(shù)據(jù)的指數(shù)型增長(zhǎng)和人類對(duì)其有限的認(rèn)識(shí)等原因，第一類方法的缺陷不僅速度較慢，而且準(zhǔn)確率不高.相較第一類來(lái)說(shuō)，第二類方法具有更堅(jiān)實(shí)的數(shù)學(xué)基礎(chǔ)，模型的物理意義也更加明顯直觀，而且，在實(shí)驗(yàn)當(dāng)中對(duì)若干基因預(yù)測(cè)軟件的測(cè)試表明，具有最高正確率的幾種基因預(yù)測(cè)軟件都屬于這一種方法[1].

在生物學(xué)、醫(yī)學(xué)、藥學(xué)等諸多方面，DNA的研究都具有重要的理論意義和實(shí)際價(jià)值[3].在面對(duì)大量、復(fù)雜的基因序列數(shù)據(jù)時(shí)，如何更好更快捷地獲取準(zhǔn)確的基因信息，如何能夠在眾多的基因序列中確定功率譜和信噪比，如何能夠?qū)γ款惢蚩焖俚氐玫狡溟撝荡_定方法，如何快速實(shí)現(xiàn)基因識(shí)別算法，是擺在我們面前的一個(gè)具有研究意義的實(shí)際課題.

本文將以現(xiàn)有的基因序列數(shù)據(jù)為研究對(duì)象研究基因識(shí)別的系統(tǒng)建模方法和相關(guān)理論算法，研究基因識(shí)別模型的建立或改進(jìn)方法，通過(guò)統(tǒng)計(jì)分析的相關(guān)方法研究基因特征的識(shí)別，并通過(guò)計(jì)算機(jī)手段進(jìn)行實(shí)現(xiàn).在基因識(shí)別問(wèn)題研究中，首先研究基因識(shí)別問(wèn)題的特征，從而研究編碼蛋白質(zhì)的基因識(shí)別方法，進(jìn)而針對(duì)生物必需基因識(shí)別的問(wèn)題進(jìn)行研究[2].本文在使用已有模型算法并輔助實(shí)現(xiàn)相應(yīng)模型的過(guò)程中，對(duì)已有的方法進(jìn)行深入分析和比較，找出各方法的不足和長(zhǎng)處，發(fā)現(xiàn)新問(wèn)題，提出更有效的解決方法，使其與已有算法互相融合，取長(zhǎng)補(bǔ)短，促進(jìn)其在實(shí)踐中的應(yīng)用[2].在對(duì)各經(jīng)典算法及相關(guān)模型掌握的基礎(chǔ)上，運(yùn)用所學(xué)的統(tǒng)計(jì)分析方法建立數(shù)學(xué)模型，并運(yùn)用數(shù)據(jù)處理軟件對(duì)數(shù)據(jù)進(jìn)行預(yù)處理、編程并求解模型，并對(duì)研究的數(shù)據(jù)進(jìn)行特征提取和識(shí)別.針對(duì)外顯子和內(nèi)含子序列片段的功率譜表現(xiàn)不同的特性，研究相應(yīng)的算法，以實(shí)現(xiàn)外顯子片段的提取和識(shí)別[4].

1 基于Voss 映射求解功率譜和信噪比方法

1.1 Voss 映射

z這是指將同一類信號(hào)從長(zhǎng)串的信號(hào)序列中單獨(dú)分離出來(lái)，并保持其在長(zhǎng)序列中的相應(yīng)位置不變.

在DNA序列研究中，首先需要把A、T、G、C4種核苷酸的符號(hào)序列，根據(jù)一定的規(guī)則映射成相應(yīng)的數(shù)值序列，以便于對(duì)其作數(shù)字處理[3].

令I(lǐng)={A,T,G,C}，長(zhǎng)度為N的任意DNA序列，可表達(dá)為

S={S[n]|S[n]I,n=0,1,2,…N-1}

即A、T、G、C的符號(hào)序列S：S[0],S[1],…,S[N-1].現(xiàn)對(duì)于任意確定的bI，

Voss映射的關(guān)系式為：

稱之為Voss映射，于是生成相應(yīng)的0-1序列(即二進(jìn)制序列){ub[n]}：ub[0]，ub[1],…,ub[N-1](bI).

1.2 傅里葉變換方法

假設(shè)給定的一段DNA序列片段為S = ATCGTACTG，則所生成的4個(gè)0-1序列分別為[4]：

{uA[n]}：{1,0,0,0,0,1,0,0,0}；{uG[n]}：{0,0,0,1,0,0,0,0,1}；

{uC[n]}：{0,0,1,0,0,0,1,0,0}；{uT[n]}：{0,1,0,0,1,0,0,1,0}。

這樣產(chǎn)生的四個(gè)數(shù)字序列又稱為DNA序列的指示序列.為研究DNA編碼序列(外顯子)的特性，對(duì)指示序列分別做離散Fourier 變換(DFT)[3].

以此可得到4個(gè)長(zhǎng)度均為N的復(fù)數(shù)序列{Ub[K]}，bI.

1.3 功率譜與信噪比

計(jì)算每個(gè)復(fù)序列{Ub[K]}的平方功率譜，并相加則得到整個(gè)DNA序列S的功率譜序列{P[K]}[4]：

P[K]=|UA[K]|2+|UT[K]|2+|UG[K]|2+|UC[K]|2,k=0,1,…N-1

筆者采用了附件中的genes6中的數(shù)據(jù)，基因長(zhǎng)度n=5800，運(yùn)行第1行第1列的基因數(shù)據(jù)，計(jì)算機(jī)主程序如下所示：

clear all;

load genes6.mat;

gene=genes(1,1).Seq;

n=length(gene);

%生成A的指示序列uA[n]

for i=1:n

if strcmp(′A′,gene(i))

ua(i)=1;

else ua(i)=0;

end

end

%生成G的指示序列uG[n]

for i=1:n

if strcmp(′G′,gene(i))

ug(i)=1;

else ug(i)=0;

end

end

%生成C的指示序列uC[n]

for i=1:n

if strcmp(′C′,gene(i))

uc(i)=1;

else uc(i)=0;

end

end

%生成T的指示序列uT[n]

for i=1:n

if strcmp('T',gene(i))

ut(i)=1;

else ut(i)=0;

end

end

fua=fft(ua,n);

fug=fft(ug,n);

fuc=fft(uc,n);

fut=fft(ut,n);

for k=1:n

p(k)=abs(fua(k))^2+abs(fug(k))^2+abs(fuc(k))^2+abs(fut(k))^2;

end

plot(p(floor(n/100):n))

e=sum(p)/n

p1=[p(round(n/3)-2),p(round(n/3)-1),p(round(n/3)),p(round(n/3)+1),p(round(n/3)+2)〗;

pn3=max(p1);

r=pn3/e

r=3.4474

通過(guò)運(yùn)行計(jì)算機(jī)程序，得到下列結(jié)果：

e=5.8000e+003

r=3.4474

所以筆者認(rèn)為：DNA序列g(shù)enes6的總功率譜的平均值為e=5.8000e+003，信噪比r=3.4474.

2 基于Z-curve 映射求解功率譜與信噪比方法

設(shè)DNA序列S的4個(gè)指示序列{ub[n]}，bI={A,C,G,T}，n=0,1,2,…N-1的累積序列bn(n=0,1，…，N-1)為.則定義3個(gè)序列 ：x[n],y[n],z[n]:

令x[-1]=0,y[-1]=0和z[-1]=0,以及Δx[n]=x[n]-x[n-1],

Δy[n]=y[n]-y[n-1]和Δz[n]=z[n]-z[n-1]，于是得到Z-curve映射：

例如，對(duì)DNA序列S(n)：ACGTTAG，則對(duì)應(yīng)的Z-curve映射序列為：

Δx[n]:1-11-1-111Δy[n]:11-1-1-11-1Δz[n]:1-1-1111-1

類似于Voss 映射那樣，定義Z-curve映射的總功率譜,

PZ[k]=|ΔX[k]|2+|ΔY[k]|2+|ΔZ[k]|2,

其中ΔX[k]，ΔY[k]和ΔZ[k]分別表示數(shù)字序列Δx[n]，Δy[n]和Δz[n]的離散傅里葉變換.

定義Z-curve映射的信噪比為：

筆者采用了附件中的genes6中的數(shù)據(jù)，基因長(zhǎng)度n=5800，運(yùn)行任意給定的的基因數(shù)據(jù)，計(jì)算機(jī)主程序運(yùn)算結(jié)果如下：

clear all;

Load genes6.mat;

gene=genes(1,1).Seq;

n=length(gene);

%生成A的指示序列uA[n]

for i=1:n

if strcmp(‘A’,gene(i))

ua(i)=1;

else ua(i)=0;

end

end

%生成G的指示序列uG[n]

for i=1:n

if strcmp(‘G’,gene(i))

ug(i)=1;

else ug(i)=0;

end

end

%生成C的指示序列uC[n]

for i=1:n

if strcmp(‘C’,gene(i))

uc(i)=1;

else uc(i)=0;

end

end

%生成T的指示序列uT[n]

for i=1:n

if strcmp(‘T’,gene(i))

ut(i)=1;

else ut(i)=0;

end

end

a=[1,-1,1,-1;1,1,-1,-1;1,-1,-1,1];

dxyz=a*[ua;uc;ug;ut];

fdx=fft(dxyz(1,:),n);

fdy=fft(dxyz(2,:),n);

fdz=fft(dxyz(3,:),n);

for k=1:n

p(k)=ads(fdx(k))^2+abs(fdy(k))^2+abs(fdz(k))^2;

end

plot(p(floor(n/100):n))

>>e=sum(p)/n

e=1.7400e+004

>>pl=[p(round(n/3)-2),p(round(n/3)-1),p(round(n/3),p(round(n/3)+1),p(round(n/3)+2)];

pn3=max(p1);

r=4.5965

通過(guò)運(yùn)行計(jì)算機(jī)程序，我們得到下列結(jié)果：

e=1.7400e+004

r=4.5965

所以筆者認(rèn)為:DNA序列g(shù)enes6的總功率譜的平均值為e=1.7400e+004，信噪比r=4.5965.

3 基因識(shí)別算法的實(shí)現(xiàn)

在此筆者采用的是基于固定長(zhǎng)度滑動(dòng)窗口上頻譜曲線的基因識(shí)別方法，對(duì)一個(gè)DNA序列S和它的指示序列{ub[n]}[3]，bI，n=0,1,2,…,N-1.取長(zhǎng)度M(通常取為3的倍數(shù))作為固定窗口長(zhǎng)度.

圖1是運(yùn)行Voss 映射程序的第1行第1列數(shù)據(jù)所得到的頻譜圖，其R值為R=3.4474.

圖1 Genes6 的第一段基因序列的頻譜圖

從圖1可以看到，該頻譜圖呈現(xiàn)較為明顯的3-周期性，大致在2000、4000和6000處出現(xiàn)峰值.

為探究某片段是否為具有三周期性的外顯子基因，筆者截取片段上的一部分，圖2為附件1第一段基因序列中3—120處的頻譜圖.

圖2 Genes6第一段基因序列中取值(3——120)

通過(guò)圖像，首先去掉開始部分未形成規(guī)律性的頻譜.然后選取最大值處，設(shè)此處為n1/3(在此我們選取48),此處為該段中的第一個(gè)峰值，再繼續(xù)尋找峰值n2/3(60)，比較前兩個(gè)峰值處的值是否相近，若相近，則筆者斷定該處形成了三周期性.利用基因的三周期性，判別處該段為外顯子，此段長(zhǎng)度為234.同理，利用此方法，尋找其他外顯子區(qū)域.兩個(gè)外顯子之間即為內(nèi)含子，由此得出外顯子的兩個(gè)端點(diǎn)(36和72).實(shí)現(xiàn)了對(duì)基因的識(shí)別.

圖2是運(yùn)行Genes6第1行第2列數(shù)據(jù)所得到的頻譜圖，得到R值為R=1.3827.筆者可以看到圖2在整段上不具有3-周期性，但是從局部上來(lái)看部分段具有3-周期性.

圖3為150—480段處的頻譜圖，利用上面的分析方法，同樣可以得到很多段外顯子和內(nèi)含子的區(qū)間，由此可以實(shí)現(xiàn)對(duì)基因的判別.

圖3 Genes6第一段基因序列中取值(150—480)

圖3是運(yùn)行Genes6第1行第3列數(shù)據(jù)所得到的頻譜圖，得到R值為R=1.2524.筆者可以看到圖3在整段上不具有3-周期性，噪聲較大，因此筆者認(rèn)為此段不具有外顯性.

圖4 Genes6 的第二段基因序列的頻譜圖

圖4是運(yùn)行Genes6第1行第4列數(shù)據(jù)所得到的頻譜圖，得到R值為R=3.5113.筆者可以看到圖4在整段上不具有3-周期性，噪聲太大，因此筆者認(rèn)為此段不具有外顯性.

圖5 Genes6 的第三段基因序列的頻譜圖

圖5是運(yùn)行Genes6第1行第5列數(shù)據(jù)所得到的頻譜圖，得到R值為R=2.4334.此段上R值數(shù)值較大，在整段上3-周期性不是很明顯，但是局部存在3-周期性的可能。因此筆者初步認(rèn)為此段可能存在外顯性.

圖6 Genes6 的第四段基因序列的頻譜圖

圖6是運(yùn)行Genes6第1行第6列數(shù)據(jù)所得到的頻譜圖，得到R值為R=0.9958.在整段上3-周期性不是很明顯，但是局部存在3-周期性的可能.因此筆者初步認(rèn)為此段可能存在外顯性。

圖7 Genes6 的第五段基因序列的頻譜圖

圖8 Genes6 的第六段基因序列的頻譜圖

4 結(jié)論

本文在了解基因存在及其構(gòu)成的基礎(chǔ)上，運(yùn)用核苷酸的非平衡性和快速Fourier 變換等方法求取功率譜和信噪比，得到理想準(zhǔn)確的基因序列.

為了檢驗(yàn)所得數(shù)據(jù)的有效性，Voss映射和Z-curve 映射進(jìn)行再次分析.通過(guò)分析與比較，筆者發(fā)現(xiàn)基于Voss映射和Z-curve 映射運(yùn)算的信噪比二者呈現(xiàn)一定的常數(shù)倍數(shù)關(guān)系，大約為1.3333.因此，筆者認(rèn)為Z-curve 映射比Voss映射更具有生物學(xué)意義.

參考文獻(xiàn)：

[1]楊莉.DNA序列4D表示及基因識(shí)別算法研究[D].長(zhǎng)沙：湖南大學(xué)，2007.11.

[2]周柚.基因識(shí)別和微陣列數(shù)據(jù)識(shí)別算法研究[D].長(zhǎng)春：吉林大學(xué)，2008.4.

[3]2012研究生數(shù)模式題[EB/OL].百度文庫(kù)，http//wenku.baidu.com/view/bde7341552d380

ebb2946d91.html.

[4]方芷君，明媚，王婷，等.基因預(yù)測(cè)中的信噪e計(jì)算問(wèn)題[C].桂林：廣西計(jì)算機(jī)學(xué)會(huì)2012年學(xué)術(shù)會(huì)議論文集.2012.