Description of ArrowheadProcess

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
ArrowheadProcess

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
