INIT / ftINIT#

MATLAB functions in RAVEN/src/INIT of the RAVEN toolbox.

Functions#

Function	Summary
`ftINIT`	Generate a model using the ftINIT algorithm.
`ftINITInternalAlg`	Run the MILP for one step of ftINIT.
`getINITModel`	Generate a model using the original tINIT algorithm.
`getINITSteps`	Define the steps used as input to ftINIT.
`groupRxnScores`	Sums the reaction scores of all reactions that were merged into one by the
`mergeLinear`	Simplifies a model by merging rxns with linear dependencies. It differs
`parseHPA`	Parse a database dump of the Human Protein Atlas (HPA).
`parseHPArna`	Parse a dump of Human Protein Atlas (HPA) RNA-Seq data.
`prepINITModel`	Precompute omics-independent data for ftINIT.
`rescaleModelForINIT`	Rescale a model to avoid large stoichiometric differences.
`runINIT`	Solve the MILP of the original INIT algorithm.
`scoreModel`	Score reactions and genes in a model with complex gene rules.

Reference#

ftINIT#

function [model, metProduction, addedRxnsForTasks, deletedRxnsInINIT, fullMipRes] = ftINIT(prepData, tissue, celltype, hpaData, varargin)

Generate a model using the ftINIT algorithm.

Generates a model from proteomics, transcriptomics, metabolomics and metabolic tasks, in any combination except metabolomics alone. The function prepINITModel must be run first on the template model (such as the generic Human-GEM). It runs once per template and precalculates the parts that are independent of the omics data. Its output, prepData, is the input to this function.

Parameters

prepData : struct: the prepData for the model, as returned by prepINITModel.
tissue : char: tissue to score for. Must exist in hpaData.tissues or transcrData.tissues.
celltype : char: cell type to score for. Must exist in hpaData.celltypes or transcrData.celltypes for this tissue. If empty, the maximum value among all the cell types of the tissue is used.
hpaData : struct: HPA data structure from parseHPA (required if transcrData is not supplied, default []).

Name-Value Arguments

transcrData : struct

gene expression data structure (required if hpaData is not supplied, default []). With fields:

genes : cell array with the unique gene names.
tissues : cell array with the tissue names. The list can contain duplicates, 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 and cell type. NaN marks a gene that was not measured.
threshold : a single value or a vector of gene expression thresholds, above which genes are considered to be "expressed". If omitted, the mean expression level of each gene across all tissues in transcrData is used as the threshold.

metabolomicsData : cell

cell array with metabolite names that the model should produce (default []).

INITSteps : cell

specifies the steps in the algorithm. For more info, see INITStepDesc and getINITSteps (default getINITSteps([],'1+1'), the standard ftINIT).

removeGenes : logical

if true, low-score genes are removed from grRules, unless they are the only gene associated with a reaction, or a subunit of an enzyme complex (see "removeLowScoreGenes" function for details). If false, grRules are not modified; genes that were associated only with removed reactions are not present in the final model (default true).

useScoresForTasks : logical

true if the calculated reaction scores are used as weights when fitting to tasks. Scores above -0.1 are set to -0.1 for this step (default true).

paramsFT : struct

solver parameters for the gap-filling step that fits the tasks, passed to fitTasks as params (default []).

verbose : logical

if true, the MILP progression is shown (default false).

resolveTies : logical

if true, the degenerate MILP of each step is solved to a canonical answer in addition to the solver Seed: the score objective is held at its optimum, then the count of removable ("on") reactions is minimised, then, among the sparsest solutions, their summed reaction-id rank is minimised. This applies to every main-extraction step and to the final gap-filling step (resolveTies of fitTasks); see raven-gecko-parity#104 and resolve_ties in raven-toolbox. In each main-extraction step except the last (ftINITInternalAlg), the fluxes set the direction and magnitude at which a reaction is forced in the next step. They are taken from the least-total-flux solution for that on/off state and, among those, from the solution with the least id-rank-weighted backward flux (a merged reaction is ranked by the smallest original id in its group). The last step skips these flux phases, as no later step reads its fluxes. A step whose primary solve stops at the time limit keeps its primary solution without any of these phases. The count and id-rank phases get half of the step's TimeLimit each. A warning names each phase that does not prove optimal or finds no solution. Measured with raven-toolbox's ftinit, whose phase handling differs from ftINIT, on one cell line (Human-GEM v2.0.0-38, seeds 1234 and 7, time limit 1800 s): two seeds with the default escalation differ by 110 reactions and 17 predicted essential genes, and with resolveTies alone by 41 reactions and 1 gene, at about 3.6 times the run time (about 14 min against 4 min) (default false).

proveAbsGap : double

if set, each step is proven to this fixed absolute MILP gap in a single solve, replacing the relative-gap escalation in stp.MILPParams. The final near-zero-objective run of the escalation accepts an arbitrary incumbent within the gap. A warning is raised when a step reaches the time limit before the proof. The figures below were measured with raven-toolbox's ftinit, whose phase handling differs from ftINIT, on one cell line (Human-GEM v2.0.0-38, time limit 1800 s) with resolveTies on and without gap-filling, in a sweep over four seeds (1234, 7, 1 and 2). The step 1 objective of the escalation ended 4.0 to 12.0 units below the proven optimum; with seeds 1234 and 7 alone it ended 7.0 and 12.0 below. A value of 5.0 stopped 2.0 below the optimum for one seed (seed 2), and values from 2.0 to 0.05 reached it for all four seeds and gave the same models. Every primary solve of the steps (48 solves) ended optimal; the flux phases phase2c-minflux (3 of 4 seeds) and phase2d-direction (4 of 4 seeds) reached the time limit and were adopted as unproven incumbents. With 1.0 and resolveTies the models of two seeds differ by 1 reaction and 0 predicted essential genes, at about 26 times the run time of the default escalation (about 100 min against 4 min); the mean of a whole build was 5766 s at 0.05 and 5641 s at 1.0. The study recommends 1.0 together with resolveTies when the variation of resolveTies alone is too large and the run time is acceptable. Its full-pipeline runs used a time limit of 600 s per step and its determinism builds 1800 s (default [], which uses the escalation).

Returns

model : struct: the resulting model structure.
metProduction : logical: array from the last MILP step that indicates which of the metabolites in metabolomicsData could be produced (true) or not (false). Metabolites with a producing reaction that is essential are not part of the MILP and are not in the array. The array is determined before the gap-filling for the defined tasks. To see which metabolites can be produced in the final model, use canExchange.
addedRxnsForTasks : cell: cell array of the reactions which were added in order to perform the tasks.
deletedRxnsInINIT : cell: cell array of the reactions of the template model that are not in the final model, excluding exchange reactions (these are all kept). It contains the reactions removed by the MILP steps that gap-filling did not add back.
fullMipRes : struct: the solver results from the last MILP step run.

Notes

The template model keeps its exchange metabolites, since checkTasks and fitTasks put specific constraints on them for each task. The INIT algorithm removes exchange metabolites if any are present. Use importModel(file,false) to import a model with exchange metabolites remaining.

ftINIT passes a prodWeight of 5 (the reward for each produced metabolite of metabolomicsData) to ftINITInternalAlg, whose own default is 0.5.

ftINITInternalAlg#

function [deletedRxns,metProduction,res,turnedOnRxns,fluxes]=ftINITInternalAlg(model,rxnScores,metData,essentialRxns,prodWeight,allowExcretion,remPosRev,params,startVals,fluxes,verbose,forceOn,resolveTies,rankKeys,pinFluxes)

Run the MILP for one step of ftINIT.

  model           a reference model structure
  rxnScores       a vector of scores for the reactions in the model.
                  Positive scores are reactions to keep and negative
                  scores are reactions to exclude. Reactions with score 0
                  are not part of the problem
  metData         boolean matrix with mets as rows and rxns as columns
                  saying which reaction produces each detected met (optional, default [])
  essentialRxns   cell array of reactions that are essential and that
                  have to be in the resulting model. ftINIT passes the
                  reactions that the tasks require (see fitTasks)
  prodWeight      a reward for each detected metabolite that has at least
                  one producing reaction turned on; net production is not
                  checked. The sum of these rewards and the scores for the
                  reactions is optimised. The reward gives a reason for
                  including reactions with a negative score for
                  connectivity reasons, which gives a more functional
                  network (default 0.5)
  allowExcretion  true if excretion of all metabolites is allowed.
                  This results in fewer reactions being considered
                  dead-ends, but not every reaction in the resulting model
                  can carry flux. If this is "false" then the
                  equality constraints are taken from model.b. If the
                  input model lacks exchange reactions, set it to true;
                  with false a large part of the model is excluded for
                  connectivity reasons
  remPosRev       If true, the positive reversible reactions are removed from the problem.
                  This is used in step 1 of ftINIT
  params          parameters for the MILP, for example MIPGap and TimeLimit
  startVals       Start values for the MILP. ftINIT passes the result of
                  the previous run when it reruns the MILP with a higher
                  MIPGap
  fluxes          Fluxes from the last run. The flux of an essential
                  irreversible reaction is forced to at least
                  min(0.99*|flux|, 0.1), clamped to the reaction's capacity
                  in the forced direction (its upper bound, as ftINIT
                  passes reactions in the forward direction)
  verbose         If true, the MILP progression is shown.
  forceOn         minimum flux magnitude enforced through a reaction once
                  it is turned on in the MILP (optional, default 0.1)
  resolveTies     if true, the MILP's degenerate optimum is solved to a
                  canonical answer in addition to the solver Seed: the
                  score objective is held at its optimum, then the count
                  of "on" removable (negative-score) reactions is
                  minimised, then, among the sparsest solutions, their
                  summed reaction-id rank is minimised. Positive-score
                  reactions have a continuous on/off variable (pinned near
                  1 by the score objective) and only the negative-score
                  ones have a binary on/off variable, so only these are
                  tied. The fluxes are then taken from the
                  least-total-flux solution for that on/off state and,
                  among those, from the one with the least id-rank-weighted
                  backward flux, so the directions and forced magnitudes
                  the next ftINIT step reads from them are fixed too.
                  This follows resolve_ties in raven-toolbox
                  (raven-gecko-parity#104). The phases are skipped when
                  the primary solve stops at params.TimeLimit, and the
                  primary solution is kept. Phases 2a and 2b get half of
                  params.TimeLimit each (rounded down, at least 1 s), the
                  flux phases the full limit. A warning names each phase
                  that does not prove optimal or finds no solution; when a
                  flux phase finds none, the fluxes of the phase before it
                  are kept (optional, default false)
  rankKeys        cell array with, for each reaction in model.rxns, the key
                  that the id-rank phases of resolveTies sort it by
                  (optional, default model.rxns). ftINIT passes the smallest
                  original reaction id of each merge group, so the order
                  does not depend on which member names the merged reaction
  pinFluxes       if false, resolveTies skips phases 2c and 2d and the
                  fluxes and metProduction are those of phase 2b
                  (optional, default true).
                  ftINIT passes false in its last step, whose fluxes no
                  later step reads

  deletedRxns     reactions which were deleted by the algorithm (only
                  rxns included in the problem)
  metProduction   logical array that indicates which of the metabolites in
                  metData could be produced (true) or not (false).
                  Metabolites with a producing reaction that is essential
                  are not part of the MILP and are not in the array
  res             The result from the MILP
  turnedOnRxns    The reactions determined to be present (only
                  rxns included in the problem)
  fluxes          The fluxes from the MILP

  This function is the implementation of the algorithm. See
  ftINIT for a higher-level function for model reconstruction.

Usage: [deletedRxns,metProduction,res,turnedOnRxns,fluxes]=ftINITInternalAlg(model,...
          rxnScores,metData,essentialRxns,prodWeight,allowExcretion,...
          remPosRev,params,startVals,fluxes,verbose,forceOn,resolveTies,...
          rankKeys,pinFluxes)

getINITModel#

function [model, metProduction, essentialRxnsForTasks, addedRxnsForTasks, deletedDeadEndRxns, deletedRxnsInINIT, taskReport]=getINITModel(refModel, tissue, varargin)

Generate a model using the original tINIT algorithm.

Generates a model using the tINIT algorithm, from proteomics, transcriptomics, metabolomics and metabolic tasks in any combination (PLoS Comput Biol. 2012;8(5):e1002518, Mol Syst Biol. 2014;10:721).

ftINIT replaces this function for context-specific model extraction. This function is kept to reproduce models built with the original tINIT. It issues a RAVEN:legacyMethod notice once per session, which warning('off','RAVEN:legacyMethod') silences.

Not all settings are available through this function; runINIT and fitTasks can be called individually to set the others.

Parameters

refModel : struct: a model structure in the closed form (no exchange reactions open), imported with importModel(filename,false). A model that is not loaded with importModel can lack the "unconstrained" field, and the function then raises RAVEN:badInput. The field can be added with model.unconstrained=false(numel(model.mets),1).
tissue : char: tissue to score for. Must exist in either hpaData.tissues or arrayData.tissues (scoreModel raises an error otherwise).

Name-Value Arguments

celltype : char

cell type to score for. Must exist in either hpaData.celltypes or arrayData.celltypes (scoreModel raises an error otherwise). The default is to use the maximum values among all the cell types for the tissue.

hpaData : struct

HPA data structure from parseHPA (optional if arrayData is supplied, default []).

arrayData : struct

gene expression data structure (optional if hpaData is supplied, default []). With fields:

genes : cell array with the unique gene names.
tissues : cell array with the tissue names. The list can contain duplicates, 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 and cell type. NaN marks a gene that was not measured.
threshold : a single value or a vector of gene expression thresholds, above which genes are considered to be "expressed". If omitted, the mean expression level of each gene across all tissues in arrayData is used as the threshold.
singleCells : binary value selecting whether to use the single-cell algorithm to identify expressed genes. If used, the cell subpopulations are specified in celltypes (default []).
plotResults : true if single cell probability distributions are plotted (default false).

metabolomicsData : cell

cell array with metabolite names that the model should produce (default []).

taskFile : char

a task list in Excel format. See parseTaskList for details (default []).

useScoresForTasks : logical

true if the calculated reaction scores are used as weights in the fitting to tasks. Scores above -0.1 are set to -0.1 (default true).

printReport : logical

true if a report is printed to the screen (default true).

taskStructure : struct

task structure as returned by parseTaskList. It defines the tasks when Excel sheets are not used and overrides taskFile (default []).

params : struct

parameter structure for use by the INIT solver (default []).

paramsFT : struct

parameter structure for the fitTasks step (default []).

seed : double

seed for the random number generator used to break ties among genes with equal uninformative (-Inf) scores, which makes the choice of the kept gene reproducible across runs (default [], which leaves the rng state unchanged).

Returns

model : struct

the resulting model structure.

metProduction : double

array that indicates which of the metabolites in metabolomicsData could be produced. The array is determined before the gap-filling for the defined tasks. To see which metabolites can be produced in the final model, use canExchange:

-2 : metabolite name not found in model.
-1 : metabolite found, but it could not be produced.
1 : metabolite could be produced.

essentialRxnsForTasks : cell

cell array of the reactions which were essential to perform the tasks.

addedRxnsForTasks : cell

cell array of the reactions which were added in order to perform the tasks.

deletedDeadEndRxns : cell

cell array of reactions deleted because they could not carry flux (INIT requires a functional input model).

deletedRxnsInINIT : cell

cell array of the reactions which were deleted by the INIT algorithm.

taskReport : struct

structure with the results for each task. With fields:

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.

Notes

The model keeps its exchange metabolites, since checkTasks and fitTasks put specific constraints on them for each task. The INIT algorithm removes exchange metabolites if any are present. Use importModel(file,false) to import a model with exchange metabolites remaining.

getINITSteps#

function steps = getINITSteps(varargin)

Define the steps used as input to ftINIT.

Defines the steps in the way ftINIT is run. Not all combinations of steps work. A new series that sets the HowToUsePrevResults option to "exclude", which no predefined series uses, may not work.

Name-Value Arguments

metsToIgnore : struct

structure describing mets that can be removed from the model before running ftINIT, such as water and other ubiquitous metabolites. With fields:

simpleMets.mets : names of metabolites to remove.
simpleMets.compsToKeep : compartments for which metabolites are kept.

(default []). Ignored for series "full".

series : char

describes the way to run ftINIT (default "1+1"):

"1+1" : standard behaviour. Step 1 and 2 described in the paper are merged into 1.
"2+1" : the 3-step procedure described in the paper. Faster and slightly less accurate than "1+1".
"1+0" : same as "1+1", but skips step 3 described in the paper. The model includes many reactions without GPRs. The series suits structural comparison, because the reactions removed in step 3 can be arbitrary and add no information. Faster than "1+1".
"2+0" : same as "2+1", but skips step 3 described in the paper. Faster and slightly less accurate than "1+0", with similar results.
"full" : 1-step run, similar to tINIT, without the model simplifications that ftINIT applies in the other series. Accurate, but slow. This series is used for testing.

Returns

steps : cell: cell array of steps, used as input to ftINIT.

groupRxnScores#

function newRxnScores = groupRxnScores(model, origRxnScores, origRxnIds, groupIds, origRxnsToZero)

Sums the reaction scores of all reactions that were merged into one by the

linear merge.

model          The model with linearly merged rxns.
origRxnScores  The rxnScores from the model before the linear merge.
origRxnIds     The rxn ids of the model before the linear merge.
groupIds       The groupIds vector output from linearMerge.
               There is one integer for each rxn in origRxnIds. 0 means
               the reaction was not merged with any other reaction. A
               non-zero integer means the reaction was merged with all
               other rxns sharing the same integer.
origRxnsToZero A logical vector saying which of the original rxns should not
               be part of the problem; the score for these reactions is set to 0.
               Since 0 means "excluded", an original score of 0 and a group
               score that sums to 0 while at least one score of the group
               is not 0 are set to 0.01.

mergeLinear#

function [reducedModel,origRxnIds,groupIds,reversedRxns]=mergeLinear(model, noMergeRxns)

Simplifies a model by merging rxns with linear dependencies. It differs

  from the grouping in simplifyModel in two ways:
  1) It returns which reactions were merged.
  2) It also merges reversible reactions.
  The reduced model has no genes and empty grRules.

  model                 a model structure
  noMergeRxns           Cell array with reaction IDs that are not allowed to be
                        merged

  reducedModel          an updated model structure
  origRxnIds            Cell array: original rxn ids, used together with the groupIds variable
  groupIds              Vector of ids: merged rxns have the same id, same length as origRxnIds
                        0 means a reaction was not merged. The groups are numbered 1, 2, ...
                        in the order of their first reaction in origRxnIds.
  reversedRxns          Vector of booleans, saying for each rxn if it is reversed in the combined rxns.

Usage: [reducedModel,origRxnIds,groupIds,reversedRxns]=mergeLinear(model,noMergeRxns)

parseHPA#

function hpaData=parseHPA(fileName, varargin)

Parse a database dump of the Human Protein Atlas (HPA).

Parameters

fileName : char: comma- or tab-separated database dump of HPA protein data. For details regarding the format, see http://www.proteinatlas.org/about/download.

Name-Value Arguments

version : double: ignored; the format is inferred from the column headers (default 19).

Returns

hpaData : struct

parsed HPA data with fields:

genes : cell array with the unique gene IDs (Ensembl in v>=18)
geneNames : cell array with the gene symbols (present in v>=18 only)
tissues : cell array with the tissue names (can repeat; one entry per tissue-cell-type combination)
celltypes : cell array with the cell type names for each tissue
levels : cell array with the unique expression level labels
types : cell array with the unique evidence types (v<18 only)
reliabilities : cell array with the unique reliability levels
gene2Level : sparse gene × tissue-cell-type matrix; value i,j is the index in hpaData.levels for gene i in tissue-cell-type j
gene2Type : sparse gene × tissue-cell-type matrix (v<18 only)
gene2Reliability : sparse gene × tissue-cell-type matrix; value i,j is the index in hpaData.reliabilities for gene i in tissue-cell-type j

Examples

hpaData = parseHPA(fileName);

parseHPArna#

function arrayData=parseHPArna(fileName, varargin)

Parse a dump of Human Protein Atlas (HPA) RNA-Seq data.

Parameters

fileName : char: tab-separated database dump of HPA RNA data. For details regarding the format, see http://www.proteinatlas.org/about/download.

Name-Value Arguments

version : double: ignored; the format is inferred from the column headers (default 19).

Returns

arrayData : struct

parsed HPA RNA data with fields:

genes : cell array with the unique Ensembl gene IDs
geneNames : cell array with the gene symbols
tissues : cell array with the unique tissue (or sample) names
levels : matrix of expression values (rows = genes, cols = tissues)

Examples

arrayData = parseHPArna(fileName);

prepINITModel#

function prepData = prepINITModel(origRefModel, taskStruct, varargin)

Precompute omics-independent data for ftINIT.

Runs the calculation steps that do not depend on the omics data, once per template model.

Parameters

origRefModel : struct: the template model, such as Human-GEM or Mouse-GEM.
taskStruct : struct: the essential tasks, for example loaded with taskStruct = parseTaskList('../data/metabolicTasks_Essential.txt'). The task ids must be character vectors: a non-text id or a missing id field raises RAVEN:badInput. Ids may repeat (for example as category codes) and may be empty; each task is one row group of the task list, and a task that cannot be performed is removed by its position.

Name-Value Arguments

spontRxnNames : cell: the spontaneous rxns (default {}).
convertGenes : logical: if true the genes are converted to gene names (from ENSEMBL). This uses translateGrRules, which is distributed with Human-GEM and not with RAVEN, so Human-GEM must be on the MATLAB path (default false).
customRxnsToIgnore : cell: reactions that are left outside the MILP when bit b7 of the RxnsToIgnoreMask of INITStepDesc is 1 (default {}).
extComp : char: name of the external compartment, such as 's' or 'e'. This is used for identifying exchange and import reactions (default 'e').
skipScaling : logical: if true the scaling step is not run on the minimal model. The scaling removes large differences between the stoichiometric coefficients within a reaction (see rescaleModelForINIT), since such differences create numerical issues in the ftINIT algorithm due to limitations in solver resolution. The scaling changes the stoichiometry of the reactions with large differences, and some reactions can become unable to carry flux, which can make the solution infeasible. If the solution is infeasible, turning the scaling off may help; the bounds then stay as in the template and the MILP caps the per-reaction big-M at 1000 (bigMCap in ftINITInternalAlg). In raven-toolbox's ftinit on Human-GEM (one cell line, big_m 1000) the pipeline without scaling was feasible but took 157 s against 83 s with scaling. The bounds of the scaled model are set to 1000 where the upper bound is above 0 and to -1000 where the lower bound is below 0; zero bounds stay 0. Only the minModel is scaled; the scaling is not present in the final model (default false).
runParallel : logical: true to run the task-essentiality check (checkTasks) across parallel workers, since each task is independent of the others (see parallelWorkersRAVEN) (default false, serial).

Returns

prepData : struct: the resulting prepData structure which is used as input to ftINIT.

rescaleModelForINIT#

function model=rescaleModelForINIT(model, varargin)

Rescale a model to avoid large stoichiometric differences.

Rescales the stoichiometric matrix for the MILP problem in ftINIT to reduce differences in flux magnitudes between reactions. Such differences cause, among other effects, difficulties with the tolerances for integer variables. In each reaction, each coefficient with an absolute value above maxStoichDiff times the smallest absolute coefficient of the reaction is set to that limit (keeping its sign). The reaction is then multiplied by a factor that sets the mean absolute coefficient to 1.

Parameters

model : struct: a model structure.

Name-Value Arguments

maxStoichDiff : double: largest allowed ratio between the largest and the smallest absolute stoichiometric coefficient in a reaction (default 25).

Returns

model : struct: the rescaled model structure.

runINIT#

function [outModel, deletedRxns, metProduction, fValue]=runINIT(model,varargin)

Solve the MILP of the original INIT algorithm.

Generates a model using the INIT algorithm, from proteomics, transcriptomics, metabolomics and metabolic tasks in any combination. This function is the implementation of the algorithm; see getINITModel for the higher-level reconstruction function, and PLoS Comput Biol. 2012;8(5):e1002518 for the formulation.

ftINIT replaces this function for context-specific model extraction and solves a different MILP in ftINITInternalAlg. This function is kept to reproduce models built with the original INIT. It issues a RAVEN:legacyMethod notice once per session, which warning('off','RAVEN:legacyMethod') silences.

Parameters

model : struct: a reference model structure.

Name-Value Arguments

rxnScores : double: a vector of scores for the reactions in the model. Positive scores are reactions to keep and negative scores are reactions to exclude (default all 0.0).
presentMets : cell: cell array with unique metabolite names that the model should produce (default []).
essentialRxns : cell: cell array of reactions that are essential and that have to be in the resulting model. fitTasks passes the reactions that the tasks require (default []).
prodWeight : double: a score for each metabolite that has net production. The sum of these weights and the scores for the reactions is optimised. Net production of metabolites gives a reason for including reactions with a negative score for connectivity reasons, which gives a more functional network (default 0.5).
allowExcretion : logical: true if excretion of all metabolites is allowed. This results in fewer reactions being considered dead-ends, but not every reaction in the resulting model can carry flux. If this is "false" then the equality constraints are taken from model.b. If the input model lacks exchange reactions, set it to true; with false a large part of the model is excluded for connectivity reasons (default false).
noRevLoops : logical: true if reversible reactions are constrained to carry flux in only one direction. This prevents reversible reactions from being assigned as connected wrongly (the forward and backward reactions can form a loop and therefore appear connected). It adds two binary variables per reversible reaction to the problem, which makes the MILP larger and more computationally intensive (default false).
params : struct: parameter structure for use by optimizeProb.
eps : double: the lower flux bound of a reaction that is turned on, which also applies to the essential reactions (capped by their upper bound) (default 1).

Returns

outModel : struct

the resulting model structure.

deletedRxns : cell

reactions which were deleted by the algorithm.

metProduction : double

array that indicates which of the metabolites in presentMets that could be produced:

-2 : metabolite name not found in model.
-1 : metabolite found, but it could not be produced.
1 : metabolite could be produced.

fValue : double

objective value (sum of (the negative of) reaction scores for the included reactions and prodWeight*number of produced metabolites).

scoreModel#

function [rxnScores, geneScores, hpaScores, arrayScores] = scoreModel(model,hpaData,arrayData,tissue,varargin)

Score reactions and genes in a model with complex gene rules.

Scores the reactions and genes in a model containing complex gene rules based on expression data from HPA, gene arrays or both. The model grRules should be standardised with standardizeGrRules before scoring.

Parameters

model : struct

a model structure.

hpaData : struct

HPA data structure from parseHPA (optional if arrayData is supplied, default []).

arrayData : struct

gene expression data structure (optional if hpaData is supplied, default []). With fields:

genes : cell array with the unique gene names.
tissues : cell array with the tissue names. The list can contain duplicates, 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 and cell type. NaN marks a gene that was not measured.
threshold : a single value or a vector of gene expression thresholds, above which genes are considered to be "expressed". If omitted, the mean expression level of each gene across all tissues in arrayData is used as the threshold.

tissue : char

tissue to score for. Must exist in either hpaData.tissues or arrayData.tissues (matched without regard to case); otherwise an error is raised.

Name-Value Arguments

celltype : char

cell type to score for. Must exist in either hpaData.celltypes or arrayData.celltypes (matched without regard to case); otherwise an error is raised (default [], which uses all cell types of the tissue, combined as set by multipleCellScoring).

noGeneScore : double

score for reactions without genes (default -2).

isozymeScoring : char

determines how scores are calculated for reactions with multiple genes joined by "OR" expression(s) ("min", "max", "median", "average") (default "max").

complexScoring : char

determines how scores are calculated for reactions with multiple genes joined by "AND" expression(s) ("min", "max", "median", "average") (default "min").

multipleCellScoring : char

determines how scores are calculated when several cell types are used ("max" or "average") (default "max").

hpaLevelScores : struct

structure with numerical scores for the expression level categories from HPA. The structure has a "names" and a "scores" field. The default names are High, Medium, Low, None, Strong, Moderate, Weak, Negative and Not detected, with the scores 20, 15, 10, -8, 20, 15, 10, -8 and -8.

dataPrecedence : char

how hpaData takes precedence over arrayData when both are supplied (default "gene"):

"gene" : per gene. A gene with HPA data is scored from it, a gene without is scored from arrayData, and a rule can mix the two.
"reaction" : per rule. As soon as any one gene in a rule has HPA data, the whole rule is scored from HPA alone and its array-only genes count as unmeasured.

Returns

rxnScores : double: scores for each of the reactions in model.
geneScores : double: scores for each of the genes in model. Genes which are not in the dataset(s) have NaN as scores. removeLowScoreGenes reads NaN as "no evidence" and leaves the gene in place, and prunes genes with a negative score.
hpaScores : double: scores for each of the genes in model if only taking hpaData into account. Genes which are not in the dataset(s) have -Inf as scores.
arrayScores : double: scores for each of the genes in model if only taking arrayData into account. Genes which are not in the dataset(s) have -Inf as scores.