Feat/replicate analysis

ComplementaryScripts/Simulation/comparativeFVA_humanModels.m

@@ -1,4 +1,4 @@
-function FVA_Dists = comparativeFVA_humanModels(cellLine)
+function results = comparativeFVA_humanModels(cellLine)
 %comparativeFVA_humanModels
 %
 % Function that runs a comparative flux variabiability analysis between a
@@ -8,36 +8,57 @@
 % cellLine  (string) cell-line name (It should be consistent with the name
 %           of the subfolder in which the model is stored.
 %
-% FVA_Dists (2x1 cell) Flux variability distributions for both GEM and
-%           ecGEM in [mmol/gDw h]
+% results   (table) Contains rxnIDs, rxn formulas, Flux variability 
+%           distributions for both GEM and ecGEM in [mmol/gDw h] and
+%           metabolic subSystems for all rxns for which a variability range
+%           in bothmodels was calculated
 %
 % Usage: FVA_Dists = comparativeFVA_humanModels(cellLine)
 %
-%   Ivan Domenzain, 2019-05-14
+%   Ivan Domenzain, 2019-12-05
 
 current    = pwd;
 %Clone GECKO and pull comparativeFVa branch
 git('clone https://github.com/SysBioChalmers/GECKO.git')
 cd GECKO
 GECKO_path = pwd;
-git ('checkout fix/comparativeFVA')
+git ('checkout fix/comparative_FVA')
 git pull
 %Load GEM and ecGEM
-load(['../../models/humanGEM_cellLines/' cellLine '/model_modified.mat'])
-load(['../../models/humanGEM_cellLines/' cellLine '/ecModel_batch.mat'])
+load(['../../../models/' cellLine '/' cellLine '.mat'])
+load(['../../../models/' cellLine '/ecModel_batch.mat'])
+eval(['model = ' cellLine ';'])
 %Set medium constraints
-glucBound   = 1;
-oxBound     = 1000;
-stdBound    = 1;
-[model,~]   = setEMEMmedium(model_modified,glucBound,oxBound,stdBound,false);
-[ecModel,~] = setEMEMmedium(ecModel_batch,glucBound,oxBound,stdBound);
+cd (current)
+model   = removeMetFields(model);
+model   = setHamsMedium(model,false,'glucose',-1);
+ecModel = setHamsMedium(ecModel_batch,true,'glucose',1);
 evalin( 'base', 'clear(''model_modified'')' )
 evalin( 'base', 'clear(''ecModel_batch'')' )
 %Use GECKO built-in function for FVA comparative analysis
 cd ([GECKO_path '/geckomat/utilities/FVA'])
-CsourceUptk       = 'HMR_9034';
-[FVA_Dists,~,~,~] = comparativeFVA(model,ecModel,CsourceUptk,false);
+CsourceUptk = 'HMR_9034';
+[FVA_Dists,indexes,~,~] = comparativeFVA(model,ecModel,CsourceUptk,false,1E-8);
 cd (current)
-cd (['../../models/humanGEM_cellLines/' cellLine])
-save('FVA_results.mat','FVA_Dists','indexes')
+%Write results in a table and save it as .txt file
+mkdir('../../Results/FVA')
+variables = {'rxns' 'formulas' 'model_ranges' 'ecModel_ranges' 'subSystems'};
+formulas  = constructEquations(model,indexes);
+results   = table(model.rxns(indexes),formulas,FVA_Dists{1},FVA_Dists{2},model.subSystems(indexes),'VariableNames',variables);
+writetable(results,['../../Results/FVA/FVA_comp_' cellLine],'Delimiter','\t','QuoteStrings',false)
+end
+%--------------------------------------------------------------------------
+function model = removeMetFields(model)
+if isfield(model,'inchis')
+    model = rmfield(model,'inchis');
+end
+if isfield(model,'metMiriams')
+    model = rmfield(model,'metMiriams');
+end
+if isfield(model,'metCharges')
+    model = rmfield(model,'metCharges');
+end
+if isfield(model,'metFrom')
+    model = rmfield(model,'metFrom');
+end
 end

ComplementaryScripts/Simulation/setHamsMedium.m

@@ -0,0 +1,152 @@
+function exchModel = setHamsMedium(model,irrev,measuredMets,fluxes)
+% setHamsMedium
+%
+% Set a Hams culture medium for a humanGEM based model. This function works
+% for either standard or enzyme constrained-GEMs
+%
+%   model           An ihuman-based GEM
+%   irrev           (logical) TRUE if model comes in an irreversible format.
+%                   Default = false
+%   measuredMets    (string) metNames for measured compounds. Optional
+%   fluxes          (vector) Measured fluxes [mmol/gDw h]. Optional
+%                   NOTE: NEGATIVE = CONSUME
+%                         POSITIVE = PRODUCE
+%
+%   exchModel       (struct) Model with Ham's media constraints
+%
+% Ivan Domenzain.      Last edited: 2019-11-29
+
+if nargin<4
+	fluxes = [];
+    if nargin<3
+    	measuredMets = [];
+        if nargin<2 || isempty(irrev)
+        	irrev = false;
+        end
+    end
+end
+%Remove unconstrained field, if available
+if isfield(model,'unconstrained')
+    model = rmfield(model,'unconstrained');
+end
+%Check if boundary metabolites are present, if so then remove them
+boundaryIndx = find(strcmpi(model.compNames,'Boundary'));
+boundary     = find(model.metComps==boundaryIndx);
+if ~isempty(boundary)
+    model = removeMets(model,boundary,false,false,false,true);
+end
+%Ham's media composition
+mediaComps ={'glucose'
+             'arginine'
+             'histidine'
+             'lysine'
+             'methionine'
+             'phenylalanine'
+             'tryptophan'
+             'tyrosine'
+             'alanine'
+             'glycine'
+             'serine'
+             'threonine'
+             'aspartate'
+             'glutamate'
+             'asparagine'
+             'glutamine'
+             'isoleucine'
+             'leucine'
+             'proline'
+             'valine'
+             'cysteine'
+             'thiamin'
+             'hypoxanthine'
+             'folate'
+             'biotin'
+             'pantothenate'
+             'choline'
+             'inositol'
+             'nicotinamide'
+             'pyridoxine'
+             'riboflavin'
+             'thymidine'
+             'aquacob(III)alamin'
+             'lipoic acid'
+             'sulfate'
+             'linoleate'
+             'linolenate'
+             'O2'
+             'H2O'
+             'retinoate'
+             'Fe2+'
+             'Pi'
+             'alpha-tocopherol'
+             'gamma-tocopherol'};
+
+%Default flux bounds [LB, UB]
+fluxBounds = [-ones(length(mediaComps),1), ones(length(mediaComps),1)]*1000;
+%Check if provided mets are part of media's formulation
+if ~isempty(measuredMets)
+    %Modify fluxBounds with the correspondent provided flux measurements
+    [iA,iB] = ismember(measuredMets,mediaComps);
+    fluxBounds(iB(iA),:) = fluxes(iA).*ones(sum(iA),2);  % force flux to be equal to measured value (LB = UB = measured val)
+    if any(~iA)
+        %If measured mets are not in media formulation, then add them
+        mediaComps = [mediaComps; measuredMets(~iA)];
+        fluxBounds = [fluxBounds; fluxes(~iA).*ones(sum(~iA),2)];
+    end
+end
+
+if ~irrev
+    modelStr = 'model';
+    %Set uptake fluxes for media mets
+    [exchModel,unusedMets] = setExchangeBounds(model,mediaComps,fluxBounds(:,1),fluxBounds(:,2),true);
+else
+    modelStr = 'ecModel';
+    [exchRxns,exchIndxs] = getExchangeRxns(model);
+    %Exclude protein pool exchange
+    exchIndxs = exchIndxs(1:end-1);
+    exchRxns  = exchRxns(1:end-1);
+    %Differentiate between uptakes and production reactions
+    uptkIndxs = exchIndxs(find(contains(exchRxns,'_REV')));
+    prodIndxs = exchIndxs(find(~contains(exchRxns,'_REV')));
+    %Open all production reactions
+    exchModel = setParam(model,'ub',prodIndxs,1000);
+    %close all uptakes
+    exchModel = setParam(exchModel,'ub',uptkIndxs,0);
+    %Open uptake of media components one by one
+    unusedMets = [];
+    for i=1:length(mediaComps)
+        %Get metabolite indx
+        metIndx = getIndexes(model,strcat(mediaComps{i},'[s]'),'metscomps');
+        %Get rxns for metabolite
+        metRxns = find(model.S(metIndx,:));
+        %Get the uptake reaction for the metabolite
+        metExchRxn = intersect(metRxns,uptkIndxs);
+        if isempty(metExchRxn)
+            unusedMets = [unusedMets; mediaComps(i)];
+        else
+            %Open it!
+            exchModel.ub(metExchRxn) = -fluxBounds(i,1);  % flip sign of LB
+            %If the metabolite was measured, also set its production rate
+            if fluxBounds(i,2) ~= 1000
+                metExchRxn = intersect(metRxns,prodIndxs);
+                exchModel.ub(metExchRxn) = fluxBounds(i,2);
+            end
+        end
+    end
+end
+
+% report unused metabolites
+if ~isempty(unusedMets)
+    fprintf('WARNING: The following metabolites are either not in the model or do not have exchange reactions:\n');
+    fprintf('\t%s\n',unusedMets{:});
+end
+
+%Check if model is feasible
+sol = solveLP(exchModel); 
+if ~isempty(sol.x)
+	disp(['Constrained ' modelStr ' is feasible'])
+else
+    disp(['Constrained ' modelStr ' is unfeasible'])
+	exchModel = [];
+end
+end

ComplementaryScripts/enhanceGEM_cellLine.m

@@ -16,50 +16,40 @@
 %
 % Usage: [ecModel,ecModel_batch] = enhanceGEM_cellLine(cellName)
 %
-% Ivan Domenzain.      Last edited: 2019-05-15
+% Ivan Domenzain.      Last edited: 2019-11-20
 
 current  = pwd;
 org_name = 'homo sapiens';
 %Load original model 
-model    = load(['../models/' cellName '/' cellName '.mat']);
+model = load(['../models/' cellName '/' cellName '.mat']);
 eval(['model = model.' cellName])
-%model = model.model;
-%%%%%%%%%%%%%%%%%%%%%%%% model preprocessing %%%%%%%%%%%%%%%%%%%%%%%%%%
-%model_modified = ravenCobraWrapper(model);
+%% model preprocessing 
 model_modified = modelModifications(model);
 % Save models
-cd (['../models/' cellName])
-save('model.mat','model')
-save('model_modified.mat','model_modified')
-%%%%%%%%%%%%%%%%%%%%%%%% GECKO modifications %%%%%%%%%%%%%%%%%%%%%%%%%%
+save(['../models/' cellName '/model_modified.mat'],'model_modified')
+%% GECKO pipeline 
 % Retrieve kcats & MWs for each rxn in the model from Uniprot database:
-cd ([GECKO_path '/geckomat/get_enzyme_data'])
+cd GECKO/geckomat/get_enzyme_data
 model_data = getEnzymeCodes(model_modified);
 %Tries to Match kinetic coefficients to every reaction with a non empty
 %grRule
 kcats = matchKcats(model_data,org_name);
 %Save ecModel data
-cd (current)
-cd (['../models/' cellName])
-mkdir Data
-cd Data
-save('enzData.mat','model_data','kcats');
-cd (current)
+newDir = ['../../../../models/' cellName '/Data'];
+mkdir(newDir)
+save([newDir '/enzData.mat'],'model_data','kcats');
 %GEt ecModel matlab structure
+cd ../../..
 model_data = removeFields(model_data);
-cd ([GECKO_path '/geckomat/change_model'])
+cd GECKO/geckomat/change_model
 ecModel = readKcatData(model_data,kcats);
 % Save output models:
-cd (current)
-cd (['../models/' cellName])
-save('ecModel.mat','ecModel')
-cd ([GECKO_path '/geckomat/limit_proteins'])
+save(['../../../../models/' cellName '/ecModel.mat'],'ecModel')
+cd ../limit_proteins
 % Constrain total protein pool
-Ptotal       = 0.609; %HepG2 total protein content [g prot/gDw]
+Ptotal       = 0.593; %HepG2 total protein content [g prot/gDw]
 protCoverage = 0.5;
 sigma        = 0.5;
 [ecModel_batch,~,~] = constrainEnzymes(ecModel,Ptotal,sigma,protCoverage);
-cd (['../../../../models/' cellName])
-save('ecModel_batch.mat','ecModel_batch')
-cd (current)
+save(['../../../../models/' cellName '/ecModel_batch.mat'],'ecModel_batch')
 end

ComplementaryScripts/generate_human_ecModels_NCI60.m

@@ -4,32 +4,34 @@
 % NCI60 cell-lines) with enzyme-constraints with the use of the GECKO
 % pipeline.
 %
-% Ivan Domenzain.      Last edited: 2019-05-15
+% Ivan Domenzain.      Last edited: 2019-11-20
 
-load('../models/humanGEM_cellLines/11models.mat')
-modelNames = who;
-current    = pwd;
 %Clone GECKO and substitute human-specific scripts
 git('clone https://github.com/SysBioChalmers/GECKO.git')
-GECKO_path =  [current '/GECKO'];
 %Replace scripts in GECKO:
 fileNames = dir('GECKO_humanFiles');
 for i = 1:length(fileNames)
     fileName = fileNames(i).name;
-    if ~strcmp(fileName,'.') && ~strcmp(fileName,'..')
+    if ~strcmp(fileName,'.') && ~strcmp(fileName,'..') && ~strcmp(fileName,'.DS_Store')
         fullName   = ['GECKO_humanFiles/' fileName];
         GECKOpath = dir(['GECKO/**/' fileName]);
         GECKOpath = GECKOpath.folder;
         copyfile(fullName,GECKOpath)
     end
 end
+clear
+clc
+load('../models/humanGEM_cellLines/11models.mat')
+modelNames = who;
+current    = pwd;
+GECKO_path =  [current '/GECKO'];
 %Generate enzyme-constrained models
 for i=1:length(modelNames)
     cd (current)
     cellName = modelNames{i};
     mkdir (['../models/' cellName])
     cd (['../models/' cellName])
-    save([cellName '.mat'])
+    save([cellName '.mat'],cellName)
     mkdir Data
     cd (current)
     %Models mat files are saved in their respective folder by enhanceGEM_cellLine

ComplementaryScripts/modelModifications.m

@@ -1,10 +1,7 @@
 function model = modelModifications(model)
-[grRules, rxnGeneMat] = standardizeGrRules(model);
-model.grRules         = grRules;
-model.rxnGeneMat      = rxnGeneMat;
-model.b               = zeros(length(model.mets),1);
-% Standardizes the metabolites names
-model  = modifyMetNames(model);
-% Substitute biomass reaction
-model  = substituteBiomassRxns(model,false);
+[GRR, RGM]       = standardizeGrRules(model);
+model.grRules    = GRR;
+model.rxnGeneMat = RGM;
+model.b          = zeros(length(model.mets),1);
+model            = setParam(model,'obj','biomass_human',1);
 end