Clustering and Graphs from steady-states correlation matrix (#105)

* initial commit for graphs

* added documentation

* added networkx dependency

* fix conflicts1

dingo/illustrations.py

@@ -11,6 +11,8 @@
 import plotly.io as pio
 import plotly.express as px
 from dingo.utils import compute_copula
+import plotly.figure_factory as ff
+from scipy.cluster import hierarchy
 
 def plot_copula(data_flux1, data_flux2, n = 5, width = 900 , height = 600, export_format = "svg"):
     """A Python function to plot the copula between two fluxes
@@ -129,3 +131,77 @@ def plot_corr_matrix(corr_matrix, reactions, removed_reactions=[], format="svg")
 
     fig_name = "CorrelationMatrix." + format
     pio.write_image(fig, fig_name, scale=2)
+
+
+
+def plot_dendrogram(dissimilarity_matrix, reactions , plot_labels=False, t=2.0, linkage="ward"):
+    """A Python function to plot the dendrogram of a dissimilarity matrix.
+
+    Keyword arguments:
+    dissimilarity_matrix -- A matrix produced from the "cluster_corr_reactions" function
+    reactions -- A list with the model's reactions
+    plot_labels -- A boolean variable that if True plots the reactions labels in the dendrogram
+    t -- A threshold that defines a threshold that cuts the dendrogram 
+         at a specific height and colors the occuring clusters accordingly
+    linkage -- linkage defines the type of linkage. 
+               Available linkage types are: single, average, complete, ward.
+    """
+
+    fig = ff.create_dendrogram(dissimilarity_matrix,
+                               labels=reactions,
+                               linkagefun=lambda x: hierarchy.linkage(x, linkage),                           
+                               color_threshold=t)
+    fig.update_layout(width=800, height=800)
+
+    if plot_labels == False:
+        fig.update_layout(
+            xaxis=dict(
+            showticklabels=False,
+            ticks="") )
+    else:
+        fig.update_layout(
+        xaxis=dict(
+            title_font=dict(size=10),
+            tickfont=dict(size=8) ),
+        yaxis=dict(
+            title_font=dict(size=10),
+            tickfont=dict(size=8) ) )
+
+    fig.show()
+
+
+
+def plot_graph(G, pos):
+    """A Python function to plot a graph created from a correlation matrix.
+
+    Keyword arguments:
+    G -- A graph produced from the "graph_corr_matrix" function.
+    pos -- A layout for the corresponding graph.
+    """
+
+    fig = go.Figure()
+
+    for u, v, data in G.edges(data=True):
+        x0, y0 = pos[u]
+        x1, y1 = pos[v]
+
+        edge_color = 'blue' if data['weight'] > 0 else 'red'
+
+        fig.add_trace(go.Scatter(x=[x0, x1], y=[y0, y1], mode='lines', 
+                                    line=dict(width=abs(data['weight']) * 1, 
+                                    color=edge_color), hoverinfo='none',
+                                    showlegend=False))
+
+    for node in G.nodes():
+        x, y = pos[node]
+        node_name = G.nodes[node].get('name', f'Node {node}')
+
+        fig.add_trace(go.Scatter(x=[x], y=[y], mode='markers', 
+                                    marker=dict(size=10),
+                                    text=[node_name],
+                                    textposition='top center',
+                                    name = node_name,
+                                    showlegend=False))
+
+    fig.update_layout(width=800, height=800)
+    fig.show()

dingo/utils.py

@@ -11,6 +11,9 @@
 from scipy.sparse import diags
 from dingo.scaling import gmscale
 from dingo.nullspace import nullspace_dense, nullspace_sparse
+from scipy.cluster import hierarchy
+from networkx.algorithms.components import connected_components
+import networkx as nx
 
 def compute_copula(flux1, flux2, n):
     """A Python function to estimate the copula between two fluxes
@@ -340,4 +343,96 @@ def correlated_reactions(steady_states, reactions=[], pearson_cutoff = 0.90, ind
         else:
             np.fill_diagonal(filtered_corr_matrix, 1)
             return filtered_corr_matrix, indicator_dict
-
+
+
+
+def cluster_corr_reactions(correlation_matrix, reactions, linkage="ward", 
+                           t = 4.0, correction=True):
+    """A Python function for hierarchical clustering of the correlation matrix
+
+    Keyword arguments:
+    correlation_matrix -- A numpy 2D array of a correlation matrix
+    reactions -- A list with the model's reactions
+    linkage -- linkage defines the type of linkage. 
+               Available linkage types are: single, average, complete, ward.
+    t -- A threshold that defines a threshold that cuts the dendrogram 
+         at a specific height and produces clusters
+    correction -- A boolean variable that if True converts the values of the
+                  the correlation matrix to absolute values.
+    """
+
+    # function to return a nested list with grouped reactions based on clustering
+    def clusters_list(reactions, labels):
+        clusters = []
+        unique_labels = np.unique(labels)
+        for label in unique_labels:
+            cluster = []
+            label_where = np.where(labels == label)[0]
+            for where in label_where:
+                cluster.append(reactions[where])
+            clusters.append(cluster)
+        return clusters
+
+    if correction == True:        
+        dissimilarity_matrix = 1 - abs(correlation_matrix)
+    else:
+        dissimilarity_matrix = 1 - correlation_matrix
+
+    Z = hierarchy.linkage(dissimilarity_matrix, linkage)
+    labels = hierarchy.fcluster(Z, t, criterion='distance')
+
+    clusters = clusters_list(reactions, labels)
+    return dissimilarity_matrix, labels, clusters
+
+
+
+def graph_corr_matrix(correlation_matrix, reactions, correction=True, 
+                      clusters=[], subgraph_nodes = 5):
+    """A Python function that creates the main graph and its subgraphs
+       from a correlation matrix.
+
+    Keyword arguments:
+    correlation_matrix -- A numpy 2D array of a correlation matrix.
+    reactions -- A list with the model's reactions.
+    correction -- A boolean variable that if True converts the values of the
+                  the correlation matrix to absolute values.
+    clusters -- A nested list with clustered reactions created from the "" function.
+    subgraph_nodes -- A variable that specifies a cutoff for a graph's nodes. 
+                      It filters subgraphs with low number of nodes..
+    """
+
+    graph_matrix = correlation_matrix.copy()
+    np.fill_diagonal(graph_matrix, 0)
+
+    if correction == True:
+        graph_matrix = abs(graph_matrix)
+
+    G = nx.from_numpy_array(graph_matrix)
+    G = nx.relabel_nodes(G, lambda x: reactions[x])
+
+    pos = nx.spring_layout(G)
+    unconnected_nodes = list(nx.isolates(G))
+    G.remove_nodes_from(unconnected_nodes)
+    G_nodes = G.nodes()
+
+    graph_list = []
+    layout_list = []
+
+    graph_list.append(G)
+    layout_list.append(pos)
+
+    subgraphs = [G.subgraph(c) for c in connected_components(G)]
+    H_nodes_list = []
+
+    for i in range(len(subgraphs)):
+        if len(subgraphs[i].nodes()) > subgraph_nodes and len(subgraphs[i].nodes()) != len(G_nodes):           
+            H = G.subgraph(subgraphs[i].nodes())
+            for cluster in clusters:
+                if H.has_node(cluster[0]) and H.nodes() not in H_nodes_list:
+                    H_nodes_list.append(H.nodes())
+
+                    pos = nx.spring_layout(H)
+                    graph_list.append(H)
+                    layout_list.append(pos)
+
+    return graph_list, layout_list

pyproject.toml

@@ -18,6 +18,7 @@ cobra = "^0.26.0"
 plotly = "^5.11.0"
 kaleido = "0.2.1"
 pyoptinterface = {version = "^0.2.7", extras = ["highs"]}
+networkx = "3.1"
 
 [tool.poetry.dev-dependencies]