Code covered by the BSD License  

Highlights from
MGraph

MGraph

by

 

12 Apr 2007 (Updated )

Probabilistic graphical models for reconstruction of genetic regulatory networks using DNA microarra

G=MGraph_add_best_edge_to_pat(d,G,sigma,alpha,L,T0,TL)
function G=MGraph_add_best_edge_to_pat(d,G,sigma,alpha,L,T0,TL)

[nr nc]=size(d);
u0=median(d);
n=nc;
L=nr;
sigma=nc+2;
alpha=nc+2;
v=ones(1,n);

current_score=(gnt_scoring_uncompleteG(G,sigma,alpha,L,T0,TL));

disp('Add edge ...');
max_score=current_score;

while max_score<=current_score 
        %add path to graph
        A=[];
        B=[];
        temp_edges={};
        idx=1;
        pscore=[];
        done=0;
        temp_graphs={};
        temp_G=[];
        %clear Z1 Z2 ZZ1 ZZ2 ZZZ
        
     
        [A B]=find(abs(triu(setdiag(G,1)))==0); %potential add edges   
 	    idx_yes=find(A<B);
	    new_A=A(idx_yes);
	    new_B=B(idx_yes);
	    A=[];B=[];
	    A=new_A;
	    B=new_B;

        if isempty(A) break; end
        
        for i=1:length(A);
          if (G(A(i),B(i))==0 & (G(B(i),A(i))==0 & A(i)<B(i)))
             tempG=G;
             %try to find unshielded  at Y
             x=A(i);
             y=B(i);
             ZZ=[];
             Z1=[];
             Z12=[];
             Z1 = find(tempG(:,y))';
             Z12=find(tempG(y,:));
             ZZ=unique([Z1,Z12]);
             ZZ1 = mysetdiff(ZZ, x);
             
             if ~isempty(ZZ1)
                len_z=length(ZZ1);
                for ii=1:len_z
                    z=ZZ1(ii);
                    tempG=G;
                    if ( (tempG(x,y)==0 & tempG(y,x)==0)  )
                        %fprintf('%d -> %d <- %d\n', x, y, z);
                        if (x<y & z<y ) 
                            tempG(x,y) = -1; tempG(y,x) = 0;
                            tempG(z,y) = -1; tempG(y,z) = 0;
                        else
                            tempG(x,y) = 1; tempG(y,x) = 1;
                            tempG(z,y) = 1; tempG(y,z) = 1;         
                        end
                        if ~isempty(find(tempG==-1))
                            %tempG,idx
                            temp_graphs{idx}=tempG;
                            %disp('at y');
                            %pause
                            pscore(idx) = (gnt_scoring_uncompleteG(tempG,sigma,alpha,L,T0,TL));
                            idx=idx+1;
                        end
                   end
                end
            end
             
             %try to find unshielded at X
             ZZ=[];
             tempG=G;
             Z2 = find(tempG(:,x))';
             Z22=find(tempG(x,:));
             ZZ=unique([Z2,Z22]);
             ZZ2 = mysetdiff(ZZ, y);
             if ~isempty(ZZ2)
                len_z=length(ZZ2);
                for ii=1:len_z
                    z=ZZ2(ii);
                    tempG=G;
                    if ( (tempG(y,x)==0 & tempG(x,y)==0)  )
                        %fprintf('%d -> %d <- %d\n', x, y, z);
                        if (y<x & z<x) 
                            tempG(y,x) = -1; tempG(x,y) = 0;
                            tempG(z,x) = -1; tempG(x,z) = 0;
                        else
                            tempG(y,x) = 1; tempG(x,y) = 1;
                            tempG(z,x) = 1; tempG(x,z) = 1;
                        end
                        if ~isempty(find(tempG==-1))
                            %tempG , idx
                            temp_graphs{idx}=tempG;
                            %disp('at x');
                            %pause
                            pscore(idx) = (gnt_scoring_uncompleteG(tempG,sigma,alpha,L,T0,TL));
                            idx=idx+1;
                        end
                    end
                end
            end
            
            %add undirected edge
            if (tempG(x,y)==0 & tempG(y,x)==0)
                tempG=G;
                tempG(x,y)=-1;
                tempG(y,x)=0;
                if ~isempty(find(tempG==-1))
                    
                    %tempG,idx
                    %disp('not x y')
                    %pause
                    temp_graphs{idx}=tempG;
                    pscore(idx) = (gnt_scoring_uncompleteG(tempG,sigma,alpha,L,T0,TL));
                    idx=idx+1;
                end
            end
         end %end if
        end %end for
        
        if ~isempty(temp_graphs)
            [best_pscore, best_p] = max(pscore);
            best_graph=temp_graphs{best_p};
           
            current_score=best_pscore;
            current_G=best_graph;
            done=1;
        else
            done=1;
        end
        if current_score>max_score
            disp('add edge')
            max_score=current_score;
            G=current_G;
            ischanged=1;
        elseif current_score==max_score & done
            break;
        end
end

Contact us