Merge pull request #6 from SysBioChalmers/develop

Rat 1.2

code/masterScriptRatGEM.m

@@ -5,79 +5,45 @@
 %          the Human-GEM as template and taking into account rat-specific
 %          pathways/reactions.
 %
+%          Annotation files in tsv format are also generated by
+%          integrating human and rat-specific annotation information.
 %
 
 
-%% Load Human-GEM as template
-load('Human-GEM.mat');
+%% Add path
+addpath(genpath('../../Rat-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 Rat-GEM using Human-GEM as template
+%% Prepare Rat ortholog pairs and species-specific network
 
 % get ortholog pairs from human to rat
 ratOrthologPairs = extractAllianceGenomeOrthologs('human2RatOrthologs.json');
-ratGEM = getModelFromOrthology(ihuman, ratOrthologPairs);
-
 
-
-%% Incorporate rat-specific reactions
-
-% get metabolic networks based on the KEGG annoation using function getKEGGModelForOrganism 
-KEGG_human=getKEGGModelForOrganism('hsa');
-KEGG_rat=getKEGGModelForOrganism('rno');
-
-% remove shared human reactions
-RatSpecificRxns=setdiff(KEGG_rat.rxns,KEGG_human.rxns);
-
-% remove reactions included in Human-GEM
-RatSpecificRxns=setdiff(RatSpecificRxns,rxnAssoc.rxnKEGGID);
-
-% get species-specific network for manual inspection
-ratSpecificNetwork=removeReactions(KEGG_rat, setdiff(KEGG_rat.rxns,RatSpecificRxns), true, true, true);
-
-% organize species-specific pathways into two tsv files:
-% "ratSpecificMets.tsv" contains new metabolites
-metsToAdd = importTsvFile('ratSpecificMets.tsv');
-
-% "ratSpecificRxns.tsv" contains new reactions
+% load species-specific rxns and mets
 rxnsToAdd = importTsvFile('ratSpecificRxns.tsv');
-rxnsToAdd.subSystems = cellfun(@(s) {{s}}, rxnsToAdd.subSystems);
+metsToAdd = importTsvFile('ratSpecificMets.tsv');
 
-% integrate rat-specific metabolic network
-[ratGEM, modelChanges] = addMetabolicNetwork(ratGEM, rxnsToAdd, metsToAdd);
 
+%% Generate Rat-GEM
+[ratGEM, speciesSpecNetwork, gapfillNetwork]=updateAnimalGEM(...
+    ratOrthologPairs,rxnsToAdd,metsToAdd,'Rat-GEM');
 
 
-%% Gap-filling for biomass formation
-[ratGEM, gapfillNetwork]=gapfill4EssentialTasks(ratGEM,ihuman);
-% Added 0 reactions for gap-filling
+%% Update annotations
+[rxnAssoc, metAssoc] = updateAnimalAnnotations(ratGEM,rxnsToAdd,metsToAdd);
 
 
+%% Save annotation into tsv files, and the model into mat, yml, and xml
 
-%% Save the generated model into Matlab, Yaml and SBML formats
+% sanity check
+if isequal(rxnAssoc.rxns, ratGEM.rxns) && isequal(metAssoc.mets, ratGEM.mets)
+	fprintf('sanity check passed.\n');
+    exportTsvFile(rxnAssoc,'../model/reactions.tsv');
+    exportTsvFile(metAssoc,'../model/metabolites.tsv');
+end
 
-ratGEM.id = 'Rat-GEM';
 save('../model/Rat-GEM.mat', 'ratGEM');
-writeHumanYaml(ratGEM, '../model/Rat-GEM.yml');
+exportYaml(ratGEM, '../model/Rat-GEM.yml');
 exportModel(ratGEM, '../model/Rat-GEM.xml');