%%Residue interaction network
% ***Please run ResNetworkAtom.m file first***

%close all hidden; clc; 
pdb1 = '3JXW'; pdb2 = '3MA3';
cres1 = load([pdb1,'_residues.dat']);
cres2 = spconvert(cres1);
cres3 = tril(cres2); %Lower diagonal only

%consider only 'long range' interactions > 3 residues
[x,y] = find(cres3);
longrangeA = zeros(size(cres3));
for i = 1:length(x)
     %ligand is always 'long range'
    if ((abs(x(i)-y(i))>3)) | (x(i) == size(cres3))
        longrangeA(x(i),y(i)) = cres3(x(i),y(i));
    end
end

%% visualising interactions using Biograph function

r = pmin:pmax;
nodesid = num2str(r'); %nodes id
bg = biograph(cres3,nodesid);
get(bg.nodes,'ID');
set(bg,'ShowArrows','off','LayoutType','equilibrium');
set(bg.nodes,'Shape','circle','Color',[0.9,0.9,1])

% colour edges representing long range interactions
[x,y] = find(longrangeA);
for i = 1:length(x)
    k = int2str((x(i)+pmin-1)); 
    l = int2str((y(i)+pmin-1));
    edgesid = getedgesbynodeid(bg,k,l);
    set(edgesid,'LineColor',[1 0 0],'LineWidth',2.5);
end

bg.Nodes(end).ID = ['L','i']; %last residue is renames as 'Ligand'
set(bg.Nodes(end),'Color',[0.5,0.5,0.7],'Shape','diamond','FontSize',15);

view (bg);

%% view long range interactions only
%bg2 = biograph(longrange,rid);
%get(bg2.nodes,'ID');
%set(bg2,'ShowArrows','off','LayoutType','equilibrium');
%set(bg2.nodes,'Shape','circle','Color',[0.9,0.9,1])
%view (bg2);

%% analyse differences in interactions

%first structure is 'longrangeA'

%similarly, find 'longrangeB'
b1 = load([pdb2,'_residues.dat']);
b2 = tril(spconvert(b1));

[x,y] = find(b2);
longrangeB = zeros(size(b2));
for i = 1:length(x)
     %ligand is always 'long range'
    if ((abs(x(i)-y(i))>3)) | (x(i) == size(cres3))
        longrangeB(x(i),y(i)) = b2(x(i),y(i));
    end
end

ab = longrangeA - longrangeB;

   %Interactions in A but not B
   [i1, k1] = find((ab==1));
   %update residue numbers
   i1 = i1+pmin-1; k1 = k1+pmin-1;
   
   %Interactions in B but not A
   [i2, k2] = find((ab==-1));
   %update residue numbers
   i2 = i2+pmin-1; k2 = k2+pmin-1;
   
    %display residue numbers
    disp (['Interactions in ', pdb1 ' but not in ', pdb2]);
    disp([i1,k1]);

    disp (['Interactions in ', pdb2, ' but not in ' pdb1]);
    disp([i2,k2]);