cy3sbml

Mapping simulation & analysis data on SBML networks.
In this example we use the species concentration from a timecourse simulation to map them to the node size.


In [16]:
from __future__ import print_function, division
%matplotlib inline

Load model

Loading Elowitz2000 - Repressilator in tellurium & roadrunner. A synthetic oscillatory network of transcriptional regulators.
The model is available from Biomodels as BIOMD0000000012.


In [17]:
import tellurium as te
biomodel = "BIOMD0000000012"
r = te.loads("{}.xml".format(biomodel))

Simulate model

Simple timecourse simulation which we want to visualize with cy3sbml.


In [18]:
r.reset()
data = r.simulate(start=0, end=600, steps=300)
r.plot();


Convert to dataframe

We store the resulting data in a pandas dataframe. roadrunner returns a NamedArray, we import this into pandas.


In [19]:
import pandas as pd
from matplotlib import pylab as plt

# sbml ids from timeCourseSelections
print(r.timeCourseSelections)
sbml_ids = [name.replace("[", "").replace("]", "") 
            for name in r.timeCourseSelections]

print(sbml_ids)

# create dataframe
df = pd.DataFrame(data[:, 1:], index=data[:,0], columns=sbml_ids[1:])
df.head(10)


['time', '[PX]', '[PY]', '[PZ]', '[X]', '[Y]', '[Z]']
['time', 'PX', 'PY', 'PZ', 'X', 'Y', 'Z']
Out[19]:
PX PY PZ X Y Z
0.0 0.000000 0.000000 0.000000 0.000000 20.000000 0.000000
2.0 218.536962 358.025299 84.590543 21.233188 23.608512 5.866284
4.0 428.931337 536.654989 131.036149 15.558712 12.481129 3.295735
6.0 542.447157 585.378463 146.531911 11.032655 6.563502 1.902593
8.0 596.545539 572.310647 146.992668 8.518252 3.532942 1.198339
10.0 621.016984 532.647739 141.004095 7.402622 1.993866 0.869033
12.0 634.134454 483.779679 133.027789 7.161553 1.214552 0.747215
14.0 646.170121 433.970110 125.381864 7.436716 0.818488 0.748140
16.0 662.168816 386.889728 119.310951 7.970864 0.613381 0.831437
18.0 683.549126 343.940842 115.549430 8.554410 0.501915 0.980585

Connect to Cytoscape

We use py2cytoscape has wrapper modules for cyREST RESTful API. Hereby we can can access Cytoscape features in more Pythonic way instead of calling raw REST API via HTTP.

http://idekerlab.github.io/cyREST/#-1099067042


In [5]:
# !!! Start Cytoscape 3 with cyREST App

In [20]:
# start client
import json
import py2cytoscape
from py2cytoscape.data.cyrest_client import CyRestClient
from py2cytoscape.data.cynetwork import CyNetwork
from py2cytoscape.data.style import StyleUtil

cy = CyRestClient()

Load SBML in cy3sbml

All the interaction with cytoscape from now on is via simple REST API! py2cytoscape is only a simple wrapper to hide some of the REST calls.

The same can be done from anywhere (different programming language, ...).


In [21]:
# Reset
cy.session.delete()

In [22]:
# Step 2: Load SBML network
networks = cy.network.create_from('BIOMD0000000012.xml')
print(networks)


[<py2cytoscape.data.cynetwork.CyNetwork object at 0x7fe0920e1c50>, <py2cytoscape.data.cynetwork.CyNetwork object at 0x7fe09207ab90>, <py2cytoscape.data.cynetwork.CyNetwork object at 0x7fe09207ad10>]

Information about SBML networks


In [23]:
# Get information about SBML networks
base_net = None
for net in networks:
    # unique id of cytoscape objects
    net_suid = net.get_id()
    net_name = net.get_network_value(column='name') 
    print(net.get_network_value(column='name'))
    if (net_name.startswith('Base')):
        base_net = net

# Selected the Base network (reaction - species graph)
print('-'*60)
print (base_net)
print('-'*60)


All: Elowitz2000 - Repressilator
Base: Elowitz2000 - Repressilator
Kinetic: Elowitz2000 - Repressilator
------------------------------------------------------------
<py2cytoscape.data.cynetwork.CyNetwork object at 0x7fe09207ab90>
------------------------------------------------------------

Information about nodes


In [24]:
# get node and edges
nodes = base_net.get_nodes()
edges = base_net.get_edges()

print('* This network has ' + str(len(nodes)) + ' nodes and ' + str(len(edges)) + ' edges\n') 

# Get a row in the node table as pandas Series object
node0 = nodes[0]
row = base_net.get_node_value(id=node0)
print(row)

# Or, pick one cell in the table
cell = base_net.get_node_value(id=node0, column='name')
print('\nThis node has name: ' + str(cell))


* This network has 18 nodes and 18 edges

SUID                            7022
chebi                    CHEBI:33699
compartmentCode                    1
cyId                         _905769
derivedUnits               substance
hasOnlySubstanceUnits           True
kegg.compound                 C00046
label                      LacI mRNA
metaId                       _905769
name                       LacI mRNA
sbml compartment                cell
sbml id                            X
sbml initial amount                0
sbml type                    species
sbml type ext                species
sbo                      SBO:0000250
selected                       False
shared name                LacI mRNA
uniprot                       P03023
value                              0
dtype: object

This node has name: LacI mRNA

Node table

All table attributes are available and can be manipulated.


In [25]:
# cell = base_net.get_node_value(id=node0, column='name')
node_table = base_net.get_node_table()
node_table.head()


Out[25]:
shared name cyId sbml type label metaId unitSid kind scale multiplier sbo ... chebi kegg.compound reversible math kineticLaw variable compartmentCode sbml type ext name selected
SUID
7024 TetR mRNA _905781 species TetR mRNA _905781 NaN NaN NaN NaN SBO:0000250 ... CHEBI:33699 C00046 NaN NaN NaN NaN 1.0 species TetR mRNA False
7026 cI mRNA _905802 species cI mRNA _905802 NaN NaN NaN NaN SBO:0000250 ... CHEBI:33699 C00046 NaN NaN NaN NaN 1.0 species cI mRNA False
7060 translation of CI _905923 reaction translation of CI _905923 NaN NaN NaN NaN SBO:0000184 ... NaN NaN False NaN k_tl*Z NaN NaN reaction irreversible translation of CI False
7028 degradation of LacI transcripts _905823 reaction degradation of LacI transcripts _905823 NaN NaN NaN NaN SBO:0000179 ... NaN NaN False NaN kd_mRNA*X NaN NaN reaction irreversible degradation of LacI transcripts False
7095 transcription of TetR _906022 reaction transcription of TetR _906022 NaN NaN NaN NaN SBO:0000183 ... NaN NaN False NaN a0_tr+a_tr*KM^n/(KM^n+PX^n) NaN NaN reaction irreversible transcription of TetR False

5 rows × 30 columns

Find nodes for sbml ids

We now find the nodes which correspond to the species we simulated in the model. The suid (unique node identifiers) are retrieved and stored for lookup.


In [26]:
node_suids = []

# logical indexing to find the nodes
for sid in sbml_ids[1:]:
    suid_index = node_table[node_table["sbml id"]==sid].index
    node_suids.append(suid_index.get_values()[0])

suid2sid = dict(zip(node_suids, sbml_ids[1:]))
suid2sid


Out[26]:
{7016: 'PX', 7018: 'PY', 7020: 'PZ', 7022: 'X', 7024: 'Y', 7026: 'Z'}

Get CyNetworkView

CyNetworkView is a reference to a network view in your current Cytoscape session. This means CyNetworkView objects do not include actual data, such as node size or edge stroke colors. Instead, they hold a location of the data and create actual REST calls when you execute get/update methods for nodes and edges.


In [27]:
# Get views for a network: Cytoscape "may" have multiple views, and that's why it returns list instead of an object.
view_ids = base_net.get_views() 
print(view_ids)

# format option specify the return type
view1 = base_net.get_view(view_ids[0], format='view')

# This is a CyNetworkView object, not dictionary
print(view1)


[7334]
<py2cytoscape.data.network_view.CyNetworkView object at 0x7fe0b37f2110>

In [28]:
from py2cytoscape.data.util_network import NetworkUtil as util
import time

# DataFrame for node views
df_vs_node = pd.DataFrame(index=node_suids, 
                  columns=['NODE_WIDTH', 'NODE_HEIGHT'])

for timepoint in df.index:
    # print(timepoint)
    for index, row in df_vs_node.iterrows():
        sbml_id = suid2sid.get(index)
        
        # Set node size from simulation results
        row['NODE_WIDTH'] = df.loc[timepoint, sbml_id]/15
        row['NODE_HEIGHT'] = df.loc[timepoint, sbml_id]/15
    
        # normalisation with max value
        row['NODE_WIDTH'] = 80 * df.loc[timepoint, sbml_id]/max(df[sbml_id])
        row['NODE_HEIGHT'] = 80 * df.loc[timepoint, sbml_id]/max(df[sbml_id])
    

    # Apply for timepoint
    view1.batch_update_node_views(df_vs_node)
    
    # wait
    time.sleep(0.03)

In [23]:
# reset to original size
for index, row in df_vs_node.iterrows():
    sbml_id = suid2sid.get(index)
    # row['NODE_FILL_COLOR'] = '#FF0000'
    row['NODE_WIDTH'] = 40
    row['NODE_HEIGHT'] = 40
view1.batch_update_node_views(df_vs_node)

In [24]:
r.plot()


Out[24]:
<module 'matplotlib.pyplot' from '/usr/local/lib/python2.7/dist-packages/matplotlib/pyplot.pyc'>

Set node images

Create images for species

In a first step we want to display images in the network. For this the individual figures are created and store in the results folder


In [ ]:
# create results folder
import os
results_dir = "./results/{}".format(biomodel)
if not os.path.exists(results_dir):
    os.mkdir()

# Create images for all species
time = df.index
for sbml_id in df.columns:
    # Create plot
    plt.figure(figsize=[2,2])
    plt.plot(time, df[sbml_id], linewidth=2, color="black")
    plt.xlabel('time')
    plt.ylabel(sbml_id)
    plt.savefig('{}/{}.png'.format(results_dir, sbml_id))

In [ ]:
## Serve images
We run a simple file server to serve the images via URL.

In [ ]:
# !!! serve the images on webserver
# python -m SimpleHTTPServer
# http://localhost:8000/results/BIOMD0000000012/

In [ ]:
The images are than reachable under  
http://localhost:8000/results/BIOMD0000000012/

In [ ]:
# DataFrame for node views
df_vs_node = pd.DataFrame(index=node_suids, 
                  columns=['NODE_CUSTOM_PAINT_1'])

# apply all visual changes for the timepoint
for index, row in df_vs_node.iterrows():
    sbml_id = suid2sid.get(index)
    # row['NODE_FILL_COLOR'] = '#FF0000'
    row['NODE_CUSTOM_PAINT_1'] = 'http://localhost:8000/results/BIOMD0000000012/{}.png'.format(sbml_id)

# Apply for timepoint
view1.batch_update_node_views(df_vs_node)

In [19]:
# apply styles
# cy.style.apply(style=minimal_style, network=base_net)
# apply layout
# cy.layout.apply(name='force-directed', network=base_net)

In [ ]: