NetDraw

 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

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