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)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