# -*- coding: utf-8 -*-
"""
Description:
Allows to visualize a metabolic network on a compounds perspectives
::
usage:
padmet visu_network -i=FILE -o=FILE [--html=FILE] [--level=STR] [--hide-currency]
options:
-h --help Show help.
-i=FILE pathname to the input file (either PADMet or SBML).
-o=FILE pathname to the output file (picture of metabolic network).
--html=FILE pathname to the output file (interactive hmtl of metabolic network).
--level=STR level of precision for the visualization (compound, reaction or pathway). By default visualization uses "compound".
--hide-currency hide currency metabolites.
"""
import docopt
import logging
import sys
try:
import igraph
except ImportError:
raise ImportError('Requires igraph with cairocffi, try:\npip install python-igraph cairocffi')
from cobra.io.sbml import read_sbml_model
from padmet.classes import PadmetRef, PadmetSpec
from padmet.utils.sbmlPlugin import convert_from_coded_id
logging.getLogger("cobra.io.sbml").setLevel(logging.CRITICAL)
[docs]
def command_help():
"""
Show help for analysis command.
"""
print(docopt.docopt(__doc__))
[docs]
def visu_network_cli(command_args):
args = docopt.docopt(__doc__, argv=command_args)
metabolic_network_file = args["-i"]
output_file = args["-o"]
html_output_file = args['--html']
visualization_level = args['--level']
hide_currency_metabolites = args["--hide-currency"]
if visualization_level is None:
visualization_level = "compound"
create_graph(metabolic_network_file, output_file, visualization_level, hide_currency_metabolites)
if html_output_file:
create_html_graph(metabolic_network_file, html_output_file, visualization_level, hide_currency_metabolites)
[docs]
def parse_compounds_padmet(padmet_file, hide_metabolites):
""" Parse padmets files to extract compounds to create edges and nodes for igraph.
Parameters
----------
padmet_file: str
pathname of the padmet file
hide_metabolites: list
list of metabolites to hide
Returns
-------
edges: list
edges between two compounds (symbolizing the reaction)
edges_label: list
for each edge the name of the reaction
weights: list
the weight associated to each edge
nodes: list
a compound
nodes_label: list
for each node the name of the compound
"""
padmetSpec = PadmetSpec(padmet_file)
all_rxns = [node for node in padmetSpec.dicOfNode.values() if node.type == "reaction"]
edges = []
edges_label = []
weights = []
nodes = {}
nodes_label = []
for rxn in all_rxns:
ins = []
outs = []
for rlt in padmetSpec.dicOfRelationIn[rxn.id]:
if rlt.type == "consumes":
ins.append(rlt.id_out)
if rlt.type == "produces":
outs.append(rlt.id_out)
for compound_in in ins:
if compound_in not in hide_metabolites:
for compound_out in outs:
if compound_out not in hide_metabolites:
if compound_in not in nodes:
new_cpd_id = len(nodes_label)
nodes_label.append(compound_in)
nodes[compound_in] = new_cpd_id
if compound_out not in nodes:
new_cpd_id = len(nodes_label)
nodes_label.append(compound_out)
nodes[compound_out] = new_cpd_id
edges.append((nodes[compound_in], nodes[compound_out]))
weights.append(1)
edges_label.append(rxn.id)
if 'REVERSIBLE' in rxn.misc['DIRECTION']:
edges.append((nodes[compound_out], nodes[compound_in]))
weights.append(1)
edges_label.append(rxn.id)
return edges, edges_label, weights, nodes, nodes_label
[docs]
def parse_compounds_sbml(sbml_file, hide_metabolites):
""" Parse sbml files to extract compounds to create edges and nodes for igraph.
Parameters
----------
sbml_file: str
pathname of the sbml file
hide_metabolites: list
list of metabolites to hide
Returns
-------
edges: list
edges between two compounds (symbolizing the reaction)
edges_label: list
for each edge the name of the reaction
weights: list
the weight associated to each edge
nodes: list
a compound
nodes_label: list
for each node the name of the compound
"""
sbml_model = read_sbml_model(sbml_file)
edges = []
edges_label = []
weights = []
nodes = {}
nodes_label = []
for reaction in sbml_model.reactions:
for reactant in reaction.reactants:
reactant = convert_from_coded_id(reactant.id)[0]
if reactant not in hide_metabolites:
for product in reaction.products:
product = convert_from_coded_id(product.id)[0]
if product not in hide_metabolites:
if reactant not in nodes:
new_cpd_id = len(nodes_label)
nodes_label.append(reactant)
nodes[reactant] = new_cpd_id
if product not in nodes:
new_cpd_id = len(nodes_label)
nodes_label.append(product)
nodes[product] = new_cpd_id
edges.append((nodes[reactant], nodes[product]))
weights.append(1)
edges_label.append(reaction.id)
if reaction.reversibility is True:
edges.append((nodes[product], nodes[reactant]))
weights.append(1)
edges_label.append(reaction.id)
return edges, edges_label, weights, nodes, nodes_label
[docs]
def parse_reactions_padmet(padmet_file):
""" Parse padmets files to extract reactions to create edges and nodes for igraph.
Parameters
----------
padmet_file: str
pathname of the padmet file
Returns
-------
edges: list
edges between two reactions
edges_label: list
for each edge the name of the reaction
weights: list
the weight associated to each edge
nodes: list
a compound
nodes_label: list
for each node the name of the compound
"""
padmetSpec = PadmetSpec(padmet_file)
edges = []
edges_label = []
weights = []
nodes = {}
nodes_label = []
all_rxns = [node for node in padmetSpec.dicOfNode.values() if node.type == "reaction"]
for rxn in all_rxns:
for rlt in padmetSpec.dicOfRelationIn[rxn.id]:
if rlt.type == "produces":
for sec_rlt in padmetSpec.dicOfRelationOut[rlt.id_out]:
if sec_rlt.type == "consumes":
sec_rxn_id = sec_rlt.id_in
if sec_rxn_id not in nodes:
new_cpd_id = len(nodes_label)
nodes_label.append(sec_rxn_id)
nodes[sec_rxn_id] = new_cpd_id
if rxn.id not in nodes:
new_cpd_id = len(nodes_label)
nodes_label.append(rxn.id)
nodes[rxn.id] = new_cpd_id
if (nodes[rxn.id], nodes[sec_rxn_id]) not in edges:
edges.append((nodes[rxn.id], nodes[sec_rxn_id]))
weights.append(1)
edges_label.append(rxn.id)
return edges, edges_label, weights, nodes, nodes_label
[docs]
def parse_pathways_padmet(padmet_file):
""" Parse padmets files to extract pathway inputs and ouputs to create edges and nodes for igraph.
Parameters
----------
padmet_file: str
pathname of the padmet file
Returns
-------
edges: list
edges between two compounds (symbolizing the pathway)
edges_label: list
for each edge the name of the pathway
weights: list
the weight associated to each edge
nodes: list
a compound
nodes_label: list
for each node the name of the compound
"""
padmetSpec = PadmetSpec(padmet_file)
# Check if the padmets and padmetref contain the INPUT-COMPOUNDS and OUTPUT-COMPOUNDS in pathway node.misc needed for this analysis.
padmet_input_compounds_in_pwys = [1 for node_pathway in padmetSpec.dicOfNode
if padmetSpec.dicOfNode[node_pathway].type == 'pathway' and 'INPUT-COMPOUNDS' in padmetSpec.dicOfNode[node_pathway].misc]
padmet_output_compounds_in_pwys = [1 for node_pathway in padmetSpec.dicOfNode
if padmetSpec.dicOfNode[node_pathway].type == 'pathway' and 'OUTPUT-COMPOUNDS' in padmetSpec.dicOfNode[node_pathway].misc]
if sum(padmet_input_compounds_in_pwys) == 0 or sum(padmet_output_compounds_in_pwys) == 0:
sys.exit("The padmet " + padmet_file + " does not contain INPUT-COMPOUNDS and OUTPUT-COMPOUNDS in the pathway node, padmet can't produce the pathway visualization without them.")
edges = []
edges_label = []
weights = []
nodes = {}
nodes_label = []
all_pwys = [node for node in padmetSpec.dicOfNode if padmetSpec.dicOfNode[node].type == "pathway"]
for pwy in all_pwys:
node_pwy = padmetSpec.dicOfNode[pwy]
if 'INPUT-COMPOUNDS' in node_pwy.misc and 'OUTPUT-COMPOUNDS' in node_pwy.misc:
reactants = node_pwy.misc['INPUT-COMPOUNDS'][0].split(',')
products = node_pwy.misc['OUTPUT-COMPOUNDS'][0].split(',')
for reactant in reactants:
if reactant not in nodes:
new_cpd_id = len(nodes_label)
nodes_label.append(reactant)
nodes[reactant] = new_cpd_id
for product in products:
if product not in nodes:
new_cpd_id = len(nodes_label)
nodes_label.append(product)
nodes[product] = new_cpd_id
reactant_product_tuples = [(reactant, product) for reactant in reactants for product in products]
for reactant, product in reactant_product_tuples:
edges.append((nodes[reactant], nodes[product]))
weights.append(1)
edges_label.append(node_pwy.id)
return edges, edges_label, weights, nodes, nodes_label
[docs]
def create_graph(metabolic_network_file, output_file, visualization_level, hide_currency_metabolites):
""" Using output of parse_compounds_padmet or parse_compounds_sbml create a network picture using igraph.
Parameters
----------
metabolic_network_file: str
pathname of the metabolic network file
output_file: str
pathname of the output picture of the metabolic network
visualization_level: str
level of visualization either compound, reaction or pathway
hide_currency_metabolites: bool
hide currency metabolites
"""
if hide_currency_metabolites:
hide_metabolites = ["PROTON", "WATER", "OXYGEN-MOLECULE", "NADP", "NADPH", "ATP",
"PPI", "CARBON-DIOXIDE", "Pi", "ADP", "CO-A", "UDP", "NAD",
"NADH", "AMP", "AMMONIA", "HYDROGEN-PEROXIDE", "Acceptor",
"Donor-H2", "3-5-ADP", "GDP", "CARBON-MONOXIDE", "GTP", "FAD"]
else:
hide_metabolites = []
if visualization_level == "compound":
if metabolic_network_file.endswith('.padmet'):
edges, edges_label, weights, nodes, nodes_label = parse_compounds_padmet(metabolic_network_file, hide_metabolites)
elif metabolic_network_file.endswith('.sbml'):
edges, edges_label, weights, nodes, nodes_label = parse_compounds_sbml(metabolic_network_file, hide_metabolites)
else:
sys.exit('No correct extension file as input. Input must be a .padmet or a .sbml file.')
elif visualization_level == "reaction":
if metabolic_network_file.endswith('.padmet'):
edges, edges_label, weights, nodes, nodes_label = parse_reactions_padmet(metabolic_network_file)
else:
sys.exit('No correct extension file as input for reaction level. Input must be a .padmet.')
elif visualization_level == "pathway":
if metabolic_network_file.endswith('.padmet'):
edges, edges_label, weights, nodes, nodes_label = parse_pathways_padmet(metabolic_network_file)
else:
sys.exit('No correct extension file as input for pathway level. Input must be a .padmet.')
n_vertices = len(nodes)
# igraph network implementation is from https://gist.github.com/Vini2/d13b12b37b01b4001cabf38b1f850d8a#file-visualise_graph_demo-ipynb
# Thanks to Vini2.
# Create graph
compounds_graph = igraph.Graph(directed=True)
# Add vertices
compounds_graph.add_vertices(n_vertices)
# Add edges to the graph
compounds_graph.add_edges(edges)
# Nodes label
compounds_graph.vs["label"] = nodes_label
# Add weights to edges in the graph
compounds_graph.es['weight'] = weights
visual_style = {}
# Define colors used for outdegree visualization
colours = ['#fecc5c', '#a31a1c']
# Set bbox and margin
visual_style["bbox"] = (5000,5000)
visual_style["margin"] = 17
# Set vertex colours
visual_style["vertex_color"] = 'grey'
# Set vertex size
visual_style["vertex_size"] = 20
# Set vertex lable size
visual_style["vertex_label_size"] = 8
# Don't curve the edges
visual_style["edge_curved"] = True
visual_style["edge_width"] = 2
visual_style["edge_arrow_size"] = 0.2
# Set the layout
my_layout = compounds_graph.layout_fruchterman_reingold()
visual_style["layout"] = my_layout
# Plot the graph
igraph.plot(compounds_graph, output_file, **visual_style)
[docs]
def create_html_graph(metabolic_network_file, output_file, visualization_level, hide_currency_metabolites):
""" Using output of parse_compounds_padmet or parse_compounds_sbml create an interactive graph in html.
Parameters
----------
metabolic_network_file: str
pathname of the metabolic network file
output_file: str
pathname of the output picture of the metabolic network
visualization_level: str
level of visualization either compound, reaction or pathway
hide_currency_metabolites: bool
hide currency metabolites
"""
try:
import pyvis
except ImportError:
raise ImportError('Requires pyvis, try:\npip install pyvis')
if hide_currency_metabolites:
hide_metabolites = ["PROTON", "WATER", "OXYGEN-MOLECULE", "NADP", "NADPH", "ATP",
"PPI", "CARBON-DIOXIDE", "Pi", "ADP", "CO-A", "UDP", "NAD",
"NADH", "AMP", "AMMONIA", "HYDROGEN-PEROXIDE", "Acceptor",
"Donor-H2", "3-5-ADP", "GDP", "CARBON-MONOXIDE", "GTP", "FAD"]
else:
hide_metabolites = []
if visualization_level == "compound":
if metabolic_network_file.endswith('.padmet'):
edges, edges_label, weights, nodes, nodes_label = parse_compounds_padmet(metabolic_network_file, hide_metabolites)
elif metabolic_network_file.endswith('.sbml'):
edges, edges_label, weights, nodes, nodes_label = parse_compounds_sbml(metabolic_network_file, hide_metabolites)
else:
sys.exit('No correct extension file as input. Input must be a .padmet or a .sbml file.')
elif visualization_level == "reaction":
if metabolic_network_file.endswith('.padmet'):
edges, edges_label, weights, nodes, nodes_label = parse_reactions_padmet(metabolic_network_file)
else:
sys.exit('No correct extension file as input for reaciton level. Input must be a .padmet.')
elif visualization_level == "pathway":
if metabolic_network_file.endswith('.padmet'):
edges, edges_label, weights, nodes, nodes_label = parse_pathways_padmet(metabolic_network_file)
else:
sys.exit('No correct extension file as input for pathway level. Input must be a .padmet.')
compounds_graph = pyvis.network.Network(height=1900, width=1900, directed=True)
compounds_graph.barnes_hut()
compounds_graph.add_nodes([node for node in nodes.values()], label=[nodes_label[node] for node in nodes.values()])
for edge in edges:
compounds_graph.add_edge(edge[0], edge[1])
compounds_graph.show(output_file)