function Graph=ReconstructionGRNfromGeneticalGenomicsData(P,T,closestMarker,tQT,dmin,TRANSWESD)
% Function for Reconstructing Gene Regulatory Networks From Genetical 
% Genomics Data (see paper of Flassig et al.)
%
% R.J. Flassig, S.Heise, K. Sundmacher and S.Klamt
% Max Planck Institute for Complex Technical Systems Magdeburg
%
%
% output:   Graph - struct defining a weighted, signed digraph 
%		Graph.adjpos - adjacency matrix storing edge weights of positive edges
%		Graph.adjneg - adjacency matrix storing edge weights of negative edges
%		Graph.List   - matrix with sorted edge list [regulator_gene target_gene weight sign]
%
% input:    P (r x m)-matrix            - matrix with genotype of marker (0's and 1's; all other values are considered as 'not available') 
%           T (r x g)-matrix            - matrix with expression of genes
%           closestMarker (g)-vector    - vector giving the closest marker to each gene 
%           tQT                         - threshold for considering Genotype-Phenotype-Correlation 
%                                         as signifikant
%           dmin                        - theshold for considering Genotype-Genotype-Correlation as significant 
%					  (for detecting genetic linkage between markers)
%           TRANSWESD                   - parameter whether to use (TRANSWESD=1)
%                                         transitive eduction (needs CellNetAnalyzer 
%					  installed) and or not (TRNASWESD=0)
%
%           r:  number of RILs
%           g: number of Genes
%           m: number of Marker



% undefined P = NaN, removed for calculation of corrMatGeno and GeneticLinkage
remove=(P~=0)+(P~=1);
P(remove~=1)=NaN;

%% 1. preprocessing and correlations
%___________________________________________________________________________________________________________________________
disp('Computing correlation matrices...');
numsEXPRESS=size(T,2);
% Linkage Analysis 
rQQ=zeros(numsEXPRESS,numsEXPRESS);          %GeneticLinkage (gxg): matrix with genetic correlation
   for i=1:numsEXPRESS 
       g1=P(:,closestMarker(i));
       for j=1:numsEXPRESS
           g2=P(:,closestMarker(j)); 
           missing=unique([find(isnan(g1)); find(isnan(g2))]);
           g1_new=g1;
           g2_new=g2;
           g1_new(missing)=[];
           g2_new(missing)=[];
           a=corrcoef(g1_new,g2_new);
           rQQ(j,i)=a(1,2);
       end
   end   
class=[];
    for i=1:numsEXPRESS
        class(i,:)=(rQQ(:,i)>=dmin);    % class with genetic linked genes
    end  

% correlation and weightmatrix
rQT=zeros(numsEXPRESS,numsEXPRESS);             %rQT (gxg): correlation between expression target and genotype marker
   for i=1:numsEXPRESS 
       target=T(:,i);
       for j=1:numsEXPRESS
           regulator=P(:,closestMarker(j)); 
           missing=find(isnan(regulator));
           target_new=target;
           target_new(missing)=[];
           regulator(missing)=[];
           a=corrcoef(target_new,regulator);
           rQT(j,i)=a(1,2);
       end
   end    
rTT=corrcoef(T);                           %rTT (gxg): correlation between expression of genes
W=(abs(rQT)+abs(rTT))/2;     %W (gxg): mean of corrMatGeno and corrMatExp 
 


%% 2. Generate Perturbationgraph
%______________________________________________________________________________________________________________________________
% edge detection -> G1 for |rQT|>tQT
disp('Generating Perturbation graph...');
G1=(abs(rQT)>tQT);                             
disp(['number of edges in G1: ',num2str(length(find(G1==1)))]) 

% sign and weights -> G2
G2adjpos=zeros(numsEXPRESS,numsEXPRESS);
G2adjneg=zeros(numsEXPRESS,numsEXPRESS);
   for i=1:numsEXPRESS 
       for j=1:numsEXPRESS
           if G1(i,j)==1
               if rTT(i,j)>0
                    G2adjpos(i,j)=W(i,j);
               elseif rTT(i,j)<0
                    G2adjneg(i,j)=W(i,j);
               end
           end
       end
   end

% identify regulator clusters (QTL) -> G3 and select candidate regulators -> G4
h= waitbar(0,'identify QTL and select candidate regulators');
 
SUM=0;
number_edge=[];
G4=zeros(numsEXPRESS,numsEXPRESS);           
    for i=1:numsEXPRESS
        QTL=zeros(numsEXPRESS,1);
        numQTL=0; 
        potentialmarker=find(G1(:,i)==1);
         waitbar(i / numsEXPRESS)
        % identify regulator QTL
        for j=potentialmarker'
            zw1=intersect(find(class(j,:)),potentialmarker);
            if  any(QTL(zw1)~=0) 
                zw2=find(QTL(zw1)>0);
                MIN=min(QTL(zw1(zw2)));
                QTL(zw1)=MIN;
            else
                numQTL=numQTL+1;
                QTL(zw1)=numQTL; 
            end
        end   
        
        %select cantidate regulators
        QTLSize=unique(QTL)';
        QTLSize(find(QTLSize==0))=[];
        for l=QTLSize        
             v=find(QTL==l)'; 
             number_edge=[number_edge length(v)];
             v2=abs(W(v,i));   
             ER=v(find(v2==max(v2)));     
             G4(ER,i)=1;
             SUM=SUM+1;
             
         end
    end
close(h)
disp(['number of QTLs: ',num2str(SUM)])
disp(['average of edges in QTL: ',num2str(mean(number_edge))])

        % remove cis-edges
        SUM2=0;
        for k=1:numsEXPRESS
            for j=1:numsEXPRESS
                if (k==j)&&(G4(j,k)==1)
                    G4(j,k)=0;
                    SUM2=SUM2+1;
                end
            end
        end
        disp(['number of removed cis effects: ',num2str(SUM2)])

% sign and weights
G4adjpos=zeros(numsEXPRESS,numsEXPRESS);
G4adjneg=zeros(numsEXPRESS,numsEXPRESS);
   for i=1:numsEXPRESS 
       for j=1:numsEXPRESS
           if G4(i,j)==1
               if rTT(i,j)>0
                    G4adjpos(i,j)=W(i,j);
               elseif rTT(i,j)<0
                    G4adjneg(i,j)=W(i,j);
               end
           end
       end
   end
Graph.adjpos=G4adjpos;
Graph.adjneg=G4adjneg;
Graph.W=W;   %just for evaluationsscript DREAM5 

%% 3. remove indirect effects
if TRANSWESD
 disp('Transitive reduction ...');
 if(exist('CNATranswesd'))
%% remove indirect path effects (e.g. with TRANSWESD) -> G5
% CellNetAnalyzer API function: 'CNATranswesd'
% --------------------------------------------
% --> Transitive reduction in weighted signed graphs (useful as false positive 
%     reduction method in reverse engineering of cellular interaction graphs).
%     This algorithm has been described in the followinhg reference:
%     Klamt S, Flassig R, Sundmacher K. 2010. TRANSWESD: inferring
%     cellular networks with transitive reduction. Bioinfomatics 26:2160-2168.
%
% input
% epsilon: confidence factor. An edge from i to k with weight w 
%    can be explained by a path from i to k with weight z
%    (and is therefore deleted) if z<w*epsilon. Default: 1.
%    Here the length of a path is computed as the maximum weight of all its edges. 
%    For further explanations see reference given above.
%  fullcheck: whether all path lengths have to be recalculated after
%    removing an edge (1) or not (0). Default: 0.95.
%  pathexact: whether path lengths have to be calculated exactly (1)
%    or approximately (0). Default: 1.
%  acyclic: if 1, only edges between nodes from different components
%    will be considered for removal by transitive reduction. Default: 0.
%  maxpathl: maximum number of edges that a path may have if it is 
%    used to explain an edge. Choose Inf if any appropriate path 
%    (see parameter epsilon above) can explain an edge. Default: Inf.
%    If 2<maxpathl<INF (and maxpathl not too large) it is recommended
%    to use this feature only in combination with pathexact=fullcheck=1. 
% 
	epsilon=0.95 ;          
	fullcheck=0;            
	pathexact=0;            
	acyclic=0;              
	maxpath=Inf;
	%transform weights for TRANSWESD
	Graph.adjpos(find(Graph.adjpos))=1-W(find(Graph.adjpos));
	Graph.adjneg(find(Graph.adjneg))=1-W(find(Graph.adjneg));

	G5=CNATranswesd(Graph,epsilon,fullcheck,pathexact,acyclic,maxpath);

	%transform to weights
	Graph.adjpos=G5.recadjpos;
	Graph.adjneg=G5.recadjneg;
	Graph.adjpos(find(Graph.adjpos))=W(find(Graph.adjpos));
	Graph.adjneg(find(Graph.adjneg))=W(find(Graph.adjneg));
 else
      disp('Could not find function CNATranswesd of CellNetAnalyzer.');
      disp('Download CellNetAnalyzer or choose parameter TRANSWESD=0.');
 end % if exist
end %TRANSWESD


%% 4. Sorted edge list
disp('Sorted edge list ...');
edges=[];
notedges=[];
h2= waitbar(0,'sorted edge list');
for g=1:numsEXPRESS 
    waitbar(g / numsEXPRESS)
    % edges
	zw=find(Graph.adjpos(g,:));  
	if(length(zw))         
        edges=[edges; g*ones(length(zw),1) zw' abs(W(g,zw))' sign(rTT(g,zw')')];                                  
    end      
	zw=find(Graph.adjneg(g,:));  
	if(length(zw)) 
        edges=[edges; g*ones(length(zw),1) zw' abs(W(g,zw))' sign(rTT(g,zw')')];
    end
    % notedges
    zw=find(~Graph.adjpos(g,:) & ~Graph.adjneg(g,:));
    zw2=find(zw==g);            
    zw(zw2)=[];                 
    if(length(zw))     
        notedges=[notedges; g*ones(length(zw),1) zw' abs(W(g,zw))' sign(rTT(g,zw')')];
    end
end
close(h2)
% Order
if(numel(edges))                 
	[zw zw1]=sort(-edges(:,3));
	edges=edges(zw1,:);    
end
	[zw zw1]=sort(-notedges(:,3));      
	notedges=notedges(zw1,:); 
    
Graph.List=[edges;notedges];
end %function