我想使用“球和棍子”绘制分子的 3D 结构图。Matlab 生物信息学工具箱展示了如何使用例如:
ubi = getpdb('1ubi');
h1 = molviewer(ubi);
有没有办法将 pdb 文件加载到 matlab 中并\或在没有生物信息学工具箱的情况下可视化 3d 结构?
以下是在 Matlab 中绘制 pdb 文件的两个示例:
来自上述 (1) 的示例代码:
function draw_protein(pdbfile)
pdb=fopen(pdbfile,'r');
l=fgets(pdb);
q=1;
resolution = 5;
while l~=-1
if (l(1:4) == 'ATOM') & (l(14) ~= 'H' & l(13) ~= 'H')
if (strncmp('PHE',l(18:20),3)==1 & strncmp('C',l(14),1)==1) | ...
(strncmp('TYR',l(18:20),3)==1 & strncmp('C',l(14),1)==1) | ...
(strncmp('TRP',l(18:20),3)==1 & strncmp('C',l(14),1)==1)
if strncmp('CG',l(14:15),2)==1 | strncmp('CD ',l(14:16),3)==1 | ...
strncmp('CE',l(14:15),2)==1 | strncmp('CZ',l(14:15),2)==1 | ...
strncmp('CD2',l(14:16),3)==1
r(q,4)=1.85;
elseif strncmp('C',l(14),1)==1
r(q,4)=2.0;
end
elseif strncmp('N',l(14),1)==1
r(q,4)=1.5;
elseif strncmp('C',l(14),1)==1
r(q,4)=2.0;
elseif strncmp('O',l(14),1)==1
r(q,4)=1.4;
elseif strncmp('S',l(14),1)==1
r(q,4)=1.85;
elseif strncmp('POT',l(14:16),3)==1
r(q,4)=1.49
%r(q,4)=3.31
elseif (l(14) == 'H' | l(13) == 'H')
if l(64) == '1'
r(q,4) = 1;
%else
% r(q,4) = 0;
end
else
%display('Unknown atom type')
l
r(q,4)=2.0; % CALL UNKNOWN ATOM A CARBON
end
r(q,1)=str2num(l(31:38)); % x
r(q,2)=str2num(l(39:46)); % y
r(q,3)=str2num(l(47:54)); % z
if strncmp('ARG',l(18:20),3)==1 | strncmp('LYS',l(18:20),3)==1 | strncmp('HSP',l(18:20),3)==1
r(q,5:7)=[0.2 0.1 0.90]; %blue for positively charged
elseif strncmp('THR',l(18:20),3)==1 | strncmp('ASN',l(18:20),3)==1 | strncmp('SER',l(18:20),3)==1 | strncmp('GLN',l(18:20),3)==1
r(q,5:7)=[0.2 0.90 0.1]; %green for uncharged polar
elseif strncmp('ASP',l(18:20),3)==1 | strncmp('GLU',l(18:20),3)==1
r(q,5:7)=[0.90 0.2 0.1]; %red for negatively charged
else
r(q,5:7)=[1 0.95 0.9]; %white for hydrophobic
end
q=q+1;
end
l=fgets(pdb);
end
display 'done loading pdb coordinates';
figure;
hold all;
for i=1:length(r(:,1))
draw_sphere(r(i,1),r(i,2),r(i,3),r(i,4),r(i,5:7),resolution);
end
light
camlight('right');
end
function draw_sphere(xd,yd,zd,rad,color,resolution)
n = resolution;
% -pi to theta to pi is a row vector.
% -pi/2 to phi to pi/2 is a column vector.
theta = (-n:2:n)/n*pi;
phi = (-n:2:n)'/n*pi/2;
cosphi = cos(phi); cosphi(1) = 0; cosphi(n+1) = 0;
sintheta = sin(theta); sintheta(1) = 0; sintheta(n+1) = 0;
x = cosphi*cos(theta);
y = cosphi*sintheta;
z = sin(phi)*ones(1,n+1);
surf(rad*x+xd,rad*y+yd,rad*z+zd,'EdgeColor','none','FaceColor',color,'FaceLighting','phong','AmbientStrength',0.1,'DiffuseStrength',0.8,'SpecularStrength',0.2);
end