gene.m

这是测定基因序列的matlab程序

clear;
clc;
fid=fopen('gene.txt');
[yi,count]=fread(fid,Inf,'char'),
y=char(yi');
fclose(fid);
%将基因序列变为四个数值序列
n=1;
for m=1:count
 if isletter(y(m))==1       
     upper(y(m));   
     if strcmp(y(m),'A')==1
        AA(n)=1;                    
     else AA(n)=0;
     end 
     if strcmp(y(m),'G')
        GG(n)=1;
     else GG(n)=0;
     end   
     if strcmp(y(m),'C')  
         CC(n)=1;
     else CC(n)=0 ;
     end    
    if strcmp(y(m),'T') 
        TT(n)=1;
    else TT(n)=0;
    end
    n=n+1;
end
end
N=n-1;
%求出四个数值序列的平均功率谱
k=1:N/2;
f=k/N;
s1=1/N*(abs(spectrum(AA,N))+abs(spectrum(GG,N))+abs(spectrum(CC,N))+abs(spectrum(TT,N)));
s=s1(:,1);

plot(f,s(k),'r');%基因功率谱
xlabel('f');
ylabel('S(f)');
%用选取窗M的长度的351，然后每隔3个序列移动窗，算出该窗内1/3处的频谱值，
%用p表示，当p大于4时，定义为编码区，当p小于4时，定义为非编码区
m=1;
M=351;
for k=1:3:(N-M+1) 
    n=351;
    X1=0;X2=0;X3=0;X4=0;
    while  n>0
    w=exp(-j*2*pi*(1/3)*n);
    X1=w*AA(n+k-1)+X1;
    X2=w*TT(n+k-1)+X2;
    X3=w*CC(n+k-1)+X3;
    X4=w*GG(n+k-1)+X4;
    n=n-1;
    end
    p(m)=1/351*(abs(X1)^2+abs(X2)^2+abs(X3)^2+abs(X4)^2);
  m=m+1;
end
for n=1:m-1
    r(n)=4;
    if p(n)>4
        s(n)=13;
    else s(n)=12;
    end
end
figure;
x=5000:7000;
plot(x,p(x));
grid on
hold on
plot(x,r(x));
hold on
plot(x,s(x));
xlabel('j');
ylabel('Pm(j)');