KeggDraw

 KeggDraw
 uses KGML files to extract KEGG rxn IDs and their
 corresponding coordinates on each map image (extracted from
 rest.kegg.jp). Based on reactions in each map, a comparison between
 current map reactions and input reaction codes is performed to choose
 those that are present in the selected pathway. Then, using flx2col and
 performing the normalization, flux values will be placed on corresponding
 rxn boxes in form of colour intensity. Model names provided in input flux
 data are used to check if data belong to one|two models. In latter, the
 box is divided in two separate and equal sections, each for one model.
 In case of dynamic data, the data are displayed on each pathway,
 consecutively. This would allow for a better understanding of the dynamic
 behaviour of target systems.
 
 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.
 flxType = An integer: 1- Entirely fills the boxes and 2- divides them into
           two identical parts.
 
 Optional inputs:
 InOpts = A struct containing the following fields:
 SVal = Check if user wants to save images : 1 (yes) or 0 (no), default:0.
 mTime = Time points corresponding to time-series values.
 Thresh = Threshold below which the flux data won't be depicted (default:
          1e-3).
 Timevalue = The time step between two consecutive images (default: 0.1 s).
 TimeUnit = A string representing time unit for time-series values
           (default: 'No unit').
 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
        file: 1 (video), 2 (GIF), 0 (no); default: 0.
 Viso = Check if user wants to adjust visual properties (text color, box
        color and font characteristics) through a GUI: 1 (activate GUI), 0
        (deactivate it); default: 0.
 
 O. Jamialahmadi
 TMU, Chem. Eng. Dept., Biotech. Group 
 Jan. 2016

0001 function KeggDraw (Outflx,RxnCds,MapChoice,flxType,InOpts)
0002 % KeggDraw
0003 % uses KGML files to extract KEGG rxn IDs and their
0004 % corresponding coordinates on each map image (extracted from
0005 % rest.kegg.jp). Based on reactions in each map, a comparison between
0006 % current map reactions and input reaction codes is performed to choose
0007 % those that are present in the selected pathway. Then, using flx2col and
0008 % performing the normalization, flux values will be placed on corresponding
0009 % rxn boxes in form of colour intensity. Model names provided in input flux
0010 % data are used to check if data belong to one|two models. In latter, the
0011 % box is divided in two separate and equal sections, each for one model.
0012 % In case of dynamic data, the data are displayed on each pathway,
0013 % consecutively. This would allow for a better understanding of the dynamic
0014 % behaviour of target systems.
0015 %
0016 % Inputs:
0017 % Outflx = A mXn matrix containing input flux data, which m corresponds to
0018 %        KEGG reaction IDs and n corresponds to time-series values.
0019 % RxnCds = A cell array of strings comprising of KEGG reaction IDs (m rows).
0020 % MapChoice = KEGG map ID.
0021 % flxType = An integer: 1- Entirely fills the boxes and 2- divides them into
0022 %           two identical parts.
0023 %
0024 % Optional inputs:
0025 % InOpts = A struct containing the following fields:
0026 % SVal = Check if user wants to save images : 1 (yes) or 0 (no), default:0.
0027 % mTime = Time points corresponding to time-series values.
0028 % Thresh = Threshold below which the flux data won't be depicted (default:
0029 %          1e-3).
0030 % Timevalue = The time step between two consecutive images (default: 0.1 s).
0031 % TimeUnit = A string representing time unit for time-series values
0032 %           (default: 'No unit').
0033 % Netstat = Cehck if user wants to use KEGG maps and KGML files online(1)
0034 %           /offline(0) (default = 1). Data in KEGGmpas folder will be
0035 %           used in case of offline mode.
0036 % Svid = Check if user wants to convert a sequence of images to an animated
0037 %        file: 1 (video), 2 (GIF), 0 (no); default: 0.
0038 % Viso = Check if user wants to adjust visual properties (text color, box
0039 %        color and font characteristics) through a GUI: 1 (activate GUI), 0
0040 %        (deactivate it); default: 0.
0041 %
0042 % O. Jamialahmadi
0043 % TMU, Chem. Eng. Dept., Biotech. Group
0044 % Jan. 2016
0045 
0046 % Optional inputs ---------------------------------------------------------
0047 if exist('InOpts','var')
0048     OptFields = fieldnames(InOpts);
0049     for count = 1:numel(OptFields)
0050         if strcmpi('SVal',OptFields{count})
0051             SVal = InOpts.(OptFields{count});
0052         elseif strcmpi('mTime',OptFields{count})
0053             mTime = InOpts.(OptFields{count});
0054         elseif strcmpi('Thresh',OptFields{count})
0055             Thresh = InOpts.(OptFields{count});
0056         elseif strcmpi('Timevalue',OptFields{count})
0057             Timevalue = InOpts.(OptFields{count});
0058         elseif strcmpi('TimeUnit',OptFields{count})
0059             TimeUnit = InOpts.(OptFields{count});
0060         elseif strcmpi('Netstat',OptFields{count})
0061             Netstat = InOpts.(OptFields{count});
0062         elseif strcmpi('Svid',OptFields{count})
0063             Svid = InOpts.(OptFields{count});
0064         elseif strcmpi('Viso',OptFields{count})
0065             Viso = InOpts.(OptFields{count});
0066         end
0067     end
0068 end
0069 % -------------------------------------------------------------------------
0070 % Check internet connection
0071 if ~exist('Netstat','var') || isempty(Netstat)
0072     Netstat = 1;
0073 end
0074 TestLink='http://rest.kegg.jp';
0075 if Netstat
0076     [~,Stat1]=urlread(TestLink);
0077     if ~Stat1
0078         msgbox({'This step requires a stable internet connection';...
0079             'Seemingly you are offline!'});
0080         return
0081     end
0082 end
0083 Pth1 = which ('Bigg2Kegg.m');
0084 tind = find(Pth1=='\',1,'last');
0085 Pth = Pth1(1:tind-1);
0086 %==========================================================================
0087 BiGG4KeggDraw = getappdata(0,'BiGG4KeggDraw'); % For further modification
0088 for mpct = 1:numel(MapChoice)
0089     h=waitbar(0,'Initializing KeggDraw...');
0090     if Netstat
0091         Link=['http://rest.kegg.jp/get/rn',MapChoice{mpct},'/kgml'];
0092         ImgLink=['http://www.kegg.jp/kegg/pathway/map/map',MapChoice{mpct},'.png'];
0093         I=imread(ImgLink);
0094         Dat=urlread(Link);
0095     else
0096         load([Pth,'\KEGGmaps\',MapChoice{mpct},'.mat']);
0097     end
0098     waitbar(20/100); % 20% of the process
0099     % Looking for KEGG rxn IDs and corresponding coordinates===============
0100     RawRNs = regexp(Dat,'(?<=reaction="rn:)[^"]*','match');
0101     IdRec = regexp(Dat,'\<rectangle" x="\d*" y="\d*" \w*="\d*" \w*="\d*"');
0102     IdLine = regexp(Dat, '(?<=type="line" coords=")[^"]*');
0103     CoordID = [IdRec,IdLine]; CoordID = sort(CoordID);
0104     Recseg = regexp(Dat,'\<rectangle" x="\d*" y="\d*" \w*="\d*" \w*="\d*"',...
0105         'match');
0106     RecsegCel = regexp(Recseg', '\d+', 'match');
0107     Reccoord = cellfun(@str2double, RecsegCel, 'UniformOutput', false);
0108     Lineseg = regexp(Dat, '(?<=type="line" coords=")[^"]*', 'match');
0109     LinesegCel = regexp(Lineseg', '\d+', 'match');
0110     Linecoord1 = cellfun(@str2double, LinesegCel, 'UniformOutput', false);
0111     Linecoord = cell(numel(Linecoord1),1);
0112     for linect = 1:numel(Linecoord1)
0113         TempLinecoord = Linecoord1{linect};
0114         TempLineX = TempLinecoord(1:2:end);
0115         TempLineY = TempLinecoord(2:2:end);
0116         Linecoord{linect} = [sum(TempLineX)./numel(TempLineX)+22,...
0117             sum(TempLineY)./numel(TempLineY)+8];
0118     end
0119     AllCoord = cell(1,numel(RawRNs));
0120     AllCoord(ismember(CoordID,IdLine)) = Linecoord;
0121     AllCoord(ismember(CoordID,IdRec)) = Reccoord;
0122     if ~isempty(Reccoord)
0123         TempW4Line = Reccoord{1}(3);
0124         TempH4Line = Reccoord{1}(4);
0125     else
0126         TempW4Line = 46;
0127         TempH4Line = 17;
0128     end
0129     idx = zeros(numel(AllCoord),1);idy = zeros(numel(AllCoord),1);
0130     idw = zeros(numel(AllCoord),1);idh = zeros(numel(AllCoord),1);
0131     for coordct = 1:numel(AllCoord)
0132         if numel(AllCoord{coordct}) < 4
0133             AllCoord{coordct}(3) = TempW4Line;
0134             AllCoord{coordct}(4) = TempH4Line;
0135         end
0136         idx(coordct) = round(AllCoord{coordct}(1));
0137         idy(coordct) = round(AllCoord{coordct}(2));
0138         idw(coordct) = AllCoord{coordct}(3);
0139         idh(coordct) = AllCoord{coordct}(4);    
0140     end
0141     LenRawRN = numel(RawRNs)+1;
0142     RN = RawRNs;
0143     if numel(RN)~=numel(idx)
0144         waitfor(msgbox(['Number of rxns and rectangles are not same ',...
0145             'for map:',MapChoice{mpct}]));
0146         close(h)
0147         continue
0148     end
0149 
0150     for rnct = 1:numel(RawRNs)
0151         if numel(RawRNs{rnct}) > 6
0152             TempCurRN = regexp(RawRNs{rnct}, 'R\d{5}', 'match');
0153             for rnct1 = 2:numel(TempCurRN)
0154                 RN{LenRawRN} = TempCurRN{rnct1};
0155                 idx(LenRawRN) = idx(rnct);
0156                 idy(LenRawRN) = idy(rnct);
0157                 idw(LenRawRN) = idw(rnct);
0158                 idh(LenRawRN) = idh(rnct);
0159                 LenRawRN = LenRawRN + 1;
0160             end
0161             RN{rnct} = RawRNs{rnct}(1:6);
0162         end
0163     end
0164     waitbar (40/100); % 40% of the process
0165     % By experience, the difference of .png images and coordinates in KGML
0166     idy=idy-8;
0167     idx=idx-22;
0168     %==========================================================================
0169     %% Analyzing flux data ====================================================
0170     % RxnCds should be compared to RN, in order to identify which reactions
0171     % from the Outflx data are present in current KEGG map
0172     if isempty(Outflx)
0173         msgbox('Flux data must be read first!','Error','error');
0174         return
0175     end
0176 
0177     % Cehck rxnCodes with HMR IDs and replace them with KEGG ones==========
0178     Hind = strfind(RxnCds,'HMR');
0179     if sum([Hind{:}])
0180         fname = 'HMR2KEGG.mat';
0181         if ~exist(fname,'file')
0182              waitfor(msgbox(['The ',fname,' Does not',...
0183                  ' exist in local directory.',...
0184                  'Open it manually'],'Warning','warn'));
0185              [pname,fname1]=uigetfile({'*.mat';'*.xlsx'},...
0186                  ['Select the ',fname,' file']);
0187              D=importdata([fname1,pname]);
0188         else
0189             D= importdata(fname);
0190         end
0191         KEGG=D.KEGG;HMR=D.HMR;
0192         [Fr,Sn]=ismember(RxnCds,HMR);
0193         for ind=1:length(Fr)
0194             if Fr(ind)
0195                 RxnCds{ind}=KEGG{Sn(ind)};
0196             end
0197         end
0198     end 
0199     %======================================================================
0200     [rxnCds,flx1] = ConsistRxns(BiGG4KeggDraw,RxnCds,Outflx,MapChoice{mpct});
0201     waitbar (60/100); % 50% of the process
0202     [f1,~]=ismember(rxnCds,RN);
0203     flx=flx1(f1,:);
0204     % Check repeated rxns on this map and let user choose one,among several
0205     % BiGG IDs, for which the flux value will be displayed on the map.
0206     [flx,rxnCds] = ChckRptRxns(flx,f1,rxnCds,MapChoice{mpct});
0207     waitbar(80/100); % 80% of the process
0208     copt2 = 1; idxt=(0);idyt=(0);idwt=(0);idht=(0);flx2=zeros(1,size(flx,2));
0209     rxnCds1={0};
0210     for copt = 1:numel(rxnCds)
0211         Temprxn = find(ismember(RN,rxnCds{copt}));
0212         for copt1 = 1:numel(Temprxn)
0213             idxt(copt2) = idx(Temprxn(copt1));
0214             idyt(copt2) = idy(Temprxn(copt1));
0215             idwt(copt2) = idw(Temprxn(copt1));
0216             idht(copt2) = idh(Temprxn(copt1));
0217             flx2(copt2,:) = flx(copt,:);
0218             rxnCds1{copt2} = rxnCds{copt};
0219             copt2 = copt2 + 1;
0220         end
0221     end
0222     flx = flx2;
0223     % Check multi-step reactions and identify if all the child rxns are
0224     % present in the current pathway, otherwise they won't be considered.
0225     % Also, there are reactions for which there is more than one BiGG ID.
0226     % Note that this case is different from ChckRptRxns function, because
0227     % these reactions are parts of multi-step reactions, and therefore, are
0228     % not equivalent.
0229     [OutFlx,IdxN,IdyN,IdwN,IdhN] = ChckMultiRxns(RN,idx,idy,idw,idh);
0230     waitbar (90/100); % 90% of the process
0231     
0232     flx2 = [flx;OutFlx];
0233     idxt1=[idxt';IdxN]; idyt1=[idyt';IdyN]; idwt1=[idwt';IdwN];idht1=[idht';IdhN];
0234     flx =flx2; idxt=idxt1;idyt=idyt1;idwt=idwt1;idht=idht1;
0235     %----------------------------------------------------------------------
0236     if isempty(flx)
0237         uiwait(msgbox('The flux data are not is current pathway!',...
0238             'Warning','warn'));
0239         close (h);
0240         return
0241     end
0242     %======================================================================
0243     %% Visualization=======================================================
0244     % Threshold below which the flux data won't be depicted:
0245     if ~exist('Thresh','var') || isempty(Thresh)
0246         Thresh=1e-6;
0247     end
0248     Cm=jet(size(idx,1)).*255; % [Warning] Changing Cm matrix is not recommended
0249 
0250     % The time step between two consecutive images:
0251     if ~exist('Timevalue','var') || isempty(Timevalue)
0252         Timevalue=0.1;
0253     end
0254     waitbar(100/100); % 100% of the process
0255     close (h);
0256 
0257     % Time points:
0258     if ~exist('mTime','var') || isempty(mTime)
0259         if flxType == 2
0260             mTime=1:size(Outflx,2)/2;
0261         else
0262             mTime=1:size(Outflx,2);
0263         end
0264     end
0265 
0266     if ~exist('TimeUnit','var') || isempty(TimeUnit)
0267         TimeUnit='No unit';
0268     end
0269     if ~exist('Svid','var') || isempty(Svid)
0270         Svid = 0; %Does not convert images to video
0271     end
0272     % Check wether user wants to save images (using 'Save images checkbox')
0273     if ~exist('SVal','var') || isempty(SVal)
0274         SVal = 0; % Does not save images.
0275     end
0276     if SVal
0277         FolderNm=[MapChoice{mpct},' images files'];
0278         Tn=uigetdir(pwd,['Folder to save images. To save in current folder:',...
0279         'Cancel']);
0280         if ~Tn
0281             if ~exist(FolderNm,'file') % Check if folder already exists
0282                 mkdir(FolderNm);
0283                 FolderPath=fullfile(pwd,FolderNm);
0284             else
0285                 FolderPath=fullfile(pwd,FolderNm);
0286             end
0287         else
0288             if ~exist(fullfile(Tn,FolderNm),'file') % Check if folder already exists
0289                 mkdir(fullfile(Tn,FolderNm));
0290                 FolderPath=fullfile(Tn,FolderNm);
0291             else
0292                 FolderPath=fullfile(Tn,FolderNm);
0293             end
0294         end
0295     end
0296     %======================================================================
0297     %Draw images on this panel
0298     figure(mpct)
0299     H1 = gcf;
0300     h1 = gca;
0301     set(H1,'position',get(0,'screensize'));
0302     set(h1,'position',[0,0,1,1],'xtick',[],'ytick',[]);
0303     set(H1,'toolbar','figure');
0304     set(H1,'menubar','figure');
0305     set(H1,'name',MapChoice{mpct},'numbertitle','off')
0306     % Arbitrary color and font properties
0307     if ~exist('Viso','var') || ~Viso
0308         TextFont.FontName = 'Arial';
0309         TextFont.FontWeight = 'normal';
0310         TextFont.FontAngle = 'normal';
0311         TextFont.FontUnits = 'points';
0312         TextFont.FontSize = 8;
0313         TextCol = [0,0,0];
0314         BoxCol = [240,248,255];
0315     end
0316     if exist('Viso','var') && Viso
0317         waitfor(VisualProp)
0318         BoxCol = getappdata(0,'BoxtCol');
0319         if ~isempty(BoxCol)
0320             BoxCol = BoxCol.*255;
0321         else
0322             BoxCol = [240,248,255];
0323         end
0324         TextFont = getappdata(0,'TextFont');
0325         if isempty(TextFont)
0326             TextFont.FontName = 'Arial';
0327             TextFont.FontWeight = 'normal';
0328             TextFont.FontAngle = 'normal';
0329             TextFont.FontUnits = 'points';
0330             TextFont.FontSize = 8;
0331         end
0332         TextCol = getappdata(0,'TextCol');
0333         if isempty(TextCol)
0334             TextCol = [0,0,0];
0335         end
0336     end
0337     switch flxType
0338         case 1
0339             for k=1:size(flx,2)
0340                 for m=1:size(flx,1)
0341                     % Normalize colors and insertion into enzyme boxes
0342                     %flxMax=max(abs(flx(:,k)));
0343                     %flxMin=min(abs(flx(:,k)));
0344                     flxMax=max(abs(flx(:)));
0345                     flxMin=min(abs(flx(:)));
0346                     if flxMax==flxMin % One row (i.e. for 1 rxn)
0347                         %flxMax=max(abs(flx(m,:)));
0348                         %flxMin=min(abs(flx(m,:)));
0349                         flxMax=max(abs(flx(:)));
0350                         flxMin=min(abs(flx(:)));
0351                     end
0352                     Tempx=(idyt(m)+1:idyt(m)+idht(m)-1);
0353                     I= flx2col(BoxCol,flx,Cm,I,k,flxMax,flxMin,m,idxt,idwt,Tempx,Thresh);
0354                     % Show the image:
0355                     imshow(I,'InitialMagnification', 'fit','Parent',h1); 
0356                     % Insert flux value on each reaction box
0357                     TextPropUDF (TextFont,TextCol,'One',idxt,idyt,idwt,...
0358                         idht,k,flx,h1,Thresh)
0359 
0360                     % Show time-series values on the corresponding image
0361                     TimePropUDF ('One',TextFont,TimeUnit,mTime,k,h1) 
0362                 end
0363                 % Show colorbar
0364                 flxstr=linspace(flxMin,flxMax,7);
0365                 flxstr1=cell(numel(flxstr),1);
0366                 for ir=1:numel(flxstr)
0367                     flxstr1{ir,1}=num2str(flxstr(ir),2);
0368                 end
0369                 ch = colorbar('peer',h1);
0370                 set(ch,'yticklabel',flxstr1)
0371 
0372                 % Save to image files within a new folder======================
0373                 if SVal
0374                     if k == 1
0375                         waitfor(SaveImgs)
0376                         ImgRes = getappdata(0,'ImgRes');
0377                         ImgFrmt = getappdata(0,'ImgFrmt');
0378                         if isempty(ImgRes)
0379                             ImgRes=1; % Corresponds to 300 dpi in SaveFlxImgs
0380                         elseif isempty(ImgFrmt)
0381                             ImgFrmt=1; % Corresponds to PNG in SaveFlxImgs
0382                         end
0383                     end
0384                     SaveFlxImgs(H1,FolderPath,k,ImgRes,ImgFrmt);
0385                 end
0386                 %==============================================================
0387                 pause(Timevalue); % Time step between consecutive images
0388             end
0389         case 2
0390             for k=1:2:size(flx,2)
0391                 for m=1:size(flx,1)
0392                     % Normalize colors and insertion into enzyme boxes
0393                     %flxMax=max(max(abs(flx(:,k))),max(abs(flx(:,k+1))));
0394                     %flxMin=min(min(abs(flx(:,k))),min(abs(flx(:,k+1))));
0395                     flxMax=max(abs(flx(:)));
0396                     flxMin=min(abs(flx(:)));
0397                     if size(flx,1)==1
0398                         %flxMax=max(abs(flx(m,:)));
0399                         %flxMin=min(abs(flx(m,:)));
0400                         flxMax=max(abs(flx(:)));
0401                         flxMin=min(abs(flx(:)));
0402                     end
0403                     Boxheight_mid=(idht(m)-2)/2;
0404                     Tempx=(idyt(m)+1:idyt(m)+idht(m)-1-Boxheight_mid);
0405                     Tempx1=round((idyt(m)+idht(m)-1-Boxheight_mid:idyt(m)...
0406                         +idht(m)-1));
0407                     % For 1st flx set:
0408                     I= flx2col(BoxCol,flx,Cm,I,k,flxMax,flxMin,m,idxt,idwt,Tempx,Thresh); 
0409                     % For 2nd flx set:
0410                     I= flx2col(BoxCol,flx,Cm,I,k+1,flxMax,flxMin,m,idxt,idwt,Tempx1,Thresh);
0411                     % Show the image:
0412                     imshow(I,'InitialMagnification', 'fit','Parent',h1); 
0413 
0414                     % Insert flux value on each reaction box for 1st set(upper)
0415                     TextPropUDF (TextFont,TextCol,'Two',idxt,idyt,...
0416                         idwt,idht,k,flx,h1,Thresh)
0417                     % Show time-series values on the corresponding image
0418                     TimePropUDF ('Two',TextFont,TimeUnit,mTime,k,h1) 
0419                 end
0420 
0421                 % Show colorbar
0422                 flxstr=linspace(flxMin,flxMax,7);
0423                 flxstr1=cell(numel(flxstr),1);
0424                 for ir=1:numel(flxstr)
0425                     flxstr1{ir,1}=num2str(flxstr(ir),2);
0426                 end
0427                 ch = colorbar('peer',h1);
0428                 set(ch,'yticklabel',flxstr1)
0429 
0430                 % Save to image files within a new folder======================
0431                 if SVal
0432                     if k == 1
0433                         waitfor(SaveImgs)
0434                         ImgRes = getappdata(0,'ImgRes');
0435                         ImgFrmt = getappdata(0,'ImgFrmt');
0436                         if isempty(ImgRes)
0437                             ImgRes=1; % Corresponds to 300 dpi in SaveFlxImgs
0438                         elseif isempty(ImgFrmt)
0439                             ImgFrmt=1; % Corresponds to PNG in SaveFlxImgs
0440                         end
0441                     end
0442                     SaveFlxImgs(H1,FolderPath,ceil(k/2),ImgRes,ImgFrmt);
0443                 end
0444                 %==============================================================
0445                 pause(Timevalue); % Time step between consecutive images
0446             end
0447     end
0448 if Svid
0449     Img2Animation(Svid)
0450 end
0451 end
0452 
0453 % Subfunctions ------------------------------------------------------------
0454 function [flxOut,rxnCdsI] = ChckRptRxns (flx,f1,rxnCds,CurrentMap)
0455 % First, rxns are pruned based on manual assessment of pathways
0456 BiGG4KeggDraw = getappdata(0,'BiGG4KeggDraw1');
0457 rxnCdsI = rxnCds(f1);
0458 BiGGI = BiGG4KeggDraw(f1);
0459 Fileid1 = fopen('rxnbiggmap.txt','r');
0460 Rxnbmap = textscan(Fileid1,'%s %s %s %d');
0461 Rtemp = Rxnbmap{1}; Btemp = Rxnbmap{2}; Mtemp = Rxnbmap{3};
0462 Dtemp = Rxnbmap{4};
0463 if any(ismember(Mtemp,CurrentMap))
0464     ChR = Rtemp(ismember(Mtemp,CurrentMap));
0465     ChB = Btemp(ismember(Mtemp,CurrentMap));
0466     ChD = Dtemp(ismember(Mtemp,CurrentMap));
0467     for ct = 1:numel(ChR)
0468         n1_chr = find(ismember(rxnCdsI,ChR{ct}));
0469         n1_chb = find(ismember(BiGGI,ChB{ct}));
0470         if ~isempty(intersect(n1_chr,n1_chb))
0471             %Rmvit = find(ismember(BiGGI,ChB{ct}));
0472             Rmvit = intersect(n1_chr,n1_chb);
0473             for ct1 = 1:numel(Rmvit)
0474                 if strcmp(BiGGI{Rmvit(ct1)},ChB{ct})
0475                     Tmpflx = flx(Rmvit(ct1),:);
0476                 end
0477             end
0478             rxnCdsI(Rmvit) = [];
0479             BiGGI(Rmvit) = [];
0480             flx(Rmvit,:) = [];
0481             if ChD(ct)
0482                 CurL1 = numel(rxnCdsI) + 1;
0483                 rxnCdsI{CurL1} = ChR{ct};
0484                 BiGGI{CurL1} = ChB{ct};
0485                 flx(CurL1,:) = Tmpflx;
0486             end
0487         end
0488     end
0489 end
0490 %------
0491 [R1,~,R3] = unique(rxnCdsI,'stable');
0492 countR = 1;RptCount=(0);
0493 while numel(R3)
0494     RptCount(countR) = numel(find(R3==R3(1)));
0495     R3(R3==R3(1)) = [];
0496     countR = countR+1;
0497 end
0498 FndRxnLoci = find(RptCount>1);
0499 RptRxns = {0}; RptBiggs = {0};RptFlx= {0};
0500 for tctr = 1:numel(FndRxnLoci)
0501     RptRxnsI = R1(FndRxnLoci(tctr));
0502     RptRxns{tctr} = rxnCdsI(ismember(rxnCdsI,RptRxnsI));
0503     RptBiggs{tctr} = BiGGI(ismember(rxnCdsI,RptRxnsI));
0504     RptFlxLoci = ismember(rxnCdsI,RptRxnsI);
0505     RptFlx{tctr} = flx(RptFlxLoci,:);
0506 end
0507 if any(RptCount>1) 
0508     setappdata(0,'RptBiggs',RptBiggs);
0509     setappdata(0,'CurrentMap',CurrentMap);
0510     setappdata(0,'RptRxns',RptRxns);
0511     setappdata(0,'RptFlx',RptFlx);
0512     waitfor(KeggDrawTable)
0513     SlctdIds = getappdata(0,'data1');
0514     if ~isempty(SlctdIds)
0515         NtSel = RptBiggs(ismember(SlctdIds,'Not slected'));
0516         OwnSel = StubSel(NtSel);
0517         SlctdIds(ismember(SlctdIds,'Not slected')) = OwnSel;
0518     else
0519         SlctdIds = StubSel(RptBiggs);
0520     end
0521     % Refine flx
0522     for count = 1:numel(SlctdIds)
0523         Fi = find(ismember(RptBiggs{count},SlctdIds{count}));
0524         for count1 = 1:size(RptFlx{count},1)
0525             RptFlx{count}(count1,:) = RptFlx{count}(Fi,:);
0526         end
0527     end
0528     for count1 = 1:size(RptRxns,2)
0529         Fi1 = find(ismember(rxnCdsI,RptRxns{count1}));
0530         flx(Fi1,:) = RptFlx{count1};
0531     end
0532     flxOut = flx;
0533     
0534 else
0535     flxOut = flx;
0536 end
0537 % Remove all appdata
0538 AT=getappdata(0);
0539 AT1 = fieldnames(AT);
0540 for co = 1:length(AT1)
0541     if ~strcmp(AT1{co},'BiGG4KeggDraw') && ~strcmp(AT1{co},'MultiB')...
0542             && ~strcmp(AT1{co},'MultiK') && ~strcmp(AT1{co},'MultiFlx')
0543         rmappdata(0,AT1{co});
0544     end
0545 end
0546 %--------------------------------------------------------------------------
0547 function OutSel = StubSel(InSel)
0548 OutSel = {0};
0549 for count1 = 1:numel(InSel)
0550     Tinder = regexp(InSel{count1},'\w*p$'); %Peroxisome
0551     Tinder = find(~cellfun('isempty', Tinder));
0552     CurSel = InSel{count1};
0553     if (numel(CurSel)-numel(Tinder))==1
0554         OutSel{count1} = CurSel{~ismember(CurSel,CurSel(Tinder))};
0555     else %Nothing found or code cannot decide
0556         OutSel{count1} = CurSel{1};
0557     end
0558 end
0559 %--------------------------------------------------------------------------
0560 function [OutFlx,IdxN,IdyN,IdwN,IdhN] = ChckMultiRxns(RN,idx,idy,idw,idh)
0561 MultiB = getappdata(0,'MultiB'); 
0562 MultiK = getappdata(0,'MultiK');
0563 MultiFlx = getappdata(0,'MultiFlx');
0564 
0565 OutFlx =[];
0566 IdxN=[];IdyN=[];IdwN=[];IdhN=[];
0567 while numel(MultiB)
0568     Bigg1Loci = ismember(MultiB,MultiB{1});
0569     KeggTemp = MultiK(Bigg1Loci);
0570     if all(ismember(KeggTemp,RN)) % All child rxns are present in map
0571         FlxTemp = MultiFlx(Bigg1Loci,:);
0572         [Tchck1,Tchck2] = ismember(RN,KeggTemp);
0573         if sum(Tchck1) > numel(KeggTemp) % Some KggTemp elements have more than
0574             % one equivalent in RN, meaning more than one set of coordinates
0575             % exists
0576             Tchck2(Tchck2==0)=[];
0577             U = unique(Tchck2);
0578             Nch = histc(Tchck2,U);
0579             RptFound = U(Nch>1);
0580             RptRxns = KeggTemp(RptFound); % RptRxns contains repeated rxns
0581             % in RN, which are shared between RN and KeggTemp.
0582             OtherRxns1 = KeggTemp(~ismember(KeggTemp,RptRxns));
0583             OtherIdx = idx(ismember(RN,OtherRxns1));
0584             OtherIdy = idy(ismember(RN,OtherRxns1));
0585             OtherIdw = idw(ismember(RN,OtherRxns1));
0586             OtherIdh = idh(ismember(RN,OtherRxns1));
0587             % Assumption: all child rxns are placed close together, so,
0588             % their coordinates are similar. Therefore, RptRxns having
0589             % similar coordinates to mean(OtherIdx & OtherIdy) wille be
0590             % selected.
0591             for count1 = 1:numel(RptRxns)
0592                 MeanIdx = mean(OtherIdx); MeanIdy = mean(OtherIdy);
0593                 RptIdx = idx(ismember(RN,RptRxns{count1}));
0594                 RptIdy = idy(ismember(RN,RptRxns{count1}));
0595                 RptIdw = idw(ismember(RN,RptRxns{count1}));
0596                 RptIdh = idh(ismember(RN,RptRxns{count1}));
0597                 RptDist = sqrt((RptIdx-MeanIdx).^2+(RptIdy-MeanIdy).^2);
0598                 [~,Id1] = min(RptDist);
0599                 CurLen = numel(OtherIdx)+1;
0600                 OtherIdx(CurLen) = RptIdx(Id1);
0601                 OtherIdy(CurLen) = RptIdy(Id1);
0602                 OtherIdw(CurLen) = RptIdw(Id1);
0603                 OtherIdh(CurLen) = RptIdh(Id1);
0604             end
0605         else
0606             OtherIdx = idx(ismember(RN,KeggTemp));
0607             OtherIdy = idy(ismember(RN,KeggTemp));
0608             OtherIdw = idw(ismember(RN,KeggTemp));
0609             OtherIdh = idh(ismember(RN,KeggTemp));
0610         end
0611         if size(OutFlx,1) == 0
0612             OutFlx = FlxTemp;
0613             IdxN = OtherIdx; IdyN = OtherIdy; IdwN = OtherIdw; 
0614             IdhN = OtherIdh;
0615         else
0616             OutFlx(size(OutFlx,1)+1:size(OutFlx,1)+size(FlxTemp,1),:)...
0617                 = FlxTemp;
0618             IdxN(numel(IdxN)+1:numel(IdxN)+numel(OtherIdx)) = OtherIdx;
0619             IdyN(numel(IdyN)+1:numel(IdyN)+numel(OtherIdy)) = OtherIdy;
0620             IdwN(numel(IdwN)+1:numel(IdwN)+numel(OtherIdw)) = OtherIdw;
0621             IdhN(numel(IdhN)+1:numel(IdhN)+numel(OtherIdh)) = OtherIdh;
0622         end
0623     end
0624     MultiB(Bigg1Loci) = [];
0625     MultiK(Bigg1Loci) = [];
0626     MultiFlx(Bigg1Loci,:) = [];
0627 end
0628 % -------------------------------------------------------------------------
0629 function [rxnOut,flxOut] = ConsistRxns(BiGG4KeggDraw,rxnIn,flxIn,MapChoice)
0630 
0631 % Read rxn2map ========================
0632 Fileid1 = fopen('rxn2map.txt','r');
0633 rxn2map = textscan(Fileid1,'%s %s');
0634 RawRxns1 = rxn2map{1};
0635 RawMaps = rxn2map{2};
0636 [Loci1,~] = regexp(RawMaps, 'path:rn');
0637 Loci2 = ~cellfun('isempty', Loci1);
0638 RawMaps(Loci2) = [];
0639 RawRxns1(Loci2) = [];
0640 Maps = strrep(RawMaps,'path:map','');
0641 Rxns1 = strrep(RawRxns1,'rn:',''); 
0642 % Read ec2rxn =========================
0643 Fileid1 = fopen('ec2rxn.txt','r');
0644 ec2rxn = textscan(Fileid1,'%s %s');
0645 RawEc= ec2rxn{1};
0646 RawRxns = ec2rxn{2};
0647 EC = strrep(RawEc,'ec:','');
0648 Rxns = strrep(RawRxns,'rn:','');
0649 FoundRxnsID = (0); RxnChk = rxnIn;
0650 % Read cpd2rxn =======================
0651 Fileid = fopen('cpd2rxn.txt','r');
0652 cpd2rxn = textscan(Fileid,'%s %s');
0653 rxnData = strrep(cpd2rxn{2},'rn:','');
0654 cpds = strrep(cpd2rxn{1},'cpd:','');
0655 % Read cpd2weight ===================
0656 Fileid = fopen('cpd2weight.txt','r');
0657 cpd2weight = textscan(Fileid,'%s %f');
0658 cpd4weight = strrep(cpd2weight{1},'cpd:','');
0659 Weight = cpd2weight{2};
0660 % Read map2cpd ===================
0661 % Fileid = fopen('map2cpd.txt','r');
0662 % map2cpd = textscan(Fileid,'%s %s');
0663 % map4cpdT = map2cpd{1};
0664 % cpd4mapT = map2cpd{2};
0665 % Tid = regexp(map4cpdT,'path:map\d*');
0666 % Tid = ~cellfun('isempty', Tid);
0667 % map4cpd = map4cpdT(Tid);
0668 % map4cpd = strrep(map4cpd,'path:map','');
0669 % cpd4map = cpd4mapT(Tid);
0670 % cpd4map = strrep(cpd4map,'cpd:','');
0671 fprintf('\nRxn consistency for map%s(of %d): ',MapChoice,numel(rxnIn))
0672 ctr = 1;
0673 for count = 1:numel(rxnIn)
0674     Rns4Map = Maps(ismember(Rxns1,rxnIn{count}));
0675       if count>1
0676           for j=0:log10(count-1)
0677               fprintf('\b'); 
0678           end
0679       end
0680       fprintf('%d', count);
0681     if ~any(ismember(Rns4Map,MapChoice)) %Otherwise there is no need to identify
0682         % reactions for this map
0683         FoundRxnsID(ctr) = count;
0684         ctr = ctr+1;
0685         TempEC=EC(ismember(Rxns,rxnIn{count}));
0686         TempRns = Rxns(ismember(EC,TempEC));
0687         TempRns = unique(TempRns);
0688         TempRns(ismember(TempRns,rxnIn{count})) = [];
0689         TempFlx = flxIn(count,:);
0690         CurLen = numel(rxnIn) + 1;
0691         cpds4curRxn = cpds(ismember(rxnData,rxnIn{count}));
0692         W4curCpds = Weight(ismember(cpd4weight,cpds4curRxn));
0693         for count1 = 1:numel(TempRns)
0694             TempRns2Map1 = Maps(ismember(Rxns1,TempRns{count1}));
0695             ChkCur = any(ismember(RxnChk,TempRns{count1}));
0696             if any(ismember(TempRns2Map1,MapChoice)) && ~ChkCur
0697                 cpds4TempRxn = cpds(ismember(rxnData,TempRns{count1}));
0698 %                 DiffCpds4Temp = cpds4TempRxn(~ismember(cpds4TempRxn,cpds4curRxn));
0699 %                 Maps4TempRxn = map4cpd(ismember(cpd4map,DiffCpds4Temp));
0700 %                 ChckTempMap1 = sum(ismember(Maps4TempRxn,MapChoice));
0701 %                 ChckTempMap = (numel(DiffCpds4Temp) == ChckTempMap1);
0702                 W4TempCpds = Weight(ismember(cpd4weight,cpds4TempRxn));
0703                 if numel(W4curCpds) == numel(W4TempCpds) &&...
0704                         isempty(setdiff(W4curCpds,W4TempCpds)) %&& ChckTempMap
0705                     rxnIn{CurLen} = TempRns{count1};
0706                     flxIn(CurLen,:) = TempFlx;
0707                     BiGG4KeggDraw{CurLen} = BiGG4KeggDraw{count};
0708                     CurLen = CurLen+1;
0709                 end
0710             end
0711         end
0712     end
0713 end
0714 fprintf('\n')
0715 if FoundRxnsID
0716     flxIn(FoundRxnsID,:) = [];
0717     rxnIn(FoundRxnsID) = [];
0718     BiGG4KeggDraw(FoundRxnsID) = [];
0719     setappdata(0,'BiGG4KeggDraw1',BiGG4KeggDraw);
0720 end
0721 rxnOut = rxnIn;
0722 flxOut = flxIn;