Matlab: как выделить остатки GLYCINE на моем графике Рамачандрана?

Я пытаюсь Matlab построить график рамачандрана без использования встроенной команды. Я тоже преуспел. Теперь я хотел найти только 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: как избежать определенных аминокислотных последовательностей из белка, чтобы построить график Рамачандрана? имеет ответ, который имеет немного другой график.

Модифицированный код выглядит следующим образом:

% 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), как показано ниже:

введите здесь описание изображения

dexterdev 10.10.2014 источник

comment

Как написано, это больше вопрос биохимии, чем вопрос программирования. Если вы хотите удалить глицин, вам придется проверить свои молекулы на некоторые свойства, уникальные для глицина (например, psi‹0 и phi›5 или что-то в этом роде). В связанном вопросе есть массив с именами остатков, и пользователь выбрал индексы, в которых есть строка GLY. Непонятно, что у вас такой массив символов. - craigim 10.10.2014

comment

@craigim: Во-первых, я хотел бы проверить код. Я изменил код с такими символами, как «GLY». Но куда мне написать этот вопрос. - dexterdev 10.10.2014

comment

Вы должны уточнить свой код и сделать это вопросом кодирования. Например, у нас нет доступа к 1UBQ.pdb, поэтому мы не можем узнать, что находится в файле. Подробные комментарии перед каждым блоком кода не только помогут нам понять ваш код и понять, что представляет собой каждая переменная и какова ее цель, но и помогут вам в будущем, когда вам понадобится вернуться к коду для какой-либо другой цели. - craigim 10.10.2014

Ответы (1)

arrow_upward
1
arrow_downward

Я бы изменил этот блок кода, чтобы использовать логическое индексирование

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')

craigim 10.10.2014

Matlab: как выделить остатки GLYCINE на моем графике Рамачандрана?

Ответы (1)

Вопросы по теме