diff --git a/src/change_model/addNewRxnsToEC.m b/src/change_model/addNewRxnsToEC.m index 7d7e99f5..2ba53466 100644 --- a/src/change_model/addNewRxnsToEC.m +++ b/src/change_model/addNewRxnsToEC.m @@ -35,7 +35,7 @@ % - rxns : cell array with unique strings that identifies each reaction. % - rxnNames : cell array with the names of each reaction. % - equations : cell array with equation strings. Decimal coefficients are expressed as "1.2". Reversibility is indicated by "<=>" or "=>". -% - grRules : cell array with the gene-reaction relationship for each reaction. For example "(A and B) or (C)" means that the reaction could be catalyzed by a complex between A & B or by C on its own. All the genes have to be present in model.genes. Add genes with addGenesRaven before calling this function if needed (opt, default ''). +% - grRules : cell array with the gene-reaction relationship for each reaction. For example "(A and B) or (C)" means that the reaction could be catalyzed by a complex between A & B or by C on its own. A rule with brackets is split into its isozymes in disjunctive normal form, so "A and (B or C)" is the two isozymes "A and B" and "A and C". All the genes have to be present in model.genes. Add genes with addGenesRaven before calling this function if needed (opt, default ''). % % The newEnzymes structure has the following fields: % @@ -136,11 +136,17 @@ % 3. Add new reactions + % The operators of a gene rule are case insensitive: write them in lower + % case, as the rest of the function and RAVEN's addRxns expect. Gene ids + % that contain the letters (Sand1, andA, ORF1) are left as they are. + newRxns.grRules = regexprep(convertCharArray(newRxns.grRules), ... + '(? 1); if ~isempty(toSplit) splitRxns = {}; splitRxnNames = {}; splitEquations = {}; splitGrRules = {}; for i=1:numel(toSplit) idx = toSplit(i); - newRules = split(newRxns.grRules(idx),'or'); - newRules = strtrim(erase(newRules, {'(' ')'})); + % One reaction per AND-clause of the gene rule in disjunctive normal form + newRules = cellfun(@(x) strjoin(x,' and '), isozymes{idx}, 'UniformOutput', false); for j=1:numel(newRules) splitRxns{end+1} = [newRxns.rxns{idx} '_EXP_' num2str(j)]; @@ -249,8 +256,8 @@ model.ec.eccodes(end+1) = {''}; % Get rules - genes = newRxns.grRules(rxnIdx); - genes = strtrim(erase(split(genes,'and'), {'(' ')'})); + genes = regexp(newRxns.grRules{rxnIdx},'[^\s()]+','match'); % Gene ids and the word 'and', without spaces and brackets + genes = genes(~ismember(lower(genes),{'and','or'})); [~,enzIdx] = ismember(genes, model.ec.genes); % Update the enzyme rxns matrix @@ -260,3 +267,59 @@ rxnsAdded = newRxns.rxns(:); end end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%Function that expands a gene rule into disjunctive normal form: a cell array +%with one cell per isozyme, each holding the gene ids that are ANDed in it. +%Gene ids are whole words, so an id that contains 'or' or 'and' is kept +%intact. An empty rule gives an empty cell array. +function dnf = gprToDnf(rule) +tokens = regexp(rule,'\(|\)|[^\s()]+','match'); +isOperator = strcmpi(tokens,'and') | strcmpi(tokens,'or'); +tokens(isOperator) = lower(tokens(isOperator)); % Operators are case insensitive +if isempty(tokens) + dnf = {}; + return +end +[dnf,pos] = parseOr(tokens,1,rule); +if pos <= numel(tokens) + error('Cannot parse the gene rule "%s"',rule) +end +end + +function [dnf,pos] = parseOr(tokens,pos,rule) +[dnf,pos] = parseAnd(tokens,pos,rule); +while pos <= numel(tokens) && strcmp(tokens{pos},'or') + [rhs,pos] = parseAnd(tokens,pos+1,rule); + dnf = [dnf,rhs]; +end +end + +function [dnf,pos] = parseAnd(tokens,pos,rule) +[dnf,pos] = parseAtom(tokens,pos,rule); +while pos <= numel(tokens) && strcmp(tokens{pos},'and') + [rhs,pos] = parseAtom(tokens,pos+1,rule); + product = cell(1,0); + for a = 1:numel(dnf) + for b = 1:numel(rhs) + product{end+1} = [dnf{a},rhs{b}]; %#ok + end + end + dnf = product; +end +end + +function [dnf,pos] = parseAtom(tokens,pos,rule) +if pos > numel(tokens) || any(strcmp(tokens{pos},{')','and','or'})) + error('Cannot parse the gene rule "%s"',rule) +elseif strcmp(tokens{pos},'(') + [dnf,pos] = parseOr(tokens,pos+1,rule); + if pos > numel(tokens) || ~strcmp(tokens{pos},')') + error('Cannot parse the gene rule "%s"',rule) + end + pos = pos + 1; +else + dnf = {tokens(pos)}; + pos = pos + 1; +end +end diff --git a/src/change_model/applyComplexData.m b/src/change_model/applyComplexData.m index b5a60944..68bcc976 100644 --- a/src/change_model/applyComplexData.m +++ b/src/change_model/applyComplexData.m @@ -37,6 +37,9 @@ % % Notes % ----- +% The subunit stoichiometry is assigned per protein ID: the proteins of a +% complex need not be listed in the order of model.ec.enzymes. +% % foundComplex and proposedComplex share the following columns: rxn % (reaction identifier), complexID and name (Complex Portal identifiers), % genes (gene names of the complex subunits), protID_model (matched protein @@ -154,8 +157,11 @@ foundComplex{foundComplexCount,7} = complexData(complexMatch).stochiometry; % Apply in rxnEnzMat + % The stoichiometry is taken per protein ID, as the proteins of the + % complex need not be listed in the order of model.ec.enzymes [~, modelProtsIdx] = ismember(modelProts, model.ec.enzymes); - model.ec.rxnEnzMat(i,modelProtsIdx) = complexData(complexMatch).stochiometry; + [~, complexPos] = ismember(modelProts, complexData(complexMatch).protID); + model.ec.rxnEnzMat(i,modelProtsIdx) = complexData(complexMatch).stochiometry(complexPos); elseif modComplexUnits > 1 % Only suggest potential complex if grRule > 1 protein moreUnitsInData = find(percMatch == 1 & totalMatch > 1); diff --git a/src/change_model/applyCustomKcats.m b/src/change_model/applyCustomKcats.m index 0cdd5e08..54f9ba67 100644 --- a/src/change_model/applyCustomKcats.m +++ b/src/change_model/applyCustomKcats.m @@ -43,9 +43,13 @@ % - gene_name : short gene name (optional, not used in matching). % - kcat : new kcat value (one per entry). % - rxns : reaction identifiers, multiple for the same kcat are separated by ',' (see further explanation below). -% - notes : will be appended to model.ec.notes (optional). +% - notes : will be appended to model.ec.notes, separated by ', ' from an existing note; an empty note changes nothing (optional). % - stoicho : complex stoichiometry, separated by ' + ' (examples: '1' or '3 + 1'), matching the order in proteins field. % +% An entry with an empty kcat is skipped. Protein and reaction identifiers +% are matched exactly, including their letter case. Light ecModels are +% supported: the rows of model.ec are matched by their ec reaction ids. +% % Matching order: % % 1. reactions are identified by .proteins and .rxns. @@ -121,125 +125,141 @@ enzInModel = cell(length(model.ec.rxns),1); ecRxnNoSuffix = regexprep(model.ec.rxns,'_EXP_\d+$',''); -% Implementation for full GECKO formulation -if ~model.ec.geckoLight - for i = 1:numel(customKcats.proteins) - if isempty(customKcats.proteins{i}) - %If only reaction ID(s) is/are specified (and no proteins), - %then apply the kcat to all isozymic reactions - rxns = strtrim(strsplit(customKcats.rxns{i}, ',')); - rxnIdxs = ismember(ecRxnNoSuffix,rxns); - rxnToUpdate(rxnIdxs) = 1; - model.ec.kcat(rxnIdxs) = customKcats.kcat(i); - - % Add note mentioning manual kcat change - model.ec.source(rxnIdxs) = {'custom'}; - if isfield(customKcats,'notes') - matchedIdxs = find(rxnIdxs); - for j = 1:numel(matchedIdxs) - if isempty(model.ec.notes{matchedIdxs(j), 1}) && ~isempty(customKcats.notes{i}) - model.ec.notes{matchedIdxs(j), 1} = customKcats.notes{i}; - else - model.ec.notes{matchedIdxs(j), 1} = [model.ec.notes{matchedIdxs(j), 1} ', ' customKcats.notes{i}]; - end - end +for i = 1:numel(customKcats.proteins) + if isnan(customKcats.kcat(i)) + printOrange(['WARNING: Entry ' num2str(i) ' has no kcat; it is skipped.']); + continue + end + if isempty(customKcats.proteins{i}) + %If only reaction ID(s) is/are specified (and no proteins), + %then apply the kcat to all isozymic reactions + rxns = strtrim(strsplit(customKcats.rxns{i}, ',')); + rxnIdxs = ismember(ecRxnNoSuffix,rxns); + rxnToUpdate(rxnIdxs) = 1; + model.ec.kcat(rxnIdxs) = customKcats.kcat(i); + + % Add note mentioning manual kcat change + model.ec.source(rxnIdxs) = {'custom'}; + if isfield(customKcats,'notes') + matchedIdxs = find(rxnIdxs); + for j = 1:numel(matchedIdxs) + model.ec.notes{matchedIdxs(j), 1} = appendNote(model.ec.notes{matchedIdxs(j), 1}, customKcats.notes{i}); end + end + else + prots = strtrim(strsplit(customKcats.proteins{i}, '+')); + + % Find the index for the enzymes which kcat will be changed. In + % case the protein is not in the model generate a warning. + try + enzIdx = cellfun(@(x) find(strcmp(model.ec.enzymes, x)), prots); + catch + enzIdx = []; + printOrange(['WARNING: Protein(s) ' customKcats.proteins{i} ' were not found in the model.']); + end + + % if not specific reactions are defined, find all the reaction + % index where the enzyme is used + if isempty(customKcats.rxns{i}) + temp_rxnIdxs = arrayfun(@(x) find(model.ec.rxnEnzMat(:, x)), enzIdx, 'UniformOutput', false); + % otherwhise, If a set of reactions if defined, only get the index for those + % but ignore any _EXP_ suffixes else - prots = strtrim(strsplit(customKcats.proteins{i}, '+')); - - % Find the index for the enzymes which kcat will be changed. In - % case the protein is not in the model generate a warning. - try - enzIdx = cellfun(@(x) find(strcmpi(model.ec.enzymes, x)), prots); - catch - enzIdx = []; - printOrange(['WARNING: Protein(s) ' customKcats.proteins{i} ' were not found in the model.']); - end + rxns = strtrim(strsplit(customKcats.rxns{i}, ',')); + temp_rxnIdxs = cellfun(@(x) find(strcmp(ecRxnNoSuffix, x)), rxns, 'UniformOutput', false); + end - % if not specific reactions are defined, find all the reaction - % index where the enzyme is used - if isempty(customKcats.rxns{i}) - temp_rxnIdxs = arrayfun(@(x) find(model.ec.rxnEnzMat(:, x)), enzIdx, 'UniformOutput', false); - % otherwhise, If a set of reactions if defined, only get the index for those - % but ignore any _EXP_ suffixes - else - rxns = strtrim(strsplit(customKcats.rxns{i}, ',')); - temp_rxnIdxs = arrayfun(@(x) find(strcmpi(ecRxnNoSuffix, x)), rxns, 'UniformOutput', false); + if ~isempty(temp_rxnIdxs) + rxnIdxs = []; + for j = 1:numel(temp_rxnIdxs) + rxnIdxs = [rxnIdxs; temp_rxnIdxs{j}]; end - if ~isempty(temp_rxnIdxs) - rxnIdxs = []; - for j = 1:numel(temp_rxnIdxs) - rxnIdxs = [rxnIdxs; temp_rxnIdxs{j}]; - end + % Check when multiple proteins are involved, since it can return same rxn n times + rxnIdxs = unique(rxnIdxs); %unique(rxnIdxs{1, :}); - % Check when multiple proteins are involved, since it can return same rxn n times - rxnIdxs = unique(rxnIdxs); %unique(rxnIdxs{1, :}); + for j = 1:numel(rxnIdxs) + % Get all the enzymes involved in the reaction + allEnzInRxn = find(model.ec.rxnEnzMat(rxnIdxs(j),:)); - for j = 1:numel(rxnIdxs) - % Get all the enzymes involved in the reaction - allEnzInRxn = find(model.ec.rxnEnzMat(rxnIdxs(j),:)); + C = intersect(enzIdx, allEnzInRxn); - C = intersect(enzIdx, allEnzInRxn); - - % Determine the match percentage bewteen the rules - if numel(C) == numel(enzIdx) && numel(C) == numel(allEnzInRxn) - match = 1; + % Determine the match percentage bewteen the rules + if numel(C) == numel(enzIdx) && numel(C) == numel(allEnzInRxn) + match = 1; + else + if numel(enzIdx) < numel(allEnzInRxn) + match = numel(C) / numel(allEnzInRxn); else - if numel(enzIdx) < numel(allEnzInRxn) - match = numel(C) / numel(allEnzInRxn); - else - match = numel(C) / numel(enzIdx); - end + match = numel(C) / numel(enzIdx); end + end - % Update the kcat value stored in the model, if match 100%, - % otherwhise if >= 0.5 inform to the user to be validated - if match == 1 - rxnToUpdate(rxnIdxs(j)) = 1; - model.ec.kcat(rxnIdxs(j)) = customKcats.kcat(i); - - % Add note mentioning manual kcat change - model.ec.source{rxnIdxs(j),1} = 'custom'; - if isfield(customKcats,'notes') - if isempty(model.ec.notes{rxnIdxs(j), 1}) && ~isempty(customKcats.notes{i}) - model.ec.notes{rxnIdxs(j), 1} = customKcats.notes{i}; - else - model.ec.notes{rxnIdxs(j), 1} = [model.ec.notes{rxnIdxs(j), 1} ', ' customKcats.notes{i}]; - end - end - elseif match >= 0.5 && match < 1 - rxnNotMatch(rxnIdxs(j)) = 1; - evaluatedRule{rxnIdxs(j), 1} = customKcats.proteins{i}; - enzInModel{rxnIdxs(j), 1} = strjoin(model.ec.enzymes(allEnzInRxn), ' + '); + % Update the kcat value stored in the model, if match 100%, + % otherwhise if >= 0.5 inform to the user to be validated + if match == 1 + rxnToUpdate(rxnIdxs(j)) = 1; + model.ec.kcat(rxnIdxs(j)) = customKcats.kcat(i); + + % Add note mentioning manual kcat change + model.ec.source{rxnIdxs(j),1} = 'custom'; + if isfield(customKcats,'notes') + model.ec.notes{rxnIdxs(j), 1} = appendNote(model.ec.notes{rxnIdxs(j), 1}, customKcats.notes{i}); end + elseif match >= 0.5 && match < 1 + rxnNotMatch(rxnIdxs(j)) = 1; + evaluatedRule{rxnIdxs(j), 1} = customKcats.proteins{i}; + enzInModel{rxnIdxs(j), 1} = strjoin(model.ec.enzymes(allEnzInRxn), ' + '); end end end end +end - % Apply the new kcat values to the model - if ~isempty(find(rxnToUpdate, 1)) - model = applyKcatConstraints(model, rxnToUpdate); +% Apply the new kcat values to the model +if ~isempty(find(rxnToUpdate, 1)) + if model.ec.geckoLight + %Light ecModels update by reaction of the model, not by ec row + updateRxns = ismember(model.rxns, extractAfter(model.ec.rxns(rxnToUpdate), 4)); else - printOrange('WARNING: No matches found. Consider checking the IDs or proteins in customKcats.'); + updateRxns = rxnToUpdate; end + model = applyKcatConstraints(model, updateRxns); +else + printOrange('WARNING: No matches found. Consider checking the IDs or proteins in customKcats.'); +end + +rxnUpdated = model.ec.rxns(find(rxnToUpdate)); - rxnUpdated = model.ec.rxns(find(rxnToUpdate)); - - % Remove from notMatch those reactions already updated - remove = and(rxnToUpdate, rxnNotMatch); - rxnNotMatch(remove) = 0; - evaluatedRule(remove) = ''; - enzInModel(remove) = ''; - - idRxns = model.ec.rxns(find(rxnNotMatch)); - fullIdx = cellfun(@(x) find(strcmpi(model.rxns, x)), idRxns); - rxnsNames = model.rxnNames(fullIdx); - evaluatedRule = evaluatedRule(~cellfun('isempty', evaluatedRule)); - enzInModel = enzInModel(~cellfun('isempty', enzInModel)); - rules = model.grRules(fullIdx); - notMatch = table(idRxns, rxnsNames, evaluatedRule, enzInModel, rules, ... - 'VariableNames',{'rxns', 'name', 'custom enzymes', 'enzymes in model', 'rules'}); +% Remove from notMatch those reactions already updated +remove = and(rxnToUpdate, rxnNotMatch); +rxnNotMatch(remove) = 0; +evaluatedRule(remove) = ''; +enzInModel(remove) = ''; + +idRxns = model.ec.rxns(find(rxnNotMatch)); +if model.ec.geckoLight + modelRxnIds = extractAfter(idRxns, 4); +else + modelRxnIds = idRxns; +end +fullIdx = cellfun(@(x) find(strcmpi(model.rxns, x), 1), modelRxnIds); +rxnsNames = model.rxnNames(fullIdx); +evaluatedRule = evaluatedRule(~cellfun('isempty', evaluatedRule)); +enzInModel = enzInModel(~cellfun('isempty', enzInModel)); +rules = model.grRules(fullIdx); +notMatch = table(idRxns, rxnsNames, evaluatedRule, enzInModel, rules, ... + 'VariableNames',{'rxns', 'name', 'custom enzymes', 'enzymes in model', 'rules'}); +end + +function notes = appendNote(notes, note) +% An empty note changes nothing; a note is separated by ', ' from an existing one. +note = strtrim(note); +if isempty(note) + return +elseif isempty(notes) + notes = note; +else + notes = [notes ', ' note]; end end diff --git a/src/change_model/applyKcatConstraints.m b/src/change_model/applyKcatConstraints.m index d8c25de2..e87f7ce2 100644 --- a/src/change_model/applyKcatConstraints.m +++ b/src/change_model/applyKcatConstraints.m @@ -3,7 +3,9 @@ % % Existing enzyme constraints are first removed (unless updateRxns is % provided), and new constraints are defined based on the content of -% model.ec.kcat. +% model.ec.kcat. When every kcat is 0, the existing constraints of a full +% ecModel are removed and no new ones are written; a light ecModel is +% returned unchanged. % % Parameters % ---------- @@ -66,8 +68,13 @@ ' not a valid GECKO3 ecModel. First run makeEcModel(model).']) end if all(model.ec.kcat==0) - printOrange('WARNING: No kcat values are provided in model.ec.kcat, model remains unchanged.\n'); - return + if model.ec.geckoLight + printOrange('WARNING: No kcat values are provided in model.ec.kcat, model remains unchanged.\n'); + return + end + % In a full model the existing enzyme constraints of the selected + % reactions are still cleared below, and no new ones are written. + printOrange('WARNING: No kcat values are provided in model.ec.kcat, existing enzyme constraints are cleared and no new ones are written.\n'); end %Clear existing incorporation of enzyme usage @@ -135,12 +142,17 @@ % summed per reaction), and divide by their kcat, to get a vector % of MW/kcat values. MWkcat = (model.ec.rxnEnzMat(ecIdx,:) * model.ec.mw) ./ model.ec.kcat(ecIdx); - % If kcat was zero, MWkcat is Inf. If no enzyme info was found, - % MWkcat is NaN. Drop those isozymes *before* selecting the - % cheapest one: correcting them to zero and including them in the - % min() below would let an isozyme with no kcat assigned always - % win over a real kcat on another isozyme of the same reaction. - validKcat = ~isinf(MWkcat) & ~isnan(MWkcat); + % If kcat was zero, MWkcat is Inf (or NaN without enzymes). An + % isozyme row without any enzyme has MWkcat zero. Drop those + % isozymes *before* selecting the cheapest one: including them in + % the min() below would let an isozyme with no kcat or no enzyme + % always win over a real kcat on another isozyme of the same reaction. + % A row with a kcat that is not positive or an enzyme mass that is not + % positive has no usable cost either. + hasEnzyme = full(sum(model.ec.rxnEnzMat(ecIdx,:) ~= 0,2)) > 0; + enzMass = full(model.ec.rxnEnzMat(ecIdx,:) * model.ec.mw); + validKcat = ~isinf(MWkcat) & ~isnan(MWkcat) & hasEnzyme & ... + model.ec.kcat(ecIdx) > 0 & enzMass > 0; if any(validKcat) % Select the lowest MW/kcat (= most efficient), and convert to /hour model.S(prot_pool_idx, hasEc(i)) = -min(MWkcat(validKcat)/3600); diff --git a/src/change_model/getKcatAcrossIsozymes.m b/src/change_model/getKcatAcrossIsozymes.m index 47c03664..71dfe485 100644 --- a/src/change_model/getKcatAcrossIsozymes.m +++ b/src/change_model/getKcatAcrossIsozymes.m @@ -4,7 +4,9 @@ % For reactions without kcat value (0 in model.ec.kcat), isozymes are found % (being based on the same reaction in the conventional GEM), that do have a % kcat value assigned. The mean kcat value of these isozymes is then used to -% fill in model.ec.kcat. +% fill in model.ec.kcat. Only isozymes with a kcat above 0 enter the mean: a +% NaN or negative kcat of an isozyme is left out, and an entry whose kcat is +% NaN is not filled. % % Parameters % ---------- @@ -34,10 +36,11 @@ end noKcats = model.ec.kcat==0; +hasKcat = model.ec.kcat>0; % NaN and negative kcats are not isozyme values rxnIDs = regexprep(model.ec.rxns,'_EXP_\d+',''); noKcatID = rxnIDs(noKcats); -yesKcatID = rxnIDs(~noKcats); -yesKcatVal = model.ec.kcat(~noKcats); +yesKcatID = rxnIDs(hasKcat); +yesKcatVal = model.ec.kcat(hasKcat); noKcatVal = cellfun(@(x) strcmp(x, yesKcatID), noKcatID, 'UniformOutput', false); noKcatVal = cell2mat(cellfun(@(x) mean(yesKcatVal(x)), noKcatVal, 'UniformOutput', false)); diff --git a/src/change_model/makeEcModel.m b/src/change_model/makeEcModel.m index f8662b7e..313e1c0d 100644 --- a/src/change_model/makeEcModel.m +++ b/src/change_model/makeEcModel.m @@ -226,10 +226,13 @@ else %Different strategy for GECKO light: Each reaction can exist multiple times in %ec.rxns and similar fields - one time per isozyme. The number of copies is - %the number of ORs in the GPR + 1 - numOrs = count(model.grRules(rxnWithGene), ' or '); - cpys = numOrs + 1; - prevNumRxns = length(numOrs); + %the number of AND-clauses of the GPR in disjunctive normal form, so + %'G1 and (G2 or G3)' gives the isozymes {G1,G2} and {G1,G3}. + isozymes = cellfun(@gprToDnf, model.grRules(rxnWithGene), 'UniformOutput', false); + % A reaction with gene entries but an empty rule keeps one isozyme without enzymes + isozymes(cellfun(@isempty,isozymes)) = {{cell(1,0)}}; + cpys = cellfun(@numel, isozymes); + prevNumRxns = numel(isozymes); cpyIndices = repelem(rxnWithGene, cpys).'; %.' only matter when number of rxns is 1 %loop through and add a prefix with an isozyme index to the rxns %we just give a fixed-length number as prefix, and assume that 999 is enough @@ -362,51 +365,18 @@ for r=1:numel(rxnWithGene) rxnGenes = model.genes(find(model.rxnGeneMat(rxnWithGene(r),:))); [~,locEnz] = ismember(rxnGenes,ec.genes); % Could also parse directly from rxnGeneMat, but some genes might be missing from Uniprot DB - if locEnz ~= 0 - ec.rxnEnzMat(r,locEnz) = 1; %Assume 1 copy per subunit or enzyme, can be modified later - end + locEnz = locEnz(locEnz ~= 0); % Subunits missing from the databases are left out; the others still couple the reaction + ec.rxnEnzMat(r,locEnz) = 1; %Assume 1 copy per subunit or enzyme, can be modified later end else %For light models, we need to split up all GPRs ec.rxnEnzMat = zeros(numel(ec.rxns),numel(ec.genes)); % Non-zeros will indicate the number of subunits nextIndex = 1; - %For full model generation, the GPRs are controlled in expandModel, but - %here we need to make an explicit format check - indexes2check = findPotentialErrors(model.grRules,model); - if ~isempty(indexes2check) - printOrange('Run standardizeGrRules(model) for a more detailed warning.\n'); - printOrange('For Human-GEM, these reactions can be corrected using simplifyGrRules.\n'); - end - for i=1:prevNumRxns - %ind is the index in the model, not to confuse with the index in the ec struct (i), - %which only contains reactions with GPRs. - ind = rxnWithGene(i); - %Get rid of all '(' and ')' since I'm not looking at complex stuff - %anyways - geneString=model.grRules{ind}; - geneString=strrep(geneString,'(',''); - geneString=strrep(geneString,')',''); - geneString=strrep(geneString,' or ',';'); - - if (numOrs(i) == 0) - geneNames = {geneString}; - else - %Split the string into gene names - geneNames=regexp(geneString,';','split'); - end - - %Now loop through the isozymes and set the rxnGeneMat - for j = 1:length(geneNames) - %Find the gene in the gene list If ' and ' relationship, first - %split the genes - fnd = strfind(geneNames{j},' and '); - if ~isempty(fnd) - andGenes=regexp(geneNames{j},' and ','split'); - ec.rxnEnzMat(nextIndex,ismember(ec.genes,andGenes)) = 1; %should be subunit stoichoimetry - else - ec.rxnEnzMat(nextIndex,ismember(ec.genes,geneNames(j)))=1;%should be subunit stoichoimetry - end + %Each isozyme is one AND-clause of the GPR and gets its own row of + %ec.rxnEnzMat, in the order of the copies of the reaction in ec.rxns. + for j = 1:numel(isozymes{i}) + ec.rxnEnzMat(nextIndex,ismember(ec.genes,isozymes{i}{j})) = 1;%should be subunit stoichoimetry nextIndex = nextIndex + 1; end end @@ -485,23 +455,57 @@ end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%Function that gets the model field grRules and returns the indexes of the -%rules in which the pattern ") and (" is present. -%Copied from standardizeGrRules -% TODO: Make this an accessible function in a separate file in RAVEN and remove this -%implementation. -function indexes2check = findPotentialErrors(grRules,model) -indxs_l = find(~cellfun(@isempty,strfind(grRules,') and ('))); -indxs_l_L = find(~cellfun(@isempty,strfind(grRules,') and'))); -indxs_l_R = find(~cellfun(@isempty,strfind(grRules,'and ('))); -indexes2check = vertcat(indxs_l,indxs_l_L,indxs_l_R); -indexes2check = unique(indexes2check); - -if ~isempty(indexes2check) - textToPrint = 'WARNING: Potentially problematic ") AND (" in the grRules for reaction(s):\n'; - for i=1:numel(indexes2check) - textToPrint=[textToPrint '\t' model.rxns{indexes2check(i)} '\n']; +%Function that expands a gene rule into disjunctive normal form: a cell array +%with one cell per isozyme, each holding the gene ids that are ANDed in it. +%Gene ids are whole words, so an id that contains 'or' or 'and' is kept +%intact. An empty rule gives an empty cell array. +function dnf = gprToDnf(rule) +tokens = regexp(rule,'\(|\)|[^\s()]+','match'); +isOperator = strcmpi(tokens,'and') | strcmpi(tokens,'or'); +tokens(isOperator) = lower(tokens(isOperator)); % Operators are case insensitive +if isempty(tokens) + dnf = {}; + return +end +[dnf,pos] = parseOr(tokens,1,rule); +if pos <= numel(tokens) + error('Cannot parse the gene rule "%s"',rule) +end +end + +function [dnf,pos] = parseOr(tokens,pos,rule) +[dnf,pos] = parseAnd(tokens,pos,rule); +while pos <= numel(tokens) && strcmp(tokens{pos},'or') + [rhs,pos] = parseAnd(tokens,pos+1,rule); + dnf = [dnf,rhs]; +end +end + +function [dnf,pos] = parseAnd(tokens,pos,rule) +[dnf,pos] = parseAtom(tokens,pos,rule); +while pos <= numel(tokens) && strcmp(tokens{pos},'and') + [rhs,pos] = parseAtom(tokens,pos+1,rule); + product = cell(1,0); + for a = 1:numel(dnf) + for b = 1:numel(rhs) + product{end+1} = [dnf{a},rhs{b}]; %#ok + end end - printOrange(textToPrint); + dnf = product; +end +end + +function [dnf,pos] = parseAtom(tokens,pos,rule) +if pos > numel(tokens) || any(strcmp(tokens{pos},{')','and','or'})) + error('Cannot parse the gene rule "%s"',rule) +elseif strcmp(tokens{pos},'(') + [dnf,pos] = parseOr(tokens,pos+1,rule); + if pos > numel(tokens) || ~strcmp(tokens{pos},')') + error('Cannot parse the gene rule "%s"',rule) + end + pos = pos + 1; +else + dnf = {tokens(pos)}; + pos = pos + 1; end end diff --git a/src/change_model/mapRxnsToConv.m b/src/change_model/mapRxnsToConv.m index e68f3dff..ee599626 100644 --- a/src/change_model/mapRxnsToConv.m +++ b/src/change_model/mapRxnsToConv.m @@ -43,11 +43,12 @@ rxnIDs = ecModel.rxns; % Invert flux of _REV reactions -revRxns = endsWith(rxnIDs,'_REV') | contains(rxnIDs,'_REV_EXP_'); +% The suffixes are added in the order _REV, _EXP_n, so strip _EXP_n first and +% then _REV, each only at the end of the id (an id like R_REVX is left intact). +rxnIDs = regexprep(rxnIDs,'_EXP_\d+$',''); +revRxns = endsWith(rxnIDs,'_REV'); fluxes(revRxns,:) = -fluxes(revRxns,:); -rxnIDs(revRxns) = replace(rxnIDs(revRxns),'_REV',''); -% Remove _EXP_. suffixes -rxnIDs = regexprep(rxnIDs,'_EXP_\d+',''); +rxnIDs(revRxns) = regexprep(rxnIDs(revRxns),'_REV$',''); % Map and sum fluxes to converted reaction IDs [rxnIDmap, convRxnID] = findgroups(rxnIDs); diff --git a/src/download_databases/loadBRENDAdata.m b/src/download_databases/loadBRENDAdata.m index 6844ae25..07a96db4 100644 --- a/src/download_databases/loadBRENDAdata.m +++ b/src/download_databases/loadBRENDAdata.m @@ -21,6 +21,16 @@ % 'median' (default: params.kcatAggregation from the model adapter, % or 'max' if the adapter does not set it). A triple's max is the % highest turnover reported for it; its median is the middle one. +% filterOutlierKcats : logical +% whether to drop a row whose kcat is implausibly low for its EC +% (default: true). A row is dropped when its kcat is more than four +% orders of magnitude below the median kcat of its EC. The reference +% median of an EC is taken over the positive kcat values of the max and +% median views of kcat.tsv and of the kcat values derived from sa.tsv +% and mw.tsv, before any row is dropped, so it does not depend on +% kcatAggregation. Wildcard ECs (containing '-') do not contribute to a +% reference. The EC level widens from w.x.y.z to w.x.y, w.x and w while +% fewer than five values are available. Set to false to keep every row. % % Returns % ------- @@ -35,6 +45,7 @@ % [KCATcell, SAcell] = loadBRENDAdata(); % [KCATcell, SAcell] = loadBRENDAdata('modelAdapter', adapter); % [KCATcell, SAcell] = loadBRENDAdata('kcatAggregation', 'median'); +% [KCATcell, SAcell] = loadBRENDAdata('filterOutlierKcats', false); % % See also % -------- @@ -47,9 +58,14 @@ 'kcatAggregation', [], @(x) assert(isempty(x) || any(strcmpi(x, ... {'max','median'})), ... 'loadBRENDAdata:invalidKcatAggregation', ... - 'kcatAggregation must be ''max'' or ''median''.')}); -modelAdapter = p.modelAdapter; -kcatAggregation = p.kcatAggregation; + 'kcatAggregation must be ''max'' or ''median''.'); ... + 'filterOutlierKcats', true, @(x) assert((islogical(x) || isnumeric(x)) && ... + isscalar(x), ... + 'loadBRENDAdata:invalidFilterOutlierKcats', ... + 'filterOutlierKcats must be a logical scalar.')}); +modelAdapter = p.modelAdapter; +kcatAggregation = p.kcatAggregation; +filterOutlierKcats = logical(p.filterOutlierKcats); if isempty(modelAdapter) modelAdapter = ModelAdapterManager.getDefault(); @@ -59,7 +75,8 @@ end kcatAggregation = resolveKcatAggregation(kcatAggregation, modelAdapter); -%Column 4 holds the per-triple maximum and column 5 the median. +%Column 4 holds the per-triple maximum and column 5 the median; both are +%read, so that the outlier filter can use the max and the median views. valueCol = 4 + strcmpi(kcatAggregation,'median'); basePath = modelAdapter.getBrendaDBFolder(); @@ -69,31 +86,31 @@ %Extract BRENDA DATA from files information. kcat.tsv/sa.tsv are seven %columns (ec_code, substrate, organism, value_max, value_median, n, -%references); valueCol picks max or median. mw.tsv is six columns (no -%substrate-level aggregation choice): ec_code, substrate, organism, value, -%n, references. -KCATcell = openDataFile(KCAT_file,1,'%q %q %q %f %f %f %q',valueCol); +%references); the data cells hold ec_code, substrate, organism, the max and +%the median. mw.tsv is six columns (no substrate-level aggregation choice): +%ec_code, substrate, organism, value, n, references. +KCAT = openDataFile(KCAT_file,1,'%q %q %q %f %f %f %q',[4 5]); scalingFactor = 1/60; %[umol/min/mg] -> [mmol/s/g] -SA = openDataFile(SA_file,scalingFactor,'%q %q %q %f %f %f %q',valueCol); +SA = openDataFile(SA_file,scalingFactor,'%q %q %q %f %f %f %q',[4 5]); scalingFactor = 1/1000; %[g/mol] -> [g/mmol] MW = openDataFile(MW_file,scalingFactor,'%q %q %q %f %f %q',4); -for i=1:4 - SAcell{i} = []; -end -previousEC = []; EC_indexes = []; - %build an index on MW{1} to speed things up a bit %first just extract the genus (i.e. the first part of the name) -MWECNum = upper(unique(MW{1})); +MWECNum = unique(upper(MW{1})); MWECNumIndices = cell(length(MWECNum),1); MWECNumHashMap = containers.Map(MWECNum,1:length(MWECNum)); for i = 1:length(MW{1}) - matchInd = cell2mat(values(MWECNumHashMap, MW{1}(i))); + matchInd = cell2mat(values(MWECNumHashMap, upper(MW{1}(i)))); MWECNumIndices{matchInd} = [MWECNumIndices{matchInd};i]; end - +%Match each specific activity to the molecular weight of its EC and +%organism. The first match is used. +previousEC = []; EC_indexes = []; +saRows = zeros(length(SA{1}),1); +mwVals = zeros(length(SA{1}),1); +nMatched = 0; for i=1:length(SA{1}) %Gets the indexes of the EC repetitions in the MW cell for every %new (different) EC @@ -111,18 +128,85 @@ % each EC# / Orgaism is reported on the file org_index = find(strcmpi(SA{3}(i),mwEC{1}),1); if ~isempty(org_index) - SAcell{1} = [SAcell{1};SA{1}(i)]; - SAcell{2} = [SAcell{2};SA{3}(i)]; - SAcell{3} = [SAcell{3}; SA{4}(i)*mwEC{2}(org_index)]; %[1/s] - SAcell{4} = [SAcell{4}; mwEC{2}(org_index)]; + nMatched = nMatched + 1; + saRows(nMatched) = i; + mwVals(nMatched) = mwEC{2}(org_index); end previousEC = SA{1}(i); end +saRows = saRows(1:nMatched); +mwVals = mwVals(1:nMatched); + +%Specific activities as kcat values [1/s], for the max and the median +saEC = SA{1}(saRows); +saMax = SA{4}(saRows).*mwVals; +saMedian = SA{5}(saRows).*mwVals; + +if valueCol == 4 + saSelected = saMax; +else + saSelected = saMedian; +end + +%Drop kcats that are implausibly low for their EC +keepKcat = true(size(KCAT{1})); +keepSA = true(size(saEC)); +if filterOutlierKcats + ecPool = [KCAT{1}; KCAT{1}; saEC; saEC]; + valPool = [KCAT{4}; KCAT{5}; saMax; saMedian]; + keepKcat = ~isOutlier(KCAT{1},KCAT{valueCol},ecPool,valPool); + keepSA = ~isOutlier(saEC,saSelected,ecPool,valPool); +end -function data_cell = openDataFile(fileName,scalingFactor,formatSpec,valueCol) +KCATcell = {KCAT{1}(keepKcat), KCAT{2}(keepKcat), KCAT{3}(keepKcat), ... + KCAT{valueCol}(keepKcat)}; +if nMatched > 0 + SAcell = {saEC(keepSA), SA{3}(saRows(keepSA)), saSelected(keepSA), ... + mwVals(keepSA)}; +else + SAcell = {[],[],[],[]}; +end +end + +function data_cell = openDataFile(fileName,scalingFactor,formatSpec,valueCols) fID = fopen(fileName); raw = textscan(fID,formatSpec,'delimiter','\t','HeaderLines',2); fclose(fID); -data_cell = {raw{1}, raw{2}, raw{3}, raw{valueCol}*scalingFactor}; +data_cell = [raw(1:3), cellfun(@(c) raw{c}*scalingFactor, num2cell(valueCols), ... + 'UniformOutput', false)]; end + +function outlier = isOutlier(ecs,vals,ecPool,valPool) +%A value is an outlier if it is more than 'decades' orders of magnitude below +%the median kcat of its own EC, taken over the pool of values. The EC widens +%from the full code to its 3-, 2- and 1-level prefix while fewer than minRef +%values are available. +decades = 4; +minRef = 5; +usePool = valPool > 0 & ~contains(ecPool,'-'); +ecPool = ecPool(usePool); +logPool = log10(valPool(usePool)); +refLog = NaN(size(ecs)); +if isempty(ecs) || isempty(logPool) + outlier = false(size(ecs)); + return +end +for level = 4:-1:1 + poolKey = ecPrefix(ecPool,level); + [grp,keys] = findgroups(poolKey); + nPerKey = accumarray(grp,1); + medPerKey = splitapply(@median,logPool,grp); + [found,loc] = ismember(ecPrefix(ecs,level),keys); + use = isnan(refLog) & found; + use(use) = nPerKey(loc(use)) >= minRef; + refLog(use) = medPerKey(loc(use)); +end +outlier = vals > 0 & ~isnan(refLog) & log10(vals) < refLog - decades; +end + +function prefix = ecPrefix(ecs,level) +%The first 'level' dot-separated parts of each EC code +prefix = regexp(ecs,['^[^.]*(?:\.[^.]*){0,' num2str(level-1) '}'],'match','once'); +prefix(cellfun(@isempty,prefix)) = {''}; +prefix = reshape(prefix,size(ecs)); end diff --git a/src/flux_analysis/ecFSEOF.m b/src/flux_analysis/ecFSEOF.m index 0ec57b93..df24c4b7 100644 --- a/src/flux_analysis/ecFSEOF.m +++ b/src/flux_analysis/ecFSEOF.m @@ -3,7 +3,9 @@ % % Runs Flux-Scanning with Enforced Objective Function (FSEOF) for a % specified production target, to identify reaction and gene targets for -% metabolic engineering. +% metabolic engineering. Reactions are classified as targets with the same +% rule as RAVEN's FSEOF: the correlation between the absolute flux of a +% reaction and the enforced production flux. % % Parameters % ---------- @@ -28,6 +30,10 @@ % modelAdapter, e.g. output/ecFSEOF_rxns.tsv). % modelAdapter : ModelAdapter % a loaded model adapter (default: the current default model adapter). +% corrThreshold : double +% minimum absolute Pearson correlation between the absolute flux of a +% reaction and the enforced production flux for the reaction to be a +% target (default 0.9). % % Returns % ------- @@ -40,25 +46,36 @@ % % - alpha : target production used for enforced objective limits (from minimum to maximum production). % - v_matrix : fluxes for each target reaction predicted and each alpha. -% - rxnTargets : a list with all reactions with fluxes that change consistently as target production increases. Contains ID, name, slope, gene rule, and equation. -% - transportTargets : a list with all transport reactions with fluxes that change consistently as target production increases. Contains ID, name, slope, gene rule, and equation. -% - geneTargets : a list with all selected targets that increase production. Contains gene, shortName, slope, action, and essentiality. +% - rxnTargets : a list with all reactions whose absolute flux is correlated with the enforced production, with an absolute correlation of at least corrThreshold. Contains ID, name, slope, gene rule, and equation. The slope is the absolute slope of the regression of the absolute flux on the enforced production flux. The list is ordered by the product of slope and correlation, highest first. +% - transportTargets : the same for transport reactions. +% - geneTargets : a list with all genes of the targets. Contains gene, shortName, slope, action, and essentiality. +% +% Only reactions with a gene rule are considered, and not those that are +% assigned to the standard enzyme. The enforced production reaction is +% included when it has a gene rule, as its flux follows the enforced +% production by definition. The reactions that draw enzyme usage +% (usage_prot_*) take part in the classification and are not reported. +% A reaction is over expressed (OE) when its flux increases with production, +% knocked out (KO) when it decreases and is at most 1e-8 at the last step, +% and knocked down (KD) when it decreases otherwise. % % Examples % -------- % % optional arguments may be given positionally or as name-value pairs: -% fseof = ecFSEOF(model,prodTargetRxn,csRxn,nSteps,outputFile,filePath,modelAdapter); +% fseof = ecFSEOF(model,prodTargetRxn,csRxn,nSteps,outputFile,filePath,modelAdapter,corrThreshold); % fseof = ecFSEOF(model,prodTargetRxn,csRxn,'nSteps',32); p = parseGECKOargs(varargin, { ... 'nSteps', 16; ... 'outputFile', false; ... 'filePath', []; ... - 'modelAdapter', []}); + 'modelAdapter', []; ... + 'corrThreshold',0.9}); nSteps = p.nSteps; outputFile = p.outputFile; filePath = p.filePath; -modelAdapter = p.modelAdapter; +modelAdapter = p.modelAdapter; +corrThreshold = p.corrThreshold; if isempty(modelAdapter) modelAdapter = ModelAdapterManager.getDefault(); @@ -133,65 +150,66 @@ v_matrix = v_matrix(withGR,:); rxnGeneM = model.rxnGeneMat(withGR,:); -% Filter out rxns that are always zero -zero_flux = ~all(abs(v_matrix(:,1:nSteps))==0,2); -target_rxns = target_rxns(zero_flux,:); -v_matrix = v_matrix(zero_flux,:); -rxnGeneM = rxnGeneM(zero_flux,:); - -% Identify those rxns that always increase or decrease, and calculate the -% slope as the difference in the flux when enforce objetive target -% production is set to 90%% of the maximum teorethical yield -% << v_matrix(i,nSteps-1) >> and the flux when the enforce objetive target -% production is set to the minimum << v_matrix(i,1) >> for the reaction i, -% divided by maxTarget-maxTarget/nSteps-1. +% Classify the reactions as in RAVEN's FSEOF. Slope and correlation are taken +% on the absolute flux against the enforced production flux, so a reversible +% reaction running backwards counts by its magnitude. A reaction is a target +% when some absolute flux exceeds 1e-10, the standard deviation of its +% absolute flux is above 0 and the absolute Pearson correlation reaches +% corrThreshold. It is suggested for over expression (OE) when the slope is +% positive, for knock out (KO) when the slope is negative and the absolute +% flux at the last step is at most 1e-8, and for knock down (KD) for any +% other negative slope. +absV = abs(v_matrix); slope_rxns = zeros(size(target_rxns)); -targets = logical(size(target_rxns)); +corr_rxns = zeros(size(target_rxns)); target_type = cell(size(target_rxns)); -for i = 1:length(target_rxns) - if issorted(abs(v_matrix(i,1:nSteps)),'strictascend') - % Those reactions that always increase while enforcing target - % production are suggested for Over Expression - targets(i) = true; - slope_rxns(i) = abs(v_matrix(i,nSteps-1)-v_matrix(i,1))/abs(maxTarget-maxTarget/nSteps-1); - target_type(i) = {'OE'}; - elseif issorted(abs(v_matrix(i,1:nSteps)),'strictdescend') - % Those reactions that always decrease while enforcing target - % production are suggested for KnockDown or KnockOut. KO are those - % reactions which have zero flux when enforcing target production - % to 90%% of the maximum theoretical yield. - targets(i) = true; - slope_rxns(i) = abs(v_matrix(i,nSteps-1)-v_matrix(i,1))/abs(maxTarget-maxTarget/nSteps-1); - if v_matrix(i,nSteps) == 0 - target_type(i) = {'KO'}; - else - target_type(i) = {'KD'}; +targets = false(size(target_rxns)); +if numel(alpha) > 1 && std(alpha) > 0 + for i = 1:length(target_rxns) + track = absV(i,:); + if any(track > 1e-10) && std(track) > 0 + fit = polyfit(alpha,track,1); + r = corrcoef(alpha,track); + if abs(r(1,2)) >= corrThreshold && fit(1) ~= 0 + targets(i) = true; + slope_rxns(i) = abs(fit(1)); + corr_rxns(i) = abs(r(1,2)); + if fit(1) > 0 + target_type(i) = {'OE'}; + elseif track(end) <= 1e-8 + target_type(i) = {'KO'}; + else + target_type(i) = {'KD'}; + end + end end end end -% Only keep those reaction that shows an increase or decrease pattern. +% Only keep those reactions that are classified as a target target_rxns = target_rxns(targets); v_matrix = v_matrix(targets,:); rxnGeneM = rxnGeneM(targets,:); slope_rxns = slope_rxns(targets); +corr_rxns = corr_rxns(targets); target_type = target_type(targets); -% Order from highest to lowest slope -[~,order] = sort(slope_rxns,'descend'); +% Order from highest to lowest score, the product of slope and correlation +[~,order] = sort(slope_rxns .* corr_rxns,'descend'); target_rxns = target_rxns(order); v_matrix = v_matrix(order,:); rxnGeneM = rxnGeneM(order,:); slope_rxns = slope_rxns(order); target_type = target_type(order); -% Filter out reactions with slope < quantile -non_zero_slope = slope_rxns > quantile(slope_rxns,0.75); -target_rxns = target_rxns(non_zero_slope); -v_matrix = v_matrix(non_zero_slope,:); -rxnGeneM = rxnGeneM(non_zero_slope,:); -slope_rxns = slope_rxns(non_zero_slope); -target_type = target_type(non_zero_slope); +% Exclude enzyme usage reactions: they are not reported, and they would otherwise +% enter the gene slopes and the choice of the action per gene +notUsage = ~startsWith(target_rxns,'usage_prot_'); +target_rxns = target_rxns(notUsage); +v_matrix = v_matrix(notUsage,:); +rxnGeneM = rxnGeneM(notUsage,:); +slope_rxns = slope_rxns(notUsage); +target_type = target_type(notUsage); % Create gene list of those connected to the remaining rxns genes = model.genes(sum(rxnGeneM,1) > 0); @@ -255,20 +273,22 @@ % Report values used to enforce objective target production fseof.alpha = alpha; -% Exclude enzyme usage reactions -toKeep = ~startsWith(target_rxns,'usage_prot_'); -rxnIdx = getIndexes(model,target_rxns(toKeep),'rxns'); +rxnIdx = getIndexes(model,target_rxns,'rxns'); % Report the fluxes at each alpha -fseof.v_matrix = array2table(v_matrix(toKeep,:), ... +fseof.v_matrix = array2table(v_matrix, ... 'VariableNames', string(alpha), 'RowNames', model.rxns(rxnIdx)); % Report reaction targets fseof.rxnTargets(:,1) = model.rxns(rxnIdx); fseof.rxnTargets(:,2) = model.rxnNames(rxnIdx); -fseof.rxnTargets(:,3) = num2cell(slope_rxns(toKeep)); +fseof.rxnTargets(:,3) = num2cell(slope_rxns); fseof.rxnTargets(:,4) = model.grRules(rxnIdx); -fseof.rxnTargets(:,5) = constructEquations(model,rxnIdx); +if isempty(rxnIdx) + fseof.rxnTargets(:,5) = cell(0,1); % constructEquations would give all reactions +else + fseof.rxnTargets(:,5) = constructEquations(model,rxnIdx); +end % Identify transport reactions; defined as reactions involving (at least) % one metabolite name in more than one compartment. diff --git a/src/flux_analysis/ecFVA.m b/src/flux_analysis/ecFVA.m index dfee0ca4..e4aee30a 100644 --- a/src/flux_analysis/ecFVA.m +++ b/src/flux_analysis/ecFVA.m @@ -42,7 +42,7 @@ 'runParallel', true}); runParallel = p.runParallel; -rxnIDs = regexprep(ecModel.rxns,'(_REV)?(_EXP_\d+)?',''); +rxnIDs = regexprep(ecModel.rxns,'(_REV)?(_EXP_\d+)?$',''); [rxnIDmap, convRxnID] = findgroups(rxnIDs); N = numel(convRxnID); diff --git a/src/gather_kcats/assignKcatValues.m b/src/gather_kcats/assignKcatValues.m index c7ec1c14..8188a3b0 100644 --- a/src/gather_kcats/assignKcatValues.m +++ b/src/gather_kcats/assignKcatValues.m @@ -21,6 +21,8 @@ % - genes : gene identifiers, matching model.ec.genes. % - substrates : substrate names, matching model.metNames. % - kcats : predicted kcat value in /sec. +% - kcatSource : optional, the source of each entry; entries with a zero +% kcat are discarded together with their source. % % Name-Value Arguments % -------------------- @@ -77,6 +79,9 @@ removeZero = kcatList.kcats == 0; kcatList.kcats(removeZero) = []; kcatList.rxns(removeZero) = []; +if isfield(kcatList,'kcatSource') + kcatList.kcatSource(removeZero) = []; +end % Map to model.ec.rxns [sanityCheck,idxInModel] = ismember(kcatList.rxns,model.ec.rxns); diff --git a/src/gather_kcats/fuzzyKcatMatching.m b/src/gather_kcats/fuzzyKcatMatching.m index 9956ee8e..d9787111 100644 --- a/src/gather_kcats/fuzzyKcatMatching.m +++ b/src/gather_kcats/fuzzyKcatMatching.m @@ -18,14 +18,20 @@ % % Name-Value Arguments % -------------------- -% ecRxns : logical +% ecRxns : logical or double % for which reactions (from model.ec.rxns) kcat values should be found, -% provided as logical vector with same length as model.ec.rxns -% (default: all reactions). +% provided as a logical vector with the same length as model.ec.rxns +% (the reactions are returned in model order), or as indices into +% model.ec.rxns (the reactions are returned in the order given, and an +% index given twice is returned twice). The literal [] selects all +% reactions, which is the default; any other empty selection, such as +% find(mask) without true entries, selects none and returns an empty +% kcatList. % modelAdapter : ModelAdapter % a loaded model adapter (default: the current default model adapter). % forceWClvl : double -% force a minimum wildcard level (default 0). +% number of extra levels, from the right, replaced by a wildcard (on top +% of any wildcard already in the code; default 0). % % Returns % ------- @@ -43,7 +49,8 @@ % - kcats : proposed kcat value in /sec. % - eccodes : as used to query BRENDA. % - wildcardLvl : which level of EC wild-card was necessary to find a match -% (0: w.x.y.z; 1: w.x.y.-; 2: w.x.-.-; 3: w.-.-.-). +% (0: w.x.y.z; 1: w.x.y.-; 2: w.x.-.-; 3: w.-.-.-; 4: -.-.-.-, which +% matches every EC number in BRENDA). % - origin : which level of specificity was necessary to find a match % (1: correct organism, correct substrate, kcat; 2: any organism, correct % substrate, kcat; 3: correct organism, any substrate, kcat; 4: any @@ -52,7 +59,9 @@ % 1, 2, 3, 5, 4, 6: a specific activity from the correct organism is % preferred over a kcat from any organism with any substrate. When a % reaction has multiple EC numbers, the EC with the fewest wildcards and -% the earliest-tried level is selected, and among those the highest kcat. +% the earliest-tried level is selected, and among those the highest kcat +% (the first one when several EC numbers share it); the wildcardLvl and +% origin are those of the selected EC. % % If a wildcard is used, origin levels 1 and 2 are ignored. The last digits % in the E.C. number indicate the substrate specificity, so if this should @@ -93,14 +102,19 @@ forceWClvl = p.forceWClvl; kcatAggregation = p.kcatAggregation; -if isempty(ecRxns) - ecRxns = true(numel(model.ec.rxns),1); -elseif isnumeric(ecRxns) - ecRxnsVec = false(numel(model.ec.rxns),1); - ecRxnsVec(ecRxns) = true; - ecRxns = ecRxnsVec; +% Get indices: [] selects all reactions, a logical vector the true entries (in +% model order) and indices are used as given (order and repeats kept). +if isnumeric(ecRxns) && isequal(size(ecRxns),[0 0]) + ecRxns = (1:numel(model.ec.rxns)).'; +elseif islogical(ecRxns) + ecRxns = find(ecRxns(:)); +else + ecRxns = ecRxns(:); + if any(ecRxns < 1 | ecRxns > numel(model.ec.rxns) | ecRxns ~= round(ecRxns)) + error('fuzzyKcatMatching:invalidEcRxns', ... + 'ecRxns must hold indices of model.ec.rxns.'); + end end -ecRxns=find(ecRxns); % Get indices instead of logical if isempty(modelAdapter) modelAdapter = ModelAdapterManager.getDefault(); @@ -118,6 +132,17 @@ error('No EC codes defined in model.ec.eccodes. First run getECfromGEM() and/or getECfromDatabase().') end eccodes = model.ec.eccodes(ecRxns); +if isempty(ecRxns) + %No reactions selected: nothing to match, and BRENDA is not read + kcatList.source = 'brenda'; + kcatList.rxns = cell(0,1); + kcatList.substrates = cell(0,1); + kcatList.kcats = zeros(0,1); + kcatList.eccodes = cell(0,1); + kcatList.wildcardLvl = NaN(0,1); + kcatList.origin = NaN(0,1); + return +end substrates = cell(numel(ecRxns),1); substrCoeffs = cell(numel(ecRxns),1); @@ -198,10 +223,7 @@ end %Apply force wildcard level -while forceWClvl > 0 - eccodes=regexprep(eccodes,'(.)*(\.\d+)(\.-)*$','$1\.-$3'); - forceWClvl = forceWClvl - 1; -end +eccodes = forceWildcards(eccodes,forceWClvl); PB = progressReport(mM, 'Gathering kcat values by fuzzy matching to BRENDA database'); %Main loop: @@ -235,6 +257,36 @@ end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%Replaces the rightmost level of each EC code that is not a wildcard yet with +%a wildcard, forceWClvl times, so that level 4 gives -.-.-.- . The codes of one +%entry are separated by ';'. A code that does not have four levels is kept. +function eccodes = forceWildcards(eccodes,forceWClvl) +if forceWClvl <= 0 + return +end +for i = 1:numel(eccodes) + if isempty(eccodes{i}) + continue + end + codes = strsplit(eccodes{i},';'); + for j = 1:numel(codes) + parts = strsplit(codes{j},'.'); + if numel(parts) ~= 4 + continue + end + for k = 1:forceWClvl + last = find(~strcmp(parts,'-'),1,'last'); + if isempty(last) + break + end + parts{last} = '-'; + end + codes{j} = strjoin(parts,'.'); + end + eccodes{i} = strjoin(codes,';'); +end +end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [kcat,dir,tot] =iterativeMatch(EC,subs,substrCoeff,i,KCATcell,dir,tot,... name,phylDist,org_index,SAcell,ECIndexIds,EcIndexIndices,kcatAggregation) %Will iteratively try to match the EC number to some registry in BRENDA, @@ -258,12 +310,20 @@ wc_num(k) = sum(EC{k}=='-'); else wild_num = sum(EC{k}=='-'); - if wild_num >= 4 + if numel(strsplit(EC{k},'.')) ~= 4 + %An EC code without four levels (for example '1.1.1') cannot be + %escalated level by level: stop with the no-match values, as + %geckopy does. + success = true; + elseif wild_num >= 4 %Already fully wildcarded ('-.-.-.-') and still no match -- %stop here instead of computing dot_pos(4-wild_num) with %wild_num==4, an out-of-range index. kcat(k)/origin(k) keep %their zero-initialised "no match" values. success = true; + elseif wild_num == 3 + %Last escalation: '-.-.-.-' matches every EC number in BRENDA + EC{k} = '-.-.-.-'; else dot_pos = [2 strfind(EC{k},'.')]; wc_text = '-.-.-.-'; @@ -284,11 +344,14 @@ %(higher) wildcard count back in whenever it happens to share best_origin, %defeating the wildcard-count filter entirely. best_pos = best_pos & (origin == best_origin); - max_pos = find(kcat == max(kcat(best_pos))); - wc_num = wc_num(max_pos(1)); - origin = origin(max_pos(1)); - matches = matches(max_pos(1)); - kcat = kcat(max_pos(1)); + %Highest kcat among those, the first one if several share it + best_idx = find(best_pos); + [~,j] = max(kcat(best_idx)); + max_pos = best_idx(j); + wc_num = wc_num(max_pos); + origin = origin(max_pos); + matches = matches(max_pos); + kcat = kcat(max_pos); %Update dir and tot: dir.org_s(i) = matches*(origin == 1); @@ -465,6 +528,8 @@ for j = 1:length(X) EC_indexes = [EC_indexes,EcIndexIndices{X(j)}]; end + elseif isempty(EC) %Not optimized, fully wildcarded: every entry matches + EC_indexes = 1:length(EC_cell); else %Not optimized for j=1:length(EC_cell) if strfind(EC_cell{j},EC)==1 diff --git a/src/gather_kcats/kcatsFromBrenda.m b/src/gather_kcats/kcatsFromBrenda.m index 64d29428..5df57419 100644 --- a/src/gather_kcats/kcatsFromBrenda.m +++ b/src/gather_kcats/kcatsFromBrenda.m @@ -19,9 +19,12 @@ % 'max' or 'median' (default: params.kcatAggregation of the model % adapter, or 'max'). % forceWClvl : double -% force a minimum wildcard level in the EC numbers (default 0). -% ecRxns : logical -% for which model.ec.rxns kcats are gathered (default: all). +% number of extra levels, from the right, replaced by a wildcard in the +% EC numbers (on top of any wildcard already in the code; default 0). +% ecRxns : logical or double +% for which model.ec.rxns kcats are gathered (default: all). A logical +% vector or indices of model.ec.rxns, see fuzzyKcatMatching; [] selects +% all reactions and any other empty selection is an error. % save : logical % write the source file (default true). % saveAs : char @@ -60,6 +63,10 @@ end end aggregation = resolveKcatAggregation(p.aggregation, modelAdapter); +if isSubset(p.ecRxns) && countSelected(p.ecRxns) == 0 + error('kcatsFromBrenda:emptySubset', ... + 'ecRxns selects no reactions; leave it out to gather for all reactions.'); +end if p.save file = fullfile(modelAdapter.params.path, 'data', 'kcats', [char(p.saveAs) '.tsv']); warnSubset(p.ecRxns, numel(model.ec.rxns), file); @@ -70,8 +77,8 @@ if p.save parameters = struct('aggregation', aggregation, 'force_wildcard_level', p.forceWClvl); - if ~isempty(p.ecRxns) - parameters.ec_rxns = countRxns(p.ecRxns); + if isSubset(p.ecRxns) + parameters.ec_rxns = numel(kcatList.rxns); end saveKcatSource(model, kcatList, p.saveAs, 'modelAdapter', modelAdapter, ... 'generator', 'kcatsFromBrenda', 'parameters', parameters, ... @@ -96,6 +103,19 @@ end end +function tf = isSubset(ecRxns) +% [] selects all reactions; any other value selects the reactions it names. +tf = ~(isnumeric(ecRxns) && isequal(size(ecRxns),[0 0])); +end + +function n = countSelected(ecRxns) +if islogical(ecRxns) + n = sum(ecRxns); +else + n = numel(ecRxns); +end +end + function n = countRxns(ecRxns) if islogical(ecRxns) n = sum(ecRxns); @@ -106,7 +126,7 @@ function warnSubset(ecRxns, nRxns, file) % A subset written under the default name replaces the full file. -if ~isempty(ecRxns) && countRxns(ecRxns) < nRxns && isfile(file) +if isSubset(ecRxns) && countRxns(ecRxns) < nRxns && isfile(file) printOrange(sprintf(['%s is replaced by kcats for %d reactions only. Give saveAs ' ... 'to keep it.\n'], file, countRxns(ecRxns))); end diff --git a/src/gather_kcats/mergeDLKcatAndFuzzyKcats.m b/src/gather_kcats/mergeDLKcatAndFuzzyKcats.m index 4bb1e0dd..364efeac 100644 --- a/src/gather_kcats/mergeDLKcatAndFuzzyKcats.m +++ b/src/gather_kcats/mergeDLKcatAndFuzzyKcats.m @@ -28,7 +28,8 @@ % -------------------- % topOriginLimit : double % origin limit for prioritized BRENDA matches. Origin is explained in -% more detail below (default 6). +% more detail below (default 4: BRENDA matches of origin 5 and 6, from +% specific activities, rank after DLKcat; use 6 to rank them above). % bottomOriginLimit : double % origin limit for low priority BRENDA matches. Origin is explained in % more detail below (default 6). @@ -74,7 +75,7 @@ % mergeKcats, readDLKcatOutput, fuzzyKcatMatching, assignKcatValues p = parseGECKOargs(varargin, { ... - 'topOriginLimit', 6; ... + 'topOriginLimit', 4; ... 'bottomOriginLimit', 6; ... 'wildcardLimit', 3}); topOriginLimit = p.topOriginLimit; diff --git a/src/gather_kcats/mergeKcats.m b/src/gather_kcats/mergeKcats.m index ad0f65c5..483c51f3 100644 --- a/src/gather_kcats/mergeKcats.m +++ b/src/gather_kcats/mergeKcats.m @@ -33,7 +33,10 @@ % Name-Value Arguments % -------------------- % topOriginLimit : double -% origin limit for database_top (default 6). +% origin limit for database_top (default 4). Origins 5 and 6 are specific +% activity times molecular weight kcats, which do not always beat a +% prediction; they rank after the other sources and are used where +% nothing else exists, or above DLKcat with topOriginLimit 6. % bottomOriginLimit : double % origin limit for database_bottom (default 6). % wildcardLimit : double @@ -87,7 +90,7 @@ % mergeDLKcatAndFuzzyKcats, readDLKcatOutput, fuzzyKcatMatching, assignKcatValues p = parseGECKOargs(varargin, { ... - 'topOriginLimit', 6; ... + 'topOriginLimit', 4; ... 'bottomOriginLimit', 6; ... 'wildcardLimit', 3; ... 'databaseSources', {'brenda','sabio_rk'}}); diff --git a/src/get_enzyme_data/findECInDB.m b/src/get_enzyme_data/findECInDB.m index f23c8dbb..50abdd05 100644 --- a/src/get_enzyme_data/findECInDB.m +++ b/src/get_enzyme_data/findECInDB.m @@ -4,10 +4,14 @@ % For each gene, finds its matching database entries and takes the EC number % of the lightest-weight matching protein when a single EC number applies to % that gene; genes with several distinct EC numbers among their matches are -% also reported as conflicts (the first match is still used for EC). When -% gene_set has more than one gene (an enzyme complex), the numbers found for -% each gene are combined by intersection where they share a match (wildcards -% such as "-" are accounted for), otherwise by union. +% also reported as conflicts (the first match is still used for EC). Genes +% without an EC number are left out. When gene_set has more than one gene +% (an enzyme complex), the numbers found for each gene are combined by +% intersection where they share a match (a wildcard such as "1.1.1.-" matches +% the specific codes of its prefix), otherwise by union. The EC set of every +% gene set, one gene or a complex, is cleaned of a wildcard code such as +% "1.1.1.-" when a more specific code of the same prefix ("1.1.1.4") is +% present; only codes with four levels take part, and empty codes are removed. % % Parameters % ---------- @@ -88,25 +92,33 @@ end end -%EC: Find union and intersection between all units (only applies for -%complexes, i.e. length(EC_set) > 1): -uni_EC = strsplit(EC_set{1},' '); -int_EC = uni_EC; -for j = 2:length(EC_set) - other_EC = strsplit(EC_set{j},' '); - if isempty(uni_EC) - uni_EC = other_EC; - int_EC = other_EC; - elseif ~isempty(other_EC) - uni_EC = compare_wild([uni_EC other_EC]); - int_EC = intersection(int_EC,other_EC); - end +%EC: Split the codes of each gene, leave out the genes without any, and +%combine the rest (union and intersection only apply for complexes, i.e. more +%than one gene with codes): +tokens = cell(size(EC_set)); +for j = 1:length(EC_set) + codes = strsplit(EC_set{j},' '); + tokens{j} = codes(~cellfun(@isempty,codes)); end -%Use the intersection found, if any. If not, use the union: -if isempty(int_EC) - EC_set = uni_EC; +tokens = tokens(~cellfun(@isempty,tokens)); +if isempty(tokens) + EC_set = {}; +elseif isscalar(tokens) + EC_set = dropWildcardParents(tokens{1}); else - EC_set = int_EC; + uni_EC = tokens{1}; + int_EC = tokens{1}; + for j = 2:length(tokens) + uni_EC = dropWildcardParents([uni_EC tokens{j}]); + int_EC = wildcardIntersection(int_EC,tokens{j}); + end + %Use the intersection found, if any. If not, use the union: + if isempty(int_EC) + EC_set = uni_EC; + else + EC_set = int_EC; + end + EC_set = dropWildcardParents(EC_set); end %Add new codes as new possible isozymes and remove the initial "EC" text: @@ -118,7 +130,11 @@ end %EC_set = cellfun(@(x) extractAfter(x,2),EC_set,'UniformOutput',false);%remove 2 first chars -EC = strjoin(EC_set,';'); +if isempty(EC_set) + EC = ''; +else + EC = strjoin(EC_set,';'); +end %Delete repeated Uniprots (and the corresponding ECs): @@ -127,58 +143,70 @@ end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -function int_EC = intersection(prev_EC,new_EC) +function int_EC = wildcardIntersection(prev_EC,new_EC) %Finds the common elements between two cell arrays, if any. Also considers %wildcards (e.g. if 'EC1.1.1.1' is in one array and 'EC1.1.1.-' is in the -%other one, then 'EC1.1.1.1' is added to the intersection). +%other one, then 'EC1.1.1.1' is added to the intersection). A pair adds its +%code only when one code of the pair is the wildcard parent of the other, or +%the two are equal four-level codes. int_EC = {}; for i = 1:length(prev_EC) for j = 1:length(new_EC) - new_int = compare_wild({prev_EC{i} new_EC{j}}); - if length(new_int) == 1 - int_EC = [int_EC new_int]; + survivors = dropWildcardParents({prev_EC{i} new_EC{j}}); + if length(survivors) == 1 + int_EC = [int_EC survivors]; end end end -%When prev_EC carries both a code and its own wildcard parent (e.g. -%'EC1.1.1.1' and 'EC1.1.1.-'), both independently pair-match the same -%new_EC entry above and the loop appends the resolved code twice; dedupe. -int_EC = unique(int_EC,'stable'); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -function int_EC = compare_wild(EC) -%Goes through a cell array of EC numbers, and erases any repetitions, -%considering also wildcards (e.g. will erase 'EC1.2.3.-' if 'EC1.2.3.4' is -%already present). - -%Trim all EC numbers of wild cards (e.g. 'EC1.2.-.-' -> 'EC1.2.'): -EC_trimmed = EC; -for i = 1:length(EC) - %Add a final dot for avoiding issues (e.g. 'EC1.1.1.1' & 'EC1.1.1.12'): - ECi = [EC{i} '.']; - pos = strfind(ECi,'-'); - if ~isempty(pos) - ECi = ECi(1:pos(1)-1); +function int_EC = dropWildcardParents(EC) +%Removes empty codes and every code that is a wildcard parent of another +%code in the cell array (e.g. 'EC1.2.3.-' when 'EC1.2.3.4' is present). Only +%codes with four levels are compared, so 'EC1.1.1' is neither a parent nor a +%child of 'EC1.1.1.4'. Equal codes collapse to the first. +EC = EC(~cellfun(@isempty,EC) & ~strcmp(EC,'EC')); +keep = true(1,numel(EC)); +for i = 1:numel(EC) + if ~keep(i) + continue; end - EC_trimmed{i} = ECi; -end - -%Compare all EC numbers between them to find if 1 fits in the other: -non_repeated = true(1,length(EC)); -for i = 1:length(EC)-1 - for j = i+1:length(EC) - ECi = EC_trimmed{i}; - ECj = EC_trimmed{j}; - %If ECj fits in ECi then ECj can be disregarded: - if strfind(ECi,ECj) - non_repeated(j) = false; - %Else, if ECi fits in ECj then ECi can be disregarded: - elseif strfind(ECj,ECi) - non_repeated(i) = false; + for j = i+1:numel(EC) + if ~keep(j) + continue; + end + iInJ = subsumes(EC{i},EC{j}); + jInI = subsumes(EC{j},EC{i}); + if iInJ && jInI + keep(j) = false; + elseif iInJ + keep(i) = false; + break; + elseif jInI + keep(j) = false; end end end - -int_EC = EC(non_repeated); - +int_EC = EC(keep); +end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function tf = subsumes(general,specific) +%True if the code general is equal to or a wildcard parent of the code +%specific; both carry the "EC" prefix and need four levels. +general = strsplit(extractAfter(general,2),'.'); +specific = strsplit(extractAfter(specific,2),'.'); +tf = false; +if numel(general) ~= 4 || numel(specific) ~= 4 + return; +end +for k = 1:4 + if strcmp(general{k},'-') + tf = true; + return; + end + if strcmp(specific{k},'-') || ~strcmp(general{k},specific{k}) + return; + end +end +tf = true; end diff --git a/src/io/loadProtData.m b/src/io/loadProtData.m index bc3ea4b2..e4e3be27 100644 --- a/src/io/loadProtData.m +++ b/src/io/loadProtData.m @@ -143,16 +143,16 @@ end % Filter by RSD remData = (std(condAbund,0,2,'omitnan')./mean(condAbund,2,'omitnan'))>maxRSD; - condAbund(remData) = nan; + condAbund(remData,:) = nan; % Add stdevs condAbund = mean(condAbund,2,'omitnan')+(addStdevs*std(condAbund,0,2,'omitnan')); % Filter by minVal remData = mean(condAbund,2,'omitnan') 0) + printOrange('WARNING: No protein usage reaction carries flux, tuning cannot continue.\n'); + break; + end [~,sel] = max(drawFluxes); % since bounds 0 to 1000 %Now get the metabolite metSel = m.S(:,sel) > 0; % positive coeff (forward usage) @@ -146,7 +158,11 @@ iteration = 1; while true res = solveLP(m,0); %skip parsimonius, only takes time - if (lastGrowth == res.f) + if isempty(res.f) + printOrange('WARNING: The model has no feasible solution with the increased kcats, tuning stops.\n'); + break; + end + if ~isnan(lastGrowth) && res.f <= lastGrowth + growthStallTol printOrange('No growth increase from increased kcats - check if the constraints on the uptake reactions are too tight.\n'); break; end @@ -162,6 +178,10 @@ %find the highest protein usage flux protPoolStoich = m.S(strcmp(m.mets, 'prot_pool'),:).'; protPoolStoich(idxToIgnore) = 0; + if ~any(res.x .* protPoolStoich < 0) + printOrange('WARNING: No reaction draws on the protein pool, tuning cannot continue.\n'); + break; + end [~,sel] = min(res.x .* protPoolStoich); %max consumption kcatList = [kcatList, sel]; rxn = m.rxns(sel.'); diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index d027bbbe..d14c6743 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -1,24 +1,24 @@ Function/File (in src) Verification Comments BrendaParser/* Not explicitly tested - The code has existed for a long time -change_model/applyComplexData.m tc0003, tc0004 -change_model/applyKcatConstraints.m tc0011 +change_model/applyComplexData.m tc0003, tc0004, tc0103 +change_model/applyKcatConstraints.m tc0011, tc0215, tc0235, tc0238, tc0250 change_model/findMetSmiles.m tc0012 change_model/getComplexData.m tc0003, tc0004 - File download not covered by test cases but is tested manually and works. -change_model/makeEcModel.m tc0001, tc0002 +change_model/makeEcModel.m tc0001, tc0002, tc0210, tc0212, tc0233, tc0234 change_model/setProtPoolSize.m tc0005 -change_model/getKcatAcrossIsoenzymes.m tc0011 -change_model/applyCustomKcats.m -gather_kcats/fuzzyKcatMatching.m tc0010 +change_model/getKcatAcrossIsoenzymes.m tc0011, tc0252 +change_model/applyCustomKcats.m tc0227, tc0228, tc0229, tc0230 +gather_kcats/fuzzyKcatMatching.m tc0010, tc0223, tc0224, tc0225, tc0226, tc0231, tc0239 gather_kcats/getStandardKcat.m -gather_kcats/mergeDlkcatAndFuzzyKcats.m tc0011 +gather_kcats/mergeDlkcatAndFuzzyKcats.m tc0011; tc0241 gather_kcats/readDLKcatOutput.m tcm0001, tcm0002 - The test doesn't really check the function, only that some output is read gather_kcats/runDLKcat.m tcm0001, tcm0002 - The test doesn't really check the function, only that we get something -gather_kcats/assignKcatValues.m tc0011 - we did not test all possible parameterizations here +gather_kcats/assignKcatValues.m tc0011 - we did not test all possible parameterizations here; also tc0104 gather_kcats/writeDLKcatInput.m partly by tc0011 gather_kcats/getEcRxnKeys.m tc0070, tc0076 gather_kcats/saveKcatSource.m tc0071, tc0073, tc0076, tc0078 gather_kcats/loadKcatSource.m tc0071, tc0072, tc0073, tc0076, tc0077, tc0078 -gather_kcats/kcatsFromBrenda.m tc0074, tc0078 +gather_kcats/kcatsFromBrenda.m tc0074, tc0078, tc0226 gather_kcats/kcatsFromDLKcat.m tc0074, tc0078 - reads an existing DLKcat output; running DLKcat needs Docker gather_kcats/kcatsFromOKP.m tc0075 - reads a stored result; submitting needs the OKP service gather_kcats/collectKcats.m tc0079, tc0080, tc0081, tc0083, tc0087, tc0088, tc0089, tc0090, tc0098 @@ -32,22 +32,28 @@ gather_kcats/kcatsFromHomologs.m tc0095, tc0096 gather_kcats/extractEnzymeSubstratePairs.m tc0096, tc0015 gather_kcats/convertCustomKcats.m tc0091, tc0092, tc0094 change_model/upgradeEcModel.m tc0093, tc0094, tc0099 +io/loadProtData.m tc0101, tc0102 +change_model/addNewRxnsToEC.m tc0211, tc0232, tc0237 +change_model/mapRxnsToConv.m tc0213 +flux_analysis/ecFVA.m tc0214 +flux_analysis/ecFSEOF.m tc0100, tc0236 utilities/upgradeGECKOproject.m tc0094, tc0099 gather_kcats/readOpenKineticsPredictorOutput.m tc0030, tc0069 gather_kcats/reviewKcatAssignment.m Has only been tested by running it manually, not tested by test cases gather_kcats/ecStatsFromBrenda.m Not verified against real BRENDA data in this environment -- see its own TODO comment +gather_kcats/mergeKcats.m tc0240 get_enzyme_data/calculateMW.m Not explicitly tested - The code has existed for a long time -get_enzyme_data/findECInDB.m Called from getECfromDatabase - tested when that function is tested. +get_enzyme_data/findECInDB.m Called from getECfromDatabase - tested when that function is tested; also tc0242, tc0243, tc0244, tc0245 get_enzyme_data/getECfromDatabase.m tc0007 - does not test download of the databases - this is tested in the manual workflows for Yeast-GEM and Human-GEM get_enzyme_data/getECfromGEM.m tc0006 get_enzyme_data/getECstring.m Called from findECInDB - tested when that function is tested. -download_databases/loadBRENDAdata.m Called from fuzzyKcatMatching - tested when that function is tested. +download_databases/loadBRENDAdata.m Called from fuzzyKcatMatching - tested when that function is tested; also tc0105, tc0220, tc0221, tc0222 download_databases/fetchBrendaAsset.m tc0085 download_databases/loadDatabases.m Called from getECfromDatabase - tested when that function is tested. kcat_tuning/evotune/* Has only been tested by running it manually, not tested by test cases kcat_tuning/fillCarbonNum.m tc0055 (deprecated-alias check only, not a real usage test) kcat_tuning/findTopLimitations.m Not explicitly tested - The code has existed for a long time -kcat_tuning/sensitivityTuning.m tcm0001 +kcat_tuning/sensitivityTuning.m tcm0001, tc0251, tc0253, tc0254 kcat_tuning/sigmaFitter.m Not explicitly tested - The code has existed for a long time kcat_tuning/topUsedEnzymes.m Not explicitly tested - The code has existed for a long time Not sure this works anymore - written for the old model structure? kcat_tuning/truncateValues.m Not explicitly tested - The code has existed for a long time diff --git a/test/unit_tests/ecModelBuildTests.m b/test/unit_tests/ecModelBuildTests.m new file mode 100644 index 00000000..523ecd7b --- /dev/null +++ b/test/unit_tests/ecModelBuildTests.m @@ -0,0 +1,206 @@ +%run this test case with the command +%results = runtests('ecModelBuildTests.m') + +function tests = ecModelBuildTests + tests = functiontests(localfunctions); +end + +function testLightNestedGeneRuleExpandsToIsozymes_tc0210(testCase) + % 'G1 and (G2 or G3)' is the two isozymes {G1,G2} and {G1,G3}; the + % brackets must not be dropped, which would make the second one {G3}. + adapter = getAdapter(); + model = getGeckoTestModel(); + model = changeGrRules(model, {'R2'}, {'G1 and (G2 or G3)'}, true); + ecModel = makeEcModel(model, true, adapter); + verifyEqual(testCase, sort(coupledGenes(ecModel, '001_R2')), {'G1';'G2'}) + verifyEqual(testCase, sort(coupledGenes(ecModel, '002_R2')), {'G1';'G3'}) + verifyEqual(testCase, sum(strcmp(ecModel.ec.rxns, '003_R2')), 0) +end + +function testAddNewRxnsToECNestedRuleAndGeneIdWithOr_tc0211(testCase) + % The gene ids of the rule are whole words, so 'orfX' is a gene and not + % the operator 'or' followed by 'fX', and the nested rule is split into + % one reaction per isozyme: 'orfX and GNEW2' and 'orfX and GNEW3'. + adapter = getAdapter(); + model = getGeckoTestModel(); + ecModel = makeEcModel(model, false, adapter); + + newRxns.rxns = {'RNEWA'}; + newRxns.rxnNames = {'RNEWA'}; + newRxns.equations = {'m1[c] => e2[e]'}; + newRxns.grRules = {'orfX and (GNEW2 or GNEW3)'}; + newEnzymes.enzymes = {'ENEW1';'ENEW2';'ENEW3'}; + newEnzymes.genes = {'orfX';'GNEW2';'GNEW3'}; + newEnzymes.mw = [15000;16000;17000]; + + ecModel = addNewRxnsToEC(ecModel, newRxns, newEnzymes, adapter); + verifyEqual(testCase, ecModel.grRules(strcmp(ecModel.rxns, 'RNEWA_EXP_1')), {'orfX and GNEW2'}) + verifyEqual(testCase, ecModel.grRules(strcmp(ecModel.rxns, 'RNEWA_EXP_2')), {'orfX and GNEW3'}) +end + +function testMakeEcModelCouplesTheKnownSubunitsOfAComplex_tc0212(testCase) + % GX is in neither database: the reaction is still coupled to the enzyme + % of G4, as the databases know it. + adapter = getAdapter(); + model = getGeckoTestModel(); + model = changeGrRules(model, {'R3'}, {'G4 and GX'}, true); + ecModel = makeEcModel(model, false, adapter); + verifyEqual(testCase, coupledGenes(ecModel, 'R3'), {'G4'}) +end + +function testMapRxnsToConvStripsOnlyTheSuffixes_tc0213(testCase) + % R_REVX is a conventional reaction whose id contains _REV; its reverse + % variant R_REVX_REV carries flux that maps back to it with the opposite sign. + adapter = getAdapter(); + model = getGeckoTestModel(); + model.rxns{strcmp(model.rxns, 'R1')} = 'R_REVX'; + ecModel = makeEcModel(model, false, adapter); + fluxVect = zeros(numel(ecModel.rxns), 1); + fluxVect(strcmp(ecModel.rxns, 'R_REVX')) = 5; + fluxVect(strcmp(ecModel.rxns, 'R_REVX_REV')) = 2; + mappedFlux = mapRxnsToConv(ecModel, model, fluxVect); + verifyEqual(testCase, mappedFlux(strcmp(model.rxns, 'R_REVX')), 3) +end + +function testEcFVAStripsOnlyTheSuffixes_tc0214(testCase) + % The same reaction id in ecFVA, which groups the variants of each + % conventional reaction by id. + adapter = getAdapter(); + model = getGeckoTestModel(); + model.rxns{strcmp(model.rxns, 'R1')} = 'R_REVX'; + ecModel = makeEcModel(model, false, adapter); + [minFlux, maxFlux] = ecFVA(ecModel, model, 'runParallel', false); + verifySize(testCase, minFlux, [numel(model.rxns) 1]) + verifySize(testCase, maxFlux, [numel(model.rxns) 1]) +end + +function testApplyKcatConstraintsLightSkipsIsozymeWithoutEnzymes_tc0215(testCase) + % R2 has two isozymes; the second one has no enzyme but a kcat. It cannot + % be the cheapest isozyme: the cost of R2 is the one of the first. + adapter = getAdapter(); + model = getGeckoTestModel(); + lecModel = makeEcModel(model, true, adapter); + lecModel.ec.kcat(:) = 10; + first = strcmp(lecModel.ec.rxns, '001_R2'); + second = strcmp(lecModel.ec.rxns, '002_R2'); + lecModel.ec.rxnEnzMat(second,:) = 0; + lecModel = applyKcatConstraints(lecModel); + expected = -(lecModel.ec.rxnEnzMat(first,:) * lecModel.ec.mw) / 10 / 3600; + actual = full(lecModel.S(strcmp(lecModel.mets, 'prot_pool'), strcmp(lecModel.rxns, 'R2'))); + verifyGreaterThan(testCase, abs(expected), 0) + verifyEqual(testCase, actual, expected, 'AbsTol', 1e-10) +end + + +function testAddNewRxnsToECGeneIdsContainingAnd_tc0232(testCase) + % Gene ids that contain 'and' (Sand1, andA, GandH) are whole words: the + % new reaction is coupled to exactly the enzymes of the genes in its rule. + adapter = getAdapter(); + model = getGeckoTestModel(); + ecModel = makeEcModel(model, false, adapter); + + newRxns.rxns = {'RNEWB'}; + newRxns.rxnNames = {'RNEWB'}; + newRxns.equations = {'m1[c] => e2[e]'}; + newRxns.grRules = {'Sand1 and (andA or GandH)'}; + newEnzymes.enzymes = {'ENEWA';'ENEWB';'ENEWC'}; + newEnzymes.genes = {'Sand1';'andA';'GandH'}; + newEnzymes.mw = [15000;16000;17000]; + + ecModel = addNewRxnsToEC(ecModel, newRxns, newEnzymes, adapter); + verifyEqual(testCase, sort(coupledGenes(ecModel, 'RNEWB_EXP_1')), {'Sand1';'andA'}) + verifyEqual(testCase, sort(coupledGenes(ecModel, 'RNEWB_EXP_2')), {'GandH';'Sand1'}) +end + +function testAddNewRxnsToECUpperCaseOperators_tc0237(testCase) + % 'AND' and 'OR' are operators like 'and' and 'or', also in a rule with a + % single clause and on a reversible reaction, where the rule is copied as + % it is to the reverse reaction. + adapter = getAdapter(); + model = getGeckoTestModel(); + ecModel = makeEcModel(model, false, adapter); + + newRxns.rxns = {'RNEWC';'RNEWD'}; + newRxns.rxnNames = {'RNEWC';'RNEWD'}; + newRxns.equations = {'m1[c] <=> e2[e]';'m1[c] => e2[e]'}; + newRxns.grRules = {'Sand1 AND andA';'Sand1 OR andA'}; + newEnzymes.enzymes = {'ENEWA';'ENEWB'}; + newEnzymes.genes = {'Sand1';'andA'}; + newEnzymes.mw = [15000;16000]; + + ecModel = addNewRxnsToEC(ecModel, newRxns, newEnzymes, adapter); + verifyEqual(testCase, ecModel.grRules(strcmp(ecModel.rxns, 'RNEWC')), {'Sand1 and andA'}) + verifyEqual(testCase, sort(coupledGenes(ecModel, 'RNEWC')), {'Sand1';'andA'}) + verifyEqual(testCase, sort(coupledGenes(ecModel, 'RNEWC_REV')), {'Sand1';'andA'}) + verifyEqual(testCase, coupledGenes(ecModel, 'RNEWD_EXP_1'), {'Sand1'}) + verifyEqual(testCase, coupledGenes(ecModel, 'RNEWD_EXP_2'), {'andA'}) +end + +function testLightGeneRuleOperatorsAreCaseInsensitive_tc0233(testCase) + % 'AND' and 'OR' are operators like 'and' and 'or'. + adapter = getAdapter(); + model = getGeckoTestModel(); + model.grRules{strcmp(model.rxns, 'R2')} = 'G1 AND (G2 OR G3)'; + ecModel = makeEcModel(model, true, adapter); + verifyEqual(testCase, sort(coupledGenes(ecModel, '001_R2')), {'G1';'G2'}) + verifyEqual(testCase, sort(coupledGenes(ecModel, '002_R2')), {'G1';'G3'}) +end + +function testLightReactionWithEmptyRuleKeepsOneRow_tc0234(testCase) + % A reaction that has gene entries but an empty gene rule keeps one row + % in ec.rxns, without enzymes. + adapter = getAdapter(); + model = getGeckoTestModel(); + model.grRules{strcmp(model.rxns, 'R2')} = ''; + ecModel = makeEcModel(model, true, adapter); + row = strcmp(ecModel.ec.rxns, '001_R2'); + verifyEqual(testCase, sum(row), 1) + verifyEmpty(testCase, coupledGenes(ecModel, '001_R2')) +end + +function testApplyKcatConstraintsLightSkipsNonPositiveKcat_tc0235(testCase) + % A negative kcat is not a usable kcat: the cost of R2 is the one of the + % isozyme with a positive kcat. + adapter = getAdapter(); + model = getGeckoTestModel(); + lecModel = makeEcModel(model, true, adapter); + lecModel.ec.kcat(:) = 10; + first = strcmp(lecModel.ec.rxns, '001_R2'); + second = strcmp(lecModel.ec.rxns, '002_R2'); + lecModel.ec.kcat(first) = -5; + lecModel = applyKcatConstraints(lecModel); + expected = -(lecModel.ec.rxnEnzMat(second,:) * lecModel.ec.mw) / 10 / 3600; + actual = full(lecModel.S(strcmp(lecModel.mets, 'prot_pool'), strcmp(lecModel.rxns, 'R2'))); + verifyGreaterThan(testCase, abs(expected), 0) + verifyEqual(testCase, actual, expected, 'AbsTol', 1e-10) +end + +function testApplyKcatConstraintsLightSkipsEnzymeWithoutMass_tc0238(testCase) + % An isozyme whose enzymes have no mass gives a cost of 0, which is not a + % usable cost: the cost of R2 is the one of the isozyme with a mass. + adapter = getAdapter(); + model = getGeckoTestModel(); + lecModel = makeEcModel(model, true, adapter); + lecModel.ec.kcat(:) = 10; + first = strcmp(lecModel.ec.rxns, '001_R2'); + second = strcmp(lecModel.ec.rxns, '002_R2'); + onlyFirst = find(full(lecModel.ec.rxnEnzMat(first,:)) > 0 & full(lecModel.ec.rxnEnzMat(second,:)) == 0); + verifyNotEmpty(testCase, onlyFirst) + lecModel.ec.mw(onlyFirst) = 0; + lecModel = applyKcatConstraints(lecModel); + expected = -full(lecModel.ec.rxnEnzMat(second,:) * lecModel.ec.mw) / 10 / 3600; + actual = full(lecModel.S(strcmp(lecModel.mets, 'prot_pool'), strcmp(lecModel.rxns, 'R2'))); + verifyGreaterThan(testCase, abs(expected), 0) + verifyEqual(testCase, actual, expected, 'AbsTol', 1e-10) +end + +function adapter = getAdapter() +geckoPath = findGECKOroot; +adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM','TestGEMAdapter.m')); +end + +function genes = coupledGenes(ecModel, ecRxn) +% The genes of the enzymes that ec reaction ecRxn is coupled to. +row = strcmp(ecModel.ec.rxns, ecRxn); +genes = ecModel.ec.genes(logical(ecModel.ec.rxnEnzMat(row,:))); +end diff --git a/test/unit_tests/geckoCoreFunctionTests.m b/test/unit_tests/geckoCoreFunctionTests.m index d0a0e7fd..308b830f 100644 --- a/test/unit_tests/geckoCoreFunctionTests.m +++ b/test/unit_tests/geckoCoreFunctionTests.m @@ -1270,28 +1270,26 @@ function testFindECInDBIntersectionDedupesWildcardPair_tc0032(testCase) function testFuzzyKcatMatchingWildcardExhaustionDoesNotCrash_tc0033(testCase) % Verifies that when an EC code has no BRENDA match at any wildcard - % level, including the fully wildcarded '-.-.-.-', fuzzyKcatMatching's - % escalation loop terminates cleanly and reports "no match" (kcat 0, - % origin NaN) instead of erroring. + % level, including the fully wildcarded '-.-.-.-' (which matches every + % EC number in BRENDA), fuzzyKcatMatching's escalation loop terminates + % cleanly and reports "no match" (kcat 0, origin NaN) instead of + % erroring. The fixture BRENDA kcat table has no rows, so nothing can match. geckoPath = findGECKOroot; adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM', 'TestGEMAdapter.m')); model = getGeckoTestModel(); ecModel = makeEcModel(model, false, adapter); ecModel = getECfromGEM(ecModel); - % An EC class unrelated, at every wildcard level, to anything in the - % fixture BRENDA file below. ecModel.ec.eccodes{strcmp(ecModel.ec.rxns, 'R2_EXP_1')} = '9.9.9.9'; scratchDir = fullfile(tempname); brendaDir = fullfile(scratchDir, 'data'); mkdir(brendaDir); cleanupDir = onCleanup(@() rmdir(scratchDir, 's')); - % Deliberately unrelated to '9.9.9.9' at every wildcard level, so the - % escalation loop runs out of levels without ever matching. + % No rows at all, so the escalation loop runs out of levels without + % ever matching. fid = fopen(fullfile(brendaDir, 'kcat.tsv'), 'w'); fprintf(fid, '# BRENDA release test-fixture - kcat in 1/s\n'); fprintf(fid, 'ec_code\tsubstrate\torganism\tkcat_max\tkcat_median\tn\treferences\n'); - fprintf(fid, '5.5.5.5\tm1\ttestus testus\t42\t42\t1\t*\n'); fclose(fid); fid = fopen(fullfile(brendaDir, 'mw.tsv'), 'w'); fprintf(fid, '# BRENDA release test-fixture - molecular weight in g/mol\n'); @@ -2192,3 +2190,76 @@ function testGECKOInstallerIsOlderVersion_tc0059(testCase) verifyFalse(testCase, older('4.0.0', '4.0.0b1')) verifyFalse(testCase, older('4.0.0b1', '4.0.0b2')) end + +function testEcFSEOFSelectsTargetsByCorrelation_tc0100(testCase) + % Verifies ecFSEOF classifies a reaction as a target when the absolute + % flux correlates with the enforced production flux (default + % corrThreshold 0.9), as RAVEN's FSEOF does, and does not cut the list to + % the top quartile by slope. With R3 enforced, the flux of R5 falls as + % the shared enzyme pool goes to R3, and R3 itself follows the enforced + % flux; both are targets, the ones without a gene rule are not. + geckoPath = findGECKOroot; + adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM', 'TestGEMAdapter.m')); + model = getGeckoTestModel(); + model.geneShortNames = model.genes; + ecModel = makeEcModel(model, false, adapter); + ecModel = getECfromGEM(ecModel); + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R2_EXP_1')) = 50; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R2_EXP_2')) = 10; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R2_REV_EXP_1')) = 50; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R2_REV_EXP_2')) = 10; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R3')) = 20; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R5')) = 100; + ecModel.ec.source(:) = {'manual'}; + ecModel = applyKcatConstraints(ecModel); + ecModel = setProtPoolSize(ecModel, [], [], [], adapter); + adapter.params.bioRxn = 'R5'; + + fseof = ecFSEOF(ecModel, 'R3', 'S1', 'nSteps', 6, 'modelAdapter', adapter); + targets = sort([fseof.rxnTargets(:,1); fseof.transportTargets(:,1)]); + verifyEqual(testCase, targets, {'R3';'R5'}) + verifyEqual(testCase, numel(fseof.alpha), 6) + + % The slope column is the absolute regression slope, highest score first + verifyTrue(testCase, all(cell2mat(fseof.rxnTargets(:,3)) > 0)) + + % The default corrThreshold is 0.9 + sameAsDefault = ecFSEOF(ecModel, 'R3', 'S1', 'nSteps', 6, 'modelAdapter', adapter, 'corrThreshold', 0.9); + verifyEqual(testCase, sameAsDefault.rxnTargets, fseof.rxnTargets) + + % The correlation of a reaction cannot exceed 1: a threshold above 1 gives no targets + none = ecFSEOF(ecModel, 'R3', 'S1', 'nSteps', 6, 'modelAdapter', adapter, 'corrThreshold', 1.5); + verifyEmpty(testCase, none.rxnTargets) + verifyEmpty(testCase, none.transportTargets) + verifyEmpty(testCase, none.geneTargets) +end + +function testEcFSEOFGeneTargetsComeFromReportedReactions_tc0236(testCase) + % The usage reactions of the enzymes are not reported, and they do not enter + % the gene targets either: the slope of a gene whose rule is the only gene of + % a reported reaction is the slope of that reaction. + geckoPath = findGECKOroot; + adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM', 'TestGEMAdapter.m')); + model = getGeckoTestModel(); + model.geneShortNames = model.genes; + ecModel = makeEcModel(model, false, adapter); + ecModel = getECfromGEM(ecModel); + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R2_EXP_1')) = 50; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R2_EXP_2')) = 10; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R2_REV_EXP_1')) = 50; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R2_REV_EXP_2')) = 10; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R3')) = 20; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns,'R5')) = 100; + ecModel.ec.source(:) = {'manual'}; + ecModel = applyKcatConstraints(ecModel); + ecModel = setProtPoolSize(ecModel, [], [], [], adapter); + adapter.params.bioRxn = 'R5'; + + fseof = ecFSEOF(ecModel, 'R3', 'S1', 'nSteps', 6, 'modelAdapter', adapter); + verifyNotEmpty(testCase, fseof.geneTargets) + for i = 1:size(fseof.geneTargets,1) + rxnRow = strcmp(fseof.rxnTargets(:,4), fseof.geneTargets{i,1}); + verifyEqual(testCase, sum(rxnRow), 1) + verifyEqual(testCase, fseof.geneTargets{i,3}, fseof.rxnTargets{rxnRow,3}, 'AbsTol', 1e-9) + end +end diff --git a/test/unit_tests/kcatAssignmentTests.m b/test/unit_tests/kcatAssignmentTests.m index 4f612967..68312c88 100644 --- a/test/unit_tests/kcatAssignmentTests.m +++ b/test/unit_tests/kcatAssignmentTests.m @@ -14,7 +14,8 @@ function testAssignKcatsReproducesChain_tc0079(testCase) adapter.params.kcat_assignment = struct('preset', 'gecko_classic'); fuzzy = fuzzyKcatMatching(base, [], adapter); dlkcat = readDLKcatOutput(base, 'modelAdapter', adapter); - merged = mergeKcats({fuzzy, dlkcat}, {'database_top','dlkcat','database_bottom'}); + % gecko_classic ranks the specific-activity matches (origin 5 and 6) above DLKcat + merged = mergeKcats({fuzzy, dlkcat}, {'database_top','dlkcat','database_bottom'}, 'topOriginLimit', 6); A = assignKcatValues(base, merged); A = applyCustomKcats(A, 'modelAdapter', adapter); A = getKcatAcrossIsozymes(A); diff --git a/test/unit_tests/kcatListTests.m b/test/unit_tests/kcatListTests.m new file mode 100644 index 00000000..db343938 --- /dev/null +++ b/test/unit_tests/kcatListTests.m @@ -0,0 +1,57 @@ +%run this test case with the command +%results = runtests('kcatListTests.m') + +function tests = kcatListTests + tests = functiontests(localfunctions); +end + +function testAssignKcatValuesSourceAfterZeroRows_tc0104(testCase) + % Rows with a zero kcat are dropped together with their source. + adapter = getTestAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + kcatList = struct(); + kcatList.source = 'test'; + kcatList.rxns = ecModel.ec.rxns(1:2); + kcatList.kcats = [0; 5]; + kcatList.kcatSource = {'a'; 'b'}; + ecModel = assignKcatValues(ecModel, kcatList); + verifyEqual(testCase, ecModel.ec.kcat(2), 5) + verifyEqual(testCase, ecModel.ec.source{2}, 'b') + verifyEqual(testCase, ecModel.ec.kcat(1), 0) +end + +function testLoadBRENDAdataLowercaseEc_tc0105(testCase) + % An EC code with a lowercase letter (provisional EC numbers end in n) + % is found in the molecular weight table. + adapter = getTestAdapter(); + scratchDir = tempname; + cleanupDir = onCleanup(@() rmdir(scratchDir, 's')); + mkdir(fullfile(scratchDir,'data')); + writeTable(fullfile(scratchDir,'data','kcat.tsv'), ... + 'ec_code\tsubstrate\torganism\tkcat_max\tkcat_median\tn\treferences\n', ... + '3.5.1.n3\tm1\ttestorg\t99\t5\t3\tPMID:1\n'); + writeTable(fullfile(scratchDir,'data','mw.tsv'), ... + 'ec_code\tsubstrate\torganism\tmw\tn\treferences\n', ... + '3.5.1.n3\t*\ttestorg\t60000\t1\t*\n'); + writeTable(fullfile(scratchDir,'data','sa.tsv'), ... + 'ec_code\tsubstrate\torganism\tsa_max\tsa_median\tn\treferences\n', ... + '3.5.1.n3\t*\ttestorg\t120\t7\t2\tPMID:2\n'); + brendaAdapter = adapter; + brendaAdapter.params.path = scratchDir; + [~, SAcell] = loadBRENDAdata('modelAdapter', brendaAdapter); + verifyEqual(testCase, SAcell{1}, {'3.5.1.n3'}) + verifyEqual(testCase, SAcell{3}, (120/60) * (60000/1000), 'AbsTol', 1e-9) +end + +function adapter = getTestAdapter() + geckoPath = findGECKOroot; + adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM','TestGEMAdapter.m')); +end + +function writeTable(file, header, row) + fid = fopen(file, 'w'); + fprintf(fid, '# BRENDA release test-fixture\n'); + fprintf(fid, header); + fprintf(fid, row); + fclose(fid); +end diff --git a/test/unit_tests/kcatMatchTests.m b/test/unit_tests/kcatMatchTests.m new file mode 100644 index 00000000..3ec73a03 --- /dev/null +++ b/test/unit_tests/kcatMatchTests.m @@ -0,0 +1,260 @@ +%run this test case with the command +%results = runtests('kcatMatchTests.m') + +function tests = kcatMatchTests + tests = functiontests(localfunctions); +end + +function testLoadBRENDAdataDropsKcatFarBelowItsEc_tc0220(testCase) + % A kcat more than four orders of magnitude below the median of its EC + % is dropped by default and kept when the filter is switched off. + adapter = brendaAdapter(testCase, ... + [kcatRows('1.1.1.1', {'o1','o2','o3','o4','o5','o6'}, 10, 10), ... + kcatRows('1.1.1.1', {'o7'}, 1e-4, 1e-4)]); + KCATcell = loadBRENDAdata('modelAdapter', adapter); + verifyEqual(testCase, KCATcell{4}, repmat(10,6,1)) + KCATcell = loadBRENDAdata('modelAdapter', adapter, 'filterOutlierKcats', false); + verifyEqual(testCase, sort(KCATcell{4}), sort([repmat(10,6,1); 1e-4])) +end + +function testLoadBRENDAdataReferenceUsesMaxAndMedian_tc0221(testCase) + % The reference of an EC is built from the max and the median values, + % whichever aggregation is read: a median that is far below its EC is an + % outlier even when the max of the same row is not. + adapter = brendaAdapter(testCase, ... + [kcatRows('2.1.1.1', {'o1','o2','o3','o4'}, 100, 100), ... + kcatRows('2.1.1.1', {'o5'}, 100, 0.001)]); + KCATcell = loadBRENDAdata('modelAdapter', adapter, 'kcatAggregation', 'median'); + verifyEqual(testCase, KCATcell{4}, repmat(100,4,1)) + KCATcell = loadBRENDAdata('modelAdapter', adapter, 'kcatAggregation', 'max'); + verifyEqual(testCase, KCATcell{4}, repmat(100,5,1)) +end + +function testLoadBRENDAdataReferenceNeedsFiveValues_tc0222(testCase) + % The values of the max and the median of a row both count towards the + % reference of its EC. With fewer than five of them at every EC level + % there is no reference and no row is dropped; with five or more at a + % shorter EC prefix, that prefix is the reference. + adapter = brendaAdapter(testCase, ... + [kcatRows('3.1.1.1', {'o1'}, 10, 10), kcatRows('3.1.1.1', {'o2'}, 1e-9, 1e-9)]); + KCATcell = loadBRENDAdata('modelAdapter', adapter); + verifyEqual(testCase, numel(KCATcell{4}), 2) + adapter = brendaAdapter(testCase, ... + [kcatRows('4.1.1.2', {'o1','o2','o3'}, 10, 10), kcatRows('4.1.1.1', {'o4'}, 1e-9, 1e-9)]); + KCATcell = loadBRENDAdata('modelAdapter', adapter); + verifyEqual(testCase, KCATcell{4}, repmat(10,3,1)) +end + +function testFuzzyKcatMatchingFullWildcard_tc0223(testCase) + % An EC number without any match at the wildcard levels 1 to 3 is matched + % against every EC number in BRENDA, at wildcard level 4. + [ecModel, fuzzyAdapter] = fuzzyModel(testCase, {'5.5.5.5\tm1\ttestus testus\t42\t42\t1\t*\n'}); + ecModel.ec.eccodes{strcmp(ecModel.ec.rxns, 'R2_EXP_1')} = '9.9.9.9'; + kcatList = fuzzyKcatMatching(ecModel, ismember(ecModel.ec.rxns,{'R2_EXP_1'}), fuzzyAdapter); + verifyEqual(testCase, kcatList.kcats, 42) + verifyEqual(testCase, kcatList.wildcardLvl, 4) +end + +function testFuzzyKcatMatchingTieReportsTheSelectedEc_tc0224(testCase) + % Two EC numbers of a reaction give the same kcat, one through a wildcard + % and one exactly: the wildcard level and origin are those of the exact one. + [ecModel, fuzzyAdapter] = fuzzyModel(testCase, ... + {'1.1.1.1\tm1\ttestus testus\t7\t7\t1\t*\n', ... + '1.1.1.9\tm1\ttestus testus\t7\t7\t1\t*\n'}); + ecModel.ec.eccodes{strcmp(ecModel.ec.rxns, 'R2_EXP_1')} = '1.1.1.5;1.1.1.1'; + kcatList = fuzzyKcatMatching(ecModel, ismember(ecModel.ec.rxns,{'R2_EXP_1'}), fuzzyAdapter); + verifyEqual(testCase, kcatList.kcats, 7) + verifyEqual(testCase, kcatList.wildcardLvl, 0) + verifyEqual(testCase, kcatList.origin, 1) +end + +function testFuzzyKcatMatchingSubsetKeepsOrderAndRepeats_tc0225(testCase) + % Indices select the reactions in the order given, a repeated index + % gives a repeated reaction, and a logical vector gives model order. + [ecModel, fuzzyAdapter] = fuzzyModel(testCase, {'1.1.1.1\tm1\ttestus testus\t7\t7\t1\t*\n'}); + ids = {'R3','R2_EXP_1','R2_EXP_1'}; + idx = cellfun(@(r) find(strcmp(ecModel.ec.rxns, r)), ids); + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', idx, 'modelAdapter', fuzzyAdapter); + verifyEqual(testCase, kcatList.rxns(:)', ids) + verifyEqual(testCase, numel(kcatList.kcats), 3) + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', ismember(ecModel.ec.rxns, ids), ... + 'modelAdapter', fuzzyAdapter); + verifyEqual(testCase, kcatList.rxns, ecModel.ec.rxns(ismember(ecModel.ec.rxns, ids))) +end + +function testFuzzyKcatMatchingEmptySubsetSelectsNone_tc0226(testCase) + % [] selects all reactions; any other empty selection selects none. + [ecModel, fuzzyAdapter] = fuzzyModel(testCase, {'1.1.1.1\tm1\ttestus testus\t7\t7\t1\t*\n'}); + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', [], 'modelAdapter', fuzzyAdapter); + verifyEqual(testCase, numel(kcatList.rxns), numel(ecModel.ec.rxns)) + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', zeros(0,1), 'modelAdapter', fuzzyAdapter); + verifyEmpty(testCase, kcatList.rxns) + verifyEmpty(testCase, kcatList.kcats) + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', false(numel(ecModel.ec.rxns),1), ... + 'modelAdapter', fuzzyAdapter); + verifyEmpty(testCase, kcatList.rxns) + adapter = getTestAdapter(); + verifyError(testCase, @() kcatsFromBrenda(ecModel, 'ecRxns', zeros(0,1), 'save', false, ... + 'modelAdapter', adapter), 'kcatsFromBrenda:emptySubset') +end + +function testFuzzyKcatMatchingForceWildcardEscalatesFromTheRight_tc0231(testCase) + % forceWClvl replaces the rightmost level that is not a wildcard yet, for + % every EC number of a reaction, also a provisional one: level 4 gives + % -.-.-.- and level 1 gives a wildcard at the last level. + [ecModel, fuzzyAdapter] = fuzzyModel(testCase, {'1.1.1.1\tm1\ttestus testus\t7\t7\t1\t*\n'}); + ecModel.ec.eccodes{strcmp(ecModel.ec.rxns, 'R2_EXP_1')} = '1.1.1.1;3.5.1.n3'; + sel = ismember(ecModel.ec.rxns,{'R2_EXP_1'}); + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', sel, 'modelAdapter', fuzzyAdapter, 'forceWClvl', 1); + verifyEqual(testCase, kcatList.eccodes, {'1.1.1.-;3.5.1.-'}) + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', sel, 'modelAdapter', fuzzyAdapter, 'forceWClvl', 4); + verifyEqual(testCase, kcatList.eccodes, {'-.-.-.-;-.-.-.-'}) + + % A code with a trailing wildcard escalates from the last level that is not + % a wildcard yet, and a level above the number of levels left stops at -.-.-.- + ecModel.ec.eccodes{strcmp(ecModel.ec.rxns, 'R2_EXP_1')} = '1.1.-.-;3.5.1.n3'; + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', sel, 'modelAdapter', fuzzyAdapter, 'forceWClvl', 1); + verifyEqual(testCase, kcatList.eccodes, {'1.-.-.-;3.5.1.-'}) + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', sel, 'modelAdapter', fuzzyAdapter, 'forceWClvl', 6); + verifyEqual(testCase, kcatList.eccodes, {'-.-.-.-;-.-.-.-'}) +end + +function testFuzzyKcatMatchingEcCodeWithoutFourLevels_tc0239(testCase) + % An EC code without four levels is not escalated level by level: it is + % left as it is by forceWClvl and, when BRENDA has no match, ends with the + % no-match values instead of an error. + [ecModel, fuzzyAdapter] = fuzzyModel(testCase, {'1.1.1.1\tm1\ttestus testus\t7\t7\t1\t*\n'}); + ecModel.ec.eccodes{strcmp(ecModel.ec.rxns, 'R2_EXP_1')} = '9.9.9'; + sel = ismember(ecModel.ec.rxns,{'R2_EXP_1'}); + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', sel, 'modelAdapter', fuzzyAdapter); + verifyEqual(testCase, kcatList.eccodes, {'9.9.9'}) + verifyEqual(testCase, kcatList.kcats, 0) + kcatList = fuzzyKcatMatching(ecModel, 'ecRxns', sel, 'modelAdapter', fuzzyAdapter, 'forceWClvl', 1); + verifyEqual(testCase, kcatList.eccodes, {'9.9.9'}) + verifyEqual(testCase, kcatList.kcats, 0) +end + +function testApplyCustomKcatsSkipsABlankKcat_tc0227(testCase) + % An entry without a kcat leaves the reaction as it is. + [ecModel, adapter] = customModel(); + isR3 = strcmp(ecModel.ec.rxns,'R3'); + isR2 = strcmp(ecModel.ec.rxns,'R2_EXP_1'); + custom = customStruct({''; ''}, [NaN; 50], {'R3'; 'R2'}, {'';''}); + ecModel = applyCustomKcats(ecModel, custom, adapter); + verifyEqual(testCase, ecModel.ec.kcat(isR3), 10) + verifyFalse(testCase, strcmp(ecModel.ec.source{isR3}, 'custom')) + verifyEqual(testCase, ecModel.ec.kcat(isR2), 50) + verifyEqual(testCase, ecModel.ec.source{isR2}, 'custom') +end + +function testApplyCustomKcatsEmptyNoteIsNotAppended_tc0228(testCase) + % An empty note does not add a separator to the note of the reaction. + [ecModel, adapter] = customModel(); + custom = customStruct({''; ''}, [50; 60], {'R3'; 'R3'}, {'first'; ''}); + ecModel = applyCustomKcats(ecModel, custom, adapter); + verifyEqual(testCase, ecModel.ec.notes{strcmp(ecModel.ec.rxns,'R3')}, 'first') + custom = customStruct({''; ''}, [50; 60], {'R3'; 'R3'}, {''; 'second'}); + ecModel = applyCustomKcats(customModelAgain(), custom, adapter); + verifyEqual(testCase, ecModel.ec.notes{strcmp(ecModel.ec.rxns,'R3')}, 'second') +end + +function testApplyCustomKcatsMatchesIdsExactly_tc0229(testCase) + % Protein and reaction ids are case sensitive. + [ecModel, adapter] = customModel(); + before = ecModel.ec.kcat; + custom = customStruct({'p3'}, 123, {''}, {''}); + updated = applyCustomKcats(ecModel, custom, adapter); + verifyEqual(testCase, updated.ec.kcat, before) + custom = customStruct({'P3'}, 123, {'r3'}, {''}); + updated = applyCustomKcats(ecModel, custom, adapter); + verifyEqual(testCase, updated.ec.kcat, before) + custom = customStruct({'P3'}, 123, {''}, {''}); + updated = applyCustomKcats(ecModel, custom, adapter); + verifyTrue(testCase, any(updated.ec.kcat == 123)) +end + +function testApplyCustomKcatsLightModel_tc0230(testCase) + % A light ecModel is updated by the proteins of its rows, and all three + % outputs are returned. + adapter = getTestAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), true, adapter); + ecModel = getECfromGEM(ecModel); + custom = customStruct({'P3'}, 123, {''}, {'light'}); + [updated, rxnUpdated, notMatch] = applyCustomKcats(ecModel, custom, adapter); + verifyNotEmpty(testCase, rxnUpdated) + verifyTrue(testCase, all(updated.ec.kcat(ismember(updated.ec.rxns, rxnUpdated)) == 123)) + verifyTrue(testCase, istable(notMatch)) +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function adapter = getTestAdapter() + geckoPath = findGECKOroot; + adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM','TestGEMAdapter.m')); +end + +function rows = kcatRows(ec, organisms, kcatMax, kcatMedian) + rows = cell(1, numel(organisms)); + for i = 1:numel(organisms) + rows{i} = sprintf('%s\tm1\t%s\t%g\t%g\t1\t*\n', ec, organisms{i}, kcatMax, kcatMedian); + end +end + +function writeBrenda(brendaDir, kcatLines) + fid = fopen(fullfile(brendaDir, 'kcat.tsv'), 'w'); + fprintf(fid, '# BRENDA release test-fixture - kcat in 1/s\n'); + fprintf(fid, 'ec_code\tsubstrate\torganism\tkcat_max\tkcat_median\tn\treferences\n'); + for i = 1:numel(kcatLines) + fprintf(fid, kcatLines{i}); + end + fclose(fid); + fid = fopen(fullfile(brendaDir, 'mw.tsv'), 'w'); + fprintf(fid, '# BRENDA release test-fixture - molecular weight in g/mol\n'); + fprintf(fid, 'ec_code\tsubstrate\torganism\tmw\tn\treferences\n'); + fprintf(fid, '0.0.0.0\t*\tplaceholder\t1\t1\t*\n'); + fclose(fid); + fid = fopen(fullfile(brendaDir, 'sa.tsv'), 'w'); + fprintf(fid, '# BRENDA release test-fixture - specific activity in umol/min/mg\n'); + fprintf(fid, 'ec_code\tsubstrate\torganism\tsa_max\tsa_median\tn\treferences\n'); + fprintf(fid, '0.0.0.0\t*\tplaceholder\t1\t1\t1\t*\n'); + fclose(fid); +end + +function adapter = brendaAdapter(testCase, kcatLines) + scratchDir = tempname; + brendaDir = fullfile(scratchDir, 'data'); + mkdir(brendaDir); + testCase.addTeardown(@() rmdir(scratchDir, 's')); + writeBrenda(brendaDir, kcatLines); + adapter = getTestAdapter(); + adapter.params.path = scratchDir; +end + +function [ecModel, fuzzyAdapter] = fuzzyModel(testCase, kcatLines) + fuzzyAdapter = brendaAdapter(testCase, kcatLines); + geckoPath = findGECKOroot; + copyfile(fullfile(geckoPath,'test','unit_tests','ecTestGEM','data','PhylDist.mat'), ... + fullfile(fuzzyAdapter.params.path, 'data', 'PhylDist.mat')); + adapter = getTestAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + ecModel = getECfromGEM(ecModel); +end + +function [ecModel, adapter] = customModel() + adapter = getTestAdapter(); + ecModel = customModelAgain(); +end + +function ecModel = customModelAgain() + adapter = getTestAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + ecModel = getECfromGEM(ecModel); + ecModel.ec.kcat(:) = 10; +end + +function custom = customStruct(proteins, kcats, rxns, notes) + custom.proteins = proteins; + custom.kcat = kcats; + custom.rxns = rxns; + custom.notes = notes; +end diff --git a/test/unit_tests/parityDecisionTests.m b/test/unit_tests/parityDecisionTests.m new file mode 100644 index 00000000..0c11fee0 --- /dev/null +++ b/test/unit_tests/parityDecisionTests.m @@ -0,0 +1,106 @@ +%run this test case with the command +%results = runtests('parityDecisionTests.m') + +function tests = parityDecisionTests + tests = functiontests(localfunctions); +end + +function testMergeKcatsDefaultTopOriginLimitIs4_tc0240(testCase) + % A BRENDA match of origin 5 (specific activity, correct organism) ranks + % after DLKcat by default and above it with topOriginLimit 6. + [fuzzy, dlkcat] = originFiveLists(); + tiers = {'database_top', 'dlkcat', 'database_bottom'}; + merged = mergeKcats({fuzzy, dlkcat}, tiers); + verifyEqual(testCase, merged.kcats(strcmp(merged.rxns, 'R1')), 50) + merged = mergeKcats({fuzzy, dlkcat}, tiers, 'topOriginLimit', 4); + verifyEqual(testCase, merged.kcats(strcmp(merged.rxns, 'R1')), 50) + merged = mergeKcats({fuzzy, dlkcat}, tiers, 'topOriginLimit', 6); + verifyEqual(testCase, merged.kcats(strcmp(merged.rxns, 'R1')), 5) +end + +function testMergeDLKcatAndFuzzyKcatsDefaultTopOriginLimitIs4_tc0241(testCase) + [fuzzy, dlkcat] = originFiveLists(); + merged = mergeDLKcatAndFuzzyKcats(dlkcat, fuzzy); + verifyEqual(testCase, merged.kcats(strcmp(merged.rxns, 'R1')), 50) + merged = mergeDLKcatAndFuzzyKcats(dlkcat, fuzzy, 'topOriginLimit', 6); + verifyEqual(testCase, merged.kcats(strcmp(merged.rxns, 'R1')), 5) +end + +function testFindECInDBSingleGeneDropsWildcardParent_tc0242(testCase) + % One gene whose protein carries '1.1.1.4 1.1.1.- 1.1.1.303' keeps the + % two specific codes and loses the wildcard of the same prefix. + geneHashMap = containers.Map({'G1'}, {1}); + geneIndex = {1}; + DBecNum = {'1.1.1.4 1.1.1.- 1.1.1.303'}; + DBMW = 10000; + EC = findECInDB({'G1'}, DBecNum, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '1.1.1.4;1.1.1.303') + % A wildcard without a specific code of its prefix stays + DBecNum = {'1.1.1.- 2.2.2.2'}; + EC = findECInDB({'G1'}, DBecNum, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '1.1.1.-;2.2.2.2') +end + +function testFindECInDBSingleGeneShortCodesAndTrailingSpace_tc0243(testCase) + % A code with fewer than four levels is neither the wildcard parent nor + % the child of a four-level code, and a trailing space leaves no empty code. + geneHashMap = containers.Map({'G1'}, {1}); + geneIndex = {1}; + DBMW = 10000; + EC = findECInDB({'G1'}, {'1.1.1 1.1.1.4'}, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '1.1.1;1.1.1.4') + EC = findECInDB({'G1'}, {'1.1.1.4 '}, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '1.1.1.4') +end + +function testFindECInDBComplexDropsWildcardParent_tc0244(testCase) + % The shared EC set of an enzyme complex loses a wildcard code next to a + % specific code of its prefix, as the set of a single gene does. + DBMW = [10000; 20000; 30000]; + geneHashMap = containers.Map({'G1','G2','G3'}, {1,2,3}); + geneIndex = {1; 2; 3}; + EC = findECInDB({'G1','G2'}, {'1.1.1.4 1.1.1.-'; '1.1.1.4 1.1.1.-'; ''}, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '1.1.1.4') + % A wildcard that only the union carries is dropped too + EC = findECInDB({'G1','G2','G3'}, {'1.1.1.4'; '1.1.1.- 2.2.2.2'; '5.5.5.5'}, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '1.1.1.4;2.2.2.2;5.5.5.5') + % A wildcard shared by every gene stays, as nothing more specific exists + EC = findECInDB({'G1','G2'}, {'1.1.1.-'; '1.1.1.-'; ''}, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '1.1.1.-') +end + +function testFindECInDBGeneWithoutECIsLeftOut_tc0245(testCase) + % A gene without an EC number does not change the codes of the others: + % the complex is treated as the genes that have codes, in either order. + DBMW = [10000; 20000]; + geneHashMap = containers.Map({'G1','G2'}, {1,2}); + geneIndex = {1; 2}; + EC = findECInDB({'G1','G2'}, {'1.1.1.4 1.1.1.-'; ''}, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '1.1.1.4') + EC = findECInDB({'G1','G2'}, {''; '1.1.1.4 1.1.1.-'}, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '1.1.1.4') + % No gene has a code + EC = findECInDB({'G1','G2'}, {''; ''}, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '') + % A gene missing from the database is left out as well + geneHashMap = containers.Map({'G2'}, {2}); + geneIndex = {[]; 2}; + EC = findECInDB({'G1','G2'}, {''; '1.1.1.4 1.1.1.-'}, DBMW, geneIndex, geneHashMap); + verifyEqual(testCase, EC, '1.1.1.4') +end + +function [fuzzy, dlkcat] = originFiveLists() + fuzzy.rxns = {'R1'}; + fuzzy.kcats = 5; + fuzzy.genes = {{}}; + fuzzy.substrates = {{}}; + fuzzy.eccodes = {'1.1.1.1'}; + fuzzy.wildcardLvl = 0; + fuzzy.origin = 5; + fuzzy.source = 'brenda'; + dlkcat.rxns = {'R1'}; + dlkcat.kcats = 50; + dlkcat.genes = {'G1'}; + dlkcat.substrates = {'m1'}; + dlkcat.source = 'DLKcat'; +end diff --git a/test/unit_tests/parityRound3Tests.m b/test/unit_tests/parityRound3Tests.m new file mode 100644 index 00000000..d230faf9 --- /dev/null +++ b/test/unit_tests/parityRound3Tests.m @@ -0,0 +1,194 @@ +%run this test case with the command +%results = runtests('parityRound3Tests.m') + +function tests = parityRound3Tests + tests = functiontests(localfunctions); +end + +function testApplyKcatConstraintsAllZeroClearsFullModel_tc0250(testCase) + % With every kcat 0 the enzyme constraints of a full ecModel are removed + % and no new ones are written. A light ecModel is returned unchanged. + adapter = getAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + ecModel.ec.kcat(:) = 10; + ecModel = applyKcatConstraints(ecModel); + protMets = startsWith(ecModel.mets, 'prot_') & ~strcmp(ecModel.mets, 'prot_pool'); + metabolic = ~startsWith(ecModel.rxns, 'usage_prot_') & ~strcmp(ecModel.rxns, 'prot_pool_exchange'); + verifyGreaterThan(testCase, nnz(ecModel.S(protMets, metabolic)), 0) + ecModel.ec.kcat(:) = 0; + cleared = applyKcatConstraints(ecModel); + verifyEqual(testCase, nnz(cleared.S(protMets, metabolic)), 0) + % the usage reactions and the pool exchange keep their fixed entries + verifyEqual(testCase, cleared.S(:, ~metabolic), ecModel.S(:, ~metabolic)) +end + +function testApplyKcatConstraintsAllZeroLightIsUnchanged_tc0250light(testCase) + % A light ecModel with every kcat 0 is returned unchanged (control: the + % behaviour is the same before and after the full-model change). + adapter = getAdapter(); + lecModel = makeEcModel(getGeckoTestModel(), true, adapter); + lecModel.ec.kcat(:) = 10; + lecModel = applyKcatConstraints(lecModel); + poolRow = strcmp(lecModel.mets, 'prot_pool'); + verifyGreaterThan(testCase, nnz(lecModel.S(poolRow, :)), 1) + lecModel.ec.kcat(:) = 0; + same = applyKcatConstraints(lecModel); + verifyEqual(testCase, same.S, lecModel.S) +end + +function testSensitivityTuningStopsOnGrowthBelowTolerance_tc0251(testCase) + % The solver reports a growth of 1 in the first iteration and 1 + 5e-10 in + % the second: the gain is below 1e-9, so the loop stops after the second + % iteration (3 solves with the one before the loop). An exact-equality + % loop would go on until two growths are identical (4 solves). + [ecModel, tuneAdapter] = tuningModel(); + cleanup = shadowSolveLP('growth', [1 1 1+5e-10]); %#ok + tunedModel = sensitivityTuning(ecModel, 'desiredGrowthRate', 50, ... + 'modelAdapter', tuneAdapter, 'protToIgnore', {'P4'}, 'verbose', false); + verifyEqual(testCase, getappdata(0, 'gecko_test_solvelp_count'), 3) + r5 = tunedModel.ec.kcat(strcmp(tunedModel.ec.rxns, 'R5')); + verifyGreaterThan(testCase, r5, 200) +end + +function testSensitivityTuningEndsWhenTheLPBecomesInfeasible_tc0254(testCase) + % The third solve (the second iteration) has no solution: RAVEN returns + % an empty objective value. The loop ends with a warning and returns + % what was tuned in the first iteration. + [ecModel, tuneAdapter] = tuningModel(); + cleanup = shadowSolveLP('infeasibleFrom', 3); %#ok + out = evalc(['tunedModel = sensitivityTuning(ecModel, ''desiredGrowthRate'', 50, ' ... + '''modelAdapter'', tuneAdapter, ''protToIgnore'', {''P4''}, ''verbose'', false);']); + verifyTrue(testCase, contains(out, 'no feasible solution')) + verifyEqual(testCase, getappdata(0, 'gecko_test_solvelp_count'), 3) + r5 = tunedModel.ec.kcat(strcmp(tunedModel.ec.rxns, 'R5')); + verifyGreaterThan(testCase, r5, 200) +end + +function testSensitivityTuningStopsWhenNothingCarriesFlux_tc0253(testCase) + % R5, the growth reaction, is blocked: growth is 0 from the start and no + % usage reaction carries flux. The loop ends with a warning instead of + % tuning an arbitrary reaction or failing, and no kcat changes. + adapter = getAdapter(); + model = loadConventionalGEM('modelAdapter', adapter); + for id = {'R2', 'R4', 'R5'} + idx = strcmp(model.rxns, id{1}); + model.lb(idx) = 0; model.ub(idx) = 0; + end + ecModel = makeEcModel(model, false, adapter); + ecModel = getECfromGEM(ecModel); + ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R3')) = 0.05; + ecModel = applyKcatConstraints(ecModel); + ecModel = setProtPoolSize(ecModel, [], [], [], adapter); + tuneAdapter = adapter; + tuneAdapter.params.bioRxn = 'R5'; + [tunedModel, tunedKcats] = sensitivityTuning(ecModel, 'desiredGrowthRate', 50, ... + 'modelAdapter', tuneAdapter, 'verbose', false); + verifyEqual(testCase, tunedModel.ec.kcat, ecModel.ec.kcat) + verifyEmpty(testCase, tunedKcats.rxns) +end + +function testGetKcatAcrossIsozymesUsesPositiveSiblingsOnly_tc0252(testCase) + % Only isozymes with a kcat above 0 enter the mean: a negative or NaN + % sibling does not fill a kcat of 0. + adapter = getAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + ecModel.ec.kcat(:) = 0; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R3')) = 8; + % R2: a zero next to a negative sibling is not filled + ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R2_EXP_2')) = -3; + % R2_REV: a zero next to a NaN sibling is not filled + ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R2_REV_EXP_2')) = NaN; + out = getKcatAcrossIsozymes(ecModel); + verifyEqual(testCase, out.ec.kcat(strcmp(out.ec.rxns, 'R2_EXP_1')), 0) + verifyEqual(testCase, out.ec.kcat(strcmp(out.ec.rxns, 'R2_REV_EXP_1')), 0) + % a positive sibling fills the zero (control: the same before and after) + ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R2_EXP_2')) = 6; + out = getKcatAcrossIsozymes(ecModel); + j = strcmp(out.ec.rxns, 'R2_EXP_1'); + verifyEqual(testCase, out.ec.kcat(j), 6) + verifyEqual(testCase, out.ec.source{j}, 'isozymes') + % a NaN sibling next to a positive one is left out of the mean: three + % isozymes of R2 with kcats 0, NaN and 6 give 6 to the first + ecModel.ec.kcat(:) = 0; + ecModel.ec.rxns{strcmp(ecModel.ec.rxns, 'R3')} = 'R2_EXP_3'; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R2_EXP_2')) = NaN; + ecModel.ec.kcat(strcmp(ecModel.ec.rxns, 'R2_EXP_3')) = 6; + out = getKcatAcrossIsozymes(ecModel); + j = strcmp(out.ec.rxns, 'R2_EXP_1'); + verifyEqual(testCase, out.ec.kcat(j), 6) + verifyEqual(testCase, out.ec.source{j}, 'isozymes') + verifyTrue(testCase, isnan(out.ec.kcat(strcmp(out.ec.rxns, 'R2_EXP_2')))) +end + +function [ecModel, tuneAdapter] = tuningModel() +% R3 limits growth, its enzyme P4 is ignored by the tests, so tuning raises R5. +adapter = getAdapter(); +model = loadConventionalGEM('modelAdapter', adapter); +for id = {'R2', 'R4'} + idx = strcmp(model.rxns, id{1}); + model.lb(idx) = 0; model.ub(idx) = 0; +end +ecModel = makeEcModel(model, false, adapter); +ecModel = getECfromGEM(ecModel); +kcats = struct('R3', 0.05, 'R5', 200); +for id = fieldnames(kcats)' + j = strcmp(ecModel.ec.rxns, id{1}); + ecModel.ec.kcat(j) = kcats.(id{1}); + ecModel.ec.source{j} = 'manual'; +end +ecModel = applyKcatConstraints(ecModel); +ecModel = setProtPoolSize(ecModel, [], [], [], adapter); +tuneAdapter = adapter; +tuneAdapter.params.bioRxn = 'R5'; +end + +function cleanup = shadowSolveLP(mode, values) +% Puts a solveLP on the path that calls RAVEN's and then reports the growth +% of solve number n as values(min(n, end)) ('growth'), or no solution from +% solve number values on ('infeasibleFrom'). The number of solves is kept in +% the root appdata gecko_test_solvelp_count. The path is restored when the +% returned object is cleared. +setappdata(0, 'gecko_test_solvelp_real', @solveLP); +setappdata(0, 'gecko_test_solvelp_count', 0); +setappdata(0, 'gecko_test_solvelp_mode', mode); +setappdata(0, 'gecko_test_solvelp_values', values); +stubDir = tempname; +mkdir(stubDir); +fid = fopen(fullfile(stubDir, 'solveLP.m'), 'w'); +fprintf(fid, '%s\n', ... + 'function [res,hs] = solveLP(varargin)', ... + 'n = getappdata(0,''gecko_test_solvelp_count'') + 1;', ... + 'setappdata(0,''gecko_test_solvelp_count'',n);', ... + 'real = getappdata(0,''gecko_test_solvelp_real'');', ... + 'mode = getappdata(0,''gecko_test_solvelp_mode'');', ... + 'values = getappdata(0,''gecko_test_solvelp_values'');', ... + '[res,hs] = real(varargin{:});', ... + 'switch mode', ... + ' case ''growth''', ... + ' res.f = values(min(n,numel(values)));', ... + ' case ''infeasibleFrom''', ... + ' if n >= values', ... + ' res.f = [];', ... + ' res.x = [];', ... + ' end', ... + 'end', ... + 'end'); +fclose(fid); +addpath(stubDir, '-begin'); +cleanup = onCleanup(@() removeShadow(stubDir)); +end + +function removeShadow(stubDir) +rmpath(stubDir); +rmdir(stubDir, 's'); +for name = {'real', 'mode', 'values'} + if isappdata(0, ['gecko_test_solvelp_' name{1}]) + rmappdata(0, ['gecko_test_solvelp_' name{1}]); + end +end +end + +function adapter = getAdapter() +geckoPath = findGECKOroot; +adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM','TestGEMAdapter.m')); +end diff --git a/test/unit_tests/protDataTests.m b/test/unit_tests/protDataTests.m new file mode 100644 index 00000000..70d1a6bb --- /dev/null +++ b/test/unit_tests/protDataTests.m @@ -0,0 +1,63 @@ +%run this test case with the command +%results = runtests('protDataTests.m') + +function tests = protDataTests + tests = functiontests(localfunctions); +end + +function testLoadProtDataRsdFilterDropsWholeRow_tc0101(testCase) + % A protein that exceeds maxRSD is discarded in every replicate, not + % only in the first one. + adapter = getTestAdapter(); + file = writeProteomics({'X1','X2','X3'}, [1 1 10; 5 5 5; 6 6 6]); + protData = loadProtData(3, 'protDataFile', file, 'modelAdapter', adapter, ... + 'maxRSD', 1, 'addStdevs', 0, 'cutLowest', 0); + verifyEqual(testCase, sort(protData.uniprotIDs), {'X2';'X3'}) + verifyEqual(testCase, sort(protData.abundances), [5;6]) +end + +function testLoadProtDataCutLowestIsUsed_tc0102(testCase) + % The lowest cutLowest percent of the values of a condition is discarded. + adapter = getTestAdapter(); + ids = strcat('P', string(1:40)); + values = repmat((1:40)', 1, 3); + file = writeProteomics(cellstr(ids), values); + p0 = loadProtData(3, 'protDataFile', file, 'modelAdapter', adapter, 'addStdevs', 0, 'cutLowest', 0); + p5 = loadProtData(3, 'protDataFile', file, 'modelAdapter', adapter, 'addStdevs', 0); + p20 = loadProtData(3, 'protDataFile', file, 'modelAdapter', adapter, 'addStdevs', 0, 'cutLowest', 20); + verifyEqual(testCase, numel(p0.uniprotIDs), 40) + verifyEqual(testCase, numel(p5.uniprotIDs), 38) + verifyEqual(testCase, numel(p20.uniprotIDs), 32) + verifyEqual(testCase, min(p20.abundances), 9) +end + +function testApplyComplexDataStoichiometryByProteinId_tc0103(testCase) + % The complex lists its proteins in another order than model.ec.enzymes; + % each protein receives its own stoichiometry. + adapter = getTestAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + complexInfo = struct('complexID', 'R2Compl', 'name', 'R2Compl', ... + 'specie', 'Testus testus; 000000', 'geneName', {{'G2';'G1'}}, ... + 'protID', {{'P2';'P1'}}, 'stochiometry', [2;1]); + ecModel = applyComplexData(ecModel, complexInfo, adapter, false); + rxn = find(strcmp(ecModel.ec.rxns, 'R2_EXP_1')); + p1 = strcmp(ecModel.ec.enzymes, 'P1'); + p2 = strcmp(ecModel.ec.enzymes, 'P2'); + verifyEqual(testCase, full(ecModel.ec.rxnEnzMat(rxn, p1)), 1) + verifyEqual(testCase, full(ecModel.ec.rxnEnzMat(rxn, p2)), 2) +end + +function adapter = getTestAdapter() + geckoPath = findGECKOroot; + adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM','TestGEMAdapter.m')); +end + +function file = writeProteomics(ids, values) + file = [tempname '.tsv']; + fid = fopen(file, 'w'); + fprintf(fid, 'uniprot\tr1\tr2\tr3\n'); + for i = 1:numel(ids) + fprintf(fid, '%s\t%g\t%g\t%g\n', ids{i}, values(i,:)); + end + fclose(fid); +end