Complex annotation

.gitignore

@@ -51,7 +51,9 @@ helpsearch*/
 *.pyc
 *.env
 
-# Non-complying tables #
-########################
-*.xls
+# Non-complying tables and files #
+##################################
+*.xls*
+*.doc*
+*.ppt*
 *.tab

README.md

@@ -33,7 +33,7 @@ This repository contains the current consensus genome-scale metabolic model of _
 
 | Taxonomy | Latest update | Version | Reactions | Metabolites | Genes |
 |:-------|:--------------|:------|:------|:----------|:-----|
-| _Saccharomyces cerevisiae_ | 14-May-2022 | develop | 4069 | 2749 | 1151 |
+| _Saccharomyces cerevisiae_ | 26-May-2022 | develop | 4069 | 2749 | 1162 |
 
 # Installation & usage
 

code/curateMetsRxnsGenes.m

@@ -36,7 +36,18 @@
 %   Output:
 %       newModel    curated RAVEN model structure
 
+if nargin==4
+    error('Provide both a ''rxnsInfo'' and a ''rxnsCoeffs'' file')
+end
+if nargin<4
+    rxnsInfo='none';
+    rxnsCoeffs='none';
+end
+if nargin<3
+    genesInfo='none';
+end
 newModel=model;
+
 %% Metabolites
 if ~strcmp(metsInfo,'none')
     fid = fopen(metsInfo);
@@ -302,6 +313,9 @@
 emptyMiriam     =   all(cellfun(@isempty,miriamValues),1);
 miriamName(emptyMiriam)     = [];
 miriamValues(:,emptyMiriam) = [];
+if ~isfield(newModel,[type 'Miriams']);
+    newModel.([type 'Miriams'])=cell(numel(newModel.([type 's'])),1);
+end
 if ~isempty(miriamName)
     for i=1:numel(miriamName)
         newModel = editMiriam(newModel,type,modelIndex,miriamName{i},miriamValues(:,i),'replace');

code/getEarlierModelVersion.m

@@ -11,7 +11,7 @@
 
 if any(regexp(version,'^\d+\.\d+\.\d+$')) % Tag is a version number
     tagpath=['git show refs/tags/v' version ':'];
-elseif any(contains(tag,{'main','develop'}))
+elseif any(contains(version,{'main','develop'}))
     tagpath=['git show ' version ':'];
 else
     error('version should be ''main'', ''develop'', or a specific release, e.g. ''8.1.0''')

code/loadYeastModel.m

@@ -36,4 +36,5 @@
         model = ravenCobraWrapper(model);
     end
 end
+disp('If there is only 1 warning and no errors, it can savely be ignored.')
 end

code/modelCuration/complexAnnotation.m

@@ -0,0 +1,28 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% complexAnnotation
+% Correct complex annotation
+% As for the reference of new GPRs, please find detailed information in:
+% data/modelCuration/complexAnnotation/complexAnnotation.tsv
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Load model
+cd ..
+model = loadYeastModel;
+
+% Add new genes
+newModel       = curateMetsRxnsGenes(model,'none','../data/modelCuration/complexAnnotation/complexAnnotationGenes.tsv');
+
+% Add gene standard name for new genes
+fid = fopen('../data/modelCuration/complexAnnotation/complexAnnotation.tsv');
+complexAnnot = textscan(fid,'%q %q %q %q %q %q %q','Delimiter','\t','HeaderLines',1);
+fclose(fid);
+newGPR.ID     = complexAnnot{1};
+newGPR.GPR    = complexAnnot{3};
+newModel=changeGrRules(newModel,newGPR.ID,newGPR.GPR);
+
+% Delete unused genes (if any)
+newModel = deleteUnusedGenes(newModel);
+
+% Save model:
+saveYeastModel(newModel)
+cd modelCuration

code/modelTests/essentialGenes.m

@@ -1,4 +1,4 @@
-function [accurancy,tp,tn,fn,fp] = essentialGenes
+function [accuracy,tp,tn,fn,fp] = essentialGenes(model)
 % essentialGenes
 %   Modify media + find essential genes in model. Adapted from:
 %   https://doi.org/10.1371/journal.pcbi.1004530
@@ -12,11 +12,11 @@
 %   Usage: [accurancy,tp,tn,fn,fp] = essentialGenes
 %
 
-initCobraToolbox
-cd ..
-model = loadYeastModel;
-model = ravenCobraWrapper(model);
-cd modelTests
+if nargin<1;
+    cd ..
+    model = loadYeastModel;
+    cd modelTests
+end
 ko_tol = 1e-6;
 
 %constraints from genotype: check the genotype of the strains used in deletion experiment
@@ -41,7 +41,7 @@
 exp_viable = intersect(exp_viable,verifiedORFs);
 
 %calculate the growth rate after the single gene deletion using the original model and update model
-grRatio = singleGeneDeletion(model);
+[~, ~, ~, ~, grRatio]=findGeneDeletions(model,'sgd','fba');
 mod_viable  = model.genes(grRatio >= ko_tol);
 mod_viable = intersect(mod_viable,verifiedORFs);
 mod_inviable = model.genes(grRatio < ko_tol );
@@ -59,7 +59,7 @@
 %compare the prediction performances of two models
 %prediction accuracy was used to evaluate the quality of model update in
 %each step
-accurancy = (n_tp+n_tn)/(n_tp+n_tn+n_fn+n_fp);
+accuracy = (n_tp+n_tn)/(n_tp+n_tn+n_fn+n_fp);
 sensitivity = (100*n_tp/(n_tp+n_fn));
 specificity = (100*n_tn/(n_tn+n_fp));
 positivePredictive = (100*n_tp/(n_tp+n_fp));
@@ -75,7 +75,7 @@
 
     % start with a clean slate: set all exchange reactions to upper bound =
     % 1000 and lower bound = 0 (ie, unconstrained excretion, no uptake)
-    exchangeRxns = findExcRxns(model);
+    [~, exchangeRxns] = getExchangeRxns(model);
     model.lb(exchangeRxns) = 0;
     model.ub(exchangeRxns) = 1000;
 
@@ -106,12 +106,9 @@
                     'r_4600'; ... % Ca(2+) exchange
                     'r_2020' };
 
-    constrainedUptakeRxnIndexes = findRxnIDs(model,constrainedUptake);
-    glucoseExchangeIndex = findRxnIDs(model,glucoseExchange);
-    unconstrainedUptakeRxnIndexes = findRxnIDs(model,unconstrainedUptake);
-    model.lb(constrainedUptakeRxnIndexes) = -0.5;
-    model.lb(glucoseExchangeIndex) = -20;
-    model.lb(unconstrainedUptakeRxnIndexes) = -1000;
+    model=setParam(model,'lb',constrainedUptake,-0.5);
+    model=setParam(model,'lb',glucoseExchange,-20);
+    model=setParam(model,'lb',unconstrainedUptake,-1000);
 end
 
 function genes = inviableORFs