getModelFromKEGG: build the global model from raven-data artefacts

.gitignore

@@ -38,7 +38,13 @@ Temporary Items
 #external/kegg/keggMets.mat
 #external/kegg/keggPhylDist.mat
 #external/kegg/keggRxns.mat
-external/kegg/keggModel.mat
+# getModelFromKEGG downloads and assembles these at reconstruction/kegg/ on
+# first use (see buildGlobalKEGGModel); none of them belong in the repo. The
+# glob (rather than a version-pinned name) keeps this working across a
+# KEGG_VERSION.md bump without needing its own update.
+reconstruction/kegg/keggModel.mat
+reconstruction/kegg/kegg*_core.tar.gz
+reconstruction/kegg/kegg*_core/
 #software
 software/scip/
 software/blast+/

installation/downloadRavenBinaries.m

@@ -29,7 +29,6 @@ function downloadRavenBinaries(tools)
 end
 
 ravenDir = findRAVENroot();
-base     = 'https://github.com/SysBioChalmers/raven-data/releases/download';
 
 % raven-data platform key for the per-platform binary ZIPs.
 if ispc
@@ -80,16 +79,7 @@ function downloadRavenBinaries(tools)
         continue;   % already provisioned
     end
 
-    url = [base '/' tag '/' asset];
-    zipPath = fullfile(ravenDir,'software',asset);
-    fprintf('Downloading %s from raven-data ...\n',tool);
-    try
-        websave(zipPath,url);
-    catch
-        error(['Failed to download %s from %s\n' ...
-               'Check your internet connection, or fetch the offline ' ...
-               '"*-binaries" RAVEN release.'],tool,url);
-    end
+    zipPath = fetchRavenDataAsset(fullfile(ravenDir,'software'),tag,asset);
     unzip(zipPath,destDir);
     delete(zipPath);
 

installation/fetchRavenDataAsset.m

@@ -0,0 +1,43 @@
+function localPath=fetchRavenDataAsset(destDir,releaseTag,assetName)
+% fetchRavenDataAsset  Download (if needed) a raven-data release asset.
+%
+% Downloads assetName from the given raven-data release tag into destDir
+% and returns its local path. An already-downloaded copy is reused as-is.
+% Generic across raven-data's release families (KEGG artefacts, HMM
+% libraries, binaries, ...): nothing here is specific to any one of
+% them, so callers supply the release tag and asset name rather than
+% this function hardcoding either.
+%
+% Parameters
+% ----------
+% destDir : char
+%     directory to download into (created if it does not exist yet).
+% releaseTag : char
+%     the raven-data release tag, e.g. 'kegg118'.
+% assetName : char
+%     the release asset file name, e.g. 'kegg118_core.tar.gz'.
+%
+% Returns
+% -------
+% localPath : char
+%     path to the downloaded file.
+
+if ~isfolder(destDir)
+    mkdir(destDir);
+end
+localPath=fullfile(destDir,assetName);
+if isfile(localPath)
+    return;
+end
+fprintf(['Downloading ' assetName '... ']);
+try
+    websave(localPath,['https://github.com/SysBioChalmers/raven-data/releases/download/' releaseTag '/' assetName]);
+catch ME
+    if strcmp(ME.identifier,'MATLAB:webservices:HTTP404StatusCodeError')
+        error('Failed to download %s, the server returned a 404 error, try again later. If the problem persists please report it on the RAVEN GitHub Issues page: https://github.com/SysBioChalmers/RAVEN/issues',assetName)
+    else
+        rethrow(ME);
+    end
+end
+fprintf('COMPLETE\n');
+end

installation/findRAVENroot.m

@@ -7,9 +7,19 @@
 
 ST=dbstack('-completenames');
 prevDir = pwd();
+% A stored preference is only trusted if it still points at a real RAVEN
+% install; otherwise fall through to walking up from the currently
+% executing copy of this file. Without this check, a stale preference
+% left over from a different RAVEN checkout on the same machine silently
+% resolves to that other copy's data, not the one actually running.
+ravenPath = '';
 if ispref('RAVEN','ravenPath')
-    ravenPath = getpref('RAVEN','ravenPath');
-else
+    prefPath = getpref('RAVEN','ravenPath');
+    if isfile(fullfile(prefPath,'RAVEN.png'))
+        ravenPath = prefPath;
+    end
+end
+if isempty(ravenPath)
     ravenPath = ST(strcmp({ST.name},'findRAVENroot')).file;
     rootFound = 0;
     while rootFound == 0

reconstruction/kegg/KEGG_VERSION.md

@@ -0,0 +1,12 @@
+kegg118
+
+This is the raven-data release tag that every KEGG-artefact download in
+this folder reads from --- getModelFromKEGG's core reference-model/table
+bundle, and getKEGGModelForOrganism's HMM libraries. It is the single
+place this version is recorded; every function that downloads a KEGG
+artefact reads it from here via keggDataVersion.m, instead of
+hardcoding a "kegg###" string of its own. Bump this file, not the
+calling code, when raven-data publishes a new KEGG release.
+
+The version is read as the file's first line; everything below this
+point is only for a human reader.

reconstruction/kegg/buildGlobalGPR.m

@@ -0,0 +1,54 @@
+function [genes,rxnGeneMat]=buildGlobalGPR(rxns,koReaction,organismGeneKO)
+% buildGlobalGPR  Join every KEGG organism's genes onto the reference
+% reactions through their shared KO (KEGG Orthology) ids.
+%
+% Builds two sparse gene-KO / KO-reaction incidence matrices sharing one
+% KO axis, and multiplies them, rather than expanding the join row by
+% row --- organismGeneKO can carry millions of rows (every gene of every
+% KEGG organism), so a per-row loop is not viable.
+%
+% Parameters
+% ----------
+% rxns : cell array
+%     the reference model's reaction ids (model.rxns), in model order.
+% koReaction : table
+%     the ko_reaction table (columns 'ko', 'reaction') from readKEGGTable.
+% organismGeneKO : table
+%     the organism_gene_ko table (columns 'organism', 'gene', 'ko') from
+%     readKEGGTable.
+%
+% Returns
+% -------
+% genes : cell array
+%     sorted unique 'organism:gene' identifiers.
+% rxnGeneMat : sparse double
+%     numel(rxns) x numel(genes); rxnGeneMat(i,j)=1 when gene j shares at
+%     least one KO with reaction i.
+
+%Map each ko_reaction row onto the reference model's own reaction order;
+%rows naming a reaction absent from rxns (should not happen for a
+%consistent artefact set, but be defensive) are dropped.
+[inModel,rxnIdx]=ismember(koReaction.reaction,rxns);
+koCol=koReaction.ko(inModel);
+rxnIdxCol=rxnIdx(inModel);
+
+%Shared KO axis for both incidence matrices below.
+[kos,~,koGroupForRxnRow]=unique(koCol);
+numKOs=numel(kos);
+numRxns=numel(rxns);
+koRxn=sparse(koGroupForRxnRow,rxnIdxCol,1,numKOs,numRxns);
+
+%'organism:gene' identifiers, one row per (organism,gene,ko) triple. A
+%gene can carry more than one KO, so genes are de-duplicated separately.
+geneId=strcat(organismGeneKO.organism,':',organismGeneKO.gene);
+[genes,~,geneRowIdx]=unique(geneId);
+numGenes=numel(genes);
+
+%Rows whose KO is not linked to any reaction cannot contribute a gene-
+%reaction edge; drop them (organism_gene_ko is published pre-restricted
+%to linked KOs, so this is normally a no-op).
+[inKOs,koIdxForGeneRow]=ismember(organismGeneKO.ko,kos);
+geneKO=sparse(geneRowIdx(inKOs),koIdxForGeneRow(inKOs),1,numGenes,numKOs);
+
+rxnGeneMat=double((geneKO*koRxn)'>0);
+end

reconstruction/kegg/buildGlobalKEGGModel.m

@@ -0,0 +1,74 @@
+function [model,KOModel,isSpontaneous,isUndefinedStoich,isIncomplete,isGeneral]=buildGlobalKEGGModel(ravenPath)
+% buildGlobalKEGGModel  Download and assemble the global KEGG model.
+%
+% Fetches the raven-data <kegg version>_core.tar.gz bundle (the
+% gene-free reference model plus the ko_reaction / organism_gene_ko /
+% rxn_flags relational tables --- see raven-toolbox's
+% docs/maintenance/kegg_data_format.md), then joins every organism's
+% genes onto the reference reactions through their shared KO ids. This
+% is the data getModelFromKEGG used to load from a pre-built
+% keggModel.mat; see RAVEN issue #704. The KEGG version is read from
+% keggDataVersion, not hardcoded here.
+%
+% Parameters
+% ----------
+% ravenPath : char
+%     the RAVEN root directory, as returned by findRAVENroot.
+%
+% Returns
+% -------
+% model : struct
+%     the full global KEGG model (all reactions/metabolites, genes and
+%     rxnGeneMat spanning every KEGG organism). Callers narrow this down
+%     to one organism (see getKEGGModelForOrganism).
+% KOModel : struct
+%     a minimal model struct whose rxns are the KO ids linked to at
+%     least one reaction (getKEGGModelForOrganism's HMM-search path uses
+%     this as a KO id lookup table; nothing else of KOModel is read).
+% isSpontaneous, isUndefinedStoich, isIncomplete, isGeneral : cell arrays
+%     reaction ids carrying the corresponding rxn_flags quality flag.
+
+kver=keggDataVersion();
+keggDir=fullfile(ravenPath,'reconstruction','kegg');
+archive=fetchRavenDataAsset(keggDir,kver,[kver '_core.tar.gz']);
+
+coreDir=fullfile(keggDir,[kver '_core']);
+if ~isfolder(coreDir)
+    fprintf('Extracting the KEGG core artefacts... ');
+    untar(archive,coreDir);
+    fprintf('COMPLETE\n');
+end
+
+fprintf('Reading the KEGG reference model... ');
+refFileGz=fullfile(coreDir,[kver '_reference_model.yml.gz']);
+refFile=refFileGz(1:end-3);
+if ~isfile(refFile)
+    gunzip(refFileGz);
+end
+model=readYAMLmodel(refFile);
+fprintf('COMPLETE\n');
+
+fprintf('Reading the KEGG relational tables... ');
+koReaction=readKEGGTable(fullfile(coreDir,[kver '_ko_reaction.tsv.gz']));
+rxnFlags=readKEGGTable(fullfile(coreDir,[kver '_rxn_flags.tsv.gz']));
+organismGeneKO=readKEGGTable(fullfile(coreDir,[kver '_organism_gene_ko.tsv.gz']));
+fprintf('COMPLETE\n');
+
+isSpontaneous=flaggedReactions(rxnFlags,'spontaneous');
+isUndefinedStoich=flaggedReactions(rxnFlags,'undefined_stoich');
+isIncomplete=flaggedReactions(rxnFlags,'incomplete');
+isGeneral=flaggedReactions(rxnFlags,'general');
+
+KOModel.id='KOModel';
+KOModel.description='KEGG Orthology ids linked to at least one reaction';
+KOModel.rxns=unique(koReaction.ko);
+
+fprintf('Joining organism genes onto the reference reactions (this can take a while for the full KEGG gene set)... ');
+[model.genes,model.rxnGeneMat]=buildGlobalGPR(model.rxns,koReaction,organismGeneKO);
+fprintf('COMPLETE\n');
+end
+
+function ids=flaggedReactions(rxnFlags,column)
+mask=strcmpi(rxnFlags.(column),'true') | strcmp(rxnFlags.(column),'1');
+ids=rxnFlags.reaction(mask);
+end

reconstruction/kegg/getKEGGModelForOrganism.m

@@ -199,8 +199,12 @@
 %gzip-compressed flatfile, queried in one hmmsearch); if it is not already
 %present it is downloaded from the corresponding raven-data release
 %(https://github.com/SysBioChalmers/raven-data).
-if ~isempty(dataDir)
-    hmmOptions={'kegg118_eukaryotes','kegg118_prokaryotes'};
+%Only needed for the protein-homology path (fastaFile supplied): the
+%annotation-only path never touches libraryFile, so skip the (100+ MB)
+%download/extraction entirely when there is no FASTA file to search.
+if ~isempty(dataDir) && ~isempty(fastaFile)
+    kver=keggDataVersion();
+    hmmOptions={[kver '_eukaryotes'],[kver '_prokaryotes']};
     if ~endsWith(dataDir,hmmOptions)
         error(['Pre-trained HMMs set is not recognised. dataDir must match one of: ' strjoin(hmmOptions,' or ')])
     end
@@ -218,15 +222,7 @@
         gunzip([libraryFile '.gz']);
         fprintf('COMPLETE\n');
     else
-        fprintf('Downloading the HMM library file... ');
-        try
-            websave([libraryFile '.gz'],['https://github.com/SysBioChalmers/raven-data/releases/download/kegg118/' hmmName '.hmm.gz']);
-        catch ME
-            if strcmp(ME.identifier,'MATLAB:webservices:HTTP404StatusCodeError')
-                error('Failed to download the HMM library file, the server returned a 404 error, try again later. If the problem persists please report it on the RAVEN GitHub Issues page: https://github.com/SysBioChalmers/RAVEN/issues')
-            end
-        end
-        fprintf('COMPLETE\n');
+        fetchRavenDataAsset(fileparts(libraryFile),kver,[hmmName '.hmm.gz']);
         fprintf('Extracting the HMM library file... ');
         gunzip([libraryFile '.gz']);
         fprintf('COMPLETE\n');

reconstruction/kegg/getModelFromKEGG.m

@@ -1,9 +1,14 @@
 function [model,KOModel]=getModelFromKEGG(varargin)
-% getModelFromKEGG  Load the pre-built global KEGG model.
+% getModelFromKEGG  Load the global KEGG model.
 %
-% Loads the pre-built global KEGG reaction/gene model from keggModel.mat.
-% The artefact is generated by the raven-toolbox Python package and
-% distributed as a raven-data release asset.
+% Loads the global KEGG reaction/gene model from keggModel.mat. On first
+% use --- when no keggModel.mat is present yet --- the underlying
+% artefacts (the gene-free reference model plus the KO/reaction/organism-
+% gene relational tables, published by the raven-toolbox Python package
+% as a raven-data release) are downloaded and assembled instead, and the
+% result is cached to keggModel.mat so later calls load instantly. The
+% first build can take a while (the organism-gene table covers every
+% KEGG organism) and needs a few hundred MB of disk.
 %
 % Name-Value Arguments
 % --------------------
@@ -45,15 +50,16 @@
 keepGeneral=p.keepGeneral;
 
 modelFile=fullfile(ravenPath,'reconstruction','kegg','keggModel.mat');
-if ~exist(modelFile,'file')
-    error('getModelFromKEGG:noModel', ...
-        ['keggModel.mat not found at ' strrep(modelFile,'\','/') '.\n' ...
-        'Generate it with the raven-toolbox Python package or download it ' ...
-        'via downloadRavenBinaries.']);
+if isfile(modelFile)
+    fprintf(['Importing the global KEGG model from ' strrep(modelFile,'\','/') '... ']);
+    load(modelFile,'model','KOModel','isSpontaneous','isUndefinedStoich','isIncomplete','isGeneral');
+    fprintf('COMPLETE\n');
+else
+    [model,KOModel,isSpontaneous,isUndefinedStoich,isIncomplete,isGeneral]=buildGlobalKEGGModel(ravenPath);
+    fprintf(['Saving the global KEGG model to ' strrep(modelFile,'\','/') ' for future use... ']);
+    save(modelFile,'model','KOModel','isSpontaneous','isUndefinedStoich','isIncomplete','isGeneral','-v7.3');
+    fprintf('COMPLETE\n');
 end
-fprintf(['Importing the global KEGG model from ' strrep(modelFile,'\','/') '... ']);
-load(modelFile);
-fprintf('COMPLETE\n');
 
 if keepSpontaneous==false
     model=removeReactions(model,intersect(isSpontaneous,model.rxns),true,true);

reconstruction/kegg/keggDataVersion.m

@@ -0,0 +1,23 @@
+function version=keggDataVersion()
+% keggDataVersion  The raven-data release tag for KEGG artefacts.
+%
+% Reads reconstruction/kegg/KEGG_VERSION.md --- the single place this
+% version is recorded --- so every function that downloads a KEGG
+% artefact (getModelFromKEGG, getKEGGModelForOrganism) stays in sync
+% when raven-data publishes a new KEGG release. The version is the
+% file's first line.
+%
+% Returns
+% -------
+% version : char
+%     the raven-data release tag, e.g. 'kegg118'.
+
+ravenPath=findRAVENroot();
+versionFile=fullfile(ravenPath,'reconstruction','kegg','KEGG_VERSION.md');
+fid=fopen(versionFile,'r');
+if fid==-1
+    error('keggDataVersion:fileNotFound','Cannot read %s.',strrep(versionFile,'\','/'))
+end
+version=strtrim(fgetl(fid));
+fclose(fid);
+end

reconstruction/kegg/readKEGGTable.m

@@ -0,0 +1,25 @@
+function tbl=readKEGGTable(gzFile)
+% readKEGGTable  Read a gzipped KEGG relational table (published by
+% raven-toolbox as gzipped TSV, see raven-toolbox's
+% docs/maintenance/kegg_data_format.md) into a table with every column
+% forced to char cellstr (matching RAVEN's own cell-array-of-char
+% convention, regardless of the caller's default text-import type).
+%
+% Parameters
+% ----------
+% gzFile : char
+%     path to the gzipped TSV file, e.g. '.../kegg118_ko_reaction.tsv.gz'.
+%
+% Returns
+% -------
+% tbl : table
+%     the table, with every column as a cellstr.
+
+plainFile=gzFile(1:end-3); %strip the trailing '.gz'
+if ~isfile(plainFile)
+    gunzip(gzFile);
+end
+opts=detectImportOptions(plainFile,'FileType','text','Delimiter',char(9));
+opts=setvartype(opts,opts.VariableNames,'char');
+tbl=readtable(plainFile,opts);
+end