GEM-PRO - List of Gene IDs

This notebook gives an example of how to run the GEM-PRO pipeline with a list of gene IDs.

Input: List of gene IDs
Output: GEM-PRO model


In [1]:
import sys
import logging
In [2]:
# Import the GEM-PRO class
from ssbio.pipeline.gempro import GEMPRO
In [3]:
# Printing multiple outputs per cell
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"


Set the logging level in logger.setLevel(logging.<LEVEL_HERE>) to specify how verbose you want the pipeline to be. Debug is most verbose.

    • Only really important messages shown
    • Major errors
    • Warnings that don’t affect running of the pipeline
  • INFO (default)
    • Info such as the number of structures mapped per gene
    • Really detailed information that will print out a lot of stuff
Warning: DEBUG mode prints out a large amount of information, especially if you have a lot of genes. This may stall your notebook!
In [4]:
# Create logger
logger = logging.getLogger()
logger.setLevel(logging.INFO)  # SET YOUR LOGGING LEVEL HERE #
In [5]:
# Other logger stuff for Jupyter notebooks
handler = logging.StreamHandler(sys.stderr)
formatter = logging.Formatter('[%(asctime)s] [%(name)s] %(levelname)s: %(message)s', datefmt="%Y-%m-%d %H:%M")
logger.handlers = [handler]

Initialization of the project

Set these three things:

    • The directory where a folder named after your PROJECT will be created
    • Your project name
    • Your list of gene IDs

A directory will be created in ROOT_DIR with your PROJECT name. The folders are organized like so:

    ├── data  # General storage for pipeline outputs
    ├── model  # SBML and GEM-PRO models are stored here
    ├── genes  # Per gene information
    │   ├── <gene_id1>  # Specific gene directory
    │   │   └── protein
    │   │       ├── sequences  # Protein sequence files, alignments, etc.
    │   │       └── structures  # Protein structure files, calculations, etc.
    │   └── <gene_id2>
    │       └── protein
    │           ├── sequences
    │           └── structures
    ├── reactions  # Per reaction information
    │   └── <reaction_id1>  # Specific reaction directory
    │       └── complex
    │           └── structures  # Protein complex files
    └── metabolites  # Per metabolite information
        └── <metabolite_id1>  # Specific metabolite directory
            └── chemical
                └── structures  # Metabolite 2D and 3D structure files
Note: Methods for protein complexes and metabolites are still in development.
In [6]:
import tempfile
ROOT_DIR = tempfile.gettempdir()

PROJECT = 'genes_GP'
LIST_OF_GENES = ['b0761', 'b0889', 'b0995', 'b1013', 'b1014', 'b1040', 'b1130', 'b1187', 'b1221', 'b1299']
PDB_FILE_TYPE = 'mmtf'
In [7]:
# Create the GEM-PRO project
my_gempro = GEMPRO(gem_name=PROJECT, root_dir=ROOT_DIR, genes_list=LIST_OF_GENES, pdb_file_type=PDB_FILE_TYPE)
[2017-09-24 03:21] [ssbio.pipeline.gempro] INFO: Creating GEM-PRO project directory in folder /tmp
[2017-09-24 03:21] [ssbio.pipeline.gempro] INFO: /tmp/genes_GP: GEM-PRO project location
[2017-09-24 03:21] [ssbio.pipeline.gempro] INFO: 10: number of genes

Mapping gene ID –> sequence

First, we need to map these IDs to their protein sequences. There are 2 ID mapping services provided to do this - through KEGG or UniProt. The end goal is to map a UniProt ID to each ID, since there is a comprehensive mapping (and some useful APIs) between UniProt and the PDB.

Note: You only need to map gene IDs using one service. However you can run both if some genes don’t map in one service and do map in another!


In [8]:
# KEGG mapping of gene ids
print('Missing KEGG mapping: ', my_gempro.missing_kegg_mapping)
[2017-09-24 03:21] [ssbio.pipeline.gempro] INFO: 10/10: number of genes mapped to KEGG
[2017-09-24 03:21] [ssbio.pipeline.gempro] INFO: Completed ID mapping --> KEGG. See the "df_kegg_metadata" attribute for a summary dataframe.

Missing KEGG mapping:  []
kegg refseq uniprot pdbs sequence_file metadata_file
b0761 eco:b0761 NP_415282 P0A9G8 1B9M;1H9S;1B9N;1O7L;1H9R eco-b0761.faa eco-b0761.kegg
b0889 eco:b0889 NP_415409 P0ACJ0 2GQQ;2L4A eco-b0889.faa eco-b0889.kegg
b0995 eco:b0995 NP_415515 P38684 1ZGZ eco-b0995.faa eco-b0995.kegg
b1013 eco:b1013 NP_415533 P0ACU2 4JYK;4XK4;4X1E;3LOC eco-b1013.faa eco-b1013.kegg
b1014 eco:b1014 NP_415534 P09546 3E2Q;4JNZ;3E2R;4JNY;2GPE;4O8A;3E2S;2FZN;1TJ1;1... eco-b1014.faa eco-b1014.kegg
In [9]:
# UniProt mapping
print('Missing UniProt mapping: ', my_gempro.missing_uniprot_mapping)
[2017-09-24 03:21] [root] INFO: getUserAgent: Begin
[2017-09-24 03:21] [root] INFO: getUserAgent: user_agent: EBI-Sample-Client/ (; Python 3.5.3; Linux) Python-requests/2.14.2
[2017-09-24 03:21] [root] INFO: getUserAgent: End
[2017-09-24 03:21] [ssbio.pipeline.gempro] INFO: 10/10: number of genes mapped to UniProt
[2017-09-24 03:21] [ssbio.pipeline.gempro] INFO: Completed ID mapping --> UniProt. See the "df_uniprot_metadata" attribute for a summary dataframe.

Missing UniProt mapping:  []
uniprot reviewed gene_name kegg refseq pdbs pfam description entry_date entry_version seq_date seq_version sequence_file metadata_file
b0761 P0A9G8 False modE ecj:JW0744;eco:b0761 NP_415282.1;WP_001147439.1 1B9M;1B9N;1H9R;1H9S;1O7L PF00126;PF03459 Transcriptional regulator ModE 2017-07-05 104 2005-07-19 1 P0A9G8.fasta P0A9G8.xml
b0889 P0ACJ0 False lrp ecj:JW0872;eco:b0889 NP_415409.1;WP_000228473.1 2GQQ;2L4A PF01037 Leucine-responsive regulatory protein 2017-08-30 105 2007-01-23 2 P0ACJ0.fasta P0ACJ0.xml
b0995 P38684 False torR ecj:JW0980;eco:b0995 NP_415515.1;WP_001120125.1 1ZGZ PF00072;PF00486 TorCAD operon transcriptional regulatory prote... 2017-08-30 147 1997-11-01 2 P38684.fasta P38684.xml
b1013 P0ACU2 False rutR ecj:JW0998;eco:b1013 NP_415533.1;WP_000191701.1 3LOC;4JYK;4X1E;4XK4 PF00440;PF08362 HTH-type transcriptional regulator RutR 2017-07-05 98 2005-11-22 1 P0ACU2.fasta P0ACU2.xml
b1014 P09546 False putA ecj:JW0999;eco:b1014 NP_415534.1;WP_001326840.1 1TIW;1TJ0;1TJ1;1TJ2;2AY0;2FZM;2FZN;2GPE;2RBF;3... PF00171;PF01619;PF14850 Bifunctional protein PutA 2017-08-30 178 1997-11-01 3 P09546.fasta P09546.xml
In [10]:
# Set representative sequences
print('Missing a representative sequence: ', my_gempro.missing_representative_sequence)
[2017-09-24 03:21] [ssbio.pipeline.gempro] INFO: 10/10: number of genes with a representative sequence
[2017-09-24 03:21] [ssbio.pipeline.gempro] INFO: See the "df_representative_sequences" attribute for a summary dataframe.

Missing a representative sequence:  []
uniprot kegg pdbs sequence_file metadata_file
b0761 P0A9G8 ecj:JW0744;eco:b0761 1B9M;1B9N;1H9R;1H9S;1O7L P0A9G8.fasta P0A9G8.xml
b0889 P0ACJ0 ecj:JW0872;eco:b0889 2GQQ;2L4A P0ACJ0.fasta P0ACJ0.xml
b0995 P38684 ecj:JW0980;eco:b0995 1ZGZ P38684.fasta P38684.xml
b1013 P0ACU2 ecj:JW0998;eco:b1013 3LOC;4JYK;4X1E;4XK4 P0ACU2.fasta P0ACU2.xml
b1014 P09546 ecj:JW0999;eco:b1014 1TIW;1TJ0;1TJ1;1TJ2;2AY0;2FZM;2FZN;2GPE;2RBF;3... P09546.fasta P09546.xml

Mapping representative sequence –> structure

These are the ways to map sequence to structure:

  1. Use the UniProt ID and their automatic mappings to the PDB
  2. BLAST the sequence to the PDB
  3. Make homology models or
  4. Map to existing homology models

You can only utilize option #1 to map to PDBs if there is a mapped UniProt ID set in the representative sequence. If not, you’ll have to BLAST your sequence to the PDB or make a homology model. You can also run both for maximum coverage.


In [11]:
# Mapping using the PDBe best_structures service
[2017-09-24 03:21] [ssbio.pipeline.gempro] INFO: Mapping UniProt IDs --> PDB IDs...
[2017-09-24 03:21] [root] INFO: getUserAgent: Begin
[2017-09-24 03:21] [root] INFO: getUserAgent: user_agent: EBI-Sample-Client/ (; Python 3.5.3; Linux) Python-requests/2.14.2
[2017-09-24 03:21] [root] INFO: getUserAgent: End
[2017-09-24 03:21] [ssbio.pipeline.gempro] INFO: 8/10: number of genes with at least one experimental structure
[2017-09-24 03:21] [ssbio.pipeline.gempro] INFO: Completed UniProt --> best PDB mapping. See the "df_pdb_ranking" attribute for a summary dataframe.

pdb_id pdb_chain_id uniprot experimental_method resolution coverage start end unp_start unp_end rank
b0761 1b9m B P0A9G8 X-ray diffraction 1.75 1.000 4 265 1 262 2
b0761 1b9m A P0A9G8 X-ray diffraction 1.75 1.000 4 265 1 262 1
b0761 1b9n A P0A9G8 X-ray diffraction 2.09 1.000 4 265 1 262 3
b0761 1h9r A P0A9G8 X-ray diffraction 1.90 0.534 1 140 123 262 11
b0761 1h9r B P0A9G8 X-ray diffraction 1.90 0.534 1 140 123 262 12
In [12]:
# Mapping using BLAST
my_gempro.blast_seqs_to_pdb(all_genes=True, seq_ident_cutoff=.9, evalue=0.00001)
[2017-09-24 03:21] [ssbio.pipeline.gempro] INFO: Completed sequence --> PDB BLAST. See the "df_pdb_blast" attribute for a summary dataframe.
[2017-09-24 03:21] [ssbio.pipeline.gempro] INFO: 1: number of genes with additional structures added from BLAST

pdb_id pdb_chain_id hit_score hit_evalue hit_percent_similar hit_percent_ident hit_num_ident hit_num_similar
b1013 4x1e A 966.0 1.382510e-104 0.910377 0.910377 193 193
b1013 4x1e B 966.0 1.382510e-104 0.910377 0.910377 193 193

Downloading and ranking structures


Warning: Downloading all PDBs takes a while, since they are also parsed for metadata. You can skip this step and just set representative structures below if you want to minimize the number of PDBs downloaded.
In [13]:
# Download all mapped PDBs and gather the metadata
[2017-09-24 03:21] [ssbio.pipeline.gempro] INFO: Updated PDB metadata dataframe. See the "df_pdb_metadata" attribute for a summary dataframe.
[2017-09-24 03:21] [ssbio.pipeline.gempro] INFO: Saved 40 structures total

chemicals description experimental_method mapped_chains pdb_id pdb_title resolution structure_file taxonomy_name
b0761 CA;CL;MOO TRANSCRIPTIONAL REGULATOR MODE X-ray diffraction A;B;C;D 1o7l Molybdate-activated form of ModE from Escheric... 2.75 1o7l.mmtf ESCHERICHIA COLI
b0761 MOO MOLYBDENUM TRANSPORT PROTEIN MODE X-ray diffraction A;B 1h9s Molybdate bound complex of Dimop domain of Mod... 1.82 1h9s.mmtf ESCHERICHIA COLI;ESCHERICHIA COLI
In [14]:
# Set representative structures
[2017-09-24 03:21] [ssbio.core.protein] WARNING: b1014: no structures meet quality checks
[2017-09-24 03:21] [ssbio.core.protein] WARNING: b0995: no structures meet quality checks
[2017-09-24 03:21] [ssbio.core.protein] WARNING: b1130: no structures meet quality checks
[2017-09-24 03:21] [ssbio.pipeline.gempro] INFO: 5/10: number of genes with a representative structure
[2017-09-24 03:21] [ssbio.pipeline.gempro] INFO: See the "df_representative_structures" attribute for a summary dataframe.

id is_experimental file_type structure_file
b0761 REP-1b9m True pdb 1b9m-A_clean.pdb
b0889 REP-2gqq True pdb 2gqq-A_clean.pdb
b1013 REP-4jyk True pdb 4jyk-A_clean.pdb
b1187 REP-1hw1 True pdb 1hw1-A_clean.pdb
b1221 REP-1a04 True pdb 1a04-A_clean.pdb
In [16]:
# Looking at the information saved within a gene
<StructProp REP-1hw1 at 0x7f03bf44e390>
{'_structure_dir': '/tmp/genes_GP/genes/b1187/b1187_protein/structures',
 'chains': [<ChainProp A at 0x7f03bee1f828>],
 'date': None,
 'file_type': 'pdb',
 'id': 'REP-1hw1',
 'is_experimental': True,
 'mapped_chains': ['A'],
 'notes': {},
 'original_pdb_id': '1hw1',
 'resolution': 1.5,
 'structure_file': '1hw1-A_clean.pdb',
 'taxonomy_name': 'Escherichia coli'}

Saving your GEM-PRO

Warning: Saving is still experimental. For a full GEM-PRO with sequences & structures, depending on the number of genes, saving can take >5 minutes.

In [17]:
import os.path as op
my_gempro.save_json(op.join(my_gempro.model_dir, '{}.json'.format(, compression=False)
[2017-09-24 03:21] [root] WARNING: json-tricks: numpy scalar serialization is experimental and may work differently in future versions
[2017-09-24 03:21] [] INFO: Saved <class 'ssbio.pipeline.gempro.GEMPRO'> (id: genes_GP) to /tmp/genes_GP/model/genes_GP.json