From 43bb829c940b78a702f8d1e956edeaeea0fa6a29 Mon Sep 17 00:00:00 2001
From: Hao-Chalmers <hao.wang@chalmers.se>
Date: Wed, 23 Jun 2021 10:38:29 +0200
Subject: [PATCH] feat: update GEM generating master script

- Use automated function for updating animal GEM
- Automated updating of tsv format annotation files
---
 code/masterScriptMouseGEM.m | 74 +++++++++++--------------------------
 1 file changed, 21 insertions(+), 53 deletions(-)

diff --git a/code/masterScriptMouseGEM.m b/code/masterScriptMouseGEM.m
index 1a2aa93..b52ec17 100644
--- a/code/masterScriptMouseGEM.m
+++ b/code/masterScriptMouseGEM.m
@@ -6,76 +6,44 @@
 %          the Human-GEM as template and taking into account mouse-specific
 %          pathways/reactions.
 %
+%          Annotation files in tsv format are also generated by
+%          integrating human and mouse-specific annotation information.
 %
 
 
-%% Load Human-GEM as template
-load('Human-GEM.mat');
+%% Add path
+addpath(genpath('../../Mouse-GEM/'));
 
 
-% convert gene identifiers from Ensembl ids to gene symbols
-[grRules,genes,rxnGeneMat] = translateGrRules(ihuman.grRules,'Name','ENSG');
-ihuman.grRules    = grRules;
-ihuman.genes      = genes;
-ihuman.rxnGeneMat = rxnGeneMat;
-
-
-
-%% Use MA reactions identifiers 
-
-% load reaction annotaiton files
-rxnAssoc = jsondecode(fileread('humanGEMRxnAssoc.JSON'));
-
-%replace reaction identifiers with MA ids if available
-ind = getNonEmptyList(rxnAssoc.rxnMAID);
-ihuman.rxns(ind) = rxnAssoc.rxnMAID(ind);
-
-
-
-%% Generate Mouse-GEM using Human-GEM as template
+%% Prepare Mouse ortholog pairs and species-specific network
 
 % get ortholog pairs from human to mouse
 mouseOrthologPairs = extractAllianceGenomeOrthologs('human2MouseOrthologs.json');
-mouseGEM = getModelFromOrthology(ihuman, mouseOrthologPairs);
-
 
-
-%% Incorporate mouse-specific reactions
-
-% get metabolic networks based on the KEGG annoation using RAVEN function
-KEGG_human=getKEGGModelForOrganism('hsa');
-KEGG_mouse=getKEGGModelForOrganism('mmu');
-
-% remove reactions shared with human
-MouseSpecificRxns=setdiff(KEGG_mouse.rxns,KEGG_human.rxns);
-
-% remove reactions included in Human-GEM
-MouseSpecificRxns=setdiff(MouseSpecificRxns,rxnAssoc.rxnKEGGID);
-
-% get species-specific network for manual inspection and then
-% organize species-specific pathways into two tsv files:
-mouseSpecificNetwork=removeReactions(KEGG_mouse, setdiff(KEGG_mouse.rxns,MouseSpecificRxns), true, true, true);
-
-% "mouseSpecificMets.tsv" contains new metabolites
+% load species-specific rxns and mets
+rxnsToAdd = importTsvFile('mouseSpecificRxns.tsv');
 metsToAdd = importTsvFile('mouseSpecificMets.tsv');
 
-% "mouseSpecificRxns.tsv" contains new reactions
-rxnsToAdd = importTsvFile('mouseSpecificRxns.tsv');
-rxnsToAdd.subSystems = cellfun(@(s) {{s}}, rxnsToAdd.subSystems);
 
-% integrate mouse-specific metabolic network
-[mouseGEM, modelChanges] = addMetabolicNetwork(mouseGEM, rxnsToAdd, metsToAdd);
+%% Generate Mouse-GEM
+[mouseGEM, speciesSpecNetwork, gapfillNetwork]=updateAnimalGEM(...
+    mouseOrthologPairs,rxnsToAdd,metsToAdd,'Mouse-GEM');
+
 
+%% Update annotation files
+[rxnAssoc, metAssoc] = updateAnimalAnnotations(mouseGEM,rxnsToAdd,metsToAdd);
 
-%% Gap-filling for biomass formation
-[mouseGEM, gapfillNetwork]=gapfill4EssentialTasks(mouseGEM,ihuman);
-% Added 0 reactions for gap-filling
 
+%% Save annotation into tsv files, and the model into mat, yml, and xml
 
-%% Save the model into mat, yml, and xml
+% sanity check
+if isequal(rxnAssoc.rxns, mouseGEM.rxns) && isequal(metAssoc.mets, mouseGEM.mets)
+	fprintf('sanity check passed.\n');
+    exportTsvFile(rxnAssoc,'../model/reactions.tsv');
+    exportTsvFile(metAssoc,'../model/metabolites.tsv');
+end
 
-mouseGEM.id = 'Mouse-GEM';
 save('../model/Mouse-GEM.mat', 'mouseGEM');
-writeHumanYaml(mouseGEM, '../model/Mouse-GEM.yml');
+writeYaml(mouseGEM, '../model/Mouse-GEM.yml');
 exportModel(mouseGEM, '../model/Mouse-GEM.xml');