acidFBA/minAminoAcidDiff.m at main · AlmaasLab/acidFBA · 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
function minAminoAcidDiff(model)
%minAminoAcidDiff Minimizes the Euclidean distance to the species-specific
%amino acid distributions across a range of growth rates. Calculates the
%absolute relative difference for each amino acid.
%
% USAGE:
%
%    minAminoAcidDiff(model)
%
% INPUTS:
%    model:             AcidFBA model structure
%
% .. Authors:
%       - Vetle Simensen 29/06/22

% Initialize
params = getParameters();
changeCobraSolver('gurobi','all');
changeCobraSolverParams('QP','feasTol',1.0000e-09);
changeCobraSolverParams('QP','optTol',1.0000e-09);

% Get amino acid mass fractions (S. cerevisiae, E. coli, B. subtilis)
currDir = cd;
cd data
aaFrac = readmatrix('aaMassFracsOrganisms.xlsx');
aaFrac = aaFrac(:,2:5);
[nAA,nOrg] = size(aaFrac);
cd(currDir);

% Wild type yeast reference conditions (experimental amino acid profile)
protIdx = findRxnIDs(model,'prot_pool_exchange');
aaYeastProfile = aaFrac(:,1);
wtModel = changeRxnBounds(model,params.aaDrains,aaYeastProfile * model.ub(protIdx),'b');
wtSol = optimizeCbModel(wtModel);

% Define range of relative growth rates
nVals = 35;
relGrowthRates = linspace(0.1,2.5,nVals);
yvals = zeros(nAA,nVals,nOrg);

% Each organism
for i = 1:nOrg
    % Select experimental organism-specific amino acid distribution
    aaProfile = model.ub(protIdx) * aaFrac(:,i);

    % Each fraction of optimal yeast growth
    for j = 1:length(relGrowthRates)
        simModel = model;

        % Set fraction of optimal growth
        simModel = changeRxnBounds(simModel,simModel.rxns(logical(simModel.c)),relGrowthRates(j) * wtSol.f,'b');

        % Construct QP problem
        QPproblem = buildLPproblemFromModel(simModel);
        [~, n] = size(simModel.S);
        objIndxs = findRxnIDs(model,params.aaDrains);

        % Formulate quadratic objective term
        quadMat = sparse(n,n);
        objIndxsVec = zeros(n,1);
        objIndxsVec(objIndxs) = 1;
        quadMat = quadMat + diag(objIndxsVec);
        QPproblem.F = quadMat;

        % Formulate linear objective term
        linVec = zeros(n,1);
        linVec(objIndxs) = -aaProfile;
        QPproblem.c = linVec;

        % Add linear constraints
        QPproblem.A = simModel.S;  % LHS matrix
        QPproblem.b = simModel.b;  % RHS vector
        QPproblem.lb = simModel.lb;    % lower flux bounds
        QPproblem.ub = simModel.ub;    % upper flux bounds
        QPproblem.osense = 1;   % objective sense linear part (minimize)
        QPproblem.csense = buildLPproblemFromModel(simModel).csense;   % constraint senses (all 'E')

        % Solve QP problem
        momaSol = solveCobraQP(QPproblem);

        % Store solution
        yvals(:,j,i) = abs(((momaSol.full(objIndxs) ./ momaSol.full(protIdx)) - aaFrac(:,i)) ./ (aaFrac(:,i)));
    end
end

% Save data
writematrix(relGrowthRates,'data/orgRelGrowthRates.csv');   % relative growth rates
fileSuffixes = {'Sce' 'Eco' 'Bsu' 'Pse'};

for i = 1:nOrg
    writematrix(yvals(:,:,i),['data/aa' fileSuffixes{i} '.csv']);
end

end