How to deal with version of gene id in modeladapter file, obj.params.kegg.geneID ?

[PkiwiBird]

I converted a cobra model to RAVEN model and then used that model with GECKO package.
The problem is that the cobra model entrez gene ids have version numbers (e.g. 9963.2).
The model adapter file should have obj.params.kegg.geneID = 'NCBI-GeneID'; but can the version of gene ids (the suffix after the dot) in the model create problems?
If so, how can we get properly rid of them?
I did this:

load 'OrganAtlas_Harvey.mat'
omc_ml=OrganCompendium_male; % organ model compendium of male
liver_ml=omc_ml.Liver.modelAllComp; % liver model of male
for i=1:length(liver_ml.genes); liver_ml.genes{i}=num2str(floor(str2num(liver_ml.genes{i}))); end
liver_ml = ravenCobraWrapper(liver_ml);

It removed the suffix in the .genes field and when converting to RAVEN it also applied the change to .grRules, but there are repeated genes in .genes (e.g. '10327'; '10327' in rows 31 and 32) and in the some rows of .grRules (e.g. '231 or 10327 or 10327' in row 4974).

When doing unique to genes after getting rid of suffixes it changed the gene ids in the .grRules, for e.g. in the same row 4974, the rule is now '259230 or 10351 or 10331':

for i=1:length(liver_ml.genes); liver_ml.genes{i}=num2str(floor(str2num(liver_ml.genes{i}))); end
liver_ml.genes=unique(liver_ml.genes);
liver_ml = ravenCobraWrapper(liver_ml);

[edkerk]

What evidence do you have that e.g. 10327.1 and 10327.2 are different genes? At NCBI there seems to be just one entry, without versioning? Identifiers.org also specifies that NCBI gene IDs contain digits only (no dots, as in this expression: ^\d+$).

If the versioning is not informative and 10327.1 and 10327.2 are the same genes, then the best approach would be to modfiy this in the model before starting in GECKO. Note that this is not trivial, as you should not only correct .genes, and .grRules, you also want to get rid of duplicate entries in .genes and you would want .rxnGeneMat to now refer to the new unique indexes in .genes.

[PkiwiBird]

I have adapted the gene rules as required but after that, when converting to RAVEN, the gene rules change completely: the gene number (which is neither the correct gene ID nor correct position in the gene list) and the gene rule format.

%% original model:
load 'OrganAtlas_Harvey.mat'
omc_ml=OrganCompendium_male; % organ model compendium of male
liver_ml=omc_ml.Liver.modelAllComp; % liver model of male
% rule for a reaction:
n=3468;
liver_ml.grRules{n}
% '( 5476.1 ) and ( 2720.1 ) and ( 2588.1 ) and ( 4758.1 ) and ( 5660.1 ) 
or ( 5476.1 ) and ( 2720.1 ) and ( 2588.1 ) and ( 4758.1 ) and ( 2760.1 ) 
or ( 5476.1 ) and ( 2720.1 ) and ( 5660.1 ) 
or ( 5476.1 ) and ( 2720.1 ) and ( 2760.1 )'

%% model after simplifying gene rules:
liver_ml = simplifyGeneGeneRules(liver_ml);
% rule for same reaction - with correct format:
n=3468;
liver_ml.grRules{n}
% '(2588 and 2720 and 2760 and 4758 and 5476) 
or (2588 and 2720 and 4758 and 5476 and 5660) 
or (2720 and 2760 and 5476) 
or (2720 and 5476 and 5660)'
% the positions in the rxnGeneMat matrix for that rule are correct:
liver_ml.rxnGeneMat(n, :) 
% (1,451)      1
% (1,512)      1
% (1,525)      1
% (1,789)      1
% (1,1007)     1
% (1,1064)     1
liver_ml.genes{451} % '2588'
liver_ml.genes{512} % '2720'
liver_ml.genes{525} % '2760'
liver_ml.genes{789} % '4758'
liver_ml.genes{1007} % '5476'
liver_ml.genes{1064} % '5660'

%% model after simplifying gene rules and converting to RAVEN:
liver_ml = ravenCobraWrapper(liver_ml);
% same rule has now incorrect gene numbers and format:
n=3468;
liver_ml.grRules{n}
% '6547) and (3157) and (2799) and (5211) and (7086) 
or (6547) and (3157) and (2799) and (5211) and (3291) 
or (6547) and (3157) and (7086) 
or (6547) and (3157) and (3291'

The format of the rule and the gene ids are different after converting to RAVEN.
Do you how to solve this?

The function simplifyGeneGeneRules:

function [model] = simplifyGeneGeneRules(model)
    %simplifyGeneRules simplifies genes and gene rules
    %   - removes gene versions of the model fields .genes and .grRules and
    %   removes duplicated genes in those fields
    %   - simplifies complicated gene rules in the formats:
    %   '(A and B) or C', 'A or B', 'A and C' where A, B and C are genes 
    %   - updates the rxnGeneMat accordingly
    %
    % Inputs:
    %   model - a model structure
    %     
    % Outputs:
    %   model - updated model structure
    %
    % Example:
    %   [model] = functionName(model);
    
    % remove the gene version from Entrez gene id:
    % in the '.genes':
    model.genes = unique(regexprep(model.genes,'\.\d+$',''));
    %model.genes = regexprep(model.genes,'\.\d+$','');
    % in the '.grRules':
    model.grRules= regexprep(model.grRules,'\.\d+\>','');
    % in '.grRules' get rid of spaces after '(' and before ')':
    model.grRules=regexprep(model.grRules, '( ', '(');
    model.grRules=regexprep(model.grRules, ' )', ')');
    
    % get rid of repeated genes in gene rules without 'and' (only with 'or' or just the gene):
    % e.g. n=11;
    for n = 1:length(model.grRules)
        rl = model.grRules{n};
        match = regexp(rl, 'and', 'match');
        % if no 'and' is present,
        % the following gets rid of ' or ' or of '(', ')' if they are present:
        if isempty(match)
            % note: if there are no occurences,
            % the following doesn't split, if they exist, it splits:
            rl = regexprep(rl,'\((.*)\)', '$1'); % the parenthesis in '(.*)' capture the group, '$1' calls that 1st group
            rl = unique(regexp(rl,'(\sor\s)','split'));                 
            gr = join(rl, ' or ');
            model.grRules(n) = gr;
        end
    end
    
    % get rid of parenthesis in single genes:
    % e.g. n= 631;
    for n = 1:length(model.grRules)
        rl = model.grRules{n};
        if any(rl) % if there is a gene rule for the reaction
            rl = regexprep(rl,'\((\d+)\)', '$1'); % the parenthesis in '(\d+)' capture the group, '$1' calls that 1st group
            model.grRules(n) = cellstr(rl);
        end
    end
    
    % simplify rules with '(( and ) or ( and )) and (( and ) or ( and ))':
    % e.g. n=195; n=3678;
    % note: not ready to deal with gene rules where 'substrings' after spliting by ')) and (('
    % don't have 'or', e.g rl = '((50617 and 533 and 527 and 8992 and 245972)) and ((51382 and 9296 and 529 and 245973 and 523 and 51606 and 9550 and 526) or (51382 and 9296 and 525 and 127124 and 51606 and 529 and 528 and 523))'
    for n = 1:length(model.grRules)
        rl = model.grRules{n};
        if any(strfind(rl, ')) and (('))
            els = regexp(rl, '\)\) and \(\(', 'split'); % split by ')) and (('
            a = cell(1, length(els));
            for i=1:length(els) % for each of resulting elements
                if contains(els{i}, ') or (')  % if element contains an 'or'
                    el = regexprep(els{i}, '\(\(|\)\)', ''); % remove '((', and also '))'
                    a{1, i} = string(regexp(el, '\) or \(', 'split')); % split by ') or ('
                end
            end
            [~, prd] = cartesianProd(a{:}); % cartesian product
            cmb = join(string(prd), ' and ');
            cmb = strcat('(', join(cmb, ') or ('), ')');
            model.grRules(n) = cellstr(cmb);
        end
    end
    
    % simplify rules with 'and ( or ))':
    % e.g. n=922
    % note: following might not work in situations where there is more than one
    % match in the same rule
    for n = 1:length(model.grRules)
        rl = model.grRules{n};
        match = regexp(rl, '\d+ and \(\d+\sor\s\d+', 'match');
        if ~isempty(match)               
            BeforeMatch = regexp(rl, '.*(?=(\<\d+ and \(\d+\sor\s\d+))', 'match');
            AfterMatch = regexp(rl, '(?<=(\d+ and \(\d+\sor\s\d+\>\))).*', 'match'); 
            els = regexp(match{:}, '\sand\s', 'split');
            a = cell(1, length(els));
            for i=1:length(els)
                el = regexprep(els{i}, '\(|\)', '');
                a{1, i} = string(regexp(el, '\sor\s', 'split'));
            end
            [~, prd] = cartesianProd(a{:}); % cartesian product
            cmb = join(string(prd), ' and ');
            cmb = strcat('(', join(cmb, ') or ('), ')');
            fn = strcat(BeforeMatch, cmb, AfterMatch);
            model.grRules(n) = cellstr(fn);
            disp(n)
        end
    end
    
    % simplify rules with '( or ) and ) ':
    % e.g. n=1250;n=1251;n=1410;n=4052
    % note: following might not work in situations where there is more than one
    % match in the same rule
    for n = 1:length(model.grRules)
        rl = model.grRules{n};
        match = regexp(rl, '\(\d+(\sor\s\d+)+\) and \d+', 'match');
        if ~isempty(match)
            BeforeMatch = regexp(rl,'.*(?=\(\d+(\sor\s\d+)+\) and \d+)', 'match');
            AfterMatch = regexp(rl, '(?<=\(\d+(\sor\s\d+)+\) and \d+\>).*', 'match');
            els = regexp(match{:}, '\sand\s', 'split');
            a = cell(1, length(els));
            for i=1:length(els)
                el = regexprep(els{i}, '\(|\)', '');
                a{1, i} = string(regexp(el, '\sor\s', 'split'));
            end
            [~, prd] = cartesianProd(a{:}); % cartesian product
            cmb = join(string(prd), ' and ');
            cmb = strcat('(', join(cmb, ') or ('), ')');
            fn = strcat(BeforeMatch, cmb, AfterMatch);
            model.grRules(n) = cellstr(fn);
            disp(n)
        end
    end
    
    % join 'and' outside with 'and' inside complexes:
    % e.g. n=456;n=863;n=864;n=865;n=904;n=4964;n=3678
    for n = 1:length(model.grRules)
        rl = model.grRules{n};
        match = regexp(rl, '(\)\sand\s\()|(\sand\s\()|(\)\sand\s)', 'match');
        if ~isempty(match)
            els = regexp(rl, '\sor\s', 'split');
            a = cell(1, length(els));
            for i=1:length(els)
                a{i} = regexprep(els{i}, '\((.*?)\)', '$1'); % this is lazy
                a{i} = regexprep(a{i}, '\((.*)\)', '$1'); % this is greedy, so the next step is sometimes also needed
                %els{i} = regexprep(els{i}, '(\)\sand\s()|(\sand\s\()|(\)\sand\s)', ' and ');
            end
            cmb = strcat('(', join(a, ') or ('), ')');
            cmb = regexprep(cmb,'\((\d+)\)', '$1'); % if parenthesis are added to single genes, remove them
            model.grRules(n) = cmb;
            disp(n)
        end
    end
    
    % simplify cases of ' or ((B and C) or (D and E))', '((A and B) or (C and
    % D)) or ', ' or ((B and C) or (D and E)) or ', ' or ((A and B) or (C and D) or (E and F))':
    % e.g. n=922; one part of n=1250;
    % rl = '67558 or ((340024 and 57393) or (340024 and 59272)) or 6541 or 55089'
    % rl = '((340024 and 57393) or (340024 and 59272)) or 6541 or 55089'
    % rl = '67558 or ((340024 and 57393) or (340024 and 59272))'
    for n = 1:length(model.grRules)
        rl = model.grRules{n};
        match = regexp(rl, '((?<=or\s)\((\(\d+\sand\s\d+\)\sor\s)+\(\d+\sand\s\d+\)\))|(\((\(\d+\sand\s\d+\)\sor\s)+\(\d+\sand\s\d+\)\)(?=\sor))', 'match');
        if ~isempty(match)
            r = regexprep(rl, '\(\((.*)\)\)', '\($1\)');
            model.grRules(n) = cellstr(r);
            disp(n)
        end
    end 
    
    % simplify cases of 'or (A or B or C) or',
    % '(A or B or C) or', 'or (A or B or C)':
    % e.g. n=1250; n=1251; n=1410; n=4052
    for n = 1:length(model.grRules)
        rl = model.grRules{n};
        match = regexp(rl, '((?<=or\s)\((\d+\sor\s)+\d+\))|(\((\d+\sor\s)+\d+\)(?=\sor))', 'match');
        if ~isempty(match)
            r = regexprep(match, '\((.*)\)', '$1'); % cannot do this in the rl directly because there are '( and )'
            literalMatch = regexptranslate('escape', match);
            rf = regexprep(rl, literalMatch, r);
            model.grRules(n) = cellstr(rf);
            disp(n)
        end
    end
    
    % correct cases of 'A and B and C or D and E or F and G' to '(A and B and C)
    % or (D and E) or (F and G)':
    % e.g. n=3468; n=5141;
    for n = 1:length(model.grRules)
        rl = model.grRules{n};
        match = regexp(rl, '((\d+\sand\s)+\d+(?=\sor))|((?<=or\s)(\d+\sand\s)+\d+)', 'match');
        if ~isempty(match)
            r = strcat('(', match, ')');
            rf = join(r, ' or ');
            model.grRules(n) = rf;
            disp(n)
        end
    end
    
    % get rid of repeated complexes in '.grRules', for cases like:
    % '(A and B and C) or (A and B and C)'
    % e.g. n=12;n=129;n=903;
    % n=863;n=864;n=865;
    % n=559; n=1041
    % n=150;n=195;n=324;n=558; ...;n=938;n=1176;n=3180;n=3344;n=3678;n=3683; 
    % rl = '221 or 435 or (70 and 897)'
    % rl = '221 and 435 and 221 and 897'
    for n = 1:length(model.grRules)
        rl = model.grRules{n};
        % if there is a 'and':
        match = regexp(rl, 'and', 'match');
        if ~isempty(match)
            % if there are no 'or' - only 'and':
            orpresent = regexp(rl, 'or', 'match');
            if isempty(orpresent)
                rlss = regexprep(rl, '\((.*)\)', '$1'); % remove parenthesis
                andsplit = regexp(rlss, ' and ', 'split');
                andjoin = join(unique(andsplit), ' and ');
                model.grRules(n) = andjoin;
                sprintf('only and: %d', n)
            % if there are 'or' - as well as 'and':
            else
                rel = regexp(rl, ' or ', 'split');
                % sometimes the genes in complexes are the same,
                % but their order in gene rule is different - so use unique and sort:
                sorted = cell(1, length(rel));
                for j=1:length(rel)
                    sb = regexprep(rel{j}, '\((.*)\)', '$1'); % remove parenthesis temporarily
                    it = strcat('(', join(sort(unique(regexp(sb, ' and ', 'split'))), ' and '), ')');
                    sorted(1, j) = regexprep(it,'\((\d+)\)', '$1'); % remove parenthesis in single genes
                end
                % if there are repeated complexes or genes, splited by 'or' unique will remove them:
                els = unique(sorted);
                gr = join(els, ' or ');
                % if there is only a complex - no 'or' after removing duplicates,
                % remove parenthesis:
                if isempty(regexp(gr{:}, 'or', 'match'))
                    gr = regexprep(gr, '\((.*)\)', '$1');
                end
                model.grRules(n) = gr;
                sprintf('and with or: %d', n)
            end
        end
    end
    [model.grRules, model.rxnGeneMat]=standardizeGrRules(model);
    % model = buildRxnGeneMat(model)
end

[edkerk]

Does the original model also have a rules field? ravenCobraWrapper currently uses this to reconstruct the grRules field. What happens if you remove the rules field before running ravenCobraWrapper, with:

liver_ml = rmfield(liver_ml,'rules');

[PkiwiBird]

It has a rules field. Doing that gives an error when running the liver_ml = ravenCobraWrapper(liver_ml);:

Converting RAVEN structure to COBRA..
Unrecognized field name "metComps".

Error in ravenCobraWrapper (line 143)
        newModel.metNames=strcat(model.metNames,' [',model.compNames(model.metComps),']');
 

[edkerk]

Right, my suggested workaround does not work as ravenCobraWrapper uses the rules field to identify it as a COBRA model. Without this field, ravenCobraWrapper thinks that it is a RAVEN model that should be converted to COBRA instead of vice versa, as you observed.

The long-term fix is implemented here. In essence, it is this change that fixes it: SysBioChalmers/RAVEN@f8b0c44. If you temporarily make this same change in your local ravenCobraWrapper, then it should work without needing to delete the rules field.

A new RAVEN version with this fix included will probably be released in 1-2 weeks time.

[PkiwiBird]

Thanks!