diff --git a/.gitignore b/.gitignore index b3c18a0f..814ebb4f 100644 --- a/.gitignore +++ b/.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+/ diff --git a/installation/downloadRavenBinaries.m b/installation/downloadRavenBinaries.m index e9fd063b..918448ea 100644 --- a/installation/downloadRavenBinaries.m +++ b/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); diff --git a/installation/fetchRavenDataAsset.m b/installation/fetchRavenDataAsset.m new file mode 100644 index 00000000..5e5266eb --- /dev/null +++ b/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 diff --git a/installation/findRAVENroot.m b/installation/findRAVENroot.m index af2336b4..8d1b1fcd 100755 --- a/installation/findRAVENroot.m +++ b/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 diff --git a/reconstruction/kegg/KEGG_VERSION.md b/reconstruction/kegg/KEGG_VERSION.md new file mode 100644 index 00000000..5095c18c --- /dev/null +++ b/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. diff --git a/reconstruction/kegg/buildGlobalGPR.m b/reconstruction/kegg/buildGlobalGPR.m new file mode 100644 index 00000000..9c88b6c6 --- /dev/null +++ b/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 diff --git a/reconstruction/kegg/buildGlobalKEGGModel.m b/reconstruction/kegg/buildGlobalKEGGModel.m new file mode 100644 index 00000000..c67e2af7 --- /dev/null +++ b/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 _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 diff --git a/reconstruction/kegg/getKEGGModelForOrganism.m b/reconstruction/kegg/getKEGGModelForOrganism.m index 72069f15..596e8cfb 100755 --- a/reconstruction/kegg/getKEGGModelForOrganism.m +++ b/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'); diff --git a/reconstruction/kegg/getModelFromKEGG.m b/reconstruction/kegg/getModelFromKEGG.m index 4ed093d7..10c17a77 100755 --- a/reconstruction/kegg/getModelFromKEGG.m +++ b/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); diff --git a/reconstruction/kegg/keggDataVersion.m b/reconstruction/kegg/keggDataVersion.m new file mode 100644 index 00000000..f8a8e9d0 --- /dev/null +++ b/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 diff --git a/reconstruction/kegg/readKEGGTable.m b/reconstruction/kegg/readKEGGTable.m new file mode 100644 index 00000000..3da37a2f --- /dev/null +++ b/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 diff --git a/testing/function_tests/tReconstruction.m b/testing/function_tests/tReconstruction.m index 39f3e27e..e85a7742 100644 --- a/testing/function_tests/tReconstruction.m +++ b/testing/function_tests/tReconstruction.m @@ -60,8 +60,23 @@ function getKEGGModelForOrganismNeedsData(testCase) function getModelFromKEGGNeedsData(testCase) matFile = fullfile(testCase.ravenRoot, 'reconstruction', 'kegg', 'keggModel.mat'); testCase.assumeFalse(exist(matFile, 'file') == 2, ... - 'keggModel.mat is present; skipping error-path test.'); - testCase.verifyError(@() getModelFromKEGG(), 'getModelFromKEGG:noModel'); + 'keggModel.mat is present; skipping the build-from-artefacts path.'); + testCase.assumeFail('Downloads and assembles the full KEGG artefact set from raven-data; not run automatically.'); + end + + function buildGlobalGPRJoinsGenesThroughSharedKO(testCase) + % Offline unit test of the join at the core of buildGlobalKEGGModel: + % two reactions sharing a KO, a KO used by two organisms' genes, and + % a reaction with no KO at all (should end up with no genes). + rxns = {'R1'; 'R2'; 'R3'}; + koReaction = table({'K1'; 'K1'; 'K2'}, {'R1'; 'R2'; 'R2'}, ... + 'VariableNames', {'ko', 'reaction'}); + organismGeneKO = table({'a'; 'a'; 'b'}, {'g1'; 'g2'; 'g1'}, {'K1'; 'K2'; 'K1'}, ... + 'VariableNames', {'organism', 'gene', 'ko'}); + [genes, rxnGeneMat] = buildGlobalGPR(rxns, koReaction, organismGeneKO); + testCase.verifyEqual(genes, {'a:g1'; 'a:g2'; 'b:g1'}); + expected = [1 0 1; 1 1 1; 0 0 0]; + testCase.verifyEqual(full(rxnGeneMat), expected); end function getPhylDistNeedsData(testCase)