calculateFfactor

 calculateFfactor
   Computes the f factor, as a proxy to the mass fraction of proteins
   accounted for in an ecModel out of the total protein content in cells.

 Input:
   model        an ecModel in GECKO 3 format (with ecModel.ec structure)
   protData     structure with proteome data, from loadProtData (Optional,
                by default it instead attempts to load data/paxDB.tsv)
   enzymes      list of enzymes (Optional, default model.ec.enzymes)
   modelAdapter a loaded model adapter (Optional, will otherwise use the
                default model adapter).

 Output:
   f            f-factor

 Usage:
   f = calculateFfactor(model, protData, enzymes, modelAdapter)

0001 function f = calculateFfactor(model, protData, enzymes, modelAdapter)
0002 % calculateFfactor
0003 %   Computes the f factor, as a proxy to the mass fraction of proteins
0004 %   accounted for in an ecModel out of the total protein content in cells.
0005 %
0006 % Input:
0007 %   model        an ecModel in GECKO 3 format (with ecModel.ec structure)
0008 %   protData     structure with proteome data, from loadProtData (Optional,
0009 %                by default it instead attempts to load data/paxDB.tsv)
0010 %   enzymes      list of enzymes (Optional, default model.ec.enzymes)
0011 %   modelAdapter a loaded model adapter (Optional, will otherwise use the
0012 %                default model adapter).
0013 %
0014 % Output:
0015 %   f            f-factor
0016 %
0017 % Usage:
0018 %   f = calculateFfactor(model, protData, enzymes, modelAdapter)
0019 
0020 if nargin < 4 || isempty(modelAdapter)
0021     modelAdapter = ModelAdapterManager.getDefault();
0022     if isempty(modelAdapter)
0023         error('Either send in a modelAdapter or set the default model adapter in the ModelAdapterManager.')
0024     end
0025 end
0026 params = modelAdapter.getParameters();
0027 
0028 if nargin < 3 || isempty(enzymes)
0029     enzymes = model.ec.enzymes;
0030 end
0031 
0032 % Gather proteome data in protData structure
0033 if nargin < 2 || isempty(protData)
0034     if exist(fullfile(params.path,'data','paxDB.tsv'),'file')
0035         protData = fullfile(params.path,'data','paxDB.tsv');
0036     else
0037         printOrange('WARNING: No proteomics data is provided or can be found. Default f value of 0.5 is returned.\n');
0038         f = 0.5;
0039     end
0040 end
0041 
0042 % Gather Uniprot database for finding MW
0043 uniprotDB = loadDatabases('uniprot', modelAdapter);
0044 uniprotDB = uniprotDB.uniprot;
0045 
0046 if ischar(protData) && endsWith(protData,'paxDB.tsv')
0047     fID         = fopen(fullfile(protData),'r');
0048     fileContent = textscan(fID,'%s','delimiter','\n');
0049     headerLines = find(startsWith(fileContent{1},'#'),1,'last');
0050     fclose(fID);
0051 
0052     %Read data file, excluding headerlines
0053     fID         = fopen(fullfile(protData),'r');
0054     fileContent = textscan(fID,'%s %s %f','delimiter','\t','HeaderLines',headerLines);
0055     genes       = fileContent{2};
0056     %Remove internal geneIDs modifiers
0057     genes       = regexprep(genes,'^\d+\.','');
0058     level       = fileContent{3};
0059     fclose(fID);
0060     [a,b]       = ismember(genes,uniprotDB.genes);
0061     uniprot     = uniprotDB.ID(b(a));
0062     level(~a)   = [];
0063     clear protData
0064     protData.uniprotIDs = uniprot;
0065     protData.level   = level;
0066     % Get MW and abundance (unit does not matter, f is fraction)
0067     [~,idx] = ismember(protData.uniprotIDs,uniprotDB.ID);
0068     protData.MW = uniprotDB.MW(idx);
0069     protData.abundances = protData.level .* protData.MW;
0070 end
0071 
0072 avgAbundances = mean(protData.abundances,2);
0073 totalProt = sum(avgAbundances,'omitnan');
0074 
0075 % Get enzymes in model
0076 enzymesInModel = ismember(protData.uniprotIDs,enzymes);
0077 totalEnz = sum(avgAbundances(enzymesInModel),'omitnan');
0078 
0079 f = totalEnz/totalProt;
0080 end