Description of ftINIT

0001 function [model, metProduction, addedRxnsForTasks, deletedRxnsInINIT, fullMipRes] = ftINIT(prepData, tissue, celltype, hpaData, transcrData, metabolomicsData, INITSteps, removeGenes, useScoresForTasks, paramsFT, verbose)
0002 % ftINIT
0003 %   Main function for generates a model using the ftINIT algorithm, based
0004 %   on proteomics and/or transcriptomics and/or metabolomics and/or metabolic
0005 %   tasks. The algorithm is not designed for running with metabolomics only.
0006 %   The function prepINITModel needs to be run first for the template
0007 %   model (such as the generic Human-GEM), but only need to be run once. This
0008 %   function precalculates things independent of the omics to speed up the model
0009 %   generation process, and outputs the prepData, which is input to this function.
0010 %
0011 %   prepData            The prepdata for the model.
0012 %   tissue              tissue to score for. Should exist in either
0013 %                       hpaData.tissues or transcrData.tissues
0014 %   celltype            cell type to score for. Should exist in either
0015 %                       hpaData.celltypes or transcrData.celltypes for this
0016 %                       tissue (optional, default is to use the max values
0017 %                       among all the cell types for the tissue.
0018 %   hpaData             HPA data structure from parseHPA (optional if transcrData
0019 %                       is supplied, default [])
0020 %   transcrData         gene expression data structure (optional if hpaData is
0021 %                       supplied, default []). Used to be called arrayData.
0022 %       genes           cell array with the unique gene names
0023 %       tissues         cell array with the tissue names. The list may not
0024 %                       be unique, as there can be multiple cell types per
0025 %                       tissue.
0026 %       celltypes       cell array with the cell type names for each tissue
0027 %       levels          GENESxTISSUES array with the expression level for
0028 %                       each gene in each tissue/celltype. NaN should be
0029 %                       used when no measurement was performed
0030 %       threshold       a single value or a vector of gene expression
0031 %                       thresholds, above which genes are considered to be
0032 %                       "expressed". default = 1(optional, by default, the mean expression
0033 %                       levels of each gene across all tissues in transcrData
0034 %                       will be used as the threshold values)
0035 %       singleCells     binary value selecting whether to use the
0036 %                       single-cell algorithm to identify expressed genes.
0037 %                       If used, specify cell subpopulations in CELLTYPES
0038 %                       (optional, default [])
0039 %       plotResults     true if single cell probability distributions
0040 %                       should be plotted (optional, default = false)
0041 %   metabolomicsData    cell array with metabolite names that the model
0042 %                       should produce (optional, default [])
0043 %   INITSteps           Specifies the steps in the algorithm. For more info,
0044 %                       see INITStepDesc and getINITSteps.
0045 %                       (optional, default getINITSteps(), which is the standard ftINIT).
0046 %   removeGenes         if true, low-abundance genes will be removed from
0047 %                       grRules, unless they are the only gene associated
0048 %                       with a reaction, or a subunit of an enzyme complex
0049 %                       (see "removeLowScoreGenes" function for details).
0050 %                       If false, grRules will not be modified; however,
0051 %                       genes that were associated only with removed
0052 %                       reactions will not be present in the final model.
0053 %                       (optional, default true).
0054 %   useScoresForTasks   true if the calculated reaction scored should be
0055 %                       used as weights when fitting to tasks (optional, default
0056 %                       true)
0057 %   paramsFT            *obsolete option*
0058 %   verbose             if true, the MILP progression will be shown.
0059 %                       (optional, default false)
0060 %
0061 %   model                   the resulting model structure
0062 %   metProduction           array that indicates which of the
0063 %                           metabolites in metabolomicsData that could be
0064 %                           produced. Note that this is before the
0065 %                           gap-filling process to enable defined tasks. To
0066 %                           see which metabolites that can be produced in
0067 %                           the final model, use canProduce.
0068 %                           -2: metabolite name not found in model
0069 %                           -1: metabolite found, but it could not be produced
0070 %                           1: metabolite could be produced
0071 %   addedRxnsForTasks       cell array of the reactions which were added in
0072 %                           order to perform the tasks
0073 %   deletedRxnsInINIT       cell array of reactions deleted because they
0074 %                           could not carry flux (INIT requires a
0075 %                           functional input model)
0076 %   fullMipRes              The solver results from the last MILP step run
0077 %
0078 %   This is the main function for automatic reconstruction of models based
0079 %   on the ftINIT algorithm ().
0080 %
0081 %   NOTE: Exchange metabolites should normally not be removed from the model
0082 %   when using this approach, since checkTasks/fitTasks rely on putting specific
0083 %   constraints for each task. The INIT algorithm will remove exchange metabolites
0084 %   if any are present. Use importModel(file,false) to import a model with
0085 %   exchange metabolites remaining.
0086 %
0087 % Usage: [model, metProduction, addedRxnsForTasks, deletedRxnsInINIT, ...
0088 %               fullMipRes] = ...
0089 %               ftINIT(prepData, tissue, celltype, hpaData, transcrData, ...
0090 %               metabolomicsData, INITSteps, removeGenes, useScoresForTasks, ...
0091 %               paramsFT);
0092 %
0093 
0094 
0095 if nargin < 5
0096     transcrData = [];
0097 end
0098 if nargin < 6
0099     metabolomicsData = [];
0100 end
0101 if nargin < 7 || isempty(INITSteps)
0102     INITSteps = getINITSteps([],'1+1');
0103 end
0104 if nargin < 8 || isempty(removeGenes)
0105     removeGenes = true;
0106 end
0107 if nargin < 9 || isempty(useScoresForTasks)
0108     useScoresForTasks = true;
0109 end
0110 if nargin < 10
0111     paramsFT = [];
0112 end
0113 
0114 if nargin < 11
0115     verbose = false;
0116 end
0117 %Handle detected mets:
0118 %Previously, this was handled by giving a bonus for secreting those metabolites,
0119 %but that doesn't work since the metabolite secretion and uptake can be lost when
0120 %we merge linearly dependent reactions.
0121 %Instead, we need to figure out which reactions either produce or take up the mets.
0122 %We then give a bonus if any of them carry flux.
0123 %To simplify things, we focus on reactions that produce the metabolite (since there must be one such reaction).
0124 %It is still a bit complicated though. In this step, we focus on identifying
0125 %producer reactions. We further reason that the direction doesn't matter -
0126 %we can force one of these reactions in any direction - if it becomes a consumer, it will
0127 %automatically force another producer on as well (otherwise we'll have a net consumption).
0128 
0129 if (~isempty(metabolomicsData))
0130     if length(unique(upper(metabolomicsData))) ~= length(metabolomicsData)
0131         dispEM('Metabolomics contains the same metabolite multiple times');
0132     end
0133     metData = false(numel(metabolomicsData), length(prepData.minModel.rxns)); %one row per metabolite that is a boolean vector
0134     for i=1:numel(metabolomicsData)
0135         %Get the matching mets
0136         metSel = ismember(upper(prepData.refModel.metNames),upper(metabolomicsData{i}));
0137         prodRxnsSel = any(prepData.refModel.S(metSel,:) > 0,1) | ... %direct producers
0138                      (any(prepData.refModel.S(metSel,:) < 0,1) & prepData.refModel.rev.'); %reversible reactions that are consumers
0139         %convert the production rxns from refModel to minModel
0140         prepData.groupIds;
0141         [~,ia,ib] = intersect(prepData.minModel.rxns,prepData.refModel.rxns);
0142         grpIdsMerged = nan(length(prepData.minModel.rxns),1);
0143         grpIdsMerged(ia) = prepData.groupIds(ib);
0144         
0145         groupIdsPos = unique(prepData.groupIds(prodRxnsSel));%gets all group ids which includes a production rxn
0146         groupIdsPos = groupIdsPos(groupIdsPos ~= 0);%remove the 0 id, it means there is no group
0147         %the other option is that there is a direct match between the rxn id in minModel and refModel:
0148         posRxns = prepData.refModel.rxns(prodRxnsSel);
0149         directMatch = ismember(prepData.minModel.rxns, posRxns).';
0150         
0151         metData(i,:) = ismember(grpIdsMerged, groupIdsPos).' | directMatch;
0152     end
0153     metData = sparse(metData);
0154 else
0155     metData = [];
0156 end
0157 
0158 % Get rxn scores and adapt them to the minimized model
0159 origRxnScores = scoreComplexModel(prepData.refModel,hpaData,transcrData,tissue,celltype);
0160 origRxnScores(origRxnScores > -0.1 & origRxnScores <= 0) = -0.1;%we don't want reaction scores that are exactly 0 (or close), this causes problems in the milp
0161 origRxnScores(origRxnScores < 0.1 & origRxnScores > 0) = 0.1;
0162 
0163 rxnsTurnedOn = false(length(prepData.minModel.rxns),1);
0164 fluxes = zeros(length(prepData.minModel.rxns),1);
0165 
0166 rxnsToIgnoreLastStep = [1;1;1;1;1;1;1;1];
0167 
0168 %We assume that all essential rxns are irrev - this is taken care of in
0169 %prepINITModel. We then use an initial flux "from last run" of 0.1 for all
0170 %reactions. This is used for knowing what flux should be forced through an
0171 %essential rxn.
0172 fluxes = ones(length(prepData.minModel.rxns), 1).*0.1;
0173 
0174 for initStep = 1:length(INITSteps)
0175     disp(['ftINIT: Running step ' num2str(initStep)])
0176     stp = INITSteps{initStep};
0177     
0178     if any ((rxnsToIgnoreLastStep - stp.RxnsToIgnoreMask) < 0)
0179         dispEM('RxnsToIgnoreMask may not cover rxns not covered in previous steps, but the other way around is fine.');
0180     end
0181     rxnsToIgnoreLastStep = stp.RxnsToIgnoreMask;
0182     
0183     mm = prepData.minModel;
0184     
0185     if (~isempty(stp.MetsToIgnore))
0186         if (~isempty(stp.MetsToIgnore.simpleMets))
0187             %Here, we remove simple metabolites that will not really affect the milp but
0188             %are very common in the S matrix. For example H2O, H+, etc.
0189             %It is also possible to leave compartments untouched, for example the i compartment in the mitochondria (for H+).
0190             metsToRem = ismember(mm.metNames,stp.MetsToIgnore.simpleMets.mets);
0191             compsToKeep = find(ismember(mm.comps, stp.MetsToIgnore.simpleMets.compsToKeep));
0192             metsToRem = metsToRem & ~ismember(mm.metComps, compsToKeep);
0193             mm.S(metsToRem,:) = 0;
0194         end    
0195     end
0196 
0197     %Set up the reaction scores and essential rxns
0198     rxnsToIgnore = getRxnsFromPattern(stp.RxnsToIgnoreMask, prepData);
0199     rxnScores = groupRxnScores(prepData.minModel, origRxnScores, prepData.refModel.rxns, prepData.groupIds, rxnsToIgnore);
0200 
0201     essentialRxns = prepData.essentialRxns;
0202     toRev = false(numel(mm.rxns),1);
0203     %Handle the results from previous steps ('ignore', 'exclude', 'essential')
0204     if strcmp(stp.HowToUsePrevResults, 'exclude')
0205         rxnScores(rxnsTurnedOn) = 0; %This is not used anymore in any step setup.
0206     elseif strcmp(stp.HowToUsePrevResults, 'essential')
0207         %Make all reversible reactions turned on in previous steps reversible
0208         %in the direction that they were previously carrying flux
0209         
0210         %first reverse the reactions that need to be reversed
0211         rev = mm.rev == 1;
0212         toRev = rxnsTurnedOn & rev & fluxes < 0;
0213         mm = reverseRxns(mm, mm.rxns(toRev));
0214         
0215         %Then make them irreversible
0216         mm.rev(rxnsTurnedOn) = 0;
0217         mm.lb(rxnsTurnedOn) = 0;
0218 
0219         essentialRxns = unique([prepData.essentialRxns;mm.rxns(rxnsTurnedOn)]);
0220     end
0221 
0222     
0223     mipGap = 1;
0224     first = true;
0225     success = false;
0226     fullMipRes = [];
0227     for rn = 1:length(stp.MILPParams)
0228         params = stp.MILPParams{rn};
0229         if ~isfield(params, 'MIPGap')
0230             params.MIPGap = 0.0004;
0231         end
0232         
0233         if ~isfield(params, 'TimeLimit')
0234             params.TimeLimit = 5000;
0235         end
0236         
0237         if ~first 
0238             %There is sometimes a problem with that the objective function becomes close to zero,
0239             %which leads to that a small percentage of that (which is the MIPGap sent in) is very small
0240             %and the MILP hence takes a lot of time to finish. We also therefore use an absolute MIP gap,
0241             %converted to a percentage using the last value of the objective function.
0242             params.MIPGap = min(max(params.MIPGap, stp.AbsMIPGaps{rn}/abs(lastObjVal)),1);
0243             params.seed = 1234;%use another seed, may work better
0244 
0245             if mipGap <= params.MIPGap
0246                 success = true;
0247                 break; %we're done - this will not happen the first time
0248             else
0249                 disp(['MipGap too high, trying with a different run. MipGap = ' num2str(mipGap) ' New MipGap Limit = ' num2str(params.MIPGap)])
0250             end
0251         end
0252         
0253         first = false;
0254         
0255         %now run the MILP
0256         try
0257             %The prodweight for metabolomics is currently set to 5 - 0.5 was default in the old version, which I deemed very small?
0258             %There could be a need to specify this somewhere in the call at some point.
0259             %This value has not been evaluated, but is assumed in the test cases - if changed, update the test case
0260             startVals = [];
0261             if ~isempty(fullMipRes)
0262                 startVals = fullMipRes.full;
0263             end
0264             [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, startVals, fluxes, verbose);
0265             %This is a bit tricky - since we reversed some reactions, those fluxes also need to be reversed
0266             fluxes1(toRev) = -fluxes1(toRev);
0267             
0268             mipGap = fullMipRes.mipgap;
0269             lastObjVal = fullMipRes.obj;
0270         catch e
0271             mipGap = Inf;
0272             lastObjVal = Inf; %we need to set something here, Inf leads to that this doesn't come into play
0273         end
0274         
0275         success = mipGap <= params.MIPGap;
0276     end
0277     
0278     if ~success
0279         dispEM(['Failed to find good enough solution within the time frame. MIPGap: ' num2str(mipGap)]);
0280     end
0281     
0282     %save the reactions turned on and their fluxes for the next step
0283     rxnsTurnedOn = rxnsTurnedOn | rxnsTurnedOn1.';
0284     %The fluxes are a bit tricky - what if they change direction between the steps?
0285     %The fluxes are used to determine the direction in which reactions are forced on
0286     %(to simplify the problem it is good if they are unidirectional).
0287     %We use the following strategy:
0288     %1. Use the fluxes from the most recent step.
0289     %2. If any flux is very low there (i.e. basically zero), use the flux from the previous steps
0290     %This could in theory cause problems, but seems to work well practically
0291     fluxesOld = fluxes;
0292     fluxes = fluxes1;
0293     %make sure that all reactions that are on actually has a flux - otherwise
0294     %things could go bad, since the flux will be set to essential in a random direction
0295     %This sometimes happens for rxns with negative score - let's just accept that.
0296     %if (sum(abs(fluxes1) < 10^-7 & rxnsTurnedOn))
0297     %    dispEM('There are rxns turned on without flux - this might cause problems');
0298     %end
0299     %fluxes(abs(fluxes1) < 10^-7) = fluxesOld(abs(fluxes1) < 10^-9);
0300 end
0301 
0302 
0303 %get the essential rxns
0304 essential = ismember(prepData.minModel.rxns,prepData.essentialRxns);
0305 %So, we only add reactions where the linearly merged scores are zero for all linearly dependent reactions
0306 % (this cannot happen by chance, taken care of in the function groupRxnScores)
0307 rxnsToIgn = rxnScores == 0; 
0308 deletedRxnsInINITSel = ~(rxnsTurnedOn | rxnsToIgn | essential);
0309 deletedRxnsInINIT = prepData.minModel.rxns(deletedRxnsInINITSel);
0310 
0311 %Here we need to figure out which original reactions (before the linear merge)
0312 %that were removed. These are all reactions with the same group ids as the removed reactions
0313 groupIdsRemoved = prepData.groupIds(ismember(prepData.refModel.rxns, deletedRxnsInINIT)); %can improve this slightly, use sel above
0314 groupIdsRemoved = groupIdsRemoved(groupIdsRemoved ~= 0);%zero means that the reaction was not grouped, all with zeros are not a group!
0315 rxnsToRem = union(prepData.refModel.rxns(ismember(prepData.groupIds,groupIdsRemoved)), deletedRxnsInINIT);%make a union here to include the ungrouped (unmerged) as well
0316 
0317 initModel = removeReactions(prepData.refModel,rxnsToRem,false,true);
0318 
0319 % remove metabolites separately to avoid removing those needed for tasks
0320 unusedMets = initModel.mets(all(initModel.S == 0,2));
0321 initModel = removeMets(initModel, setdiff(unusedMets, prepData.essentialMetsForTasks));
0322 
0323 %if printReport == true
0324 %    printScores(initModel,'INIT model statistics',hpaData,transcrData,tissue,celltype);
0325 %    printScores(removeReactions(cModel,setdiff(cModel.rxns,rxnsToRem),true,true),'Reactions deleted by INIT',hpaData,transcrData,tissue,celltype);
0326 %end
0327 
0328 %The full model has exchange reactions in it. ftINITFillGapsForAllTasks calls
0329 %ftINITFillGaps, which automatically removes exchange metabolites (because it
0330 %assumes that the reactions are constrained when appropriate). In this case the
0331 %uptakes/outputs are retrieved from the task sheet instead. To prevent
0332 %exchange reactions being used to fill gaps, they are deleted from the
0333 %reference model here.
0334 initModel.id = 'INITModel';
0335 
0336 %If gaps in the model should be filled using a task list
0337 if ~isempty(prepData.taskStruct)
0338     %Remove exchange reactions and reactions already included in the INIT
0339     %model
0340     %We changed strategy and instead include all rxns except the exchange rxns in the ref model
0341     %But we do keep the exchange rxns that are essential.
0342     %Let's test to remove all, that should work
0343     
0344     %At this stage the model is fully connected and most of the genes with
0345     %good scores should have been included. The final gap-filling should
0346     %take the scores of the genes into account, so that "rather bad"
0347     %reactions are preferred to "very bad" reactions. However, reactions
0348     %with positive scores will be included even if they are not connected
0349     %in the current formulation. Therefore, such reactions will have to be
0350     %assigned a small negative score instead.
0351     exchRxns = getExchangeRxns(prepData.refModel);
0352     refModelNoExc = removeReactions(prepData.refModelWithBM,exchRxns,false,true);
0353     exchRxns = getExchangeRxns(initModel);
0354     initModelNoExc = removeReactions(closeModel(initModel),exchRxns,false,true);
0355     
0356     if useScoresForTasks == true
0357         %map the rxn scores to the model without exchange rxns
0358         [~,ia,ib] = intersect(refModelNoExc.rxns,prepData.refModel.rxns);
0359         rxnScores2nd = NaN(length(refModelNoExc.rxns),1);
0360         rxnScores2nd(ia) = origRxnScores(ib);
0361         %all(rxnScores2nd == refRxnScores);%should be the same, ok!
0362         [outModel,addedRxnMat] = ftINITFillGapsForAllTasks(initModelNoExc,refModelNoExc,[],true,min(rxnScores2nd,-0.1),prepData.taskStruct,paramsFT,verbose);
0363     else
0364         [outModel,addedRxnMat] = ftINITFillGapsForAllTasks(initModelNoExc,refModelNoExc,[],true,[],prepData.taskStruct,paramsFT,verbose);
0365     end
0366     %if printReport == true
0367     %    printScores(outModel,'Functional model statistics',hpaData,transcrData,tissue,celltype);
0368     %    printScores(removeReactions(outModel,intersect(outModel.rxns,initModel.rxns),true,true),'Reactions added to perform the tasks',hpaData,transcrData,tissue,celltype);
0369     %end
0370     
0371     addedRxnsForTasks = refModelNoExc.rxns(any(addedRxnMat,2));
0372 else
0373     outModel = initModel;
0374     addedRxnMat = [];
0375     addedRxnsForTasks = {};
0376 end
0377 
0378 % The model can now perform all the tasks defined in the task list.
0379 model = outModel;
0380 
0381 
0382 % At this stage the model will contain some exchange reactions but probably
0383 % not all (and maybe zero). This can be inconvenient, so all exchange
0384 % reactions from the reference model are added, except for those which
0385 % involve metabolites that are not in the model.
0386 
0387 %Start from the original model, and just remove the reactions that are no longer there (and keep exchange rxns). The model we got out
0388 %from the problem is not complete, it doesn't have GRPs etc.
0389 %The logic below is a bit complicated. We identify the reactions that should be removed from the full model as
0390 %reactions that have been removed in the init model except the ones that were added back. In addition, we make
0391 %sure that no exchange rxns are removed - they can be removed in the init model if they were linearly merged with other
0392 %reactions that were decided to be removed from the model. We want to keep all exchange rxns to make sure the tasks can
0393 %be performed also without manipulating the b vector in the model (which is what is done in the gap-filling).
0394 exchRxns = getExchangeRxns(prepData.refModel);
0395 deletedRxnsInINIT = setdiff(prepData.refModel.rxns,union(union(initModel.rxns, addedRxnsForTasks), exchRxns));
0396 outModel = removeReactions(prepData.refModel, deletedRxnsInINIT, true); %we skip removing the genes for now, I'm not sure it is desirable
0397 
0398 % If requested, attempt to remove negative-score genes from the model,
0399 % depending on their role (isozyme or complex subunit) in each grRule.
0400 % See the "removeLowScoreGenes" function more more details, and to adjust
0401 % any default parameters therein.
0402 if ( removeGenes )
0403     [~, geneScores] = scoreComplexModel(outModel,hpaData,transcrData,tissue,celltype);
0404     outModel = removeLowScoreGenes(outModel,geneScores);
0405 end
0406 
0407 
0408 model = outModel;
0409 
0410 end
0411 
0412 %This is for printing a summary of a model
0413 function [rxnS, geneS] = printScores(model,name,hpaData,transcrData,tissue,celltype)
0414     [a, b] = scoreComplexModel(model,hpaData,transcrData,tissue,celltype);
0415     rxnS = mean(a);
0416     geneS = mean(b,'omitnan');
0417     fprintf([name ':\n']);
0418     fprintf(['\t' num2str(numel(model.rxns)) ' reactions, ' num2str(numel(model.genes)) ' genes\n']);
0419     fprintf(['\tMean reaction score: ' num2str(rxnS) '\n']);
0420     fprintf(['\tMean gene score: ' num2str(geneS) '\n']);
0421     fprintf(['\tReactions with positive scores: ' num2str(100*sum(a>0)/numel(a)) '%%\n\n']);
0422 end
0423 
0424 function rxnsToIgnore = getRxnsFromPattern(rxnsToIgnorePattern, prepData)
0425     rxnsToIgnore = false(length(prepData.toIgnoreExch),1);
0426     if rxnsToIgnorePattern(1) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreExch; end;
0427     if rxnsToIgnorePattern(2) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreImportRxns; end;
0428     if rxnsToIgnorePattern(3) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreSimpleTransp; end;
0429     if rxnsToIgnorePattern(4) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreAdvTransp; end;
0430     if rxnsToIgnorePattern(5) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreSpont; end;
0431     if rxnsToIgnorePattern(6) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreS; end;
0432     if rxnsToIgnorePattern(7) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreCustomRxns; end;
0433     if rxnsToIgnorePattern(8) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreAllWithoutGPRs; end;
0434 end
0435
ftINIT

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
