Description of NetDraw

Home > BiKEGG > NetDraw.m

# NetDraw

## PURPOSE

NetDraw

## SYNOPSIS

function NetDraw (Outflx,RxnCds,MapChoice,InOpts)

## DESCRIPTION

     NetDraw
     uses KGML file of global metabolic pathway of KEGG (map01100)
     to extract KEGG rxn IDs and their corresponding coordinates.
     The overall structure of NetDraw is similar to that of KeggDraw
     with minor modifications on subfunctions. Further details can be found in
     the BiKEGG user manual

     Inputs:
     Outflx = A mXn matrix containing input flux data, which m corresponds to
            KEGG reaction IDs and n corresponds to time-series values.
     RxnCds = A cell array of strings comprising of KEGG reaction IDs (m rows).
     MapChoice = KEGG map ID.

     Optional inputs:
     InOpts = A struct containing the following fields:
     SVal = Check if user wants to save images : 1 (yes) or 0 (no), default:0.
     Timevalue = The time step between two consecutive images (default: 0.1 s).
     Netstat = Cehck if user wants to use KEGG maps and KGML files online(1)
               /offline(0) (default = 1). Data in KEGGmpas folder will be
               used in case of offline mode.
     Svid = Check if user wants to convert a sequence of images to an animated
             GIF file: 2 (GIF), 0 (no); default: 0 :: Video option is not valid
             for NetDraw

     O. Jamialahmadi
     TMU, Chem. Eng. Dept., Biotech. Group
     June. 2016

## CROSS-REFERENCE INFORMATION

This function calls:

- ArrowheadProcess ArrowheadProcess
- BaseMapReader BaseMapReader
- Img2Animation Img2Animation
- KeggDrawTable KeggDrawTable(GUI)
- MapAdjuster MapAdjuster(GUI)
- MapTrimmer MapTrimmer
- MapTrimmerF MapTrimmerF
- OverlapChecker OverlapChecker
- SaveFlxImgs SavFlxImgs
- SaveImgs SaveImgs(GUI)
- arrowhead arrowhead
- pixFill pixFill
- pixHover pixHover

This function is called by:

## SUBFUNCTIONS

- function [flxOut,rxnCdsI] = ChckRptRxns (flx,rxnCds,CurrentMap,BiGGI)
- function OutSel = StubSel(InSel)
- function [OutFlx,OutRxn,OutBigg] = ChckMultiRxns(BiGG4KeggDraw,RxnCds,rxns4thismap)
- function [rxnOut,flxOut] = ConsistRxns(BiGG4KeggDraw,rxnIn,flxIn,MapChoice)

## SOURCE CODE

    0001 function NetDraw (Outflx,RxnCds,MapChoice,InOpts)
    0002 % NetDraw
    0003 % uses KGML file of global metabolic pathway of KEGG (map01100)
    0004 % to extract KEGG rxn IDs and their corresponding coordinates.
    0005 % The overall structure of NetDraw is similar to that of KeggDraw
    0006 % with minor modifications on subfunctions. Further details can be found in
    0007 % the BiKEGG user manual
    0008 %
    0009 % Inputs:
    0010 % Outflx = A mXn matrix containing input flux data, which m corresponds to
    0011 %        KEGG reaction IDs and n corresponds to time-series values.
    0012 % RxnCds = A cell array of strings comprising of KEGG reaction IDs (m rows).
    0013 % MapChoice = KEGG map ID.
    0014 %
    0015 % Optional inputs:
    0016 % InOpts = A struct containing the following fields:
    0017 % SVal = Check if user wants to save images : 1 (yes) or 0 (no), default:0.
    0018 % Timevalue = The time step between two consecutive images (default: 0.1 s).
    0019 % Netstat = Cehck if user wants to use KEGG maps and KGML files online(1)
    0020 %           /offline(0) (default = 1). Data in KEGGmpas folder will be
    0021 %           used in case of offline mode.
    0022 % Svid = Check if user wants to convert a sequence of images to an animated
    0023 %         GIF file: 2 (GIF), 0 (no); default: 0 :: Video option is not valid
    0024 %         for NetDraw
    0025 %
    0026 % O. Jamialahmadi
    0027 % TMU, Chem. Eng. Dept., Biotech. Group
    0028 % June. 2016
    0029
    0030 % Optional inputs ---------------------------------------------------------
    0031 if exist('InOpts','var')
    0032     OptFields = fieldnames(InOpts);
    0033     for count = 1:numel(OptFields)
    0034         if strcmpi('SVal',OptFields{count})
    0035             SVal = InOpts.(OptFields{count});
    0036         elseif strcmpi('Timevalue',OptFields{count})
    0037             Timevalue = InOpts.(OptFields{count});
    0038         elseif strcmpi('Netstat',OptFields{count})
    0039             Netstat = InOpts.(OptFields{count});
    0040         elseif strcmpi('Svid',OptFields{count})
    0041             Svid = InOpts.(OptFields{count});
    0042         end
    0043     end
    0044 end
    0045 % -------------------------------------------------------------------------
    0046 % Check internet connection
    0047 if ~exist('Netstat','var') || isempty(Netstat)
    0048     Netstat = 1;
    0049 end
    0050 TestLink='http://rest.kegg.jp';
    0051 if Netstat
    0052     [~,Stat1]=urlread(TestLink);
    0053     if ~Stat1
    0054         msgbox({'This step requires a stable internet connection';...
    0055             'Seemingly you are offline!'});
    0056         return
    0057     end
    0058 end
    0059 Pth1 = which ('Bigg2Kegg.m');
    0060 tind = find(Pth1=='\',1,'last');
    0061 Pth = Pth1(1:tind-1);
    0062 %==========================================================================
    0063 BiGG4KeggDraw = getappdata(0,'BiGG4KeggDraw'); % For further modification
    0064 % Read basemap data ==================
    0065 basemap = BaseMapReader(Netstat);
    0066 % Read rxn2map ========================
    0067 Fileid1 = fopen('rxn2map.txt','r');
    0068 rxn2map = textscan(Fileid1,'%s %s');
    0069 RawRxns1 = rxn2map{1};
    0070 RawMaps = rxn2map{2};
    0071 [Loci1,~] = regexp(RawMaps, 'path:rn');
    0072 Loci2 = ~cellfun('isempty', Loci1);
    0073 RawMaps(Loci2) = [];
    0074 RawRxns1(Loci2) = [];
    0075 Maps = strrep(RawMaps,'path:map','');
    0076 Rxns1 = strrep(RawRxns1,'rn:','');
    0077 fclose(Fileid1);
    0078 AllRN = ({});
    0079 AllRNids = [];
    0080 AllRNcol = ({});
    0081 Allflx = [];
    0082 AllrxnCds = ({});
    0083 AllBiGGI = ({});
    0084 for mpct = 1:numel(MapChoice)
    0085     rxns4thismap = Rxns1(ismember(Maps,MapChoice{mpct})); % rxns in the current map
    0086     wherxn = find(ismember(basemap.r.rxn,rxns4thismap));
    0087     RN = basemap.r.rxn(wherxn); % rxns of current map in basemap
    0088     RNids = basemap.r.rxnid(wherxn);
    0089     RNcol = basemap.r.col(wherxn);
    0090     %==========================================================================
    0091     % Analyzing flux data ====================================================
    0092     % RxnCds should be compared to RN, in order to identify which reactions
    0093     % from the Outflx data are present in current KEGG map
    0094     if isempty(Outflx)
    0095         msgbox('Flux data must be read first!','Error','error');
    0096         return
    0097     end
    0098
    0099     % Cehck rxnCodes with HMR IDs and replace them with KEGG ones==========
    0100     Hind = strfind(RxnCds,'HMR');
    0101     if sum([Hind{:}])
    0102         fname = 'HMR2KEGG.mat';
    0103         if ~exist(fname,'file')
    0104              waitfor(msgbox(['The ',fname,' Does not',...
    0105                  ' exist in local directory.',...
    0106                  'Open it manually'],'Warning','warn'));
    0107              [pname,fname1]=uigetfile({'*.mat';'*.xlsx'},...
    0108                  ['Select the ',fname,' file']);
    0109              D=importdata([fname1,pname]);
    0110         else
    0111             D= importdata(fname);
    0112         end
    0113         KEGG=D.KEGG;HMR=D.HMR;
    0114         [Fr,Sn]=ismember(RxnCds,HMR);
    0115         for ind=1:length(Fr)
    0116             if Fr(ind)
    0117                 RxnCds{ind}=KEGG{Sn(ind)};
    0118             end
    0119         end
    0120     end
    0121     %======================================================================
    0122     [MultiFlx,MultiRxn,MultiBigg] = ChckMultiRxns(BiGG4KeggDraw,RxnCds,rxns4thismap);
    0123     MultiRxn1 = []; MultiFlx1 = []; MultiRN = []; MultiRNids = []; MultiRNcol = [];
    0124     MultiBigg1 = [];
    0125     for m1 = 1:numel(MultiRxn)
    0126         where_multi = find(ismember(basemap.r.rxn,MultiRxn{m1}));
    0127         if ~isempty(where_multi)
    0128             MultiRxn1 = [MultiRxn1;MultiRxn(m1)];
    0129             MultiFlx1 = [MultiFlx1;MultiFlx(m1,:)];
    0130             MultiRN = [MultiRN,basemap.r.rxn(where_multi)];
    0131             MultiRNids = [MultiRNids,basemap.r.rxnid(where_multi)];
    0132             MultiRNcol = [MultiRNcol,basemap.r.col(where_multi)];
    0133             MultiBigg1 = [MultiBigg1;MultiBigg(m1)];
    0134         end
    0135     end
    0136     MultiRxn = MultiRxn1;
    0137     MultiFlx = MultiFlx1;
    0138     MultiBigg = MultiBigg1;
    0139     if size(RN,2) > size(RN,1)
    0140         RN = RN';
    0141     end
    0142     if size(MultiRN,2) > size(MultiRN,1)
    0143         MultiRN = MultiRN';
    0144     end
    0145     if size(RNids,2) > size(RNids,1)
    0146         RNids = RNids';
    0147     end
    0148     if size(MultiRNids,2) > size(MultiRNids,1)
    0149         MultiRNids = MultiRNids';
    0150     end
    0151     if size(RNcol,2) > size(RNcol,1)
    0152         RNcol = RNcol';
    0153     end
    0154     if size(MultiRNcol,2) > size(MultiRNcol,1)
    0155         MultiRNcol = MultiRNcol';
    0156     end
    0157     if size(MultiRxn,2) > size(MultiRxn,1)
    0158         MultiRxn = MultiRxn';
    0159     end
    0160     if size(MultiBigg,2) > size(MultiBigg,1)
    0161         MultiBigg = MultiBigg';
    0162     end
    0163     AllrxnCds = [AllrxnCds;MultiRxn];
    0164     AllRN = [AllRN;RN;MultiRN];
    0165     AllRNids = [AllRNids;RNids;MultiRNids];
    0166     Allflx = [Allflx;MultiFlx];
    0167     AllBiGGI = [AllBiGGI;MultiBigg];
    0168     AllRNcol = [AllRNcol;RNcol;MultiRNcol];
    0169
    0170 end
    0171 [rxnCds,flx1] = ConsistRxns(BiGG4KeggDraw,RxnCds,Outflx,MapChoice);
    0172 if size(rxnCds,2) > size(rxnCds,1)
    0173     rxnCds = rxnCds';
    0174 end
    0175 [f1,~]=ismember(rxnCds,AllRN);
    0176 flx=flx1(f1,:);
    0177 rxnCds = rxnCds(f1);
    0178 BiGG4KeggDraw1 = getappdata(0,'BiGG4KeggDraw1');
    0179 BiGGI = BiGG4KeggDraw1(f1);
    0180 if size(BiGGI,2) > size(BiGGI,1)
    0181     BiGGI = BiGGI';
    0182 end
    0183 AllrxnCds = [AllrxnCds;rxnCds];
    0184 Allflx = [Allflx;flx];
    0185 AllBiGGI = [AllBiGGI;BiGGI];
    0186
    0187 MixBiggRxn = strcat(AllrxnCds,AllBiGGI);
    0188 [~, nflx] = unique(MixBiggRxn);
    0189 Allflx = Allflx(nflx,:);
    0190 AllrxnCds = AllrxnCds(nflx); AllBiGGI = AllBiGGI(nflx);
    0191 % Check repeated rxns on this map and let user choose one,among several
    0192 % BiGG IDs, for which the flux value will be displayed on the map.
    0193 [flx,rxnCds] = ChckRptRxns(Allflx,AllrxnCds,MapChoice,AllBiGGI);
    0194
    0195 [rxnCds,nflx1] = unique(rxnCds);
    0196 flx = flx(nflx1,:);
    0197 InrxnCds = rxnCds; Inflx = flx;
    0198 %----------------------------------------------------------------------
    0199 if isempty(flx)
    0200     uiwait(msgbox('The flux data are not is current pathway!',...
    0201         'Warning','warn'));
    0202     close (h);
    0203     return
    0204 end
    0205 %======================================================================
    0206 load(which('Pixdata.mat'))
    0207 fnames = fieldnames(pix_x);
    0208 copt2 = 1;
    0209 rxnCds1={0};rxnPlusid = ({}); rxnCol = ({}); these_ids1 = (0);
    0210 for copt = 1:numel(AllRN)
    0211     Temprxn = find(ismember(AllRN,AllRN{copt}));
    0212     for copt1 = 1:numel(Temprxn)
    0213         these_ids1(copt2) = AllRNids(Temprxn(copt1));
    0214         rxnPlusid{copt2} = [AllRN{Temprxn(copt1)},num2str(these_ids1(copt2))];
    0215         rxnCol{copt2} = AllRNcol{Temprxn(copt1)}; % Identify rxns based on their color to remove rxns in other pathways
    0216         rxnCds1{copt2} = AllRN{copt};
    0217         copt2 = copt2 + 1;
    0218     end
    0219 end
    0220 [rxnPlusid,plusn] = unique(rxnPlusid);
    0221 rxnCds = rxnCds1(plusn);
    0222 rxnCol = rxnCol(plusn);
    0223 these_ids1 = these_ids1(plusn);
    0224 rcid = basemap.rc.rxnid;
    0225 rc_subid = basemap.rc.subid;
    0226 rc_proid = basemap.rc.prodid;
    0227 [shared_ids,~,rcid_n] = intersect(these_ids1,rcid);
    0228 S_rxn = ({}); s3 = 1; S_mat = (0); S_rxnp = ({}); S_cpdcum = []; S_color = ({});
    0229 for s11 = 1:numel(shared_ids)
    0230     S_sub_temp = rc_subid{rcid_n(s11)};
    0231     S_pro_temp = rc_proid{rcid_n(s11)};
    0232     S_cpdcum = [S_cpdcum,S_sub_temp,S_pro_temp];
    0233 end
    0234 S_cpdcum = unique(S_cpdcum); % Unique all compounds for these rxns
    0235
    0236 % Generate unsigned-stoichiometric matrix and subsequent rxn-rxn graph
    0237 for s1 = 1:numel(shared_ids)
    0238     where_rxn = find(these_ids1==shared_ids(s1));
    0239     S_sub_temp = rc_subid{rcid_n(s1)};
    0240     S_pro_temp = rc_proid{rcid_n(s1)};
    0241     for s2 = 1:numel(where_rxn)
    0242         S_rxnp{s3} = rxnPlusid{where_rxn(s2)};
    0243         S_rxn{s3} = rxnCds{where_rxn(s2)}; % Rxns for
    0244         S_color{s3} = rxnCol{where_rxn(s2)};
    0245         [~,S_sub_temp_idx,~] = intersect(S_cpdcum,S_sub_temp);
    0246         [~,S_pro_temp_idx,~] = intersect(S_cpdcum,S_pro_temp);
    0247         if size(S_sub_temp_idx,1) > size(S_sub_temp_idx,2)
    0248             S_sub_temp_idx = S_sub_temp_idx';
    0249         end
    0250         if size(S_pro_temp_idx,1) > size(S_pro_temp_idx,2)
    0251             S_pro_temp_idx = S_pro_temp_idx';
    0252         end
    0253         S_mat([S_sub_temp_idx,S_pro_temp_idx],s3) = 1;
    0254         s3 = s3 + 1;
    0255     end
    0256 end
    0257 S_mat = sparse(S_mat);
    0258 % Cpd-cpd graph------------------------------------------------------------
    0259 % S_graph = zeros(size(S_mat,1)); t_graph = ({});
    0260 % for i1=1:size(S_mat,1)
    0261 %     nr = find(S_mat(i1,:));
    0262 %     if ~isempty(nr)
    0263 %         for i2 = 1:numel(nr)
    0264 %             t_graph{i2} = find(S_mat(:,nr(i2)));
    0265 %             if size(t_graph{i2},1) > size(t_graph{i2},2)
    0266 %                 t_graph{i2} = t_graph{i2}';
    0267 %             end
    0268 %         end
    0269 %         t_con=[t_graph{:}]; t_con=unique(t_con);
    0270 %         clear t_graph
    0271 %         S_graph(i1,setdiff(t_con,i1)) = 1;
    0272 %     end
    0273 % end
    0274 % Rxn-rxn graph -----------------------------------------------------------
    0275 S_graph = zeros(size(S_mat,2)); t_graph = ({});
    0276 for i1=1:size(S_mat,2)
    0277     nr = find(S_mat(:,i1));
    0278     if ~isempty(nr)
    0279         for i2 = 1:numel(nr)
    0280             t_graph{i2} = find(S_mat(nr(i2),:));
    0281             if size(t_graph{i2},1) > size(t_graph{i2},2)
    0282                 t_graph{i2} = t_graph{i2}';
    0283             end
    0284         end
    0285         t_con=[t_graph{:}]; t_con=unique(t_con);
    0286         clear t_graph
    0287         S_graph(i1,setdiff(t_con,i1)) = 1;
    0288     end
    0289 end
    0290 b_graph = biograph(S_graph);
    0291 [S_bins, where_bins] = conncomp(b_graph); % Connected parts of rxn-rxn graph
    0292 freq_bins = histc(where_bins,1:S_bins);
    0293 [~,max_freq_bin] = max(freq_bins);
    0294 Secplace_freq_bin = find(freq_bins>2);
    0295  % All cpds present in the strongest (most connected) part of cpd-cpd
    0296  % graph:
    0297 % core_cpds = S_cpdcum(where_bins == max_freq_bin);
    0298 core_rxns= S_rxnp(where_bins == max_freq_bin); % All cpds present in the strongest (most connected) part of rxn-rxn graph
    0299 Secids = [];
    0300 for sctr = 1:numel(Secplace_freq_bin)
    0301     Secids = [Secids,find(where_bins == Secplace_freq_bin(sctr))];
    0302 end
    0303 Secondcore_rxns= S_rxnp(Secids);
    0304 [~,orph_id] = intersect(where_bins,find(freq_bins == 1));
    0305 % [~,double_id] = intersect(where_bins,find(freq_bins == 2));
    0306 % double_rxns = S_rxnp(double_id);
    0307 orphan_rxns = S_rxnp(orph_id); % Single rxns which are not connected to any part: These rxns should be removed.
    0308 % Identify repeated rxns on the final map and remove redundancies based on
    0309 % core_cpds (from graph connectivity assessment): Among
    0310 % different version of same rxn (same rxn name and different rxn id), only
    0311 % those with high connectedness are retained.
    0312 % Dump color check ////////////////////////////////////////////////////////
    0313 Prominent_col = S_color(where_bins == max_freq_bin);
    0314 Prominent_col = unique(Prominent_col);
    0315 %//////////////////////////////////////////////////////////////////////////
    0316 [~,S_rxn_idx,~] = unique(S_rxn);
    0317 S_rxn_rptd = S_rxn(setdiff(1:numel(S_rxn),S_rxn_idx));
    0318 rmv_rxnPlus_id = []; plus_id = 1;
    0319 for c1 = 1:numel(S_rxn_rptd)
    0320     where_rptd = find(ismember(S_rxn,S_rxn_rptd{c1}));
    0321     rmv_mat = []; rmv_id = 1;
    0322     for c2 = 1:numel(where_rptd)
    0323 %         cpds4rptd = find(S_mat(:,where_rptd(c2))); % Only for cpd-cpd graph
    0324         if ~isempty(setdiff(S_rxnp(where_rptd(c2)),core_rxns)) && ...
    0325               ~isempty(setdiff(S_color(where_rptd(c2)),Prominent_col))
    0326             rmv_mat(rmv_id) = where_rptd(c2); % Those who failed the connectivity test
    0327             rmv_id = rmv_id + 1;
    0328         end
    0329     end
    0330     for c3 = 1:numel(rmv_mat)
    0331         rmv_rxnPlus_id(plus_id) = rmv_mat(c3);
    0332         plus_id = plus_id + 1;
    0333     end
    0334 end
    0335 % Final check for orphan rxns after removal of rmv_rxnPlus_id ids
    0336 S_graph1 = S_graph;
    0337 S_graph1(rmv_rxnPlus_id,:) = [];
    0338 S_graph1(:,rmv_rxnPlus_id) = [];
    0339 b_graph1 = biograph(S_graph1);
    0340 [S_bins1, where_bins1] = conncomp(b_graph1);
    0341 freq_bins1 = histc(where_bins1,1:S_bins1);
    0342 [~,orph_id1] = intersect(where_bins1,find(freq_bins1 == 1));
    0343 S_rxnp1 = S_rxnp;
    0344 S_rxnp1(rmv_rxnPlus_id) = [];
    0345 orphan_rxns1 = S_rxnp1(orph_id1);
    0346 %--------------------------------------------------------------------------
    0347 % Check for orphan multiple rxns or their children : These rxns do not have
    0348 % any cpd on basemap, so, do not appear in rxn-rxn graph and therefore,
    0349 % cannot be identified. To identify these rxns, cpd2rxn is used and through
    0350 % comparison with cpds in core rxns, its orphanage state can be stablished.
    0351 multirxns_id = setdiff(these_ids1,shared_ids);
    0352 [~, multirxns_id] = intersect(these_ids1,multirxns_id);
    0353 multirxns = rxnCds(multirxns_id);
    0354 % Read cpd2rxn
    0355 Fileid = fopen('cpd2rxn.txt','r');
    0356 cpd2rxn = textscan(Fileid,'%s %s');
    0357 rxnData = strrep(cpd2rxn{2},'rn:','');
    0358 cpds = strrep(cpd2rxn{1},'cpd:','');
    0359 fclose(Fileid);
    0360 % Identify cpds of rxns with good connectedness (>2)
    0361 SecCore_ids = these_ids1(ismember(S_rxnp,Secondcore_rxns));
    0362 [~,SecCore_idx] = intersect(basemap.rc.rxnid,SecCore_ids);
    0363 core2_subs = basemap.rc.sub(SecCore_idx); core2_subs = [core2_subs{:}];
    0364 core2_pros = basemap.rc.prod(SecCore_idx); core2_pros = [core2_pros{:}];
    0365 core2_cpds = [core2_subs,core2_pros]; core2_cpds = unique(core2_cpds);
    0366 m2 = 1; rmv_mult_id = [];
    0367 for m1 = 1:numel(multirxns)
    0368     curr_cpds = cpds(ismember(rxnData,multirxns(m1)));
    0369     if isempty(intersect(curr_cpds,core2_cpds)) || ...
    0370             (~any(strcmp(MapChoice,'00640')) && strcmp(multirxns(m1),'R00928'))
    0371         rmv_mult_id(m2) = multirxns_id(m1);
    0372         m2 = m2 + 1;
    0373     end
    0374 end
    0375 rxnCds(rmv_mult_id) = [];
    0376 rxnPlusid(rmv_mult_id) = [];
    0377 these_ids1(rmv_mult_id) = [];
    0378 %--------------------------------------------------------------------------
    0379 % Modify final rxns and remove orphans and redundants
    0380 rmv_rxnPlus = S_rxnp(rmv_rxnPlus_id); % Rxns to be removed from final map
    0381 rxnCds(ismember(rxnPlusid,rmv_rxnPlus)) = [];
    0382 these_ids1(ismember(rxnPlusid,rmv_rxnPlus)) = [];
    0383 rxnPlusid(ismember(rxnPlusid,rmv_rxnPlus)) = [];
    0384 rxnCds(ismember(rxnPlusid,orphan_rxns)) = [];
    0385 these_ids1(ismember(rxnPlusid,orphan_rxns)) = [];
    0386 rxnPlusid(ismember(rxnPlusid,orphan_rxns)) = []; % Remove orphan rxns
    0387 rxnCds(ismember(rxnPlusid,orphan_rxns1)) = [];
    0388 these_ids1(ismember(rxnPlusid,orphan_rxns1)) = [];
    0389 rxnPlusid(ismember(rxnPlusid,orphan_rxns1)) = []; % Remove orphan rxns after removal of redundant rxns
    0390 %Flux data arrangement-----------------------------------------------------
    0391 [On2,~] = ismember(rxnPlusid,fnames); % Idx of all rxns in these maps present in basemap
    0392 field_intsect = rxnPlusid(On2); % Reference to all rxns in these maps for reading coordinates
    0393 rxnCds = rxnCds(On2); % Rearrange all rxn IDs in these maps
    0394
    0395 overlay_rxnPlus = rxnPlusid(ismember(rxnCds,InrxnCds)); % Idx of rxns for which flux rates are provided
    0396 [~,o2] = ismember(rxnCds,InrxnCds); % Rearrange rxn IDs and fluxes based on Idx
    0397 F_flx = [];
    0398 o2 = o2(find(o2)); OverlayRxnsIDs = InrxnCds(o2);
    0399 for m1 = 1:numel(o2)
    0400     F_flx(m1,:) = Inflx(o2(m1),:);
    0401 end
    0402 [On1,~] = (ismember(overlay_rxnPlus,fnames)); % Find idxs of rxns for which flux rates are provided and present in basemap
    0403 field_intsect_overlay = overlay_rxnPlus(On1); % Reference to all rxns with flux rates
    0404 F_flx = F_flx(On1,:); % All flux rates
    0405 OverlayRxnsIDs = OverlayRxnsIDs(On1); % KEGG rxn IDs corresponding to field_intsect_overlay
    0406 flx = F_flx;
    0407 %--------------------------------------------------------------------------
    0408 LineStrength = 7;
    0409 pixx = []; pixy =[];
    0410 for i4 = 1:numel(field_intsect)
    0411     pixx = [pixx;pix_x.(field_intsect{i4})];
    0412     pixy = [pixy;pix_y.(field_intsect{i4})];
    0413 end
    0414 min_x = min(pixx); max_x = max(pixx); min_y = min(pixy); max_y = max(pixy);
    0415 pixx = []; pixy = []; Allx = ({}); Ally = ({}); Allx4Arrow = ({});Ally4Arrow = ({});
    0416 RefOverlapData = cell(numel(field_intsect),1);
    0417 I = ones(max_y-min_y+50,max_x-min_x+50,3);
    0418 I(:,:,1) = I(:,:,1).*0.8; I(:,:,2) = I(:,:,2).*0.992; I(:,:,3) = I(:,:,3).*0.8;
    0419
    0420 for i4 = 1:numel(field_intsect)
    0421     pixx = pix_x.(field_intsect{i4})-min_x+10;
    0422     pixy = pix_y.(field_intsect{i4})-min_y+10;
    0423     RefOverlapData{i4} = [num2str(numel(pixx)),',',...
    0424         num2str(mean(pixx)),',',num2str(numel(pixy)),',',num2str(mean(pixy))];
    0425     Allx4Arrow{i4} = pixx; Ally4Arrow{i4} = pixy;
    0426     [I,pixx_temp,pixy_temp,pixx_temp1,pixy_temp1] = pixFill(I,pixx,pixy,LineStrength,[0 0 0]);
    0427     if numel(pixx_temp)>20; pixx_temp(1:10)=[]; pixx_temp(end-10:end)=[];end;
    0428     if numel(pixy_temp)>20; pixy_temp(1:10)=[]; pixy_temp(end-10:end)=[];end;
    0429     if numel(pixx_temp1)>20; pixx_temp1(1:10)=[]; pixx_temp1(end-10:end)=[];end;
    0430     if numel(pixy_temp1)>20; pixy_temp1(1:10)=[]; pixy_temp1(end-10:end)=[];end;
    0431     Allx{i4} = [pixx_temp;pixx_temp1];
    0432     Ally{i4} = [pixy_temp;pixy_temp1];
    0433     clear pixx pixy
    0434 end
    0435
    0436 % Compounds ---------------------------------------------------------------
    0437 [~,n1] = intersect(basemap.rc.rxnid,these_ids1);
    0438 a1 = basemap.rc.subid(n1); a1 = [a1{:}];
    0439 b1 = basemap.rc.prodid(n1); b1 = [b1{:}];
    0440 a=[a1,b1];
    0441 a=unique(a);
    0442 [~,n21]=intersect(basemap.c.cpdid,a);
    0443 [~,n22]=intersect(basemap.c.glid,a);
    0444 cpdnames1 = basemap.c.cpd(n21);
    0445 cpdnames2 = basemap.c.gl(n22);
    0446 cpdnames = [cpdnames1,cpdnames2];
    0447 cpdxy1 = basemap.c.cpdxy(n21,:);
    0448 cpdxy2 = basemap.c.glxy(n22,:);
    0449 cpdxy = [cpdxy1;cpdxy2];
    0450 x_cpd = cpdxy(:,1) - min_x - 18;
    0451 y_cpd = cpdxy(:,2) - min_y;
    0452
    0453 % Flux overlay ============================================================
    0454 % The time step between two consecutive images:
    0455 if ~exist('Timevalue','var') || isempty(Timevalue)
    0456     Timevalue=0.1;
    0457 end
    0458
    0459 % Time points: % Disabled in current settings /////////////////////////////
    0460 % if ~exist('mTime','var') || isempty(mTime)
    0461 %     mTime=1:size(Outflx,2);
    0462 % end
    0463 %
    0464 % if ~exist('TimeUnit','var') || isempty(TimeUnit)
    0465 %     TimeUnit=' ';
    0466 % end
    0467 %//////////////////////////////////////////////////////////////////////////
    0468 if ~exist('Svid','var') || isempty(Svid)
    0469     Svid = 0; %Does not convert images to video
    0470 end
    0471 % Check wether user wants to save images (using 'Save images checkbox')
    0472 if ~exist('SVal','var') || isempty(SVal)
    0473     SVal = 0; % Does not save images.
    0474 end
    0475 if SVal
    0476     FolderNm='NetDraw output images';
    0477     Tn=uigetdir(pwd,['Folder to save images. To save in current folder:',...
    0478     'Cancel']);
    0479     if ~Tn
    0480         if ~exist(FolderNm,'file') % Check if folder already exists
    0481             mkdir(FolderNm);
    0482             FolderPath=fullfile(pwd,FolderNm);
    0483         else
    0484             FolderPath=fullfile(pwd,FolderNm);
    0485         end
    0486     else
    0487         if ~exist(fullfile(Tn,FolderNm),'file') % Check if folder already exists
    0488             mkdir(fullfile(Tn,FolderNm));
    0489             FolderPath=fullfile(Tn,FolderNm);
    0490         else
    0491             FolderPath=fullfile(Tn,FolderNm);
    0492         end
    0493     end
    0494 end
    0495 Cm=jet(numel(field_intsect_overlay)); % [Warning] Changing Cm matrix is not recommended
    0496 if numel(MapChoice) > 4
    0497     TempTitle = MapChoice(1:4);
    0498     TempTitle = strjoin(TempTitle,',');
    0499     TempTitle = [TempTitle,',...'];
    0500 else
    0501     TempTitle = strjoin(MapChoice,',');
    0502 end
    0503 %Draw images on this panel
    0504 H1 = figure;
    0505 h1 = axes;
    0506 set(H1,'position',get(0,'screensize'));
    0507 set(h1,'position',[0,0,1,1],'xtick',[],'ytick',[]);
    0508 set(H1,'toolbar','figure');
    0509 set(H1,'menubar','figure');
    0510 set(H1,'name',TempTitle,'numbertitle','off')
    0511 % Arbitrary color and font properties % Disabled in current settings///////
    0512 % if ~exist('Viso','var') || ~Viso
    0513 %     TextFont.FontName = 'Arial';
    0514 %     TextFont.FontWeight = 'normal';
    0515 %     TextFont.FontAngle = 'normal';
    0516 %     TextFont.FontUnits = 'points';
    0517 %     TextFont.FontSize = 8;
    0518 % end
    0519 % if exist('Viso','var') && Viso
    0520 %     waitfor(VisualProp)
    0521 %     TextFont = getappdata(0,'TextFont');
    0522 %     if isempty(TextFont)
    0523 %         TextFont.FontName = 'Arial';
    0524 %         TextFont.FontWeight = 'normal';
    0525 %         TextFont.FontAngle = 'normal';
    0526 %         TextFont.FontUnits = 'points';
    0527 %         TextFont.FontSize = 8;
    0528 %     end
    0529 % end
    0530 %//////////////////////////////////////////////////////////////////////////
    0531 % Visualization -----------------------------------------------------------
    0532 flx_abs = abs(flx);
    0533 for k=1:size(flx,2)
    0534     ArX = ({}); ArY = ArX; McmCollect = zeros(size(flx,1),3);
    0535     OverlapData = cell(size(flx,1),1);
    0536     for m=1:size(flx,1)
    0537         % Normalize colors and insertion into enzyme boxes
    0538         flxMax=max(abs(flx(:)));
    0539         flxMin=min(abs(flx(:)));
    0540         if flxMax==flxMin % One row (i.e. for 1 rxn)
    0541             flxMax=max(abs(flx(:)));
    0542             flxMin=min(abs(flx(:)));
    0543         end
    0544         flxSl=(size(Cm,1)-1)./(flxMax-flxMin);
    0545         if isinf(flxSl) % in case of flxMax == flxMin
    0546             flxSl = 1e15.*(size(Cm,1)-1);
    0547         end
    0548         Mcm=ceil((flx_abs(m,k)-flxMax).*flxSl+size(Cm,1));
    0549
    0550         % Show the image:
    0551         pixx = pix_x.(field_intsect_overlay{m})-min_x+10;
    0552         pixy = pix_y.(field_intsect_overlay{m})-min_y+10;
    0553         OverlapData{m} = [num2str(numel(pixx)),',',...
    0554             num2str(mean(pixx)),',',num2str(numel(pixy)),',',num2str(mean(pixy))];
    0555         I = pixFill(I,pixx,pixy,LineStrength,[Cm(Mcm,1) Cm(Mcm,2) Cm(Mcm,3)]);
    0556         ArX{m} = pixx; ArY{m} = pixy;
    0557         McmCollect(m,:) = [Cm(Mcm,1) Cm(Mcm,2) Cm(Mcm,3)];
    0558         clear pixx pixy
    0559     end
    0560      %---------------------------------------------------------------------
    0561     [I,ArXo,ArYo,ArrowChek,RxnId] = OverlapChecker(I,ArX,ArY,OverlapData,...
    0562         flx(:,k),McmCollect,field_intsect_overlay,LineStrength);
    0563     [ArrowX,ArrowY,Cntrs] = ArrowheadProcess(Allx4Arrow,Ally4Arrow,...
    0564         min_x,min_y,field_intsect,basemap,RxnId,RefOverlapData,'bigg',field_intsect_overlay);
    0565     I = MapTrimmer(I,ArrowX,ArrowY,Cntrs,Allx4Arrow,Ally4Arrow,LineStrength,[0.8 0.992 0.8]);
    0566     [ArrowX1,ArrowY1,Cntrs1,Director] = ArrowheadProcess(ArXo,...
    0567         ArYo,min_x,min_y,field_intsect_overlay,basemap,[],[],'bigg',[]);
    0568     I = MapTrimmerF(I,ArrowX1,ArrowY1,Cntrs1,ArXo,ArYo,RxnId,...
    0569         field_intsect_overlay,LineStrength,[0.8 0.992 0.8]);
    0570     % Compounds overlay ---------------------------------------------------
    0571     for i1 = 1:size(cpdxy,1)
    0572         cc1 = x_cpd(i1)-0.7*cpdxy(i1,4):x_cpd(i1)+cpdxy(i1,4)*0.7;
    0573         cc2 = y_cpd(i1)-cpdxy(i1,4)*0.7:y_cpd(i1)+cpdxy(i1,4)*0.7;
    0574         [x,y] = meshgrid(cc1,cc2);
    0575         f1= ((x-x_cpd(i1)).^2+(y-y_cpd(i1)).^2) <=(cpdxy(i1,4)*0.7)^2;
    0576         f2= ((x-x_cpd(i1)).^2+(y-y_cpd(i1)).^2) <=(cpdxy(i1,4)*0.55)^2;
    0577         x1 = x(f1); y1= y(f1); x1=round(x1); y1=round(y1);
    0578         x2 = x(f2); y2= y(f2); x2=round(x2); y2=round(y2);
    0579         x1(~y1) = []; y1(~y1) = [];
    0580         x2(~y2) = []; y2(~y2) = [];
    0581         for i2 = 1:numel(x1)
    0582             I(y1(i2),x1(i2),1) = 0.502;
    0583             I(y1(i2),x1(i2),2) = 0.502;
    0584             I(y1(i2),x1(i2),3) = 0;
    0585         end
    0586         for i2 = 1:numel(x2)
    0587             I(y2(i2),x2(i2),1) = 0;
    0588             I(y2(i2),x2(i2),2) = 1;
    0589             I(y2(i2),x2(i2),3) = 0.498;
    0590         end
    0591         clear f1 x1 x2 y1 y2 x y cc1 cc2
    0592     end
    0593     imshow(I,'InitialMagnification', 'fit','Parent',h1);
    0594     % Disabled in current settings ////////////////////////////////////////
    0595     % Show time-series values on the corresponding image
    0596 %     TimePropUDF ('One',TextFont,TimeUnit,mTime,k,h1)
    0597     %//////////////////////////////////////////////////////////////////////
    0598     arrowhead(ArrowX,ArrowY,[],0,4,[]);
    0599     flxDirect.D = Director; flxDirect.F = flx(:,k);
    0600     arrowhead(ArrowX1,ArrowY1,ArrowChek,McmCollect,4,flxDirect);
    0601
    0602     % Process pixhover settings--------------------------------------------
    0603     Mdl = getappdata(0,'Mdl');
    0604     Uni = load(which('UniModelKEGG.mat'));
    0605     CRxns = getappdata(0,'CRxns'); % Map of consistent rxns
    0606     if ~isempty(CRxns) && k == 1
    0607         ConOverlayRxnID = OverlayRxnsIDs;
    0608         for cnter = 1:numel(CRxns.C)
    0609             ConOverlayRxnID(ismember(OverlayRxnsIDs,CRxns.C{cnter})) = ...
    0610                 {[CRxns.C{cnter},'(',CRxns.O{cnter},')']};
    0611             OverlayRxnsIDs(ismember(OverlayRxnsIDs,CRxns.C{cnter})) = ...
    0612                 CRxns.O(cnter);
    0613         end
    0614     elseif isempty(CRxns)
    0615         ConOverlayRxnID = OverlayRxnsIDs;
    0616     end
    0617     K = Uni.B2Kegg.K; K1 = Mdl.B2Kegg.K;
    0618     B = Uni.B2Kegg.B; B1 = Mdl.B2Kegg.B;
    0619     clear Uni
    0620     HoverUs.K = K1; HoverUs.B = B1;
    0621     HoverUs.R = OverlayRxnsIDs; HoverUs.P = field_intsect_overlay';
    0622     HoverUs.CR = ConOverlayRxnID;
    0623     HoverAll.R = rxnCds; HoverAll.P = field_intsect;
    0624     HoverAll.K = K; HoverAll.B = B;
    0625     pixHover(x_cpd,y_cpd,0.7*cpdxy(1,4),cpdnames,Allx,Ally,...
    0626         HoverAll,HoverUs,abs(flx(:,k)),h1,H1)
    0627
    0628     % Show colorbar--------------------------------------------------------
    0629     flxstr=linspace(flxMin,flxMax,7);
    0630     flxstr1=cell(numel(flxstr),1);
    0631     for ir=1:numel(flxstr)
    0632         flxstr1{ir,1}=num2str(flxstr(ir),2);
    0633     end
    0634     ch = colorbar('peer',h1);
    0635     set(ch,'yticklabel',flxstr1)
    0636     % Check for further settings ------------------------------------------
    0637     IsWanting = input('Do you want to post-process the image? (y/n):','s');
    0638     if isempty(IsWanting) || strcmp(IsWanting,'y')
    0639         setappdata(0,'Img',I);
    0640         setappdata(0,'ParentAx',h1);
    0641         setappdata(0,'ParentFig',H1);
    0642         PostData.field_intsect = field_intsect;PostData.Allx=Allx;
    0643         PostData.Ally=Ally;PostData.min_x=min_x;PostData.min_y=min_y;
    0644         PostData.RefOverlapData=RefOverlapData;PostData.Allx4Arrow=Allx4Arrow;
    0645         PostData.Ally4Arrow=Ally4Arrow;PostData.ArX=ArX;PostData.ArY=ArY;
    0646         PostData.OverlapData=OverlapData;PostData.field_intsect_overlay=field_intsect_overlay;
    0647         PostData.ArXo=ArXo;PostData.ArYo=ArYo;PostData.ArrowChek=ArrowChek;
    0648         PostData.RxnId=RxnId;PostData.ArrowX=ArrowX;PostData.ArrowY=ArrowY;
    0649         PostData.Cntrs=Cntrs;PostData.ArrowX1=ArrowX1;PostData.ArrowY1=ArrowY1;
    0650         PostData.Cntrs1=Cntrs1;PostData.McmCollect=McmCollect;
    0651         PostData.cpdnames=cpdnames;PostData.cpdxy=cpdxy;PostData.x_cpd=x_cpd;
    0652         PostData.y_cpd=y_cpd;PostData.flx=flx;PostData.rxnCds=rxnCds;
    0653         PostData.max_x=max_x;PostData.max_y=max_y;PostData.basemap=basemap;
    0654         PostData.OverlayRxnsIDs = OverlayRxnsIDs; PostData.k = k;
    0655         PostData.ConOverlayRxnID = ConOverlayRxnID;
    0656 %         PostData.TimeUnit=TimeUnit;PostData.TextFont=TextFont;
    0657 %         PostData.mTime=mTime;
    0658         setappdata(0,'PostData',PostData);
    0659         uiwait(MapAdjuster);
    0660     end
    0661     % Save to image files within a new folder------------------------------
    0662     if SVal
    0663         if k == 1
    0664             waitfor(SaveImgs)
    0665             ImgRes = getappdata(0,'ImgRes');
    0666             ImgFrmt = getappdata(0,'ImgFrmt');
    0667             if isempty(ImgRes)
    0668                 ImgRes=1; % Corresponds to 300 dpi in SaveFlxImgs
    0669             elseif isempty(ImgFrmt)
    0670                 ImgFrmt=1; % Corresponds to PNG in SaveFlxImgs
    0671             end
    0672         end
    0673         SaveFlxImgs(getappdata(0,'ParentFig'),FolderPath,k,ImgRes,ImgFrmt);
    0674     end
    0675     %----------------------------------------------------------------------
    0676     pause(Timevalue); % Time step between consecutive images
    0677 end
    0678 if Svid
    0679     Img2Animation(Svid)
    0680 end
    0681
    0682 %==========================================================================
    0683 % Subfunctions ------------------------------------------------------------
    0684 %==========================================================================
    0685 function [flxOut,rxnCdsI] = ChckRptRxns (flx,rxnCds,CurrentMap,BiGGI) % NetDraw
    0686 % First, rxns are pruned based on manual assessment of pathways
    0687 rxnCdsI = rxnCds;
    0688
    0689 Fileid1 = fopen('rxnbiggmap.txt','r');
    0690 Rxnbmap = textscan(Fileid1,'%s %s %s %d');
    0691 Rtemp = Rxnbmap{1}; Btemp = Rxnbmap{2}; Mtemp = Rxnbmap{3};
    0692 Dtemp = Rxnbmap{4};
    0693 fclose(Fileid1);
    0694 AllrxnCdsI = ({}); AllBiGGI = ({}); Allflx = [];
    0695 for i1 = 1:numel(CurrentMap)
    0696     if any(ismember(Mtemp,CurrentMap{i1}))
    0697         ChR = Rtemp(ismember(Mtemp,CurrentMap{i1}));
    0698         ChB = Btemp(ismember(Mtemp,CurrentMap{i1}));
    0699         ChD = Dtemp(ismember(Mtemp,CurrentMap{i1}));
    0700         for ct = 1:numel(ChR)
    0701             n1_chr = find(ismember(rxnCdsI,ChR{ct}));
    0702             n1_chb = find(ismember(BiGGI,ChB{ct}));
    0703             if ~isempty(intersect(n1_chr,n1_chb))
    0704                 Rmvit = intersect(n1_chr,n1_chb);
    0705                 for ct1 = 1:numel(Rmvit)
    0706                     if strcmp(BiGGI{Rmvit(ct1)},ChB{ct})
    0707                         Tmpflx = flx(Rmvit(ct1),:);
    0708                     end
    0709                 end
    0710                 rxnCdsI(Rmvit) = [];
    0711                 BiGGI(Rmvit) = [];
    0712                 flx(Rmvit,:) = [];
    0713                 if ChD(ct)
    0714                     CurL1 = numel(rxnCdsI) + 1;
    0715                     rxnCdsI{CurL1} = ChR{ct};
    0716                     BiGGI{CurL1} = ChB{ct};
    0717                     flx(CurL1,:) = Tmpflx;
    0718                 end
    0719             end
    0720         end
    0721     end
    0722     if size(rxnCdsI,2) > size(rxnCdsI,1)
    0723         rxnCdsI = rxnCdsI';
    0724     end
    0725     if size(BiGGI,2) > size(BiGGI,2)
    0726         BiGGI = BiGGI';
    0727     end
    0728     AllrxnCdsI = [AllrxnCdsI;rxnCdsI];
    0729     AllBiGGI = [AllBiGGI;BiGGI];
    0730     Allflx = [Allflx;flx];
    0731 end
    0732 MixBiggRxn = strcat(AllrxnCdsI,AllBiGGI);
    0733 [~, nflx] = unique(MixBiggRxn);
    0734 Allflx = Allflx(nflx,:);
    0735 AllrxnCdsI = AllrxnCdsI(nflx); AllBiGGI = AllBiGGI(nflx);
    0736 rxnCdsI = AllrxnCdsI;
    0737 BiGGI = AllBiGGI;
    0738 flx = Allflx;
    0739 %------
    0740 [R1,~,R3] = unique(rxnCdsI,'stable');
    0741 countR = 1;RptCount=(0);
    0742 while numel(R3)
    0743     RptCount(countR) = numel(find(R3==R3(1)));
    0744     R3(R3==R3(1)) = [];
    0745     countR = countR+1;
    0746 end
    0747 FndRxnLoci = find(RptCount>1);
    0748 RptRxns = {0}; RptBiggs = {0};RptFlx= {0};
    0749 for tctr = 1:numel(FndRxnLoci)
    0750     RptRxnsI = R1(FndRxnLoci(tctr));
    0751     RptRxns{tctr} = rxnCdsI(ismember(rxnCdsI,RptRxnsI));
    0752     RptBiggs{tctr} = BiGGI(ismember(rxnCdsI,RptRxnsI));
    0753     RptFlxLoci = ismember(rxnCdsI,RptRxnsI);
    0754     RptFlx{tctr} = flx(RptFlxLoci,:);
    0755 end
    0756 if any(RptCount>1)
    0757     setappdata(0,'RptBiggs',RptBiggs);
    0758     setappdata(0,'CurrentMap','All maps');
    0759     setappdata(0,'RptRxns',RptRxns);
    0760     setappdata(0,'RptFlx',RptFlx);
    0761     waitfor(KeggDrawTable)
    0762     SlctdIds = getappdata(0,'data1');
    0763     if ~isempty(SlctdIds)
    0764         NtSel = RptBiggs(ismember(SlctdIds,'Not slected'));
    0765         OwnSel = StubSel(NtSel);
    0766         SlctdIds(ismember(SlctdIds,'Not slected')) = OwnSel;
    0767     else
    0768         SlctdIds = StubSel(RptBiggs);
    0769     end
    0770     % Refine flx
    0771     for count = 1:numel(SlctdIds)
    0772         Fi = find(ismember(RptBiggs{count},SlctdIds{count}));
    0773         for count1 = 1:size(RptFlx{count},1)
    0774             RptFlx{count}(count1,:) = RptFlx{count}(Fi,:);
    0775         end
    0776     end
    0777     for count1 = 1:size(RptRxns,2)
    0778         Fi1 = find(ismember(rxnCdsI,RptRxns{count1}));
    0779         flx(Fi1,:) = RptFlx{count1};
    0780     end
    0781     flxOut = flx;
    0782
    0783 else
    0784     flxOut = flx;
    0785 end
    0786 % Remove all appdata
    0787 AT=getappdata(0);
    0788 AT1 = fieldnames(AT);
    0789 for co = 1:length(AT1)
    0790     if ~strcmp(AT1{co},'BiGG4KeggDraw') && ~strcmp(AT1{co},'MultiB')...
    0791             && ~strcmp(AT1{co},'MultiK') && ~strcmp(AT1{co},'MultiFlx') ...
    0792             && ~strcmp(AT1{co},'BiGGModel') && ~strcmp(AT1{co},'Mdl') ...
    0793             && ~strcmp(AT1{co},'CRxns')
    0794         rmappdata(0,AT1{co});
    0795     end
    0796 end
    0797 %--------------------------------------------------------------------------
    0798 function OutSel = StubSel(InSel)
    0799 OutSel = {0};
    0800 for count1 = 1:numel(InSel)
    0801     Tinder = regexp(InSel{count1},'\w*p$'); %Peroxisome
    0802     Tinder = find(~cellfun('isempty', Tinder));
    0803     CurSel = InSel{count1};
    0804     if (numel(CurSel)-numel(Tinder))==1
    0805         OutSel{count1} = CurSel{~ismember(CurSel,CurSel(Tinder))};
    0806     else %Nothing found or code cannot decide
    0807         OutSel{count1} = CurSel{1};
    0808     end
    0809 end
    0810 %--------------------------------------------------------------------------
    0811 function [OutFlx,OutRxn,OutBigg] = ChckMultiRxns(BiGG4KeggDraw,RxnCds,rxns4thismap)
    0812 MultiB = getappdata(0,'MultiB');
    0813 MultiK = getappdata(0,'MultiK');
    0814 MultiFlx = getappdata(0,'MultiFlx');
    0815
    0816 OutFlx =[];
    0817 OutBigg = ({});
    0818 OutRxn = ({});
    0819 while numel(MultiB)
    0820     Bigg1Loci = find(ismember(MultiB,MultiB{1}));
    0821     KeggTemp = MultiK(Bigg1Loci);
    0822     if all(ismember(KeggTemp,rxns4thismap)) % All child rxns are present in map
    0823         FlxTemp = MultiFlx(Bigg1Loci(1),:);
    0824         OutFlx = [OutFlx;FlxTemp];
    0825         OutRxn = [OutRxn;RxnCds(ismember(BiGG4KeggDraw,MultiB{1}))];
    0826         OutBigg = [OutBigg;MultiB{1}];
    0827     end
    0828     MultiB(Bigg1Loci) = [];
    0829     MultiK(Bigg1Loci) = [];
    0830     MultiFlx(Bigg1Loci,:) = [];
    0831 end
    0832 % -------------------------------------------------------------------------
    0833 function [rxnOut,flxOut] = ConsistRxns(BiGG4KeggDraw,rxnIn,flxIn,MapChoice)
    0834
    0835 % Read rxn2map ========================
    0836 Fileid1 = fopen('rxn2map.txt','r');
    0837 rxn2map = textscan(Fileid1,'%s %s');
    0838 RawRxns1 = rxn2map{1};
    0839 RawMaps = rxn2map{2};
    0840 [Loci1,~] = regexp(RawMaps, 'path:rn');
    0841 Loci2 = ~cellfun('isempty', Loci1);
    0842 RawMaps(Loci2) = [];
    0843 RawRxns1(Loci2) = [];
    0844 Maps = strrep(RawMaps,'path:map','');
    0845 Rxns1 = strrep(RawRxns1,'rn:','');
    0846 % Read ec2rxn =========================
    0847 Fileid1 = fopen('ec2rxn.txt','r');
    0848 ec2rxn = textscan(Fileid1,'%s %s');
    0849 RawEc= ec2rxn{1};
    0850 RawRxns = ec2rxn{2};
    0851 EC = strrep(RawEc,'ec:','');
    0852 Rxns = strrep(RawRxns,'rn:','');
    0853 FoundRxnsID = (0); RxnChk = rxnIn;
    0854 fclose(Fileid1);
    0855 % Read cpd2rxn =======================
    0856 Fileid = fopen('cpd2rxn.txt','r');
    0857 cpd2rxn = textscan(Fileid,'%s %s');
    0858 rxnData = strrep(cpd2rxn{2},'rn:','');
    0859 cpds = strrep(cpd2rxn{1},'cpd:','');
    0860 % Read cpd2weight ===================
    0861 Fileid = fopen('cpd2weight.txt','r');
    0862 cpd2weight = textscan(Fileid,'%s %f');
    0863 cpd4weight = strrep(cpd2weight{1},'cpd:','');
    0864 Weight = cpd2weight{2};
    0865 fclose(Fileid);
    0866
    0867 fprintf('\nRxn consistency for %s(of %d): ','All pathways',numel(rxnIn))
    0868 ctr = 1; cnter = 1; OriginalRxns = ({}); ConsRxns = ({});
    0869 for count = 1:numel(rxnIn)
    0870     Rns4Map = Maps(ismember(Rxns1,rxnIn{count}));
    0871       if count>1
    0872           for j=0:log10(count-1)
    0873               fprintf('\b');
    0874           end
    0875       end
    0876       fprintf('%d', count);
    0877     if ~any(ismember(Rns4Map,MapChoice)) %Otherwise there is no need to identify
    0878         % reactions for this map
    0879         FoundRxnsID(ctr) = count;
    0880         ctr = ctr+1;
    0881         TempEC=EC(ismember(Rxns,rxnIn{count}));
    0882         TempRns = Rxns(ismember(EC,TempEC));
    0883         TempRns = unique(TempRns);
    0884         TempRns(ismember(TempRns,rxnIn{count})) = [];
    0885         TempFlx = flxIn(count,:);
    0886         CurLen = numel(rxnIn) + 1;
    0887         cpds4curRxn = cpds(ismember(rxnData,rxnIn{count}));
    0888         W4curCpds = Weight(ismember(cpd4weight,cpds4curRxn));
    0889         for count1 = 1:numel(TempRns)
    0890             TempRns2Map1 = Maps(ismember(Rxns1,TempRns{count1}));
    0891             ChkCur = any(ismember(RxnChk,TempRns{count1}));
    0892             if any(ismember(TempRns2Map1,MapChoice)) && ~ChkCur
    0893                 cpds4TempRxn = cpds(ismember(rxnData,TempRns{count1}));
    0894                 W4TempCpds = Weight(ismember(cpd4weight,cpds4TempRxn));
    0895                 if numel(W4curCpds) == numel(W4TempCpds) &&...
    0896                         isempty(setdiff(W4curCpds,W4TempCpds)) && ...
    0897                         ~any(strcmp(TempRns{count1},{'R02630','R03232'})) % Exceptional condition
    0898                     rxnIn{CurLen} = TempRns{count1};
    0899                     flxIn(CurLen,:) = TempFlx;
    0900                     BiGG4KeggDraw{CurLen} = BiGG4KeggDraw{count};
    0901                     CurLen = CurLen+1;
    0902                     OriginalRxns{cnter} = rxnIn{count};
    0903                     ConsRxns{cnter} = TempRns{count1};
    0904                     cnter = cnter + 1;
    0905                 end
    0906             end
    0907         end
    0908     end
    0909 end
    0910 fprintf('\n')
    0911 if FoundRxnsID
    0912     flxIn(FoundRxnsID,:) = [];
    0913     rxnIn(FoundRxnsID) = [];
    0914     BiGG4KeggDraw(FoundRxnsID) = [];
    0915     setappdata(0,'BiGG4KeggDraw1',BiGG4KeggDraw);
    0916 end
    0917 rxnOut = rxnIn;
    0918 flxOut = flxIn;
    0919 if ~isempty(OriginalRxns) % A map between consisten rxns and those in model
    0920     CRxns.O = OriginalRxns; CRxns.C = ConsRxns;
    0921     setappdata(0,'CRxns',CRxns)
    0922 end