Skip to content

Commit

Permalink
Improved visualization of incomplete disorder-order transitions
Browse files Browse the repository at this point in the history
  • Loading branch information
@gjeschke
gjeschke committed Sep 26, 2023
1 parent 1cde59f commit f57018e
Show file tree
Hide file tree
Showing 7 changed files with 208 additions and 74 deletions.
143 changes: 96 additions & 47 deletions source/cluster_transition.m
View file
Expand Up @@ -83,6 +83,10 @@
options.mode = 'drms';
end

if ~isfield(options,'graphics') || isempty(options.graphics) || isempty(options.graphics{1})
options.graphics = {};
end

populations = [entity1.populations;entity2.populations];
[D,~,Rg] = pair_drms_matrix(entity1,options.chain1,[],entity2,options.chain2,[]);

Expand Down Expand Up @@ -388,61 +392,106 @@
entity2.name = 'ORDE';
put_pdb(entity1,options.fname1);
put_pdb(entity2,options.fname2);
[pname,bname,ext] = fileparts(options.visualize);
graphics_x = fullfile(pname,sprintf('%s_view_x.png',bname));
graphics_y = fullfile(pname,sprintf('%s_view_y.png',bname));
graphics_mx = fullfile(pname,sprintf('%s_view_-x.png',bname));
graphics_my = fullfile(pname,sprintf('%s_view_-y.png',bname));
graphics_z = fullfile(pname,sprintf('%s_view_z.png',bname));
graphics_mz = fullfile(pname,sprintf('%s_view_-z.png',bname));
if isempty(ext)
options.visualize = fullfile(pname,sprintf('%s.mmm',bname));
end
ofid = fopen(options.visualize,'wt');
fprintf(ofid,'%% MMMx transition visualization script\n');
fprintf(ofid,'new !\n');
fprintf(ofid,'pdbload %s\n',options.fname1);
fprintf(ofid,'pdbload %s\n',options.fname2);
[pname,bname] = fileparts(options.visualize);
visualize_cs = fullfile(pname,sprintf('%s_conformational_selection.mmm',bname));
ofid_cs = fopen(visualize_cs,'wt');
visualize_dp = fullfile(pname,sprintf('%s_depopulated.mmm',bname));
ofid_dp = fopen(visualize_dp,'wt');
visualize_if = fullfile(pname,sprintf('%s_induced_fit.mmm',bname));
ofid_if = fopen(visualize_if,'wt');
fprintf(ofid_cs,'%% MMMx transition visualization script for conformational selection\n');
fprintf(ofid_cs,'new !\n');
fprintf(ofid_cs,'pdbload %s\n',options.fname1);
fprintf(ofid_cs,'pdbload %s\n',options.fname2);
fprintf(ofid_dp,'%% MMMx transition visualization script for depopulated conformers\n');
fprintf(ofid_dp,'new !\n');
fprintf(ofid_dp,'pdbload %s\n',options.fname1);
fprintf(ofid_dp,'pdbload %s\n',options.fname2);
fprintf(ofid_if,'%% MMMx transition visualization script for induced fit\n');
fprintf(ofid_if,'new !\n');
fprintf(ofid_if,'pdbload %s\n',options.fname1);
fprintf(ofid_if,'pdbload %s\n',options.fname2);
for clust = 1:clusters.nc
color_index = 1 + round(100*rescaled(clust));
col = colors(color_index,:);
switch clusters.type(clust)
case 0
ofid = ofid_cs;
col = [0,0.6,0];
case 1
ofid = ofid_dp;
col = [0.75,0,0];
case 2
ofid = ofid_if;
col = [0,0,0.75];
end
% color_index = 1 + round(100*rescaled(clust));
% col = colors(color_index,:);
conformers = clusters.members{clust,1};
for c = 1:length(conformers)
pop = entity1.populations(conformers(c));
fprintf(ofid,'show [DISO]{%i} coil %6.3f\n',conformers(c),sqrt(pop));
fprintf(ofid,'color [DISO]{%i} %6.3f %6.3f %6.3f\n',conformers(c),col);
pop = entity1.populations(conformers(c))/max(entity1.populations);
if isempty(options.graphics)
fprintf(ofid,'show [DISO]{%i} coil %6.3f\n',conformers(c),sqrt(pop));
fprintf(ofid,'color [DISO]{%i} %6.3f %6.3f %6.3f\n',conformers(c),col);
else
[ngcmd,~] = size(options.graphics);
for gcmd = 1:ngcmd
fprintf(ofid,'%s [DISO]{%i}%s %s\n',options.graphics{gcmd,1},conformers(c),options.graphics{gcmd,2},options.graphics{gcmd,3});
end
fprintf(ofid,'transparency [DISO]{%i} %5.3f\n',conformers(c),pop);
end
end
conformers = clusters.members{clust,2};
for c = 1:length(conformers)
pop = entity2.populations(conformers(c));
fprintf(ofid,'show [ORDE]{%i} coil %6.3f\n',conformers(c),sqrt(pop));
fprintf(ofid,'color [ORDE]{%i} %6.3f %6.3f %6.3f\n',conformers(c),col);
pop = entity2.populations(conformers(c))/max(entity2.populations);
if isempty(options.graphics)
fprintf(ofid,'show [ORDE]{%i} coil %6.3f\n',conformers(c),sqrt(pop));
fprintf(ofid,'color [ORDE]{%i} %6.3f %6.3f %6.3f\n',conformers(c),col);
else
[ngcmd,~] = size(options.graphics);
for gcmd = 1:ngcmd
fprintf(ofid,'%s [ORDE]{%i}%s %s\n',options.graphics{gcmd,1},conformers(c),options.graphics{gcmd,2},options.graphics{gcmd,3});
end
fprintf(ofid,'transparency [ORDE]{%i} %5.3f\n',conformers(c),pop);
end
end
end
fprintf(ofid,'color [DISO]{:}%s darkgrey\n',options.superimposed);
fprintf(ofid,'color [ORDE]{:}%s darkgrey\n',options.superimposed);
fprintf(ofid,'view x\n');
fprintf(ofid,'detach\n');
fprintf(ofid,'zoom out\n');
fprintf(ofid,'copy %s png\n',graphics_x);
fprintf(ofid,'view y\n');
fprintf(ofid,'zoom out\n');
fprintf(ofid,'copy %s png\n',graphics_y);
fprintf(ofid,'view -x\n');
fprintf(ofid,'detach\n');
fprintf(ofid,'zoom out\n');
fprintf(ofid,'copy %s png\n',graphics_mx);
fprintf(ofid,'view -y\n');
fprintf(ofid,'zoom out\n');
fprintf(ofid,'copy %s png\n',graphics_my);
fprintf(ofid,'view z\n');
fprintf(ofid,'detach\n');
fprintf(ofid,'zoom out\n');
fprintf(ofid,'copy %s png\n',graphics_z);
fprintf(ofid,'view -z\n');
fprintf(ofid,'zoom out\n');
fprintf(ofid,'copy %s png\n',graphics_mz);
fclose(ofid);
if isempty(options.graphics)
fprintf(ofid_cs,'color [DISO]{:}%s darkgrey\n',options.superimposed);
fprintf(ofid_cs,'color [ORDE]{:}%s darkgrey\n',options.superimposed);
fprintf(ofid_dp,'color [DISO]{:}%s darkgrey\n',options.superimposed);
fprintf(ofid_dp,'color [ORDE]{:}%s darkgrey\n',options.superimposed);
fprintf(ofid_if,'color [DISO]{:}%s darkgrey\n',options.superimposed);
fprintf(ofid_if,'color [ORDE]{:}%s darkgrey\n',options.superimposed);
end
types = {'conformational_selection' 'depopulated' 'induced_fit'};
ofids = [ofid_cs,ofid_dp,ofid_if];
for tp = 1:3
ofid = ofids(tp);
my_type = types{tp};
fprintf(ofid,'view x\n');
fprintf(ofid,'detach\n');
fprintf(ofid,'zoom out\n');
fprintf(ofid,'copy %s png\n',fullfile(pname,sprintf('%s_%s_view_x.png',bname,my_type)));
fprintf(ofid,'view y\n');
fprintf(ofid,'zoom out\n');
fprintf(ofid,'copy %s png\n',fullfile(pname,sprintf('%s_%s_view_y.png',bname,my_type)));
fprintf(ofid,'view -x\n');
fprintf(ofid,'detach\n');
fprintf(ofid,'zoom out\n');
fprintf(ofid,'copy %s png\n',fullfile(pname,sprintf('%s_%s_view_-x.png',bname,my_type)));
fprintf(ofid,'view -y\n');
fprintf(ofid,'zoom out\n');
fprintf(ofid,'copy %s png\n',fullfile(pname,sprintf('%s_%s_view_-y.png',bname,my_type)));
fprintf(ofid,'view z\n');
fprintf(ofid,'detach\n');
fprintf(ofid,'zoom out\n');
fprintf(ofid,'copy %s png\n',fullfile(pname,sprintf('%s_%s_view_z.png',bname,my_type)));
fprintf(ofid,'view -z\n');
fprintf(ofid,'zoom out\n');
fprintf(ofid,'copy %s png\n',fullfile(pname,sprintf('%s_%s_view_-z.png',bname,my_type)));
end
fclose(ofid_cs);
fclose(ofid_dp);
fclose(ofid_if);
end

function my_map = cluster_color_map
Expand Down
28 changes: 15 additions & 13 deletions source/domain_partitioning.m
View file
Expand Up @@ -72,6 +72,8 @@
% This file is a part of MMMx. License is MIT (see LICENSE.md).
% Copyright(c) 2023: Gunnar Jeschke

PAE2stddev = 1; % factor for converting PAE to standrad deviation

my_base = load('deer_kernel.mat','base','t','r');
Pake_base.kernel = my_base.base-ones(size(my_base.base)); % kernel format for pcf2deer
Pake_base.t = my_base.t;
Expand Down Expand Up @@ -432,8 +434,8 @@
end
site2 = ref2(kr2).site(poi+1:end);
rmean = norm(coor1{1}-coor2{1});
sigr = max([entity.pae(res1,res2),entity.pae(res2,res1)])/2; % assuming CI95 = 2 sigma
if sigr <= 15
sigr = max([entity.pae(res1,res2),entity.pae(res2,res1)])*PAE2stddev;
if sigr <= 30*PAE2stddev
fprintf(fid,' %s %s %4.1f %4.1f %% ref %i.%i to %i.%i\n',site1,site2,rmean,sigr,nr1,kr1,nr2,kr2);
sim_deer_file = sprintf('sim_deer_CA_CA_%s_%s.csv',site1,site2);
[t,ff] = get_ff(rmean,sigr,Pake_base);
Expand Down Expand Up @@ -461,10 +463,10 @@
id_pae_lower = entity.pae(domains(d2,1):domains(d2,2),domains(d1,1):domains(d1,2))';
id_sigr = id_pae_upper;
id_sigr(id_sigr < id_pae_lower) = id_pae_lower(id_sigr < id_pae_lower);
id_sigr = id_sigr/2; % assuming tha PAE is a CI95
id_sigr = id_sigr*PAE2stddev;
min_err = min(min(id_sigr));
aux = 0;
while min_err <= 15 && aux < 20
while min_err <= 30*PAE2stddev && aux < 20
[colmin,rows] = min(id_sigr);
[min_err,col] = min(colmin);
row = rows(col);
Expand Down Expand Up @@ -532,10 +534,10 @@
id_pae_lower = entity.pae(domains(refd,1):domains(refd,2),res)';
id_sigr = id_pae_upper;
id_sigr(id_sigr < id_pae_lower) = id_pae_lower(id_sigr < id_pae_lower);
id_sigr = id_sigr/2; % assuming tha PAE is a CI95
id_sigr = id_sigr*PAE2stddev;
[min_err,refpoi] = min(id_sigr);
if min_err <= 15
if min_err/unrestrained <= 5
if min_err <= 30*PAE2stddev
if min_err/unrestrained <= 10*PAE2stddev
coor1 = get_label(entity,'atom.CA','coor',sprintf('(A)%i',refpoi + domains(refd,1)-1));
coor2 = get_label(entity,'atom.CA','coor',sprintf('(A)%i',res));
rmean = norm(coor1{1}-coor2{1});
Expand Down Expand Up @@ -590,10 +592,10 @@
id_pae_lower = entity.pae(domains(refd,1):domains(refd,2),res)';
id_sigr = id_pae_upper;
id_sigr(id_sigr < id_pae_lower) = id_pae_lower(id_sigr < id_pae_lower);
id_sigr = id_sigr/2; % assuming tha PAE is a CI95
id_sigr = id_sigr*PAE2stddev;
[min_err,refpoi] = min(id_sigr);
if min_err <= 15
if min_err/unrestrained <= 5
if min_err <= 30*PAE2stddev
if min_err/unrestrained <= 10*PAE2stddev
coor1 = get_label(entity,'atom.CA','coor',sprintf('(A)%i',refpoi + domains(refd,1)-1));
coor2 = get_label(entity,'atom.CA','coor',sprintf('(A)%i',res));
rmean = norm(coor1{1}-coor2{1});
Expand Down Expand Up @@ -644,10 +646,10 @@
id_pae_lower = entity.pae(domains(refd,1):domains(refd,2),res)';
id_sigr = id_pae_upper;
id_sigr(id_sigr < id_pae_lower) = id_pae_lower(id_sigr < id_pae_lower);
id_sigr = id_sigr/2; % assuming tha PAE is a CI95
id_sigr = id_sigr*PAE2stddev;
[min_err,refpoi] = min(id_sigr);
if min_err <= 15
if min_err/unrestrained <= 5
if min_err <= 30*PAE2stddev
if min_err/unrestrained <= 10*PAE2stddev
coor1 = get_label(entity,'atom.CA','coor',sprintf('(A)%i',refpoi + domains(refd,1)-1));
coor2 = get_label(entity,'atom.CA','coor',sprintf('(A)%i',res));
rmean = norm(coor1{1}-coor2{1});
Expand Down
26 changes: 24 additions & 2 deletions source/get_backbone.m
View file
@@ -1,8 +1,8 @@
function backbone = get_backbone(entity,chain,conformer,range)
function backbone = get_backbone(entity,chain,conformer,range,options)
%
% GET_BACKBONE Retrieve backbone coordinates for a chain (one conformer)
%
% backbone = GET_BACKBONE(entity,chain,conformer,range)
% backbone = GET_BACKBONE(entity,chain,conformer,range,options)
% Provides coordinates and corresponding residue numbers for all CA and
% C4' atoms in one conformer of a chain, optionally, a residue range can
% be specified
Expand All @@ -14,6 +14,8 @@
% supply this argument
% conformer conformer (model) number, defaults to 1
% range double(1,2), optional residue range, defaults to all residues
% options .full flag, if true, N and C coordinates for amino acid
% residues are also extracted
%
% OUTPUT
% backbone struct with fields, all fields are empty if the chain is
Expand All @@ -24,6 +26,9 @@
% .nt (1,nnt) double, numbers of residues with a C4' atom
% .C4p (nnt,3) coordinate array for C4' atoms
% .nuc string with length ntt, nucleotide single-letter codes
% if options.full is true
% .N (naa,3) coordinate array for N atoms
% .C (naa,3) coordinate array for C atoms
%

% This file is a part of MMMx. License is MIT (see LICENSE.md).
Expand All @@ -32,6 +37,10 @@
% load monomer definitions
defs = load('monomers.mat');

if ~exist('options','var') || isempty(options) || ~isfield(options,'full')
options.full = false;
end


% initialize empty output
backbone.aa = [];
Expand All @@ -55,6 +64,10 @@
maxres = 10000; % maximum number of residues
aa = zeros(1,maxres);
CA = zeros(maxres,3);
if options.full
N = zeros(maxres,3);
C = zeros(maxres,3);
end
res = blanks(maxres);
aa_poi = 0; % counter for detected CA atoms
nt = zeros(1,maxres);
Expand All @@ -78,6 +91,10 @@
else
res(aa_poi) = '?';
end
if options.full
N(aa_poi,:) = entity.xyz(entity.(chain).(residue).N.tab_indices(conformer),:);
C(aa_poi,:) = entity.xyz(entity.(chain).(residue).C.tab_indices(conformer),:);
end
end
if isfield(entity.(chain).(residue),'C4_') % is there a C4' atom?
nt_poi = nt_poi + 1;
Expand All @@ -100,3 +117,8 @@
backbone.nt = nt(1:nt_poi);
backbone.C4p = C4p(1:nt_poi,:);
backbone.nuc = nuc(1:nt_poi);

if options.full
backbone.N = N(1:aa_poi,:);
backbone.C = C(1:aa_poi,:);
end

0 comments on commit f57018e

Please sign in to comment.