markPathwayWithExpression

 markPathwayWithExpression
   Adds gene expression data to a pathway object

   pathway           map structure as generated by constructPathwayFromCelldesigner
   model             a model structure
   experiment        experiment structure as generated by getExpressionStructure
   experimentOrder   cell array with two experiment names as defined in
                     'experiment'

   pathway           updated pathway object

 Usage: pathway=markPathwayWithExpression(pathway,model,experiment,experimentOrder)

0001 function pathway=markPathwayWithExpression(pathway,model,experiment,experimentOrder)
0002 % markPathwayWithExpression
0003 %   Adds gene expression data to a pathway object
0004 %
0005 %   pathway           map structure as generated by constructPathwayFromCelldesigner
0006 %   model             a model structure
0007 %   experiment        experiment structure as generated by getExpressionStructure
0008 %   experimentOrder   cell array with two experiment names as defined in
0009 %                     'experiment'
0010 %
0011 %   pathway           updated pathway object
0012 %
0013 % Usage: pathway=markPathwayWithExpression(pathway,model,experiment,experimentOrder)
0014 
0015 if numel(experimentOrder)~=2
0016     EM='This can only be done for two cases at the moment (experimentOrder must be two elements)';
0017     dispEM(EM);
0018 end
0019 
0020 %Check that experiment fit with experimentOrder
0021 [present, expIds]=ismember(experimentOrder,experiment.experiments);
0022 
0023 if ~all(present)
0024     EM='Not all experiments could be found in the experiment structure';
0025     dispEM(EM);
0026 end
0027 
0028 experiment.data=experiment.data(:,expIds);
0029 experiment.experiments=experiment.experiments(:,expIds);
0030 
0031 %Go through each species in pathway and if it's a reaction, get it's genes
0032 %and log2-fold change in expression from model and experiment
0033 for i=1:numel(pathway.listOfSpecies)
0034     if strcmp('PROTEIN', pathway.listOfSpecies(i).type)
0035         if isfield(pathway.listOfSpecies(i),'note')
0036             if ~isempty(pathway.listOfSpecies(i).note)
0037                 %Get the reaction if present in model
0038                 [present, index]=ismember(pathway.listOfSpecies(i).note,model.rxns);
0039                 
0040                 %If present, then get the genes
0041                 if any(present)
0042                     [~, genes]=find(model.rxnGeneMat(index,:));
0043                     
0044                     %If it was associated with genes match them to the ORFs
0045                     if any(genes)
0046                         [present, experimentIndexes]=ismember(model.genes(genes),experiment.orfs);
0047                         
0048                         %Add annotation to pathway structure
0049                         if any(present)
0050                             pathway.listOfSpecies(i).orfs=experiment.orfs(experimentIndexes(present));
0051                             pathway.listOfSpecies(i).expA=experiment.data(experimentIndexes(present),1);
0052                             pathway.listOfSpecies(i).expB=experiment.data(experimentIndexes(present),2);
0053                         end
0054                     end
0055                 end
0056             end
0057         end
0058     end
0059 end
0060 end