Description of ArrowheadProcess

Home > BiKEGG > ArrowheadProcess.m

# ArrowheadProcess

## PURPOSE

ArrowheadProcess

## SYNOPSIS

function [ArX,ArY,Cntrs,Director] = ArrowheadProcess(Xin1,Yin1,min_x,min_y,idx,basemap,JumpData,RefOverlapData,ID,flx_idx)

## DESCRIPTION

     ArrowheadProcess
     is a subfunction of NetDraw for identifiying the
     appropriate coordinates for drawing arrowheads on metabolic maps.
     Inputs:
     Xin1, Yin1: Reaction coordinates
     min_x, min_y: Minimum of coordinates on the current customized map
     basemap: KGML data of KEGG global metabolic pathway
     JumpData: Details (ID) for overlapping flux carrying reactions
     RefOverlapData: Details (ID) for all reactions in the current customized
     map.
     ID: Reaction directionality based on either KEGG or BiGG (string)
     flx_idx: Details (ID) for all flux carrying reactions

     Outputs:
     ArX,ArY: Appropriate coordinates for overlaying arrowheads on map.
     Cntrs: Centers' coordinates corresponding to ArX and ArY for MapTrimmer
     or MapTrimmerF functions.
     Director: Reaction directionality for arrowhead function.

## CROSS-REFERENCE INFORMATION

This function calls:

- KEGG2BiGGMets KEGG2BiGGMets

This function is called by:

- NetDraw NetDraw
- NetPoster NetPoster

## SUBFUNCTIONS

- function AllCoords = getBiGGrxns(idx,basemap)
- function metCoords = getCoords(basemap,indx,Stat)
- function [ProdCoords,SubCoords] = getCoords4KEGG(basemap,indx)
- function [ArX,ArY,Cntrs,Director] = PostMod(ArX,ArY,Cntrs,Director,

## SOURCE CODE

    0001 function [ArX,ArY,Cntrs,Director] = ArrowheadProcess(Xin1,Yin1,min_x,min_y,...
    0002     idx,basemap,JumpData,RefOverlapData,ID,flx_idx)
    0003 % ArrowheadProcess
    0004 % is a subfunction of NetDraw for identifiying the
    0005 % appropriate coordinates for drawing arrowheads on metabolic maps.
    0006 % Inputs:
    0007 % Xin1, Yin1: Reaction coordinates
    0008 % min_x, min_y: Minimum of coordinates on the current customized map
    0009 % basemap: KGML data of KEGG global metabolic pathway
    0010 % JumpData: Details (ID) for overlapping flux carrying reactions
    0011 % RefOverlapData: Details (ID) for all reactions in the current customized
    0012 % map.
    0013 % ID: Reaction directionality based on either KEGG or BiGG (string)
    0014 % flx_idx: Details (ID) for all flux carrying reactions
    0015 %
    0016 % Outputs:
    0017 % ArX,ArY: Appropriate coordinates for overlaying arrowheads on map.
    0018 % Cntrs: Centers' coordinates corresponding to ArX and ArY for MapTrimmer
    0019 % or MapTrimmerF functions.
    0020 % Director: Reaction directionality for arrowhead function.
    0021
    0022 % O. Jamialahmadi
    0023 % TMU, Chem. Eng. Dept., Biotech. Group
    0024 % July 2016
    0025 %--------------------------------------------------------------------------
    0026
    0027 ArX = ({}); ArY = ArX; Cntrs = ArX; Jumpidx = []; Director = ArX;
    0028
    0029 if ~isempty(RefOverlapData)
    0030     [~,Idx1,~] = unique(RefOverlapData);
    0031     Chck1 = setdiff(1:size(RefOverlapData,1),Idx1);
    0032     Jumpidx = JumpData.Id;
    0033     RefJumpidx = Jumpidx;
    0034     for m2 = 1:numel(Chck1)
    0035         WhrC = find(ismember(RefOverlapData,RefOverlapData{Chck1(m2)}));
    0036         RedundantIdx = idx(WhrC);
    0037         if any(ismember(RedundantIdx,RefJumpidx))
    0038             Addme = setdiff(RedundantIdx,RefJumpidx);
    0039             for m3 = 1:numel(Addme)
    0040                 CurrLen = numel(Jumpidx);
    0041                 Jumpidx{CurrLen+1} = Addme{m3};
    0042             end
    0043         end
    0044     end
    0045     Jumpidx = unique(Jumpidx);
    0046 end
    0047 if strcmp(ID,'bigg') % Rxns directionality based on BiGG-------------------
    0048     Bgcoord = getBiGGrxns(idx,basemap);
    0049 end
    0050 % -------------------------------------------------------------------------
    0051 for m1 = 1:numel(idx)
    0052     if ~isempty(Jumpidx)
    0053         if any(ismember(idx{m1},Jumpidx))
    0054             ArX{m1} = ''; ArY{m1} = ''; Cntrs{m1} = ''; Director{m1} = '';
    0055             continue
    0056         end
    0057     end
    0058     [~,in1] = intersect(basemap.rc.rxnid,str2double(idx{m1}(7:end)));
    0059     Xin = Xin1{m1};
    0060     Yin = Yin1{m1};
    0061     if isempty(in1)
    0062         ArX{m1} = ''; ArY{m1} = ''; Cntrs{m1} = ''; Director{m1} = '';
    0063         continue
    0064     end
    0065     if strcmp(ID,'bigg')
    0066         if strcmp(Bgcoord{m1},'None')
    0067             continue
    0068         end
    0069         if ~isempty(Bgcoord{m1})
    0070             ProdCoords = Bgcoord{m1}.P; SubCoords = Bgcoord{m1}.S;
    0071         else
    0072             [ProdCoords,SubCoords] = getCoords4KEGG(basemap,in1);
    0073         end
    0074     else
    0075         [ProdCoords,SubCoords] = getCoords4KEGG(basemap,in1);
    0076     end
    0077
    0078     if isempty(SubCoords)
    0079         MixCoords = ProdCoords;
    0080         Director{m1} = ones(size(MixCoords,1),1);
    0081     else
    0082         MixCoords = [ProdCoords;SubCoords];
    0083         Director{m1} = [ones(size(ProdCoords,1),1);-1.*ones(size(SubCoords,1),1)];
    0084     end
    0085     MixCoords(:,1) = MixCoords(:,1) - min_x - 18;
    0086     MixCoords(:,2) = MixCoords(:,2) - min_y;
    0087     X = (0); Y =(0); Cntr =(0);
    0088     for cnt = 1:size(MixCoords,1)
    0089          DistancesToCpd = hypot(Xin - MixCoords(cnt,1), Yin - MixCoords(cnt,2));
    0090          ind3 = find(DistancesToCpd >= 0.7*MixCoords(cnt,4)-1 & DistancesToCpd <= MixCoords(cnt,4));
    0091          if isempty(ind3)
    0092              [~,ind3] = min(abs(DistancesToCpd));
    0093          end
    0094          InrangeX = Xin(ind3); InrangeY = Yin(ind3);
    0095          DistancesToCpd1 = hypot(InrangeX - MixCoords(cnt,1), InrangeY - MixCoords(cnt,2));
    0096          [~,ind4] = min(abs(DistancesToCpd1));
    0097           X(cnt,2) = InrangeX(ind4); Y(cnt,2) = InrangeY(ind4);
    0098          ind5 = find(DistancesToCpd >= MixCoords(cnt,4)+7);
    0099          OutrangeX = Xin(ind5); OutrangeY = Yin(ind5);
    0100          DistancesToPint1_temp = hypot(OutrangeX - X(cnt,2), OutrangeY - Y(cnt,2));
    0101          [~,ind6] = min(abs(DistancesToPint1_temp));
    0102          X(cnt,1) = OutrangeX(ind6); Y(cnt,1) = OutrangeY(ind6);
    0103          Cntr(cnt,1) = MixCoords(cnt,1); Cntr(cnt,2) = MixCoords(cnt,2);
    0104     end
    0105     ArX{m1} = X; ArY{m1} = Y; Cntrs{m1} = Cntr;
    0106 end
    0107 if ~isempty(flx_idx)
    0108     [ArX,ArY,Cntrs,Director] = PostMod(ArX,ArY,Cntrs,Director,...
    0109         idx,RefOverlapData,flx_idx);
    0110 end
    0111 % Subfunction =============================================================
    0112 function AllCoords = getBiGGrxns(idx,basemap)
    0113 % Read BiGG model ---------------------------------------------------------
    0114 BiGGModel = getappdata(0,'BiGGModel'); Mdl = getappdata(0,'Mdl');
    0115 if isempty(BiGGModel) || isempty(Mdl)
    0116     [pname,fname1]=uigetfile({'*.*'},'Select the model file (SBML)');
    0117     Testname=pname(strfind(pname,'.')+1:end);
    0118     if ~strcmp(Testname,'xml')
    0119         msgbox('Only xml file format!','Error','error');
    0120         return
    0121     else
    0122         Idf = fopen([fname1,pname]);
    0123         D = textscan(Idf,'%s');
    0124         fclose(Idf);
    0125         CbModel = readCbModel([fname1,pname]);
    0126         [Metkegg,Metbigg] = KEGG2BiGGMets(D,CbModel);
    0127         BiGGModel.Metkegg = Metkegg;
    0128         BiGGModel.CbModel = CbModel;
    0129         BiGGModel.Metbigg = Metbigg;
    0130         setappdata(0,'BiGGModel',BiGGModel);
    0131     end
    0132     Chk = which([pname(1:strfind(pname,'.')-1),'KEGG.mat']);
    0133     if isempty(Chk)
    0134         msgbox({[pname(1:strfind(pname,'.')-1),'KEGG.mat',' cannot be found!'];...
    0135             'Make sure the file is present in BiGG2KEGG folder'},'Error','error');
    0136         return
    0137     end
    0138     Mdl = load(which([pname(1:strfind(pname,'.')-1),'KEGG.mat']));
    0139     setappdata(0,'Mdl',Mdl);
    0140 end
    0141
    0142 Metbigg = BiGGModel.Metbigg;
    0143 Metkegg = BiGGModel.Metkegg;
    0144 CbModel = BiGGModel.CbModel;
    0145 % -------------------------------------------------------------------------
    0146 modelmets = CbModel.mets;
    0147 modelmets = strrep(modelmets,'[','_'); modelmets = strrep(modelmets,']','');
    0148 modelmets = regexp(modelmets,'\w*(?=_\w{1}\>)','match'); % Remove compartments
    0149 modelmets = [modelmets{:}];
    0150 modelrxns = CbModel.rxns; S = CbModel.S; K = Mdl.B2Kegg.K; B = Mdl.B2Kegg.B;
    0151 Rev=CbModel.rev;
    0152
    0153 for m1 = 1:numel(idx)
    0154     Isrxn = B(ismember(K,idx{m1}(1:6)));
    0155     if ~isempty(Isrxn)% If rxn has a corresponding BiGG ID
    0156         IsMdlRxn = find(ismember(modelrxns,Isrxn(1)));
    0157         RevState = Rev(IsMdlRxn);
    0158         if ~isempty(IsMdlRxn) % If corresponding BiGG rxn exists in the model
    0159             [~,in1] = intersect(basemap.rc.rxnid,str2double(idx{m1}(7:end)));
    0160             if isempty(in1)
    0161                 AllCoords{m1} = 'None';
    0162                 continue
    0163             end
    0164             MetIdx = find(S(:,IsMdlRxn));
    0165             Metsign = S(MetIdx,IsMdlRxn); Metsign(Metsign>0)=1;Metsign(Metsign<0)=-1;
    0166             ctr = 1; foundMets = ({}); MetsignEx = (0);
    0167             for m2 = 1:numel(MetIdx)
    0168                 Whrmt = find(ismember(Metbigg,modelmets(MetIdx(m2))));
    0169                 for m3 = 1:numel(Whrmt)
    0170                     foundMets{ctr} = Metkegg{Whrmt(m3)};
    0171                     MetsignEx(ctr) = Metsign(m2);
    0172                     ctr = ctr + 1;
    0173                 end
    0174             end
    0175             rxnSubs = foundMets(MetsignEx<0);
    0176             rxnProds = foundMets(MetsignEx>0);
    0177             ProdNames = basemap.rc.prod{in1};
    0178             SubNames = basemap.rc.sub{in1};
    0179             if RevState %Reversible rxn
    0180                 if any(ismember(rxnProds,SubNames)) || any(ismember(rxnSubs,ProdNames))
    0181                     metCoords = getCoords(basemap,in1,'Rev_back');
    0182                 else
    0183                     metCoords = getCoords(basemap,in1,'Rev_forward');
    0184                 end
    0185             else % Irreversible rxn
    0186                 if any(ismember(rxnProds,SubNames)) || any(ismember(rxnSubs,ProdNames))
    0187                     metCoords = getCoords(basemap,in1,'Irr_back');
    0188                 elseif any(ismember(rxnSubs,SubNames)) || any(ismember(rxnProds,ProdNames))
    0189                     metCoords = getCoords(basemap,in1,'Irr_forward');
    0190                 else % All KEGG rxns are reversible!
    0191                     metCoords = getCoords(basemap,in1,'Rev_forward');
    0192                 end
    0193             end
    0194             AllCoords{m1}.P = metCoords.P; AllCoords{m1}.S = metCoords.S;
    0195         end
    0196     else
    0197         AllCoords{m1} = [];
    0198     end
    0199 end
    0200 % -------------------------------------------------------------------------
    0201 function metCoords = getCoords(basemap,indx,Stat)
    0202 switch Stat
    0203     case 'Rev_forward' %////////////////////////////////
    0204         ProdId = basemap.rc.prodid{indx};
    0205         SubsId = basemap.rc.subid{indx};
    0206         [~,in1] = intersect(basemap.c.cpdid,ProdId);
    0207         [~,in2] = intersect(basemap.c.cpdid,SubsId);
    0208         ProdCoords = basemap.c.cpdxy(in1,:);
    0209         SubCoords = basemap.c.cpdxy(in2,:);
    0210         ProdCoords1 = []; SubCoords1 = [];
    0211         if isempty(in1)
    0212             [~,in1] = intersect(basemap.c.glid,ProdId);
    0213             ProdCoords = basemap.c.glxy(in1,:);
    0214         elseif numel(in1) < numel(ProdId)
    0215             [~,in1] = intersect(basemap.c.glid,ProdId);
    0216             ProdCoords1 = basemap.c.glxy(in1,:);
    0217         end
    0218         if isempty(in2)
    0219             [~,in2] = intersect(basemap.c.glid,SubsId);
    0220             SubCoords = basemap.c.glxy(in2,:);
    0221         elseif numel(in2) < numel(SubsId)
    0222             [~,in2] = intersect(basemap.c.glid,SubsId);
    0223             SubCoords1 = basemap.c.glxy(in2,:);
    0224         end
    0225         ProdCoords = [ProdCoords;ProdCoords1];
    0226         SubCoords = [SubCoords;SubCoords1];
    0227     case 'Rev_back' %////////////////////////////////
    0228         ProdId = basemap.rc.prodid{indx};
    0229         SubsId = basemap.rc.subid{indx};
    0230         [~,in1] = intersect(basemap.c.cpdid,ProdId);
    0231         [~,in2] = intersect(basemap.c.cpdid,SubsId);
    0232         ProdCoords = basemap.c.cpdxy(in1,:);
    0233         SubCoords = basemap.c.cpdxy(in2,:);
    0234         ProdCoords1 = []; SubCoords1 = [];
    0235         if isempty(in1)
    0236             [~,in1] = intersect(basemap.c.glid,ProdId);
    0237             ProdCoords = basemap.c.glxy(in1,:);
    0238         elseif numel(in1) < numel(ProdId)
    0239             [~,in1] = intersect(basemap.c.glid,ProdId);
    0240             ProdCoords1 = basemap.c.glxy(in1,:);
    0241         end
    0242         if isempty(in2)
    0243             [~,in2] = intersect(basemap.c.glid,SubsId);
    0244             SubCoords = basemap.c.glxy(in2,:);
    0245         elseif numel(in2) < numel(SubsId)
    0246             [~,in2] = intersect(basemap.c.glid,SubsId);
    0247             SubCoords1 = basemap.c.glxy(in2,:);
    0248         end
    0249         ProdCoordsR = [ProdCoords;ProdCoords1];
    0250         SubCoordsR = [SubCoords;SubCoords1];
    0251         SubCoords = ProdCoordsR; ProdCoords = SubCoordsR; % Switch coords
    0252     case 'Irr_forward' %////////////////////////////
    0253         ProdId = basemap.rc.prodid{indx};
    0254         [~,in2] = intersect(basemap.c.cpdid,ProdId);
    0255         ProdCoords = basemap.c.cpdxy(in2,:);
    0256         ProdCoords1 = [];
    0257         if isempty(in2)
    0258             [~,in2] = intersect(basemap.c.glid,ProdId);
    0259             ProdCoords = basemap.c.glxy(in2,:);
    0260         elseif numel(in2) < numel(ProdId)
    0261             [~,in2] = intersect(basemap.c.glid,ProdId);
    0262             ProdCoords1 = basemap.c.glxy(in2,:);
    0263         end
    0264         ProdCoords = [ProdCoords;ProdCoords1];
    0265         SubCoords = [];
    0266     case 'Irr_back' %//////////////////////////////
    0267         ProdId = basemap.rc.subid{indx}; % Switch sub/prod
    0268         [~,in2] = intersect(basemap.c.cpdid,ProdId);
    0269         ProdCoords = basemap.c.cpdxy(in2,:);
    0270         ProdCoords1 = [];
    0271         if isempty(in2)
    0272             [~,in2] = intersect(basemap.c.glid,ProdId);
    0273             ProdCoords = basemap.c.glxy(in2,:);
    0274         elseif numel(in2) < numel(ProdId)
    0275             [~,in2] = intersect(basemap.c.glid,ProdId);
    0276             ProdCoords1 = basemap.c.glxy(in2,:);
    0277         end
    0278         ProdCoords = [ProdCoords;ProdCoords1]; SubCoords = [];
    0279 end
    0280 metCoords.P = ProdCoords; metCoords.S = SubCoords;
    0281 %--------------------------------------------------------------------------
    0282 function [ProdCoords,SubCoords] = getCoords4KEGG(basemap,indx)
    0283 ProdId = basemap.rc.prodid{indx};
    0284 SubsId = basemap.rc.subid{indx};
    0285 [~,in1] = intersect(basemap.c.cpdid,ProdId);
    0286 [~,in2] = intersect(basemap.c.cpdid,SubsId);
    0287 ProdCoords = basemap.c.cpdxy(in1,:);
    0288 SubCoords = basemap.c.cpdxy(in2,:);
    0289 ProdCoords1 = []; SubCoords1 = [];
    0290 if isempty(in1)
    0291     [~,in1] = intersect(basemap.c.glid,ProdId);
    0292     ProdCoords = basemap.c.glxy(in1,:);
    0293 elseif numel(in1) < numel(ProdId)
    0294     [~,in1] = intersect(basemap.c.glid,ProdId);
    0295     ProdCoords1 = basemap.c.glxy(in1,:);
    0296 end
    0297 if isempty(in2)
    0298     [~,in2] = intersect(basemap.c.glid,SubsId);
    0299     SubCoords = basemap.c.glxy(in2,:);
    0300 elseif numel(in2) < numel(SubsId)
    0301     [~,in2] = intersect(basemap.c.glid,SubsId);
    0302     SubCoords1 = basemap.c.glxy(in2,:);
    0303 end
    0304 ProdCoords = [ProdCoords;ProdCoords1];
    0305 SubCoords = [SubCoords;SubCoords1];
    0306 %--------------------------------------------------------------------------
    0307 function [ArX,ArY,Cntrs,Director] = PostMod(ArX,ArY,Cntrs,Director,...
    0308     idx,RefOverlapData,flx_idx)
    0309
    0310 [~,Idx1,~] = unique(RefOverlapData);
    0311 Chck1 = setdiff(1:size(RefOverlapData,1),Idx1);
    0312 for ctr = 1:numel(Chck1)
    0313      WhrC = find(ismember(RefOverlapData,RefOverlapData{Chck1(ctr)}));
    0314      if any(ismember(idx(WhrC),flx_idx)) % Overlap with flux carrying rxns
    0315          WhchrxnID = find(ismember(idx(WhrC),flx_idx));
    0316          WhchX = ArX(WhrC(WhchrxnID)); sizeX = (0);
    0317          WhchY = ArY(WhrC(WhchrxnID));
    0318          WhchC = Cntrs(WhrC(WhchrxnID));
    0319          WhchD = Director(WhrC(WhchrxnID));
    0320          for i1 = 1:numel(WhchX)
    0321              sizeX(i1) = size(WhchX{i1},1);
    0322          end
    0323          [~,maxId] = max(sizeX);
    0324          for i1 = 1:numel(WhrC)
    0325              ArX{WhrC(i1)} = WhchX{maxId};
    0326              ArY{WhrC(i1)} = WhchY{maxId};
    0327              Cntrs{WhrC(i1)} = WhchC{maxId};
    0328              Director{WhrC(i1)} = WhchD{maxId};
    0329          end
    0330      end
    0331 end