Home > INIT > getINITModel.m

getINITModel

PURPOSE ^

getINITModel_legacy

SYNOPSIS ^

function [model, metProduction, essentialRxnsForTasks, addedRxnsForTasks, deletedDeadEndRxns, deletedRxnsInINIT, taskReport]=getINITModel(refModel, tissue, celltype, hpaData, arrayData, metabolomicsData, taskFile, useScoresForTasks, printReport, taskStructure)

DESCRIPTION ^

 getINITModel_legacy
   Generates a model using the INIT algorithm, based on proteomics and/or
   transcriptomics and/or metabolomics and/or metabolic tasks. This is the original 
   implementation of tINIT, which is replaced by ftINIT.

 Input:
   refModel            a model structure. The model should be in the
                       closed form (no exchange reactions open). Import
                       using import(filename,false). If the model is not
                       loaded using importModel, it might be that there
                       is no "unconstrained" field. In that case,
                       manually add the field like:
                       model.unconstrained=false(numel(model.mets),1);
   tissue              tissue to score for. Should exist in either
                       hpaData.tissues or arrayData.tissues
   celltype            cell type to score for. Should exist in either
                       hpaData.celltypes or arrayData.celltypes for this
                       tissue (optional, default is to use the best values
                       among all the cell types for the tissue. Use [] if
                       you want to supply more arguments)
   hpaData             HPA data structure from parseHPA (optional if arrayData is
                       supplied, default [])
   arrayData           gene expression data structure (optional if hpaData is
                       supplied, default [])
       genes           cell array with the unique gene names
       tissues         cell array with the tissue names. The list may not be
                       unique, as there can be multiple cell types per tissue
       celltypes       cell array with the cell type names for each tissue
       levels          GENESxTISSUES array with the expression level for
                       each gene in each tissue/celltype. NaN should be
                       used when no measurement was performed
       threshold       a single value or a vector of gene expression 
                       thresholds, above which genes are considered to be
                       "expressed". (optional, by default, the mean expression
                       levels of each gene across all tissues in arrayData
                       will be used as the threshold values)
       singleCells     binary value selecting whether to use the
                       single-cell algorithm to identify expressed genes.
                       If used, specify cell subpopulations in CELLTYPES
                       (optional, default [])
       plotResults     true if single cell probability distributions
                       should be plotted (optional, default = False)
   metabolomicsData    cell array with metabolite names that the model
                       should produce (optional, default [])
   taskFile            a task list in Excel format. See parseTaskList for
                       details (optional, default [])
   useScoresForTasks   true if the calculated reaction scored should be used as
                       weights in the fitting to tasks (optional, default true)
   printReport         true if a report should be printed to the screen
                       (optional, default true)
   taskStructure       task structure as from parseTaskList. Can be used
                       as an alternative way to define tasks when Excel
                       sheets are not suitable. Overrides taskFile (optional,
                       default [])

 Output:
   model                   the resulting model structure
   metProduction           array that indicates which of the
                           metabolites in metabolomicsData that could be
                           produced. Note that this is before the
                           gap-filling process to enable defined tasks. To
                           see which metabolites that can be produced in
                           the final model, use canProduce.
                           -2: metabolite name not found in model
                           -1: metabolite found, but it could not be produced
                           1: metabolite could be produced
   essentialRxnsForTasks   cell array of the reactions which were
                           essential to perform the tasks
   addedRxnsForTasks       cell array of the reactions which were added in
                           order to perform the tasks
   deletedDeadEndRxns      cell array of reactions deleted because they
                           could not carry flux (INIT requires a
                           functional input model)
   deletedRxnsInINIT       cell array of the reactions which were deleted by
                           the INIT algorithm
   taskReport              structure with the results for each task
       id                  cell array with the id of the task
       description         cell array with the description of the task
       ok                  boolean array with true if the task was successful
       essential           cell array with cell arrays of essential
                           reactions for the task
       gapfill             cell array of cell arrays of reactions included
                           in the gap-filling for the task

   This is the main function for automatic reconstruction of models based
   on the (t)INIT algorithm (PLoS Comput Biol. 2012;8(5):e1002518, 
   Mol Syst Biol. 2014;10:721). Not all settings are possible using this
   function, and you may want to call the functions scoreModel, runINIT
   and fitTasks individually instead.

   NOTE: Exchange metabolites should normally not be removed from the model
   when using this approach, since checkTasks/fitTasks rely on putting specific
   constraints for each task. The INIT algorithm will remove exchange metabolites
   if any are present. Use importModel(file,false) to import a model with
   exchange metabolites remaining.

 Usage: [model, metProduction, essentialRxnsForTasks, addedRxnsForTasks,...
               deletedDeadEndRxns, deletedRxnsInINIT, taskReport]=...
               getINITModel(refModel, tissue, celltype, hpaData, arrayData,...
               metabolomicsData, taskFile, useScoresForTasks, printReport,...
               taskStructure, params, paramsFT)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [model, metProduction, essentialRxnsForTasks, addedRxnsForTasks, deletedDeadEndRxns, deletedRxnsInINIT, taskReport]=getINITModel(refModel, tissue, celltype, hpaData, arrayData, metabolomicsData, taskFile, useScoresForTasks, printReport, taskStructure)
0002 % getINITModel_legacy
0003 %   Generates a model using the INIT algorithm, based on proteomics and/or
0004 %   transcriptomics and/or metabolomics and/or metabolic tasks. This is the original
0005 %   implementation of tINIT, which is replaced by ftINIT.
0006 %
0007 % Input:
0008 %   refModel            a model structure. The model should be in the
0009 %                       closed form (no exchange reactions open). Import
0010 %                       using import(filename,false). If the model is not
0011 %                       loaded using importModel, it might be that there
0012 %                       is no "unconstrained" field. In that case,
0013 %                       manually add the field like:
0014 %                       model.unconstrained=false(numel(model.mets),1);
0015 %   tissue              tissue to score for. Should exist in either
0016 %                       hpaData.tissues or arrayData.tissues
0017 %   celltype            cell type to score for. Should exist in either
0018 %                       hpaData.celltypes or arrayData.celltypes for this
0019 %                       tissue (optional, default is to use the best values
0020 %                       among all the cell types for the tissue. Use [] if
0021 %                       you want to supply more arguments)
0022 %   hpaData             HPA data structure from parseHPA (optional if arrayData is
0023 %                       supplied, default [])
0024 %   arrayData           gene expression data structure (optional if hpaData is
0025 %                       supplied, default [])
0026 %       genes           cell array with the unique gene names
0027 %       tissues         cell array with the tissue names. The list may not be
0028 %                       unique, as there can be multiple cell types per tissue
0029 %       celltypes       cell array with the cell type names for each tissue
0030 %       levels          GENESxTISSUES array with the expression level for
0031 %                       each gene in each tissue/celltype. NaN should be
0032 %                       used when no measurement was performed
0033 %       threshold       a single value or a vector of gene expression
0034 %                       thresholds, above which genes are considered to be
0035 %                       "expressed". (optional, by default, the mean expression
0036 %                       levels of each gene across all tissues in arrayData
0037 %                       will be used as the threshold values)
0038 %       singleCells     binary value selecting whether to use the
0039 %                       single-cell algorithm to identify expressed genes.
0040 %                       If used, specify cell subpopulations in CELLTYPES
0041 %                       (optional, default [])
0042 %       plotResults     true if single cell probability distributions
0043 %                       should be plotted (optional, default = False)
0044 %   metabolomicsData    cell array with metabolite names that the model
0045 %                       should produce (optional, default [])
0046 %   taskFile            a task list in Excel format. See parseTaskList for
0047 %                       details (optional, default [])
0048 %   useScoresForTasks   true if the calculated reaction scored should be used as
0049 %                       weights in the fitting to tasks (optional, default true)
0050 %   printReport         true if a report should be printed to the screen
0051 %                       (optional, default true)
0052 %   taskStructure       task structure as from parseTaskList. Can be used
0053 %                       as an alternative way to define tasks when Excel
0054 %                       sheets are not suitable. Overrides taskFile (optional,
0055 %                       default [])
0056 %
0057 % Output:
0058 %   model                   the resulting model structure
0059 %   metProduction           array that indicates which of the
0060 %                           metabolites in metabolomicsData that could be
0061 %                           produced. Note that this is before the
0062 %                           gap-filling process to enable defined tasks. To
0063 %                           see which metabolites that can be produced in
0064 %                           the final model, use canProduce.
0065 %                           -2: metabolite name not found in model
0066 %                           -1: metabolite found, but it could not be produced
0067 %                           1: metabolite could be produced
0068 %   essentialRxnsForTasks   cell array of the reactions which were
0069 %                           essential to perform the tasks
0070 %   addedRxnsForTasks       cell array of the reactions which were added in
0071 %                           order to perform the tasks
0072 %   deletedDeadEndRxns      cell array of reactions deleted because they
0073 %                           could not carry flux (INIT requires a
0074 %                           functional input model)
0075 %   deletedRxnsInINIT       cell array of the reactions which were deleted by
0076 %                           the INIT algorithm
0077 %   taskReport              structure with the results for each task
0078 %       id                  cell array with the id of the task
0079 %       description         cell array with the description of the task
0080 %       ok                  boolean array with true if the task was successful
0081 %       essential           cell array with cell arrays of essential
0082 %                           reactions for the task
0083 %       gapfill             cell array of cell arrays of reactions included
0084 %                           in the gap-filling for the task
0085 %
0086 %   This is the main function for automatic reconstruction of models based
0087 %   on the (t)INIT algorithm (PLoS Comput Biol. 2012;8(5):e1002518,
0088 %   Mol Syst Biol. 2014;10:721). Not all settings are possible using this
0089 %   function, and you may want to call the functions scoreModel, runINIT
0090 %   and fitTasks individually instead.
0091 %
0092 %   NOTE: Exchange metabolites should normally not be removed from the model
0093 %   when using this approach, since checkTasks/fitTasks rely on putting specific
0094 %   constraints for each task. The INIT algorithm will remove exchange metabolites
0095 %   if any are present. Use importModel(file,false) to import a model with
0096 %   exchange metabolites remaining.
0097 %
0098 % Usage: [model, metProduction, essentialRxnsForTasks, addedRxnsForTasks,...
0099 %               deletedDeadEndRxns, deletedRxnsInINIT, taskReport]=...
0100 %               getINITModel(refModel, tissue, celltype, hpaData, arrayData,...
0101 %               metabolomicsData, taskFile, useScoresForTasks, printReport,...
0102 %               taskStructure, params, paramsFT)
0103 
0104 if nargin<3
0105     celltype=[];
0106 else
0107     celltype=char(celltype);
0108 end
0109 if nargin<4
0110     hpaData=[];
0111 end
0112 if nargin<5
0113     arrayData=[];
0114 end
0115 if nargin<6
0116     metabolomicsData=[];
0117 end
0118 if nargin<7
0119     taskFile=[];
0120 else
0121     taskFile=char(taskFile);
0122 end
0123 if nargin<8 || isempty(useScoresForTasks)
0124     useScoresForTasks=true;
0125 end
0126 if nargin<9 || isempty(printReport)
0127     printReport=true;
0128 end
0129 if nargin<10
0130     taskStructure=[];
0131 end
0132 if nargin<11
0133     params=[];
0134 end
0135 if nargin<12
0136     paramsFT=[];
0137 end
0138 
0139 %Check that the model is in the closed form
0140 if ~isfield(refModel,'unconstrained')
0141     EM='Exchange metabolites should normally not be removed from the model when using getINITModel. Use importModel(file,false) to import a model with exchange metabolites remaining (see the documentation for details)';
0142     dispEM(EM);
0143 end
0144 
0145 %Create the task structure if not supplied
0146 if any(taskFile) && isempty(taskStructure)
0147     taskStructure=parseTaskList(taskFile);
0148 end
0149 
0150 
0151 % sc-tINIT to identify confidence levels of gene expression
0152 if ~isempty(arrayData) && isfield(arrayData,'singleCells')
0153     if arrayData.singleCells == 1
0154         % Check to ensure cell type is defined
0155         if ~isfield(arrayData,'celltypes')
0156             dispEM('arrayData must contain cell type information if sc-tINIT is to be used','false');   
0157         end
0158         if ~ismember(upper(celltype),upper(arrayData.celltypes))
0159             dispEM('The cell type name does not match');   
0160         end
0161         
0162         % Analyze only cell type of interest
0163         J= strcmpi(arrayData.celltypes,celltype);
0164         
0165         % Analyze only genes included in the reference model
0166         I=ismember(arrayData.genes,refModel.genes);
0167         
0168         % Convert expression data to population fractions
0169         binary_levels = arrayData.levels(I,J)~=0; % Classify each gene as detected (1) or not (0) in each cell
0170         cell_count_levels = sum(binary_levels,2); % Number of cells expressing each transcript
0171         cell_frac_levels = cell_count_levels/size(binary_levels,2); % Number of cells expressing each transcript
0172 
0173         % Bin cell_frac_counts manually
0174         x = 0:.01:1;
0175         for(i = 1:length(x))
0176             cell_frac_count(i) = sum(round(cell_frac_levels,2)==x(i));
0177         end
0178         
0179         % Fit four beta distributions
0180         cell_frac_count(cell_frac_count==0) = NaN; % Remove zeros from optimization
0181         cell_frac_count(1) = NaN; % Remove non-expressed genes from optimization
0182         x_lim = 1; % Somewhat arbitrary parameter to fit left tail of distr.
0183         myfun = @(par) nansum((cell_frac_count(1:find(x>=x_lim,1)) - ...
0184             abs(par(1))*betapdf(x(1:find(x>=x_lim,1)),abs(par(2)),abs(par(3))) - ...
0185             abs(par(4))*betapdf(x(1:find(x>=x_lim,1)),abs(par(5)),abs(par(6))) - ...
0186             abs(par(7))*betapdf(x(1:find(x>=x_lim,1)),abs(par(8)),abs(par(9))) - ...
0187             abs(par(10))*betapdf(x(1:find(x>=x_lim,1)),abs(par(11)),abs(par(12)))).^2);
0188         
0189         par0 = [4,2,100,7,2,30,7,5,20,5,15,20];
0190         opts = optimset('Display','off');
0191         [par,f_val] = fminsearch(myfun,par0,opts);
0192         par = abs(par);
0193         
0194         % Plot results
0195         if (isfield(arrayData,'plotResults'))
0196             if arrayData.plotResults == true
0197                 figure(); hold on; plot(x,cell_frac_count,'ko','MarkerSize',5);
0198                 plot(x,abs(par(1))*betapdf(x,abs(par(2)),abs(par(3))),'b-','LineWidth',1)
0199                 plot(x,abs(par(4))*betapdf(x,abs(par(5)),abs(par(6))),'b-','LineWidth',1)
0200                 plot(x,abs(par(7))*betapdf(x,abs(par(8)),abs(par(9))),'b-','LineWidth',1)
0201                 plot(x,abs(par(10))*betapdf(x,abs(par(11)),abs(par(12))),'b-','LineWidth',1)
0202                 plot(x,abs(par(1))*betapdf(x,abs(par(2)),abs(par(3))) + ...
0203                     abs(par(4))*betapdf(x,abs(par(5)),abs(par(6))) + ...
0204                     abs(par(7))*betapdf(x,abs(par(8)),abs(par(9))) + ...
0205                     abs(par(10))*betapdf(x,abs(par(11)),abs(par(12))),'-','Color',[.5 .5 .5],'LineWidth',2)
0206                 xlabel('Expression Probability');ylabel('# of genes');set(gca,'FontSize',14,'LineWidth',1.25);
0207                 title('Expression prediction','FontSize',18,'FontWeight','bold')
0208             end
0209         end
0210         
0211         % Score genes based on population expression (p = .05)
0212         exprs_cutoff_1 = find(cdf('beta',x,par(2),par(3)) >.95,1)-1; % Find index of no confidence genes
0213         exprs_cutoff_2 = find(cdf('beta',x,par(5),par(6)) >.95,1)-1; % Find index of low confidence genes
0214         exprs_cutoff_3 = find(cdf('beta',x,par(8),par(9)) >.95,1)-1; % Find index of low confidence genes
0215         exprs_cutoffs = sort([exprs_cutoff_1,exprs_cutoff_2,exprs_cutoff_3]);
0216         gene_scores = cell_frac_levels*0;
0217         gene_scores(cell_frac_levels <= x(exprs_cutoffs(1))) = 4; % Not detected
0218         gene_scores(logical((cell_frac_levels >= x(exprs_cutoffs(1))).*(cell_frac_levels < x(exprs_cutoffs(2))))) = 3; % Low detection
0219         gene_scores(logical((cell_frac_levels >= x(exprs_cutoffs(2))).*(cell_frac_levels < x(exprs_cutoffs(3))))) = 2; % Medium detection
0220         gene_scores(cell_frac_levels > x(exprs_cutoffs(3))) = 1; % High detection
0221 
0222         % Replace hpaData with singleCellData
0223         if printReport==true
0224             dispEM('Single cell data is not currently compatible with HPA data. \n         Replacing hpaData with single cell-based scoring.',false);
0225         end
0226         hpaData.genes = arrayData.genes;
0227         hpaData.tissues = arrayData.tissues;
0228         hpaData.celltypes = arrayData.celltypes;
0229         hpaData.levels = [{'High'},{'Medium'},{'Low'},{'None'}];
0230         hpaData.gene2Level = zeros(length(arrayData.genes),length(arrayData.celltypes));
0231         for i = 1:length(find(J))
0232             find_var = find(J,i);
0233             hpaData.gene2Level(I,find_var(end)) = gene_scores;
0234         end
0235         
0236         % Remove arrayData from the analysis (Might be a bad idea)
0237         clear arrayData
0238         arrayData=[];
0239     end
0240 end
0241 
0242 
0243 if printReport==true
0244     if any(celltype)
0245         fprintf(['***Generating model for: ' tissue ' - ' celltype '\n']);
0246     else
0247         fprintf(['***Generating model for: ' tissue '\n']);
0248     end
0249     if ~isempty(hpaData)
0250         fprintf('-Using HPA data\n');
0251     end
0252     if ~isempty(arrayData)
0253         fprintf('-Using array data\n');
0254     end
0255     if ~isempty(metabolomicsData)
0256         fprintf('-Using metabolomics data\n');
0257     end
0258     if ~isempty(taskFile) || ~isempty(taskStructure)
0259         fprintf('-Using metabolic tasks\n');
0260     end
0261     fprintf('\n');
0262     
0263     printScores(refModel,'Reference model statistics',hpaData,arrayData,tissue,celltype);
0264 end
0265 
0266 %Remove dead-end reactions to speed up the optimization and to
0267 %differentiate between reactions removed by INIT and those that are
0268 %dead-end
0269 [~, deletedDeadEndRxns]=simplifyModel(refModel,true,false,true,true,true);
0270 cModel=removeReactions(refModel,deletedDeadEndRxns,false,true);
0271 
0272 %Store the connected model like this to keep track of stuff
0273 if printReport==true
0274     printScores(cModel,'Pruned model statistics',hpaData,arrayData,tissue,celltype);
0275 end
0276 
0277 %If tasks have been defined, then go through them and get essential
0278 %reactions
0279 if ~isempty(taskStructure)
0280     [taskReport, essentialRxnMat]=checkTasks(cModel,[],printReport,true,true,taskStructure);
0281     
0282     essentialRxnsForTasks=cModel.rxns(any(essentialRxnMat,2));
0283     
0284     %Remove tasks that cannot be performed
0285     taskStructure(taskReport.ok==false)=[];
0286     if printReport==true
0287         printScores(removeReactions(cModel,setdiff(cModel.rxns,essentialRxnsForTasks),true,true),'Reactions essential for tasks',hpaData,arrayData,tissue,celltype);
0288     end
0289 else
0290     essentialRxnsForTasks={};
0291 end
0292 
0293 %Score the connected model
0294 [rxnScores, geneScores]=scoreModel(cModel,hpaData,arrayData,tissue,celltype);
0295 
0296 %Run the INIT algorithm. The exchange reactions that are used in the final
0297 %reactions will be open, which doesn't fit with the last step. Therefore
0298 %delete reactions from the original model instead of taking the output. The
0299 %default implementation does not constrain reversible reactions to only
0300 %carry flux in one direction. Runs without the constraints on reversibility
0301 %and with all output allowed. This is to reduce the complexity of the
0302 %problem.
0303 [~, deletedRxnsInINIT, metProduction]=runINIT(simplifyModel(cModel),rxnScores,metabolomicsData,essentialRxnsForTasks,0,true,false,params);
0304 initModel=removeReactions(cModel,deletedRxnsInINIT,true,true);
0305 if printReport==true
0306     printScores(initModel,'INIT model statistics',hpaData,arrayData,tissue,celltype);
0307     printScores(removeReactions(cModel,setdiff(cModel.rxns,deletedRxnsInINIT),true,true),'Reactions deleted by INIT',hpaData,arrayData,tissue,celltype);
0308 end
0309 
0310 %The full model has exchange reactions in it. fitTasks calls on fillGaps,
0311 %which automatically removes exchange metabolites (because it assumes that
0312 %the reactions are constrained when appropriate). In this case the
0313 %uptakes/outputs are retrieved from the task sheet instead. To prevent
0314 %exchange reactions being used to fill gaps, they are delete from the
0315 %reference model here.
0316 initModel.id='INITModel';
0317 
0318 %If gaps in the model should be filled using a task list
0319 if ~isempty(taskStructure)
0320     %Remove exchange reactions and reactions already included in the INIT
0321     %model
0322     refModelNoExc=removeReactions(refModel,union(initModel.rxns,getExchangeRxns(refModel)),true,true);
0323     
0324     %At this stage the model is fully connected and most of the genes with
0325     %good scores should have been included. The final gap-filling should
0326     %take the scores of the genes into account, so that "rather bad"
0327     %reactions are preferred to "very bad" reactions. However, reactions
0328     %with positive scores will be included even if they are not connected
0329     %in the current formulation. Therefore, such reactions will have to be
0330     %assigned a small negative score instead.
0331     if useScoresForTasks==true
0332         refRxnScores=scoreModel(refModelNoExc,hpaData,arrayData,tissue,celltype);
0333         [outModel, addedRxnMat]=fitTasks(initModel,refModelNoExc,[],true,min(refRxnScores,-0.1),taskStructure,paramsFT);
0334     else
0335         [outModel, addedRxnMat]=fitTasks(initModel,refModelNoExc,[],true,[],taskStructure,paramsFT);
0336     end
0337     if printReport==true
0338         printScores(outModel,'Functional model statistics',hpaData,arrayData,tissue,celltype);
0339         printScores(removeReactions(outModel,intersect(outModel.rxns,initModel.rxns),true,true),'Reactions added to perform the tasks',hpaData,arrayData,tissue,celltype);
0340     end
0341     
0342     addedRxnsForTasks=refModelNoExc.rxns(any(addedRxnMat,2));
0343 else
0344     outModel=initModel;
0345     addedRxnMat=[];
0346     addedRxnsForTasks={};
0347 end
0348 
0349 %The model can now perform all the tasks defined in the task list. The
0350 %algorithm cannot deal with gene-complexes at the moment. It is therefore
0351 %ok to remove bad genes from a reaction (as long as at least one gene is
0352 %kept)
0353 model=outModel;
0354 
0355 [~, I]=ismember(model.genes,cModel.genes); %All should be found
0356 %This is a little weird way to make sure that only one bad gene is included
0357 %if there are no good ones (since all -Inf==max(-Inf))
0358 geneScores(isinf(geneScores))=-1000+rand(sum(isinf(geneScores)),1);
0359 
0360 model.grRules(:)={''};
0361 for i=1:numel(model.rxns)
0362     ids=find(model.rxnGeneMat(i,:));
0363     if numel(ids)>1
0364         scores=geneScores(I(ids));
0365         %Only keep the positive ones if possible
0366         model.rxnGeneMat(i,ids(~(scores>0 | scores==max(scores))))=0;
0367     end
0368     %Rewrite the grRules to be only OR
0369     if isfield(model,'grRules')
0370         J=find(model.rxnGeneMat(i,:));
0371         for j=1:numel(J)
0372             model.grRules{i}=[model.grRules{i} '(' model.genes{J(j)} ')'];
0373             if j<numel(J)
0374                 model.grRules{i}=[model.grRules{i} ' or '];
0375             end
0376         end
0377     end
0378 end
0379 
0380 %Find all genes that are not used and delete them
0381 I=sum(model.rxnGeneMat)==0;
0382 model.genes(I)=[];
0383 model.rxnGeneMat(:,I)=[];
0384 if isfield(model,'geneShortNames')
0385     model.geneShortNames(I)=[];
0386 end
0387 if isfield(model,'proteins')
0388     model.proteins(I)=[];
0389 end
0390 if isfield(model,'geneMiriams')
0391     model.geneMiriams(I)=[];
0392 end
0393 if isfield(model,'geneFrom')
0394     model.geneFrom(I)=[];
0395 end
0396 if isfield(model,'geneComps')
0397     model.geneComps(I)=[];
0398 end
0399 
0400 %At this stage the model will contain some exchange reactions but probably
0401 %not all (and maybe zero). This can be inconvenient, so all exchange
0402 %reactions from the reference model are added, except for those which
0403 %involve metabolites that are not in the model.
0404 
0405 %First delete and included exchange reactions in order to prevent the order
0406 %from changing
0407 model=removeReactions(model,getExchangeRxns(model));
0408 
0409 %Create a model with only the exchange reactions in refModel
0410 excModel=removeReactions(refModel,setdiff(refModel.rxns,getExchangeRxns(refModel)),true,true);
0411 
0412 %Find the metabolites there which are not exchange metabolites and which do
0413 %not exist in the output model
0414 I=~ismember(excModel.mets,model.mets) & excModel.unconstrained==0;
0415 
0416 %Then find those reactions and delete them
0417 [~, J]=find(excModel.S(I,:));
0418 excModel=removeReactions(excModel,J,true,true);
0419 
0420 %Merge with the output model
0421 model=mergeModels({model;excModel},'metNames');
0422 model.id='INITModel';
0423 model.name=['Automatically generated model for ' tissue];
0424 if any(celltype)
0425     model.name=[model.name ' - ' celltype];
0426 end
0427 
0428 if printReport==true
0429     printScores(model,'Final model statistics',hpaData,arrayData,tissue,celltype);
0430 end
0431 
0432 %Add information about essential reactions and reactions included for
0433 %gap-filling and return a taskReport
0434 if ~isempty(taskStructure)
0435     I=find(taskReport.ok); %Ignore failed tasks
0436     for i=1:numel(I)
0437         taskReport.essential{I(i),1}=cModel.rxns(essentialRxnMat(:,I(i)));
0438         taskReport.gapfill{I(i),1}=refModelNoExc.rxns(addedRxnMat(:,i));
0439     end
0440 else
0441     taskReport=[];
0442 end
0443 
0444 %Fix grRules and reconstruct rxnGeneMat
0445 [grRules,rxnGeneMat] = standardizeGrRules(model,true);
0446 model.grRules = grRules;
0447 model.rxnGeneMat = rxnGeneMat;
0448 end
0449 
0450 %This is for printing a summary of a model
0451 function [rxnS, geneS]=printScores(model,name,hpaData,arrayData,tissue,celltype)
0452 [a, b]=scoreModel(model,hpaData,arrayData,tissue,celltype);
0453 rxnS=mean(a);
0454 geneS=mean(b(~isinf(b)));
0455 fprintf([name ':\n']);
0456 fprintf(['\t' num2str(numel(model.rxns)) ' reactions, ' num2str(numel(model.genes)) ' genes\n']);
0457 fprintf(['\tMean reaction score: ' num2str(rxnS) '\n']);
0458 fprintf(['\tMean gene score: ' num2str(geneS) '\n']);
0459 fprintf(['\tReactions with positive scores: ' num2str(100*sum(a>0)/numel(a)) '%%\n\n']);
0460 end

Generated by m2html © 2005