Human 1.13

README.md

@@ -67,8 +67,15 @@ Detailed instructions on the installation and use of the Human-GEM model and rep
 
 ### Installation Instructions
 * Clone the [main branch](https://github.com/SysBioChalmers/Human-GEM/tree/main) of this repository, or [download the latest release](https://github.com/SysBioChalmers/Human-GEM/releases/latest).
-* Add the directory to your MATLAB path (instructions [here](https://se.mathworks.com/help/matlab/ref/addpath.html?requestedDomain=www.mathworks.com)).
-
+* Add the directory to your MATLAB path either by using the lines below or manually (instructions [here](https://se.mathworks.com/help/matlab/ref/addpath.html?requestedDomain=www.mathworks.com)).
+```matlab
+% Replace "/my/path/" with the actual path to the Human-GEM folder
+cd /my/path/Human-GEM/code
+% This will add the relevant paths to the path variable in MATLAB
+HumanGEMInstaller.install
+% It is also possible to remove Human-GEM from the MATLAB path using
+HumanGEMInstaller.uninstall
+```
 
 ## Model Files
 
@@ -149,4 +156,4 @@ A collection of manually curated 2D metabolic maps associated with Human-GEM are
 <!-- markdownlint-restore -->
 <!-- prettier-ignore-end -->
 
-<!-- ALL-CONTRIBUTORS-LIST:END -->
+<!-- ALL-CONTRIBUTORS-LIST:END -->

code/HumanGEMInstaller.m

@@ -0,0 +1,44 @@
+classdef HumanGEMInstaller
+% HumanGEMInstaller
+%   Support for installing and uninstalling
+%   Run HumanGEMInstaller.install to install (will set up the path in MATLAB)
+%   Run HumanGEMInstaller.uninstall to clear the path from MATLAB
+%   To install, you first need to cd to the HumanGEM root.
+
+    methods (Static)
+        function install
+            sourceDir = fileparts(which(mfilename));
+            paths = HumanGEMInstaller.GetFilteredSubPaths(sourceDir);
+            addpath(paths);
+            savepath;
+        end
+        function uninstall
+            sourceDir = fileparts(which(mfilename));
+            paths = HumanGEMInstaller.GetFilteredSubPaths(sourceDir);
+            rmpath(paths);
+            savepath;
+        end
+        function path = getHumanGEMMainPath()
+			path = fileparts(which(mfilename));
+			path = strrep(path, '\', '/'); %get rid of backslashes in Windows
+			if ~endsWith(path, '/')
+				path = strcat(path,'/');
+			end
+		end
+
+        function newPaths = GetFilteredSubPaths(path_)
+            pathSep = pathsep();
+			%Check that there are no separators in the path - that will cause 
+            %problems since the separator is used to separate paths in a string
+			if contains(path_, pathSep)
+				error('The path in which Human-GEM resides may not contain path separator chars such as semicolon for this installation to work!');
+			end
+            subpaths = arrayfun(@(subf) sprintf('%s', genpath([path_ filesep subf{:}])), HumanGEMInstaller.SUBFOLDERS, 'UniformOutput', false);
+            newPaths = strjoin(subpaths, '');
+        end
+    end
+
+    properties (Constant)
+      SUBFOLDERS = {'../code', '../data', '../model'};
+   end
+end

code/README.md

@@ -1,6 +1,10 @@
 # Human-GEM code
 
-This directory contains functions and scripts that facilitate work with the Human-GEM model and other content in the repository. This directory is organized as follows:
+A handy installer script is provided via _HumanGEMInstaller.m_, which can add and remove, via _HumanGEMInstaller.install_ and 
+_HumanGEMInstaller.uninstall_, all necessary subfolders of Human-GEM to the Matlab path.
+
+This directory contains functions and scripts that facilitate work with the Human-GEM model and other content in the repository. 
+This directory is organized as follows:
 
 ```
 code

code/io/importYaml.m

@@ -303,6 +303,17 @@
 [~, metIdx] = ismember(model.mets, newMets);
 model.S = S(metIdx, :);
 
+% remove empty fields
+modelFields = fieldnames(model);
+emptyFields = modelFields(structfun(@isempty, model));
+keepFields = {'description', 'version'};  % some fields can be empty
+removeFields = setdiff(emptyFields, keepFields);
+model = rmfield(model, removeFields);
+if ~silentMode && (numel(removeFields) > 0)
+    fprintf('\nThe following empty fields were removed from the model:\n');
+    fprintf('\t%s\n', removeFields{:});
+end
+
 if ~silentMode
     fprintf(' Done!\n');
 end

code/modelCuration/README.md

@@ -2,4 +2,4 @@
 
 This directory contains curation-related scripts and functions used to make changes to the Human-GEM repository. These curation scripts help to improve transparency of changes made to the model when the number of changes is too large to view practically.
 
-
+- `getCompFromUniprotCellAtlas.py`: Code for collecting subcellular localization information for existing metabolic enzymes from **Swissprot** and **Cell Atlas** ([HPA](https://www.proteinatlas.org/search/has_protein_data_in%3ACell)), and then incorporating the compartment info from both sources into `genes.tsv`.

code/modelCuration/addRxnACOD1_20221102.m

@@ -0,0 +1,33 @@
+% load model and new reaction info
+ihuman = importYaml('../../model/Human-GEM.yml');
+rxnsToAdd = importTsvFile('../../data/modelCuration/addRxnACOD1_20221102.tsv');
+
+% add new genes to Human-GEM
+newGEM = ihuman;
+newGEM.genes = [newGEM.genes; rxnsToAdd.grRules];
+newGEM.rxnGeneMat(:, end+1) = 0;
+
+% reformat subsystem 
+if ~iscell(rxnsToAdd.subSystems{1})
+    rxnsToAdd.subSystems = cellfun(@(s) {{s}}, rxnsToAdd.subSystems);
+end
+% add new reaction
+newGEM = addRxns(newGEM, rxnsToAdd, 3);
+
+% add reaction annotation
+rxnAssoc = importTsvFile('../../model/reactions.tsv');
+annoNames = fieldnames(rxnAssoc);
+for i=1:length(annoNames)
+    if ismember(annoNames{i}, fieldnames(rxnsToAdd))
+        rxnAssoc.(annoNames{i}) = [rxnAssoc.(annoNames{i}); rxnsToAdd.(annoNames{i})];
+    elseif isequal(annoNames{i},'spontaneous')
+        rxnAssoc.(annoNames{i}) = [rxnAssoc.(annoNames{i}); 0];
+    else
+        rxnAssoc.(annoNames{i}) = [rxnAssoc.(annoNames{i}); {''}];
+    end
+end
+
+% update yaml model and reaction association file
+exportYaml(newGEM, '../../model/Human-GEM.yml');
+exportTsvFile(rxnAssoc,'../../model/reactions.tsv');
+

code/modelCuration/getCompFromDeepLoc2.m

@@ -0,0 +1,60 @@
+% This function is to extract the compartment prediction from DeepLoc2 tool
+% to expand the gene.tsv file.
+% Only when there is not such info from cellAtlas or SwissPort, we use the
+% compartment prediction from the DeepLoc2.
+% DeepLoc2 is used for the localization prediction, the protein fasta file
+% for all metabolic proteins in the model was used as the input file, the
+% default parameter was used in the prediction, output file is stored as
+% DeepLoc2_compartment.tsv in ../../data/modelCuration
+
+% Read the predicted compartment info
+fileName = '../../data/modelCuration/DeepLoc2_compartment.csv';
+fid  = fopen(fileName);
+DeeplocData = textscan(fid,[repmat('%s ',[1,12]) '%s'],'Delimiter',',');
+DeeplocData = [DeeplocData{1:end}];
+DeeplocData = DeeplocData(2:end,:);
+DeeplocData(:,[1,14]) = split(DeeplocData(:,1),';');
+fclose(fid);
+
+% Read gene.tsv
+fileName = '../../model/genes.tsv';
+fid  = fopen(fileName);
+geneData = textscan(fid,[repmat('%s ',[1,9]) '%s'],'Delimiter','\t');
+geneData = [geneData{1:end}];
+fclose(fid);
+geneData = strrep(geneData,'"','');
+
+% index the DeeplocDdata by geneData order
+DeeplocDdata_sorted(1:length(geneData(:,1)),1:length(DeeplocData(1,:))) = {''};
+[~,idx] = ismember(geneData(:,1),DeeplocData(:,1));
+DeeplocDdata_sorted(idx~=0,:) = DeeplocData(idx(idx~=0),:);
+
+
+% build a dict for key word mapping of compartment assignment
+DeepLoc_keywords = {'Cytoplasm','Cytosol';
+'Nucleus','Nucleus';
+'Extracellular','Extracellular';
+'Cell membrane','Cell membrane';
+'Mitochondrion','Mitochondria';
+'Plastid','';
+'Endoplasmic reticulum','Endoplasmic reticulum';
+'Lysosome/Vacuole','Lysosome';
+'Golgi apparatus','Golgi apparatus';
+'Peroxisome','Peroxisome'
+'|',';'};
+for i = 1:length(DeepLoc_keywords(:,1))
+    DeeplocDdata_sorted(:,2) = strrep(DeeplocDdata_sorted(:,2),DeepLoc_keywords{i,1},DeepLoc_keywords{i,2});
+end
+% fill the compartment for protein without localization info from CellAtlas
+%and SwissPort
+idx = cellfun(@isempty,geneData(:,9)) & ~(cellfun(@isempty,DeeplocDdata_sorted(:,1)));
+geneData(idx,9) = DeeplocDdata_sorted(idx,2);
+geneData(idx,10) = {'DeepLoc2'};
+
+% write out gene.tsv
+out = geneData';
+fID = fopen('../../model/genes.tsv','w');
+fprintf(fID,[repmat('%s\t',[1,9]) '%s\n'],out{1:10});
+fprintf(fID,[repmat('"%s"\t',[1,9]) '"%s"\n'],out{11:end});
+fclose(fID);
+

code/modelCuration/getCompFromUniprotCellAtlas.py

@@ -0,0 +1,141 @@
+# collect subcellular localization for existing enzymes from Uniprot/Swissprot and Cell Atlas
+# compartment info from both sources were combined and integrated into genes.tsv
+
+import pandas as pd
+import csv
+
+# import subcellular location from SwissProt annotation
+SwissProt_tsv = pd.read_table("~/Downloads/SwissProt_20221115.tsv")
+SwissProt_proteins = SwissProt_tsv['Entry'].to_list()
+SwissProt_subCellularlocation = SwissProt_tsv['Subcellular location [CC]'].to_list()
+
+
+# build a dict for key word mapping of compartment assignment
+swissprot_keywords = {
+    #'Extracellular': ['',''],   # e
+    'Peroxisome': ['peroxisome'],   # x
+    'Mitochondria': ['mitochondrion {','mitochondrion matrix'],   # m
+    'Cytosol': ['cytoplasm','cytosol','cytoskeleton'],   # c
+    'Lysosome': ['lysosome','lysosome lumen','lysosome membrane'],   # l
+    'Endoplasmic reticulum': ['endoplasmic reticulum','Endoplasmic reticulum lumen','endoplasmic reticulum membrane'],   # r
+    'Golgi apparatus': ['golgi apparatus','golgi apparatus lumen','golgi apparatus membrane'],   # g
+    'Nucleus': ['nucleus','nucleolus'],   # n
+    'Inner mitochondria': ['mitochondrion inner membrane','mitochondrion intermembrane space']    # i
+}
+
+
+# store swissport compartment info as a dict with a key of uniprot id
+comps_from_swissprot = {}
+for i, text in enumerate(SwissProt_subCellularlocation):
+    out_list = []  # output compartmetns as a list
+    comps_from_swissprot[SwissProt_proteins[i]] = out_list
+    if not pd.isna(text):
+        for key in swissprot_keywords:
+            # check if cellular location text contains any keyword from a compartment
+            if any(word in text.lower() for word in swissprot_keywords[key]):
+                out_list.append(key)
+        comps_from_swissprot[SwissProt_proteins[i]] = out_list
+
+
+# load Human-GEM gene annotation file
+HumanGenes_tsv = pd.read_table("../../model/genes.tsv")
+Human_genes = HumanGenes_tsv['genes'].to_list()
+Human_proteins = HumanGenes_tsv['geneUniProtID'].to_list()
+
+
+# save SwissProt compartment information to tsv file
+swissprot_compList = ['']*len(Human_genes)
+for i in range(len(Human_genes)):
+    if not pd.isna(Human_proteins[i]):
+        if Human_proteins[i] in SwissProt_proteins:
+            swissprot_compList[i] = ';'.join(comps_from_swissprot[Human_proteins[i]])
+swissprot_compartments = pd.DataFrame()
+swissprot_compartments['genes'] = Human_genes
+swissprot_compartments['geneUniProtID'] = Human_proteins
+swissprot_compartments['compartments'] = swissprot_compList
+swissprot_compartments.to_csv('../../data/modelCuration/Swissprot_compartments.tsv', sep="\t", index=False)
+
+
+# import Cell Atlas compartment info
+cell_atlas_compartments = pd.read_table("../../data/modelCuration/CellAtlasCompartments_science_2017.tsv")
+cellAtlas_comps = cell_atlas_compartments['Subcellular location'].to_list()
+cellAtlas_ensembl_id = cell_atlas_compartments['Ensembl'].to_list()
+cellAtlas_uniprot_id = cell_atlas_compartments['Uniprot'].to_list()
+
+# investigate key words
+compartment_names = []
+for comps in cellAtlas_comps:
+    compartment_names.extend(comps.split(','))
+unique_compartments = set(compartment_names)
+unique_compartments
+# Cell Atlas compartments were sorted and the following ones were discarded:
+# 'Microtubules', 'Midbody ring', 'Midbody', 'Cytokinetic bridge', 'Microtubule ends', 'Mitotic spindle'
+# 'Intermediate filaments', 'Actin filaments', 'Focal adhesion sites', 'Cleavage furrow', 
+# 'Lipid droplets', 'Vesicles', 'Endosomes', 'Plasma membrane', 'Cell Junctions'
+# 'Centrosome', 'Centriolar satellite'
+
+# construct a dict for key word mapping to Cell Atlas compartments
+cellatlas_keywords = {
+    #'Extracellular': [''],   # e
+    'Peroxisome': ['Peroxisomes'],   # x
+    'Mitochondria': ['Mitochondria'],   # m
+    'Cytosol': ['Cytosol', 'Cytoplasmic bodies', 'Rods & Rings', 'Aggresome'],   # c
+    'Lysosome': ['Lysosomes'],   # l
+    'Endoplasmic reticulum': ['Endoplasmic reticulum'],   # r
+    'Golgi apparatus': ['Golgi apparatus'],   # g
+    'Nucleus': ['Nuclear membrane', 'Nucleoli rim', 'Nucleoli', 'Nuclear bodies', 'Nucleoli fibrillar center', 'Nucleoplasm', 'Kinetochore', 'Mitotic chromosome', 'Nuclear speckles'],   # n
+    #'Inner mitochondria': ['']    # i, Cell Atlas does provide such location info
+}
+
+
+# store CellAtlas compartment info to a dict with keys as ensembl ids
+geneComps_from_cell_atlas = {}
+for gene in Human_genes:
+    out_list = []  # store compartments as a list
+    if gene not in cellAtlas_ensembl_id:
+        geneComps_from_cell_atlas[gene] = out_list
+    else:
+        gene_ind = cellAtlas_ensembl_id.index(gene)
+        for key in cellatlas_keywords:
+            # check if cellular location text contains any keyword from a compartment
+            if any(word in cellAtlas_comps[gene_ind] for word in cellatlas_keywords[key]):
+                out_list.append(key)
+        geneComps_from_cell_atlas[gene] = out_list
+
+
+# integrate compartment info from cellAtlas and SwissProt with following rules:
+# 1. output two columns: "compartments" and "compDataSource"
+# 2. use one source if another has no assigned compartments
+# 3. union compartments if both source are provided
+swissprot_comps = pd.read_table("../../data/modelCuration/Swissprot_compartments.tsv")
+geneComps_from_swissprot = swissprot_comps['compartments'].to_list()
+
+# combine compartment info from two data sources
+geneComps_combined = []
+source = ['']*len(Human_genes)
+for i, gene in enumerate(Human_genes):
+    out_list = []  # save compartmetns as a list
+    if pd.isna(geneComps_from_swissprot[i]) and geneComps_from_cell_atlas[gene] != []:
+        out_list = geneComps_from_cell_atlas[gene]
+        source[i] = 'CellAtlas'
+    elif not pd.isna(geneComps_from_swissprot[i]) and geneComps_from_cell_atlas[gene] == []:
+        out_list = geneComps_from_swissprot[i].split(';')
+        source[i] = 'SwissProt'
+    elif not pd.isna(geneComps_from_swissprot[i]) and geneComps_from_cell_atlas[gene] != []:
+        union = set(geneComps_from_swissprot[i].split(';')+ geneComps_from_cell_atlas[gene])
+        out_list = list(union)
+        source[i] = 'SwissProt;CellAtlas'
+        if 'Mitochondria' not in geneComps_from_swissprot[i].split(';') and 'Inner mitochondria' in geneComps_from_swissprot[i].split(';'):
+            if 'Mitochondria' in geneComps_from_cell_atlas[gene]:
+                out_list.remove('Mitochondria')
+        # manual inspection of integration
+        print(gene+': '+geneComps_from_swissprot[i]+'\t'+';'.join(geneComps_from_cell_atlas[gene])+'\t'+';'.join(out_list))
+    geneComps_combined.append(';'.join(out_list))
+
+
+# append columns to genes.tsv
+HumanGenes_tsv['compartments'] = geneComps_combined
+HumanGenes_tsv['compDataSource'] = source
+HumanGenes_tsv[:0].to_csv('../../model/genes.tsv', sep="\t", index=False)
+HumanGenes_tsv.to_csv('../../model/genes.tsv', sep="\t", index=False, quoting=csv.QUOTE_NONNUMERIC, header=False, mode="a")
+