ecFactory/code/constructMinimalMutant_random.m at c92ede331e7ce755f238c91ff8d652f759353c30 · SysBioChalmers/ecFactory · 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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
function [optMutant,remaining] = constructMinimalMutant_random(model,candidates,modelParam)
tol =-1E-12;%-1E-6;
%get mutant with all modifications
optMutant=model;
toRemove = [];
[candidates,~] = sortrows(candidates,{'priority' 'k_scores'},{'ascend' 'ascend'});
%candidates = candidates(randperm(height(candidates)),:);
for i=1:height(candidates)
    gene   = candidates.genes(i);
    action = candidates.actions(i);
    mutF   = 1;%candidates.OE(i);
    enzIdx = candidates.enz_pos(i);
    tempModel = optMutant;

    if enzIdx>0
        tempModel.ub(enzIdx) = 1.01*candidates.maxUsage(i);
        tempModel.lb(enzIdx) = 0;
        if tempModel.ub(enzIdx) <=tempModel.lb(enzIdx)
            tempModel.lb(enzIdx) = 0.95*tempModel.ub(enzIdx);
        end
       optMutant = tempModel;
        %for reactions without enzymatic reaction

    else
        if strcmpi(action,'OE')
            enzUsage = 1.01*candidates.maxUsage(i);
            if enzUsage <= 1E-15
                enzUsage = 1.01*candidates.maxUsage(i);
            end
        else
            enzUsage = candidates.pUsage(i);
            if strcmpi(action,'KO')
                action = {'KD'};
                enzUsage = 0;
            end
        end
        modifications = {gene action mutF};

        [tempModel,success] = getMutantModel(optMutant,modifications,enzUsage);
         if success
            optMutant = tempModel;
        else
            toRemove = [toRemove; i];
        end
    end

end

if ~isempty(toRemove)
    %disp('The following gene modifications are not compatible with the rest of remaining candidate targets')
    %disp(candidates(toRemove,[1 2 3 6]))
    candidates(toRemove,:) = [];
end

optMutant = setParam(optMutant,'obj',modelParam.targetIndx,1);
%obtain optimal production rate and yield
[mutSol_r,~] = solveECmodel(optMutant,model,'pFBA',modelParam.prot_indx);
[mutSol_y,~] = solveECmodel(optMutant,model,'pFBA',modelParam.CUR_indx);
OptprodR =  mutSol_y(modelParam.targetIndx);%mutSol_r(modelParam.targetIndx);
Optyield = mutSol_y(modelParam.targetIndx)/(mutSol_y(modelParam.CUR_indx));
%Randomize order of targets
levelCandidates = candidates(randperm(height(candidates)),:);
counter   = 0;
remaining = table();
for i=1:height(levelCandidates)
    %reverse modifications
    gene   = levelCandidates.genes(i);
    action = levelCandidates.actions(i);
    enzIdx = levelCandidates.enz_pos(i);
    short  = levelCandidates.shortNames{i};
    mutF = 1;
    tempMutant = optMutant;
    %revert mutation
    if enzIdx>0
        saturationOpt         =  mutSol_r(enzIdx)/(optMutant.ub(enzIdx)+1E-15);
        tempMutant.ub(enzIdx) = 1.01*model.ub(enzIdx);
        tempMutant.lb(enzIdx) = 0.99*model.lb(enzIdx);
        if tempMutant.ub(enzIdx) <=tempMutant.lb(enzIdx)
            tempMutant.lb(enzIdx) = 0.95*tempMutant.ub(enzIdx);
        end

    %for reactions without enzymatic reaction
    else
        enzUsage = 1E-12;
        if strcmpi(action,'KO')
            reversal = {'OE'};
        else
            if strcmpi(action,'KD')
                reversal = {'OE'};
            else
                reversal = {'KD'};
            end
        end
        modifications = {gene reversal mutF};
        tempMutant    = getMutantModel(optMutant,modifications,enzUsage);
        saturationOpt = NaN;
    end
    %Get max WT production rate
    mutsol_prod = solveECmodel(tempMutant,tempMutant,'pFBA',modelParam.prot_indx);
    %Get max WT production yield
    mutsol_yield = solveECmodel(tempMutant,tempMutant,'pFBA',modelParam.CUR_indx);
    mutprodR     = mutsol_prod(modelParam.targetIndx);
    mutyield     = mutsol_yield(modelParam.targetIndx)/(mutsol_yield(modelParam.CUR_indx));
    flag = true;
    if enzIdx>0 && numel(enzIdx)==1
        saturationM =  mutsol_yield(enzIdx)/(tempMutant.ub(enzIdx)+1E-15);
        if (strcmpi(action,'OE') & saturationM <= (model.ub(enzIdx)/levelCandidates.maxUsage(i)))
            flag = false;
        end

    else
        saturationM = NaN;
    end
    FC_y  = mutyield/Optyield;
    FC_p  = mutprodR/OptprodR;
    score = mean([FC_y,FC_p]);
    %Discard genes that don't affect the optimal phenotype
    %discard OE targets that show low saturation after reversing the
    %modification

    if isnan(score)
        score = 0;
    end

    if flag && ...
       (...
        (score<=1+tol) || ...
        (strcmpi(action,'OE') && saturationOpt >= 0.99) ...%%|| ((~strcmpi(action,'OE') && saturationOpt <= saturationM)) ...
       )

        remaining = [remaining;levelCandidates(i,:)];
        counter = counter+1;
    else
    	optMutant = tempMutant;
    end
end
remaining = sortrows(remaining,{'priority' 'k_scores'},{'ascend' 'descend'});
remaining= removevars(remaining,{'enz_pos' 'OE' 'minUsage' 'maxUsage' 'pUsage' 'minUsageBio' 'maxUsageBio' 'pUsageBio'});
end