%% last modified 2/12/2020 10:05 AM by manu
% DOCA : Domain Correlation Analysis. you can define domain regions @ Domain_Idx
% COCA : Community correlation Analysis. you can define community regions @ Domain_Idx
% First Run the Chesca Script and give the correlation matrix path here.
% Both COCA and DOCA are based on evaluation of Correlation score.
% Correlation score = [number pf correlation > cutoff]/[maximum number of correlations for assigned residues]


clear all
clc
set(0,'defaulttextInterpreter','latex') 

CorrMatFile = 'C:\Users\manu\Dropbox\M\Projects\190603_Chesca_Concise\Experiments\200204_DOmain_Correlation_Analysis_DOCA\CorrMat_out\200210_wt_human_PKI_full.mat';
CorrMatFile = 'C:\Users\manu\Dropbox\M\Projects\190603_Chesca_Concise\Experiments\200204_DOmain_Correlation_Analysis_DOCA\CorrMat_out\200210_wt_human_PKI_5-24.mat';
CorrMatFile = 'C:\Users\manu\Dropbox\M\Projects\190603_Chesca_Concise\Experiments\200204_DOmain_Correlation_Analysis_DOCA\CorrMat_out\200210_wt_mouse_PKI_5-24.mat';

ColorMap_Range = [0.2  0.5];

Correlation_CutOff  = 0.95;
colorscheme         = 'white_purple' ; %  'white_red' or 'blue_yellow_red' 'white_purple' 'blue_yellow_red'


% for DOCA
% Domain_Idx = {1:42;  43:64; 65:83; 84:113; 114:138; 139:162; 163:178; 179:193; 194:211; 212:240; 241:260; 261:297; 298:349};
% Domain_names = {'1' '2' '3' '4' '5' '6' '7' '8' '9' '10' '11' '12' '13'};
 
% for DOCA
Domain_Idx = {1:42;  43:64; 65:83;  84:98;  99:113; 114:138; 139:162; 163:178; 179:193; 194:211; 212:240; 241:260; 261:297; 298:349};
Domain_names = {' '   'I'    'II'   'III'    'IV'     'V'     'VIA'    'VIB'    'VII'    'VIII'    'IX'     'X'      'XI'     'XII'};
 
% for DOCA
% Domain_Idx = {1:32; 43:49; 50:55; 57:62; 67:77; 78:83; 84:98; 107:113; 114:119; 129:138; 139:162; 163:170; 171:175; 179:183; 184:186; 187:193; 194:201; 202:205; 206:208; 217:232; 244:253; 261:272; 289:294; 302:308};
% Domain_names = {'aA' 'b1' 'Gly-rich' 'b2' 'b3' 'aB' 'aC' 'b4' 'b5' 'aD' 'aE' 'Catalytic' 'B7' 'B8' 'DFG loop' 'B9' 'Activation ' 'P+1 loop' 'APE' 'aF' 'aG' 'aH' 'aI' 'aJ' };
% Domain_names = {'\alphaA' '\beta1' 'Gly-rich' '\beta2' '\beta3' '\alphaB' '\alphaC' '\beta4' '\beta5' '\alphaD' '\alphaE' 'Catalytic' '\beta7' '\beta8' 'DFG loop' '\beta9' 'Act' 'P+1' 'APE' '\alphaF' '\alphaG' '\alphaH' '\alphaI' '\alphaJ' };


xlabel_String = 'Domains'; 



% for COCA
% Domain_Idx = { [1 35:52  56:73  107:110  117:123  319:322  328:334]; [53:55  74:85  88  111:116  335:349];[14:34  86:87  89:101  103:106  150:164  185:190 ] ; [102  124:138  170:184  314:318  323:327];  [139:149  298:313]; [165:169  191:228  255  269  271:283]; [229:237  256:261]; [238:254] ; [261:268  284:297]};
% Domain_names = {'A' 'B' 'C' 'D' 'E' 'F' 'F1' 'G' 'H'};
% xlabel_String = 'Community'; 


%%   
load(CorrMatFile);
CorrMat     = s(1).CorrMat; 


op.MinCorrelation = Correlation_CutOff;
op.Figure = figure;


for k = 1:length(Domain_Idx)
    DomainLengths(k) = length(Domain_Idx{k}); 
end
% bar( (    DomainLengths))
%%
ND = length(Domain_Idx);

CM = CorrMat;
CM = CM - diag(diag(CM));
CDM_{ ND, ND } =[];
%% 
for k1 = 1:max([Domain_Idx{:}])
    disp(k1)
    for k2 = 1:max([Domain_Idx{:}])
        
        for kc =1:length(Domain_Idx) 
            if ~isempty(find(Domain_Idx{kc} ==  k1, 1 ))
                id1 = kc; end 
            if ~isempty(find(Domain_Idx{kc} ==  k2, 1 ))
                id2 = kc; end  
        end 
        
        CDM_{id1, id2} = [CDM_{id1, id2} CM(k1, k2)];
    end
end


for k1 = 1:size(CDM_, 1)
    for k2 = 1:size(CDM_, 2)
        temp1 = CDM_{k1, k2};

        DCOA_(k1, k2) = mean(temp1(temp1>0));  
        
        if sum(temp1 ) == 0
            DCOA_(k1, k2) = 0;
        else
            temp1(temp1  == 0) = nan;
            temp1(temp1  <  op.MinCorrelation) = 0;
            temp1(temp1  >= op.MinCorrelation) = 1;
            DCOA_(k1, k2) = nanmean(temp1 );  
        end
        DCOA_sizemap(k1, k2) = length(temp1(~isnan(temp1)));
%         length( temp1(temp1 == 1))
%         bar(temp1)
    end
end

%%

% CorrMat(Domain_Idx{3}, Domain_Idx{2})
op.MarkerSize = 22;
op.colors =  colorscheme;
% DCOA_ = DCOA_ - diag(nan*diag(DCOA_));


op.MinCorrelation = ColorMap_Range(1);
op.MaxCorrelation = ColorMap_Range(2);


Plot_CorrMat_2(DCOA_, op  )
% Plot_CorrMat_2(CorrMat, op  )

if exist('Domain_names', 'var') && length(Domain_names) == length(Domain_Idx)
    xticks(1:length(Domain_names)); xticklabels(Domain_names);
    yticks(1:length(Domain_names)); yticklabels(Domain_names);
end

xlabel(xlabel_String); ylabel(xlabel_String)  

set(gcf,'units','inches')
set(gcf, 'position', [ 1    1    op.MarkerSize/22*[3.4438    3.3271]*length(Domain_Idx)/9 ]);

[~,title_name ,~] = fileparts(CorrMatFile);
title(title_name)

ytickangle(0)
xtickangle(90)

Domain.DCorrMat      = DCOA_;
Domain.Domain_Idx   = Domain_Idx;
Domain.Domain_names = Domain_names;

if ~exist('DCorrMat_out', 'dir')
    mkdir('DCorrMat_out'); end
save(['DCorrMat_out\' title_name], 'Domain') 

[~, fileAtt ] = fileattrib(['DCorrMat_out\' title_name '.mat']); 
fprintf('DCorrMat Saved at \n %s', fileAtt.Name)
clipboard('copy',fileAtt.Name)
disp('Copied to Clipboard');
%% Functions 
function Plot_CorrMat_2(CM, varargin  )
    if ~isempty(varargin)
        op = varargin{1}; 
    else
        op.colors = 'blue_cyan_green_yellow_red';
        op.MinCorrelation = 0.95; 
    end
 

    [~, CM_Sorted_i ]= sort(CM(:), 'ascend');
    %% 
    set(op.Figure,'color','white')
    Z = CM;
    N = length(Z);
    for k = 1:(N*N) 
        ki = CM_Sorted_i(k);
%             if abs(Z(ki)) > op.MinCorrelation

                if Z(ki) < op.MinCorrelation
                    color_meter = 0;
                elseif  Z(ki) >= op.MaxCorrelation
                    color_meter = 1;
                else 
                    color_meter = (Z(ki)- op.MinCorrelation)/(op.MaxCorrelation- op.MinCorrelation); 
                end
                
                
                clr_rgb = gen_rgb_color_1(color_meter, op.colors); 
                [k1, k2] = ind2sub(size(Z), ki); 
                
               
                    if k1 == k2; clr_rgb = [0 0 0 ]+0.5;  end 
                    
                    plot(k1,k2,'s', 'MarkerFaceColor', clr_rgb, 'MarkerEdgeColor', [0 0 0 ]+0.5, 'MarkerSize', op.MarkerSize );
                    hold on
                
                
%             end 
    end
    %%

    set(gca, 'xlim',[0 N+1]  )
    set(gca, 'ylim',[0 N+1] ) 


    ytickangle(90)
    xlabel('Residue Number');ylabel('Residue Number');
    % set(gca, 'position', [0.1300    0.1100    0.7750    0.745])
end





function clr_rgb = gen_rgb_color_1(color_meter, clr_scheme)

%% Available colorschemes
% 1. blue_yellow_red    
% 2. blue_green_red
% 3. white_red
% 4. blue_cyan_green_yellow_red
% 5. green_yellow_red

%%
    switch clr_scheme 
        case 'blue_yellow_red' 
            if color_meter < 0.5 
                clr_rgb = [2*color_meter, 2*color_meter, 1-2*color_meter];
            else 
                clr_rgb = [1, 1-2*(color_meter-0.5), 0];
            end
            
        case 'blue_green_red'
            if color_meter < 0.5 
                clr_rgb = [0, 2*color_meter, 1-2*color_meter];
            else 
                clr_rgb = [2*(color_meter-0.5), 1-2*(color_meter-0.5), 0];
            end
            
         case 'blue_cyan_green_yellow_red'
            if color_meter < 0.25 
                clr_rgb = [0, 4*color_meter, 1];
            elseif color_meter < 0.5 
                clr_rgb = [0, 1, 1-4*(color_meter - 0.25)];
            elseif color_meter < 0.75
                clr_rgb = [4*(color_meter - 0.5), 1, 0];
            else 
                clr_rgb = [1, 1-4*(color_meter - 0.75), 0];
            end
            
        case 'white_red' 
            clr_rgb = [1, 1-color_meter, 1-color_meter];
            
        case 'green_yellow_red' 
            if color_meter < 0.5 
                clr_rgb = [2*color_meter, 1, 0];
            else 
                clr_rgb = [1, 1-2*(color_meter-0.5), 0];
            end
            
        case 'white_purple'
               clr_rgb = [1-color_meter*(1-0.75), 1-color_meter*(1-0), 1-color_meter*(1-0.75)];
        otherwise
            clr_rgb = [0, 0, 1];
    end

end