cobratoolbox/changeGeneAssociation.m at master · wbryant/cobratoolbox · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
function model = changeGeneAssociation(model,rxnName,grRule,geneNameList,systNameList,addRxnGeneMat)
% Change gene associations in a model
%
% model = changeGeneAssociation(model,rxnName,grRule,geneName,systName)
%
%INPUTS
% model             COBRA Toolbox model structure
% rxnName           Name of the new reaction
% grRule            Gene-reaction rule in boolean format (and/or allowed)
%
%OPTIONAL INPUTS
% geneNameList      List of gene names (used only for translation from
%                   common gene names to systematic gene names)
% systNameList      List of systematic names
% addRxnGeneMat     adds rxnGeneMat to model structure (default = true)
%
%OUTPUT
% model             COBRA Toolbox model structure with new gene reaction
%                   associations
%

% Markus Herrgard 1/12/07
% Ines Thiele 08/03/2015, made rxnGeneMat optional

% IT: updated the nargin statement to accommodate the additional option
if exist('geneNameList','var') && exist('systNameList','var')
    translateNamesFlag = true;
else
    translateNamesFlag = false;
end

if ~exist('addRxnGeneMat','var')
    addRxnGeneMat = 1;
end

[isInModel,rxnID] = ismember(rxnName,model.rxns);

if (~isInModel)
    error(['Reaction ' rxnName ' not in the model']);
end

if ~isfield(model,'genes')
    model.genes = {};
end
nGenes = length(model.genes);
model.rules{rxnID} = '';
% IT 01/2010 - this line caused problems for xls2model.m
if addRxnGeneMat ==1
    model.rxnGeneMat(rxnID,:) = zeros(1,nGenes);
end
% Remove extra white space
grRule = regexprep(grRule,'\s{2,}',' ');
grRule = regexprep(grRule,'( ','(');
grRule = regexprep(grRule,' )',')');


if (~isempty(grRule))
    % Ronan & Stefan 13/9/2011 - moved this inside check if empty
    % Remove extra white space
    grRule = regexprep(grRule,'\s{2,}',' ');
    grRule = regexprep(grRule,'( ','(');
    grRule = regexprep(grRule,' )',')');
    [genes,rule] = parseBoolean(grRule);

    for i = 1:length(genes)
        if (translateNamesFlag)
            % Translate gene names to systematic names
            [isInList,translID] = ismember(genes{i},geneNameList);
            if isInList
                newGene = systNameList{translID};
                grRule = regexprep(grRule,[genes{i} '$'],newGene);
                grRule = regexprep(grRule,[genes{i} '\s'],[newGene ' ']);
                grRule = regexprep(grRule,[genes{i} ')'],[newGene ')']);
                genes{i} = newGene;
            else
                warning(['Gene name ' genes{i} ' not in translation list']);
            end
        end
        geneID = find(strcmp(model.genes,genes{i}));
        if (isempty(geneID))
            warning(['New gene ' genes{i} ' added to model']);
            % Append gene
            model.genes = [model.genes; genes(i)];
            nGenes = length(model.genes);
            if addRxnGeneMat == 1
                model.rxnGeneMat(rxnID,end+1) = 1;
            end
            rule = strrep(rule,['x(' num2str(i) ')'],['x(' num2str(nGenes) ')']);
        else
            if addRxnGeneMat == 1
                model.rxnGeneMat(rxnID,geneID) = 1;
            end
            rule = strrep(rule,['x(' num2str(i) ')'],['x(' num2str(geneID) ')']);
        end
    end
    model.rules{rxnID} = rule;
end

model.grRules{rxnID} = grRule;

%make sure variables are column vectors
model.rules = columnVector(model.rules);
model.grRules = columnVector(model.grRules);