闫世军 neu
2008-03-30, 13:48
%输入n个蛋白质二级结构序列，我们首先用2D图形表示方法得到n个蛋白质的M-曲线，从
%曲线中我们提取一些结构特征，同时考虑原有序列的一些结构特征，得到这n个蛋白质的
%结构特征向量CV，在利用欧几里德度量公式得到一个n*n距离矩阵D
clear;
clc;
n=input('the number of the proteins:'); % 输入蛋白质的个数
%输入每个蛋白质二级结构序列SSS(Secondary Structure Sequence)
for i=1:n
fprintf('the sequence %d:\n',i);
SSS{1,i}=input(' ','s');
end
%输入每个蛋白质二级结构序列的名称NAME
for i=1:n
fprintf('the name of the protein %d:\n',i);
NAME{1,i}=input(' ','s');
end
%构造横坐标X的元胞数组X{1,n},纵坐标Y的元胞数Y{1,n}
for i=1:n
m=length(SSS{1,i});
X{1,i}=zeros(1,m+1);Y{1,i}=zeros(1,m+1);
end
%构造蛋白质二级结构序列中H和E的个数向量nH和nE,并给它们赋初值
nH=zeros(1,n);nE=zeros(1,n);
%构造蛋白质二级结构序列对应的类三角形总数n3和类梯形总数n4,以及它们的总数N,并给它们赋初值
n3=zeros(1,n);n4=zeros(1,n);N=zeros(1,n);
%构造蛋白质二级结构的特征向量CV(Characteristic Vector)的元胞数组CV{1,n},并给其赋初值
for i=1:n
CV{1,i}=zeros(1,4);
end
%构造一个n*n的距离矩阵D,并给其赋初值
D=zeros(n,n);
%每个蛋白质二级结构序列中H和E的含量 cH,cE
for i=1:n
m=length(SSS{1,i});
for j=1:m
if SSS{1,i}(j)=='H'
nH(i)=nH(i)+1;
elseif SSS{1,i}(j)=='E'
nE(i)=nE(i)+1;
end
end
cH(i)=nH(i)/m;
cE(i)=nE(i)/m;
end
%用2D图形表示方法得到n个蛋白质二级结构序列的n条M-曲线
for i=1:n
m=length(SSS{1,i});
for j=2:(m+1)
X{1,i}(j)=j-1;
if SSS{1,i}(j-1)=='H'
Y{1,i}(j)=Y{1,i}(j-1)+1;
elseif SSS{1,i}(j-1)=='E'
Y{1,i}(j)=Y{1,i}(j-1)-1;
else Y{1,i}(j)=Y{1,i}(j-1);
end
end
figure(i)
plot(X{1,i},Y{1,i},'m');
title(NAME{1,i});
end
%计算每个蛋白质二级结构序列对应的所有类三角形和类梯形总数
for i=1:n
m=length(SSS{1,i});
j=1;
while j<m
if Y{1,i}(j+1)==Y{1,i}(j)
for k=1:(m-1-j)
if Y{1,i}(j+k+1)==Y{1,i}(j+k)
continue;
else break;
end
end
j=j+k;
elseif Y{1,i}(j+1)>Y{1,i}(j)
for s=1:(m-j-1)
if Y{1,i}(j+s+1)>Y{1,i}(j+s)
continue;
else break;
end
end
j=j+s;
if Y{1,i}(j+1)==Y{1,i}(j)
for t=1:(m-j-1)
if Y{1,i}(j+t+1)==Y{1,i}(j+t)
continue;
elseif Y{1,i}(j+t+1)>Y{1,i}(j+t)
break;
else
n4(i)=n4(i)+1;break;
end
j=j+t;
end
elseif Y{1,i}(j+1)<Y{1,i}(j)
n3(i)=n3(i)+1;
end
else
for p=1:(m-j-1)
if Y{1,i}(j+p+1)<Y{1,i}(j+p)
continue;
else break;
end
end
j=j+p;
if Y{1,i}(j+1)==Y{1,i}(j)
for q=1:(m-j-1)
if Y{1,i}(j+q+1)==Y{1,i}(j+q)
continue;
elseif Y{1,i}(j+q+1)<Y{1,i}(j+q)
break;
else
n4(i)=n4(i)+1;break;
end
j=j+q;
end
elseif Y{1,i}(j+1)>Y{1,i}(j)
n3(i)=n3(i)+1;
end
end
end
N(i)=n3(i)+n4(i);
end
%计算每个蛋白质二级结构序列被H分割的E个数分布的向量pE,它们构成了一个元胞数组pE{1,n}
%同时计算每个蛋白质二级结构序列被E分割的H个数分布的向量pH,它们构成了一个元胞数组pH{1,n}
for i=1:n
m=length(SSS{1,i});
%计算每个蛋白质二级结构序列被H分割的E个数分布的向量pE,它们构成了一个元胞数组pE{1,n}
for j=1:m
if SSS{1,i}(j)=='E'
k=j;break;
end
end
s=1;pE{1,i}(s)=0;
while k<=m
if SSS{1,i}(k)=='C'
k=k+1;
elseif SSS{1,i}(k)=='E'
pE{1,i}(s)=pE{1,i}(s)+1;k=k+1;
else
for n=1:(m-k)
if (SSS{1,i}(k+n)=='H')|(SSS{1,i}(k+n)=='C')
continue;
else break;
end
end
k=k+n; s=s+1;pE{1,i}(s)=0;
end
end
if pE{1,i}(s)==0
pE{1,i}(s)=[ ];
end
%同时计算每个蛋白质二级结构序列被E分割的H个数分布的向量pH,它们构成了一个元胞数组pH{1,n}
for j=1:m
if SSS{1,i}(j)=='H'
k=j;break;
end
end
s=1;pH{1,i}(s)=0;
while k<=m
if SSS{1,i}(k)=='C'
k=k+1;
elseif SSS{1,i}(k)=='H'
pH{1,i}(s)=pH{1,i}(s)+1;k=k+1;
else
for n=1:(m-k)
if (SSS{1,i}(k+n)=='E')|(SSS{1,i}(k+n)=='C')
continue;
else break;
end
end
k=k+n; s=s+1;pH{1,i}(s)=0;
end
end
if pH{1,i}(s)==0
pH{1,i}(s)=[ ];
end
end
%得到每个蛋白质二级结构序列对应的特征向量 CV(Characteristic Vector)
for i=1:n
if cE(i)<1/20----------------------------------------------------------------------------------------179行
CV{1,i}(1)=cE(i);
else CV{1,i}(1)=-1;
end
if cH(i)<1/20
CV{1,i}(2)=cH(i);
else CV{1,i}(2)=0;
end
CV{1,i}(3)=5*N(i)/length(SSS{1,i});
if length(pH{1,i})==0
CV{1,i}(4)=sqrt(sum((pE{1,i}-mean(pE{1,i})*ones(1, length(pE{1,i}))).^2)/ length(pE{1,i}));
elseif length(pE{1,i})==0
CV{1,i}(4)=sqrt(sum((pH{1,i}-mean(pH{1,i})*ones(1,length(pH{1,i}))).^2)/length(pH{1,i}));
else
CV{1,i}(4)=max(sqrt(sum((pH{1,i}-mean(pH{1,i})*ones(1,length(pH{1,i}))).^2)/length(pH{1,i})),sqrt(sum((pE{1,i}-mean(pE{1,i})*ones(1, length(pE{1,i}))).^2)/ length(pE{1,i})));
end
end
%得到 n个蛋白质二级结构序列之间的距离矩阵D(n*n)
for i=1:n
for j=1:n
D(i,j)=sqrt(sum((CV{1,i}-CV{1,j}).^2))%利用了欧几里德距离公式
end
end
出现了如下错误,不知如何修改,请高手指教,多谢!!!
??? Index exceeds matrix dimensions.
Error in ==> graphmatrix at 179
if cE(i)<1/20
%曲线中我们提取一些结构特征，同时考虑原有序列的一些结构特征，得到这n个蛋白质的
%结构特征向量CV，在利用欧几里德度量公式得到一个n*n距离矩阵D
clear;
clc;
n=input('the number of the proteins:'); % 输入蛋白质的个数
%输入每个蛋白质二级结构序列SSS(Secondary Structure Sequence)
for i=1:n
fprintf('the sequence %d:\n',i);
SSS{1,i}=input(' ','s');
end
%输入每个蛋白质二级结构序列的名称NAME
for i=1:n
fprintf('the name of the protein %d:\n',i);
NAME{1,i}=input(' ','s');
end
%构造横坐标X的元胞数组X{1,n},纵坐标Y的元胞数Y{1,n}
for i=1:n
m=length(SSS{1,i});
X{1,i}=zeros(1,m+1);Y{1,i}=zeros(1,m+1);
end
%构造蛋白质二级结构序列中H和E的个数向量nH和nE,并给它们赋初值
nH=zeros(1,n);nE=zeros(1,n);
%构造蛋白质二级结构序列对应的类三角形总数n3和类梯形总数n4,以及它们的总数N,并给它们赋初值
n3=zeros(1,n);n4=zeros(1,n);N=zeros(1,n);
%构造蛋白质二级结构的特征向量CV(Characteristic Vector)的元胞数组CV{1,n},并给其赋初值
for i=1:n
CV{1,i}=zeros(1,4);
end
%构造一个n*n的距离矩阵D,并给其赋初值
D=zeros(n,n);
%每个蛋白质二级结构序列中H和E的含量 cH,cE
for i=1:n
m=length(SSS{1,i});
for j=1:m
if SSS{1,i}(j)=='H'
nH(i)=nH(i)+1;
elseif SSS{1,i}(j)=='E'
nE(i)=nE(i)+1;
end
end
cH(i)=nH(i)/m;
cE(i)=nE(i)/m;
end
%用2D图形表示方法得到n个蛋白质二级结构序列的n条M-曲线
for i=1:n
m=length(SSS{1,i});
for j=2:(m+1)
X{1,i}(j)=j-1;
if SSS{1,i}(j-1)=='H'
Y{1,i}(j)=Y{1,i}(j-1)+1;
elseif SSS{1,i}(j-1)=='E'
Y{1,i}(j)=Y{1,i}(j-1)-1;
else Y{1,i}(j)=Y{1,i}(j-1);
end
end
figure(i)
plot(X{1,i},Y{1,i},'m');
title(NAME{1,i});
end
%计算每个蛋白质二级结构序列对应的所有类三角形和类梯形总数
for i=1:n
m=length(SSS{1,i});
j=1;
while j<m
if Y{1,i}(j+1)==Y{1,i}(j)
for k=1:(m-1-j)
if Y{1,i}(j+k+1)==Y{1,i}(j+k)
continue;
else break;
end
end
j=j+k;
elseif Y{1,i}(j+1)>Y{1,i}(j)
for s=1:(m-j-1)
if Y{1,i}(j+s+1)>Y{1,i}(j+s)
continue;
else break;
end
end
j=j+s;
if Y{1,i}(j+1)==Y{1,i}(j)
for t=1:(m-j-1)
if Y{1,i}(j+t+1)==Y{1,i}(j+t)
continue;
elseif Y{1,i}(j+t+1)>Y{1,i}(j+t)
break;
else
n4(i)=n4(i)+1;break;
end
j=j+t;
end
elseif Y{1,i}(j+1)<Y{1,i}(j)
n3(i)=n3(i)+1;
end
else
for p=1:(m-j-1)
if Y{1,i}(j+p+1)<Y{1,i}(j+p)
continue;
else break;
end
end
j=j+p;
if Y{1,i}(j+1)==Y{1,i}(j)
for q=1:(m-j-1)
if Y{1,i}(j+q+1)==Y{1,i}(j+q)
continue;
elseif Y{1,i}(j+q+1)<Y{1,i}(j+q)
break;
else
n4(i)=n4(i)+1;break;
end
j=j+q;
end
elseif Y{1,i}(j+1)>Y{1,i}(j)
n3(i)=n3(i)+1;
end
end
end
N(i)=n3(i)+n4(i);
end
%计算每个蛋白质二级结构序列被H分割的E个数分布的向量pE,它们构成了一个元胞数组pE{1,n}
%同时计算每个蛋白质二级结构序列被E分割的H个数分布的向量pH,它们构成了一个元胞数组pH{1,n}
for i=1:n
m=length(SSS{1,i});
%计算每个蛋白质二级结构序列被H分割的E个数分布的向量pE,它们构成了一个元胞数组pE{1,n}
for j=1:m
if SSS{1,i}(j)=='E'
k=j;break;
end
end
s=1;pE{1,i}(s)=0;
while k<=m
if SSS{1,i}(k)=='C'
k=k+1;
elseif SSS{1,i}(k)=='E'
pE{1,i}(s)=pE{1,i}(s)+1;k=k+1;
else
for n=1:(m-k)
if (SSS{1,i}(k+n)=='H')|(SSS{1,i}(k+n)=='C')
continue;
else break;
end
end
k=k+n; s=s+1;pE{1,i}(s)=0;
end
end
if pE{1,i}(s)==0
pE{1,i}(s)=[ ];
end
%同时计算每个蛋白质二级结构序列被E分割的H个数分布的向量pH,它们构成了一个元胞数组pH{1,n}
for j=1:m
if SSS{1,i}(j)=='H'
k=j;break;
end
end
s=1;pH{1,i}(s)=0;
while k<=m
if SSS{1,i}(k)=='C'
k=k+1;
elseif SSS{1,i}(k)=='H'
pH{1,i}(s)=pH{1,i}(s)+1;k=k+1;
else
for n=1:(m-k)
if (SSS{1,i}(k+n)=='E')|(SSS{1,i}(k+n)=='C')
continue;
else break;
end
end
k=k+n; s=s+1;pH{1,i}(s)=0;
end
end
if pH{1,i}(s)==0
pH{1,i}(s)=[ ];
end
end
%得到每个蛋白质二级结构序列对应的特征向量 CV(Characteristic Vector)
for i=1:n
if cE(i)<1/20----------------------------------------------------------------------------------------179行
CV{1,i}(1)=cE(i);
else CV{1,i}(1)=-1;
end
if cH(i)<1/20
CV{1,i}(2)=cH(i);
else CV{1,i}(2)=0;
end
CV{1,i}(3)=5*N(i)/length(SSS{1,i});
if length(pH{1,i})==0
CV{1,i}(4)=sqrt(sum((pE{1,i}-mean(pE{1,i})*ones(1, length(pE{1,i}))).^2)/ length(pE{1,i}));
elseif length(pE{1,i})==0
CV{1,i}(4)=sqrt(sum((pH{1,i}-mean(pH{1,i})*ones(1,length(pH{1,i}))).^2)/length(pH{1,i}));
else
CV{1,i}(4)=max(sqrt(sum((pH{1,i}-mean(pH{1,i})*ones(1,length(pH{1,i}))).^2)/length(pH{1,i})),sqrt(sum((pE{1,i}-mean(pE{1,i})*ones(1, length(pE{1,i}))).^2)/ length(pE{1,i})));
end
end
%得到 n个蛋白质二级结构序列之间的距离矩阵D(n*n)
for i=1:n
for j=1:n
D(i,j)=sqrt(sum((CV{1,i}-CV{1,j}).^2))%利用了欧几里德距离公式
end
end
出现了如下错误,不知如何修改,请高手指教,多谢!!!
??? Index exceeds matrix dimensions.
Error in ==> graphmatrix at 179
if cE(i)<1/20