Skip to content

Commit

Permalink
feat: update GEM generating master script
Browse files Browse the repository at this point in the history
- Use automated function for updating animal GEM
- Automated updating of tsv format annotation files
  • Loading branch information
haowang-bioinfo committed Jun 23, 2021
1 parent c423f78 commit 43bb829
Showing 1 changed file with 21 additions and 53 deletions.
74 changes: 21 additions & 53 deletions code/masterScriptMouseGEM.m
Original file line number Diff line number Diff line change
Expand Up @@ -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');

0 comments on commit 43bb829

Please sign in to comment.