Description of KeggDraw

Home > BiKEGG > KeggDraw.m

# KeggDraw

## PURPOSE

KeggDraw

## SYNOPSIS

function KeggDraw (Outflx,RxnCds,MapChoice,flxType,InOpts)

## DESCRIPTION

     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

## CROSS-REFERENCE INFORMATION

This function calls:

- Img2Animation Img2Animation
- KeggDrawTable KeggDrawTable(GUI)
- SaveFlxImgs SavFlxImgs
- SaveImgs SaveImgs(GUI)
- TextPropUDF TextPropUDF
- TimePropUDF TimePropUDF
- TimeUnit TimeUnit(GUI)
- VisualProp VisualProp(GUI)
- flx2col flx2col

This function is called by:

## SUBFUNCTIONS

- function [flxOut,rxnCdsI] = ChckRptRxns (flx,f1,rxnCds,CurrentMap)
- function OutSel = StubSel(InSel)
- function [OutFlx,IdxN,IdyN,IdwN,IdhN] = ChckMultiRxns(RN,idx,idy,idw,idh)
- function [rxnOut,flxOut] = ConsistRxns(BiGG4KeggDraw,rxnIn,flxIn,MapChoice)

## SOURCE CODE

    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;