|
| 1 | + function [model,enzUsages,modifications,GAM,massCoverage] = constrainEnzymes(model,f,GAM,Ptot,pIDs,data,gRate,c_UptakeExp,parameters) |
| 2 | +% constrainEnzymes |
| 3 | +% |
| 4 | +% Main function for overlaying proteomics data on an enzyme-constrained |
| 5 | +% model. If chosen, also scales the protein content, optimizes GAM, and |
| 6 | +% flexibilizes the proteomics data. |
| 7 | +% |
| 8 | +% model ecModel. |
| 9 | +% f (Opt) Estimated mass fraction of enzymes in model. |
| 10 | +% GAM (Opt) Growth-associated maintenance value. If not |
| 11 | +% provided, it will be fitted to chemostat data. |
| 12 | +% Ptot (Opt) Total protein content, provide if desired content |
| 13 | +% is different from the one reported in getModelParameters [gProt/gDw] |
| 14 | +% pIDs (Opt) Protein IDs from proteomics data. |
| 15 | +% data (Opt) Protein abundances from proteomics data [mmol/gDW]. |
| 16 | +% gRate (Opt) Experimental growth rate at which the proteomics |
| 17 | +% data were obtained [1/h] |
| 18 | +% c_UptakeExp (Opt) Experimentally measured glucose uptake rate |
| 19 | +% [mmol/gDW h]. |
| 20 | +% |
| 21 | +% model ecModel with calibrated enzyme usage upper bounds |
| 22 | +% enzUsages Calculated enzyme usages after final calibration |
| 23 | +% (enzyme_i demand/enzyme_i upper bound) |
| 24 | +% modifications Table with all the modified values |
| 25 | +% (Protein ID/old value/Flexibilized value) |
| 26 | +% GAM Fitted GAM value for the ecModel |
| 27 | +% massCoverage Ratio between measured and total mass of protein in the model |
| 28 | +% |
| 29 | +% Usage: [model,enzUsages,modifications, GAM,massCoverage] = constrainEnzymes(model,f,GAM,Ptot,pIDs,data,gRate,c_UptakeExp) |
| 30 | +% |
| 31 | +% Benjamin J. Sanchez. Last update 2018-12-11 |
| 32 | +% Ivan Domenzain. Last update 2020-03-02 |
| 33 | +% |
| 34 | + |
| 35 | +%get model parameters |
| 36 | +if nargin<9 |
| 37 | + cd .. |
| 38 | + parameters = getModelParameters; |
| 39 | + cd limit_proteins |
| 40 | +end |
| 41 | +sigma = parameters.sigma; |
| 42 | +c_source = parameters.c_source; |
| 43 | + |
| 44 | +%Compute f if not provided: |
| 45 | +if nargin < 2 |
| 46 | + [f,~] = measureAbundance(model.enzymes); |
| 47 | +else |
| 48 | + if isempty(f) |
| 49 | + [f,~] = measureAbundance(model.enzymes); |
| 50 | + end |
| 51 | +end |
| 52 | +%Leave GAM empty if not provided (will be fitted later): |
| 53 | +if nargin < 3 |
| 54 | + GAM = []; |
| 55 | +end |
| 56 | +%Load Ptot if not provided: |
| 57 | +if nargin < 4 |
| 58 | + Ptot = parameters.Ptot; |
| 59 | +end |
| 60 | +%No UB will be changed if no data is available -> pool = all enzymes(FBAwMC) |
| 61 | +if nargin < 5 |
| 62 | + pIDs = cell(0,1); |
| 63 | + data = zeros(0,1); |
| 64 | +end |
| 65 | +%Remove zeros or negative values |
| 66 | +data = cleanDataset(data); |
| 67 | +%Assign concentrations as UBs [mmol/gDW]: |
| 68 | +model.concs = nan(size(model.enzymes)); %OBS: min value is zero!! |
| 69 | +disp('Matching data to enzymes in model...') |
| 70 | +for i = 1:length(model.enzymes) |
| 71 | + match = false; |
| 72 | + for j = 1:length(pIDs) |
| 73 | + if strcmpi(pIDs{j},model.enzymes{i}) && ~match |
| 74 | + model.concs(i) = data(j)*model.MWs(i); %g/gDW |
| 75 | + rxn_name = ['prot_' model.enzymes{i} '_exchange']; |
| 76 | + pos = strcmpi(rxn_name,model.rxns); |
| 77 | + model.ub(pos) = data(j); |
| 78 | + match = true; |
| 79 | + end |
| 80 | + end |
| 81 | +end |
| 82 | +%Count mass of non-measured enzymes: |
| 83 | +measured = ~isnan(model.concs); |
| 84 | +concs_measured = model.concs(measured); |
| 85 | +Pmeasured = sum(concs_measured); |
| 86 | +%Get protein content in biomass pseudoreaction: |
| 87 | +Pbase = sumProtein(model); |
| 88 | +if Pmeasured > 0 |
| 89 | + %Expected total enzyme concentration |
| 90 | + enzymeConc=Ptot*f; |
| 91 | + %Non-measured part will be pooled |
| 92 | + Ppool=enzymeConc-Pmeasured; |
| 93 | + fs=Ppool/Pbase; |
| 94 | +else |
| 95 | + fs = f*sigma; |
| 96 | +end |
| 97 | +%Constrain the rest of enzymes with the pool assumption: |
| 98 | +if sum(strcmp(model.rxns,'prot_pool_exchange')) == 0 |
| 99 | + model = constrainPool(model,~measured,full(fs*Pbase)); |
| 100 | +end |
| 101 | +if sum(data)==0 |
| 102 | + %Modify protein/carb content and GAM: |
| 103 | + [model,GAM] = scaleBioMass(model,Ptot,GAM); |
| 104 | +end |
| 105 | +%Display some metrics: |
| 106 | +disp(['Total protein amount measured = ' num2str(Pmeasured) ' g/gDW']) |
| 107 | +disp(['Total enzymes measured = ' num2str(sum(measured)) ' enzymes']) |
| 108 | +disp(['Enzymes in model with 0 g/gDW = ' num2str(sum(concs_measured==0)) ' enzymes']) |
| 109 | +disp(['Total protein amount not measured = ' num2str(Ptot - Pmeasured) ' g/gDW']) |
| 110 | +disp(['Total enzymes not measured = ' num2str(sum(~measured)) ' enzymes']) |
| 111 | +disp(['Total protein in model = ' num2str(Ptot) ' g/gDW']) |
| 112 | +enzUsages = []; |
| 113 | +if nargin > 7 |
| 114 | + model = updateProtPool(model,Ptot,f*sigma); |
| 115 | + [tempModel,enzUsages,modifications] = flexibilizeProteins(model,gRate,c_UptakeExp,c_source); |
| 116 | + model = updateProtPool(tempModel,Ptot,f*sigma); |
| 117 | +end |
| 118 | +massCoverage = Pmeasured/Ptot; |
| 119 | +if isempty(enzUsages) |
| 120 | + enzUsages = table({},zeros(0,1),'VariableNames',{'prot_IDs' 'usage'}); |
| 121 | + modifications = table({},zeros(0,1),zeros(0,1),zeros(0,1),'VariableNames',{'protein_IDs' 'previous_values' 'modified_values' 'flex_mass'}); |
| 122 | +else |
| 123 | + plotHistogram(enzUsages.usage,'Enzyme usage [-]',[0,1],'Enzyme usages','usages') |
| 124 | +end |
| 125 | +%Plot histogram (if there are measurements): |
| 126 | +%plotHistogram(concs_measured,'Protein amount [mg/gDW]',[1e-3,1e3],'Modelled Protein abundances','abundances') |
| 127 | +end |
| 128 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 129 | +function plotHistogram(variable,xlabelStr,xlimits,titleStr,option) |
| 130 | +if iscell(variable) |
| 131 | + cell2mat(variable); |
| 132 | +end |
| 133 | +if sum(variable) > 0 |
| 134 | + variable(variable==0) = 1E-15; |
| 135 | + figure |
| 136 | + if strcmpi(option,'abundances') |
| 137 | + hist(variable*1e3,10.^(-3:0.5:3)) |
| 138 | + set(gca,'xscale','log') |
| 139 | + else |
| 140 | + hist(variable,(0:0.05:1)) |
| 141 | + end |
| 142 | + xlim(xlimits) |
| 143 | + xlabel(xlabelStr) |
| 144 | + ylabel('Frequency'); |
| 145 | + title(titleStr) |
| 146 | +end |
| 147 | +end |
| 148 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 149 | +function data = cleanDataset(data) |
| 150 | +for i=1:length(data) |
| 151 | + if data(i)<=0 |
| 152 | + data(i) = NaN; |
| 153 | + end |
| 154 | +end |
| 155 | +end |
0 commit comments