0

我正在尝试 matlab 绘制 ramachandran 图,而不使用内置命令。我也成功了。现在我想在 scatter 数组中单独发现 GLYCINE。任何想法如何做到这一点?(链接到 1UBQ.pdb 文件:http ://www.rcsb.org/pdb/download/downloadFile.do?fileFormat=pdb&compression=NO&structureId=1UBQ )

% Program to plot Ramanchandran plot of Ubiquitin
close all; clear ; clc; % close all figure windows, clear variables, clear screen

pdb1 ='/home/devanandt/Documents/VMD/1UBQ.pdb';
p=pdbread(pdb1); % read pdb file corresponding to ubiquitin protein
atom={p.Model.Atom.AtomName};

n_i=find(strcmp(atom,'N')); % Find indices of atoms
ca_i=find(strcmp(atom,'CA'));
c_i=find(strcmp(atom,'C'));

X = [p.Model.Atom.X];
Y = [p.Model.Atom.Y];
Z = [p.Model.Atom.Z];

X_n = X(n_i(2:end)); % X Y Z coordinates of atoms
Y_n = Y(n_i(2:end));
Z_n = Z(n_i(2:end));

X_ca = X(ca_i(2:end));
Y_ca = Y(ca_i(2:end));
Z_ca = Z(ca_i(2:end));

X_c = X(c_i(2:end));
Y_c = Y(c_i(2:end));
Z_c = Z(c_i(2:end));

X_c_ = X(c_i(1:end-1)); % the n-1 th C (C of cabonyl)
Y_c_ = Y(c_i(1:end-1));
Z_c_ = Z(c_i(1:end-1));

V_c_ = [X_c_' Y_c_' Z_c_'];
V_n = [X_n' Y_n' Z_n'];
V_ca = [X_ca' Y_ca' Z_ca'];
V_c = [X_c' Y_c' Z_c'];

V_ab = V_n - V_c_;
V_bc = V_ca - V_n;
V_cd = V_c - V_ca;

phi=0;
for k=1:numel(X_c)
    n1=cross(V_ab(k,:),V_bc(k,:))/norm(cross(V_ab(k,:),V_bc(k,:)));
    n2=cross(V_bc(k,:),V_cd(k,:))/norm(cross(V_bc(k,:),V_cd(k,:)));
    x=dot(n1,n2);
    m1=cross(n1,(V_bc(k,:)/norm(V_bc(k,:))));
    y=dot(m1,n2);
    phi=cat(2,phi,-atan2d(y,x));

end

phi=phi(1,2:end);

X_n_ = X(n_i(2:end)); %  (n+1) nitrogens
Y_n_ = Y(n_i(2:end));
Z_n_ = Z(n_i(2:end));

X_ca = X(ca_i(1:end-1));
Y_ca = Y(ca_i(1:end-1));
Z_ca = Z(ca_i(1:end-1));

X_n = X(n_i(1:end-1));
Y_n = Y(n_i(1:end-1));
Z_n = Z(n_i(1:end-1));

X_c = X(c_i(1:end-1)); 
Y_c = Y(c_i(1:end-1));
Z_c = Z(c_i(1:end-1));

V_n_ = [X_n_' Y_n_' Z_n_'];
V_n = [X_n' Y_n' Z_n'];
V_ca = [X_ca' Y_ca' Z_ca'];
V_c = [X_c' Y_c' Z_c'];

V_ab = V_ca - V_n;
V_bc = V_c - V_ca;
V_cd = V_n_ - V_c;

psi=0;
for k=1:numel(X_c)
    n1=cross(V_ab(k,:),V_bc(k,:))/norm(cross(V_ab(k,:),V_bc(k,:)));
    n2=cross(V_bc(k,:),V_cd(k,:))/norm(cross(V_bc(k,:),V_cd(k,:)));
    x=dot(n1,n2);
    m1=cross(n1,(V_bc(k,:)/norm(V_bc(k,:))));
    y=dot(m1,n2);
    psi=cat(2,psi,-atan2d(y,x));

end

psi=psi(1,2:end);

scatter(phi,psi)
box on
axis([-180 180 -180 180])
title('Ramachandran Plot for Ubiquitn Protein','FontSize',16)
xlabel('\Phi^o','FontSize',20)
ylabel('\Psi^o','FontSize',20)
grid

输出是:

在此处输入图像描述


编辑:我的情节正确吗?Biopython:如何避免蛋白质中的特定氨基酸序列以绘制 Ramachandran 图?有一个情节略有不同的答案。

修改后的代码如下:

% Program to plot Ramanchandran plot of Ubiquitin with no glycines
close all; clear ; clc; % close all figure windows, clear variables, clear screen

pdb1 ='/home/devanandt/Documents/VMD/1UBQ.pdb';
p=pdbread(pdb1); % read pdb file corresponding to ubiquitin protein
atom={p.Model.Atom.AtomName};

n_i=find(strcmp(atom,'N')); % Find indices of atoms
ca_i=find(strcmp(atom,'CA'));
c_i=find(strcmp(atom,'C'));

X = [p.Model.Atom.X];
Y = [p.Model.Atom.Y];
Z = [p.Model.Atom.Z];

X_n = X(n_i(2:end)); % X Y Z coordinates of atoms
Y_n = Y(n_i(2:end));
Z_n = Z(n_i(2:end));

X_ca = X(ca_i(2:end));
Y_ca = Y(ca_i(2:end));
Z_ca = Z(ca_i(2:end));

X_c = X(c_i(2:end));
Y_c = Y(c_i(2:end));
Z_c = Z(c_i(2:end));

X_c_ = X(c_i(1:end-1)); % the n-1 th C (C of cabonyl)
Y_c_ = Y(c_i(1:end-1));
Z_c_ = Z(c_i(1:end-1));

V_c_ = [X_c_' Y_c_' Z_c_'];
V_n = [X_n' Y_n' Z_n'];
V_ca = [X_ca' Y_ca' Z_ca'];
V_c = [X_c' Y_c' Z_c'];

V_ab = V_n - V_c_;
V_bc = V_ca - V_n;
V_cd = V_c - V_ca;

phi=0;
for k=1:numel(X_c)
    n1=cross(V_ab(k,:),V_bc(k,:))/norm(cross(V_ab(k,:),V_bc(k,:)));
    n2=cross(V_bc(k,:),V_cd(k,:))/norm(cross(V_bc(k,:),V_cd(k,:)));
    x=dot(n1,n2);
    m1=cross(n1,(V_bc(k,:)/norm(V_bc(k,:))));
    y=dot(m1,n2);
    phi=cat(2,phi,-atan2d(y,x));

end

phi=phi(1,2:end);

X_n_ = X(n_i(2:end)); %  (n+1) nitrogens
Y_n_ = Y(n_i(2:end));
Z_n_ = Z(n_i(2:end));

X_ca = X(ca_i(1:end-1));
Y_ca = Y(ca_i(1:end-1));
Z_ca = Z(ca_i(1:end-1));

X_n = X(n_i(1:end-1));
Y_n = Y(n_i(1:end-1));
Z_n = Z(n_i(1:end-1));

X_c = X(c_i(1:end-1)); 
Y_c = Y(c_i(1:end-1));
Z_c = Z(c_i(1:end-1));

V_n_ = [X_n_' Y_n_' Z_n_'];
V_n = [X_n' Y_n' Z_n'];
V_ca = [X_ca' Y_ca' Z_ca'];
V_c = [X_c' Y_c' Z_c'];

V_ab = V_ca - V_n;
V_bc = V_c - V_ca;
V_cd = V_n_ - V_c;

psi=0;
for k=1:numel(X_c)
    n1=cross(V_ab(k,:),V_bc(k,:))/norm(cross(V_ab(k,:),V_bc(k,:)));
    n2=cross(V_bc(k,:),V_cd(k,:))/norm(cross(V_bc(k,:),V_cd(k,:)));
    x=dot(n1,n2);
    m1=cross(n1,(V_bc(k,:)/norm(V_bc(k,:))));
    y=dot(m1,n2);
    psi=cat(2,psi,-atan2d(y,x));

end

psi=psi(1,2:end);

res=strsplit(p.Sequence.ResidueNames,' ');
angle =[phi;psi];
angle(:,find(strcmp(res,'GLY'))-1)=[];


scatter(angle(1,:),angle(2,:))
box on
axis([-180 180 -180 180])
title('Ramachandran Plot for Ubiquitn Protein','FontSize',16)
xlabel('\Phi^o','FontSize',20)
ylabel('\Psi^o','FontSize',20)
grid

它给出如下输出(没有 GLY):

在此处输入图像描述

4

1 回答 1

1

我会更改此代码块以使用逻辑索引

res=strsplit(p.Sequence.ResidueNames,' ');
angle =[phi;psi];
angle(:,find(strcmp(res,'GLY'))-1)=[];

反而:

residues = strsplit(p.Sequency.ResidueNames,' ');
glycine = ismember(residues,'GLY');

angle = [phi;psi];
angleNoGLY= angle(:,~glycine);

这样做,如果你想突出甘氨酸(或任何其他残基),你可以很容易地把它叫出来:

angleGLY = angle(:,glycine);

plot(angleNoGLY(1,:),angleNoGLY(2,:),'ob')
line(angleGLY(1,:),angleGLY(2,:),'Marker','o','Color','r','LineStyle','none')
于 2014-10-10T15:41:33.007 回答