feat: update GEM generating master script

- Use automated function for updating animal GEM
- Automated updating of tsv format annotation files

code/masterScriptFruitflyGEM.m

@@ -7,77 +7,44 @@
 %          pathways/reactions.
 %
 %
+%          Annotation files in tsv format are also generated by
+%          integrating human and fruitfly-specific annotation information.
+%
 
 
-%% Load Human-GEM as template
-load('Human-GEM.mat');
-
-
-% 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);
-
+%% Add path
+addpath(genpath('../../Fruitfly-GEM/'));
 
 
-%% Generate Fruitfly-GEM by using Human-GEM as template
+%% Prepare Fruitfly ortholog pairs and species-specific network
 
-% get ortholog pairs from human to worm
+% get ortholog pairs from human to fruitfly
 fruitflyOrthologPairs = extractAllianceGenomeOrthologs('human2FruitflyOrthologs.json');
-fruitflyGEM = getModelFromOrthology(ihuman, fruitflyOrthologPairs);
-
 
-
-%% Incorporate species-specific reactions
-
-% get metabolic networks based on the KEGG annoation using RAVEN function
-KEGG_human=getKEGGModelForOrganism('hsa');
-KEGG_fruitfly=getKEGGModelForOrganism('dme');
-
-% remove reactions shared with human
-FruitflySpecificRxns=setdiff(KEGG_fruitfly.rxns,KEGG_human.rxns);
-
-% remove reactions included in Human-GEM
-FruitflySpecificRxns=setdiff(FruitflySpecificRxns,rxnAssoc.rxnKEGGID);
-
-
-% get species-specific network for manual inspection and then
-% organize species-specific pathways into two tsv files:
-fruitflySpecificNetwork=removeReactions(KEGG_fruitfly,...
-    setdiff(KEGG_fruitfly.rxns,FruitflySpecificRxns), true, true, true);
-
-% "fruitflySpecificMets.tsv" contains species-specific metabolites
+% load species-specific rxns and mets
+rxnsToAdd = importTsvFile('fruitflySpecificRxns.tsv');
 metsToAdd = importTsvFile('fruitflySpecificMets.tsv');
 
-% "fruitflySpecificRxns.tsv" contains species-specific reactions
-rxnsToAdd = importTsvFile('fruitflySpecificRxns.tsv');
-rxnsToAdd.subSystems = cellfun(@(s) {{s}}, rxnsToAdd.subSystems);
 
-% integrate fruitfly-specific metabolic network
-[fruitflyGEM, modelChanges] = addMetabolicNetwork(fruitflyGEM, rxnsToAdd, metsToAdd);
+%% Generate Fruitfly-GEM
+[fruitflyGEM, speciesSpecNetwork, gapfillNetwork]=updateAnimalGEM(...
+    fruitflyOrthologPairs,rxnsToAdd,metsToAdd,'Fruitfly-GEM');
+
 
+%% Update annotations
+[rxnAssoc, metAssoc] = updateAnimalAnnotations(fruitflyGEM,rxnsToAdd,metsToAdd);
 
-%% Gap-filling for biomass formation
-[fruitflyGEM, gapfillNetwork]=gapfill4EssentialTasks(fruitflyGEM,ihuman);
-% Added 39 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, fruitflyGEM.rxns) && isequal(metAssoc.mets, fruitflyGEM.mets)
+	fprintf('sanity check passed.\n');
+    exportTsvFile(rxnAssoc,'../model/reactions.tsv');
+    exportTsvFile(metAssoc,'../model/metabolites.tsv');
+end
 
-fruitflyGEM.id = 'Fruitfly-GEM';
 save('../model/Fruitfly-GEM.mat', 'fruitflyGEM');
-writeHumanYaml(fruitflyGEM, '../model/Fruitfly-GEM.mat');
+exportYaml(fruitflyGEM, '../model/Fruitfly-GEM.yml');
 exportModel(fruitflyGEM, '../model/Fruitfly-GEM.xml');