import numpy as np

import json
import pandas as pd
import networkx as nx
import [Link] as plt
import seaborn as sns
from [Link] import pearsonr
from mpl_toolkits.mplot3d import Axes3D
from [Link] import Normalize
import [Link] as cm
from collections import defaultdict
import random
from [Link] import PCA
import gseapy as gp
import os

class TriRNSC:
"""
Implementation of the TriRNSC algorithm for triclustering gene expression data.

TriRNSC (Triclustering using Restricted Neighborhood Search Clustering) is

a graph-based approach that considers genes, conditions, and time points to
discover meaningful triclusters in gene expression profiles.
"""

def __init__(self, threshold=0.7, tabu_length_ratio=0.02,

naive_stopping_tolerance=15,
scaled_stopping_tolerance=15, diversification_frequency=50,
diversification_length=10, max_experiments=30):
"""
Initialize the TriRNSC algorithm with the given parameters.

Parameters:
-----------
threshold : float
Correlation threshold for constructing the gene co-expression network.
tabu_length_ratio : float
Ratio of the number of vertices to determine the tabu list length.
naive_stopping_tolerance : int
Number of moves without improvement before switching to scaled cost.
scaled_stopping_tolerance : int
Number of moves without improvement before terminating an experiment.
diversification_frequency : int
Period for shuffling diversification.
diversification_length : int
Number of random moves in a diversification phase.
max_experiments : int
Maximum number of experiments to run.
"""

[Link] = threshold
self.tabu_length_ratio = tabu_length_ratio
self.naive_stopping_tolerance = naive_stopping_tolerance
self.scaled_stopping_tolerance = scaled_stopping_tolerance
self.diversification_frequency = diversification_frequency
self.diversification_length = diversification_length
self.max_experiments = max_experiments
 [Link] = None
self.gene_names = None
self.condition_names = None
self.time_points = None
[Link] = None
self.adjacency_list = None
self.best_triclusters = None
self.correlation_matrix = None

def load_data(self, data_matrix, gene_names=None, condition_names=None,

time_points=None):
"""
Load the 3D gene expression data.

Parameters:
-----------
data_matrix : [Link]
3D array with dimensions (genes, conditions, time_points).
gene_names : list, optional
List of gene names.
condition_names : list, optional
List of condition names.
time_points : list, optional
List of time points.
"""

[Link] = data_matrix

n_genes, n_conditions, n_time_points = data_matrix.shape

self.gene_names = gene_names if gene_names is not None else [f"Gene_{i}"

for i in range(n_genes)]
self.condition_names = condition_names if condition_names is not None else
[f"Condition_{i}" for i in range(n_conditions)]
self.time_points = time_points if time_points is not None else [f"Time_{i}"
for i in range(n_time_points)]

print(f"Loaded data with {n_genes} genes, {n_conditions} conditions, and

{n_time_points} time points")

def construct_gene_co_expression_network(self):
"""
Construct a gene co-expression network (GCN) based on Pearson correlation.

Returns:
--------
[Link]
The gene co-expression network.
"""

print("Constructing gene co-expression network...")

n_genes = [Link][0]

G = [Link]()
G.add_nodes_from(range(n_genes))

correlation_matrix = [Link]((n_genes, n_genes))

edges_added = []
 for i in range(n_genes):
for j in range(i + 1, n_genes):

gene_i_profile = [Link][i].reshape(-1)
gene_j_profile = [Link][j].reshape(-1)

corr, _ = pearsonr(gene_i_profile, gene_j_profile)

correlation_matrix[i, j] = corr
correlation_matrix[j, i] = corr

if abs(corr) >= [Link]:

G.add_edge(i, j, weight=abs(corr))
edges_added.append((i, j, abs(corr)))

[Link] = G
self.correlation_matrix = correlation_matrix
self.adjacency_list = {node: list([Link](node)) for node in [Link]()}

[Link]['expression_data'] = [Link]

print(f"GCN constructed with {G.number_of_nodes()} nodes and

{G.number_of_edges()} edges.")

if G.number_of_edges() == 0:
print(f"WARNING: No edges in the network. The correlation threshold
({[Link]}) may be too high.")
print("Try lowering the threshold or check your data.")

return G, correlation_matrix, edges_added

def initialize_random_clusters(self, n_clusters=None):

"""
Initialize random clusters from the graph.

Parameters:
-----------
n_clusters : int, optional
Number of clusters to initialize. If None, a random number is chosen.

Returns:
--------
dict
A dictionary mapping cluster IDs to sets of nodes.
"""

nodes = list([Link]())

if n_clusters is None:
n_clusters = [Link](2, int([Link](len(nodes))))

clusters = {i: set() for i in range(n_clusters)}

for node in nodes:

cluster_id = [Link](0, n_clusters - 1)
clusters[cluster_id].add(node)

clusters = {k: v for k, v in [Link]() if v}

 return clusters

def naive_cost(self, clusters):

"""
Calculate the naive cost function for a clustering.

Parameters:
-----------
clusters : dict
A dictionary mapping cluster IDs to sets of nodes.

Returns:
--------
float
The naive cost value.
"""

cost = 0

for node in [Link]():

cluster_id = self._get_cluster_id(node, clusters)
if cluster_id is None:
continue

cluster = clusters[cluster_id]

cross_edges = sum(1 for neighbor in self.adjacency_list[node] if

neighbor not in cluster)

missing_edges = sum(1 for other_node in cluster if other_node != node

and other_node not in self.adjacency_list[node])

cost += cross_edges + missing_edges

return cost / 2

def scaled_cost(self, clusters):

"""
Calculate the scaled cost function for a clustering.

Parameters:
-----------
clusters : dict
A dictionary mapping cluster IDs to sets of nodes.

Returns:
--------
float
The scaled cost value.
"""

n = len([Link]())
cost = 0

for node in [Link]():

cluster_id = self._get_cluster_id(node, clusters)
if cluster_id is None:
continue
 cluster = clusters[cluster_id]

cross_edges = sum(1 for neighbor in self.adjacency_list[node] if

neighbor not in cluster)

missing_edges = sum(1 for other_node in cluster if other_node != node

and other_node not in self.adjacency_list[node])

neighborhood = set(self.adjacency_list[node]).union(cluster)
if not neighborhood:
continue

cost += ((n - 1) / 3) * (cross_edges + missing_edges) /

len(neighborhood)

return cost

def _get_cluster_id(self, node, clusters):

"""
Get the cluster ID for a node.

Parameters:
-----------
node : int
The node to find.
clusters : dict
A dictionary mapping cluster IDs to sets of nodes.

Returns:
--------
int or None
The cluster ID or None if the node is not in any cluster.
"""

for cluster_id, cluster_nodes in [Link]():

if node in cluster_nodes:
return cluster_id
return None

def _move_node(self, node, from_cluster_id, to_cluster_id, clusters):

"""
Move a node from one cluster to another.

Parameters:
-----------
node : int
The node to move.
from_cluster_id : int
The source cluster ID.
to_cluster_id : int
The destination cluster ID.
clusters : dict
A dictionary mapping cluster IDs to sets of nodes.

Returns:
--------
dict
The updated clusters.
 """

new_clusters = {k: set(v) for k, v in [Link]()}

new_clusters[from_cluster_id].remove(node)

if to_cluster_id not in new_clusters:

new_clusters[to_cluster_id] = set()

new_clusters[to_cluster_id].add(node)

new_clusters = {k: v for k, v in new_clusters.items() if v}

return new_clusters

def run_naive_phase(self, initial_clusters):

"""
Run the naive phase of the RNSC algorithm.

Parameters:
-----------
initial_clusters : dict
Initial clustering.

Returns:
--------
dict
The best clusters found during the naive phase.
"""

print("Running naive phase...")

current_clusters = initial_clusters
best_clusters = current_clusters
best_cost = self.naive_cost(current_clusters)

tabu_length = max(int(len([Link]()) * self.tabu_length_ratio), 1)

tabu_list = []

no_improvement_count = 0

iteration = 0

while no_improvement_count < self.naive_stopping_tolerance:

iteration += 1

if iteration % self.diversification_frequency == 0:
for _ in range(self.diversification_length):

node = [Link](list([Link]()))

from_cluster_id = self._get_cluster_id(node, current_clusters)

if from_cluster_id is None:
continue

if [Link]() < 0.8 and current_clusters:

to_cluster_id =
[Link](list(current_clusters.keys()))
while to_cluster_id == from_cluster_id and
len(current_clusters) > 1:
to_cluster_id =
[Link](list(current_clusters.keys()))
else:
to_cluster_id = max(current_clusters.keys()) + 1 if
current_clusters else 0

current_clusters = self._move_node(node, from_cluster_id,

to_cluster_id, current_clusters)

node = [Link](list([Link]()))

from_cluster_id = self._get_cluster_id(node, current_clusters)

if from_cluster_id is None:
continue

possible_clusters = list(current_clusters.keys())
possible_clusters = [cid for cid in possible_clusters if cid !=
from_cluster_id]
if not possible_clusters:
continue # no other cluster to move to
to_cluster_id = [Link](possible_clusters)

if (node, to_cluster_id) in tabu_list:

continue

new_clusters = self._move_node(node, from_cluster_id, to_cluster_id,

current_clusters)
new_cost = self.naive_cost(new_clusters)

if new_cost < best_cost:

best_cost = new_cost
best_clusters = new_clusters
no_improvement_count = 0
else:
no_improvement_count += 1

current_clusters = new_clusters

tabu_list.append((node, to_cluster_id))
if len(tabu_list) > tabu_length:
tabu_list.pop(0)

if iteration % 1000 == 0:
print(f"Naive phase iteration {iteration}, current cost:
{new_cost:.4f}, best cost: {best_cost:.4f}")

print(f"Naive phase completed after {iteration} iterations with best cost:

{best_cost:.4f}")
return best_clusters

def run_scaled_phase(self, initial_clusters):

"""
Run the scaled phase of the RNSC algorithm.

Parameters:
-----------
initial_clusters : dict
Initial clustering from the naive phase.
 Returns:
--------
dict
The best clusters found during the scaled phase.
"""

print("Running scaled phase...")

current_clusters = initial_clusters
best_clusters = current_clusters
best_cost = self.scaled_cost(current_clusters)

tabu_length = max(int(len([Link]()) * self.tabu_length_ratio), 1)

tabu_list = []

no_improvement_count = 0

iteration = 0

while no_improvement_count < self.scaled_stopping_tolerance:

iteration += 1

if iteration % self.diversification_frequency == 0:
for _ in range(self.diversification_length):

node = [Link](list([Link]()))

from_cluster_id = self._get_cluster_id(node, current_clusters)

if from_cluster_id is None:
continue

if [Link]() < 0.8 and current_clusters:

to_cluster_id =
[Link](list(current_clusters.keys()))
while to_cluster_id == from_cluster_id and
len(current_clusters) > 1:
to_cluster_id =
[Link](list(current_clusters.keys()))
else:
to_cluster_id = max(current_clusters.keys()) + 1 if
current_clusters else 0

current_clusters = self._move_node(node, from_cluster_id,

to_cluster_id, current_clusters)

node = [Link](list([Link]()))

from_cluster_id = self._get_cluster_id(node, current_clusters)

if from_cluster_id is None:
continue

possible_clusters = list(current_clusters.keys())
possible_clusters = [cid for cid in possible_clusters if cid !=
from_cluster_id]
if not possible_clusters:
continue # no other cluster to move to
to_cluster_id = [Link](possible_clusters)
 if (node, to_cluster_id) in tabu_list:
continue

new_clusters = self._move_node(node, from_cluster_id, to_cluster_id,

current_clusters)
new_cost = self.scaled_cost(new_clusters)

if new_cost < best_cost:

best_cost = new_cost
best_clusters = new_clusters
no_improvement_count = 0
else:
no_improvement_count += 1

current_clusters = new_clusters

tabu_list.append((node, to_cluster_id))
if len(tabu_list) > tabu_length:
tabu_list.pop(0)

if iteration % 1000 == 0:
print(f"Scaled phase iteration {iteration}, current cost:
{new_cost:.4f}, best cost: {best_cost:.4f}")

print(f"Scaled phase completed after {iteration} iterations with best cost:

{best_cost:.4f}")
return best_clusters

def run_rnsc(self):
"""
Run the RNSC algorithm to find triclusters.

Returns:
--------
dict
The best triclusters found.
"""

if [Link] is None:
self.construct_gene_co_expression_network()

best_overall_clusters = None
best_overall_cost = float('inf')

for exp in range(self.max_experiments):

print(f"\nStarting experiment {exp + 1}/{self.max_experiments}")

random_clusters = self.initialize_random_clusters(n_clusters=2)

naive_clusters = self.run_naive_phase(random_clusters)

scaled_clusters = self.run_scaled_phase(naive_clusters)

current_cost = self.scaled_cost(scaled_clusters)

print(f"Experiment {exp + 1} cost: {current_cost:.4f}")

if current_cost < best_overall_cost:

 best_overall_cost = current_cost
best_overall_clusters = scaled_clusters
print(f"New best overall cost: {best_overall_cost:.4f}")

self.best_triclusters = best_overall_clusters

print(f"\nRNSC completed with best overall cost: {best_overall_cost:.4f}")

print(f"Found {len(best_overall_clusters)} triclusters")

return self.best_triclusters

def get_tricluster_info(self):
"""
Get information about the found triclusters.

Returns:
--------
list of dict
List of triclusters with gene names.
"""

if self.best_triclusters is None:
print("No triclusters found. Run the RNSC algorithm first.")
return None

tricluster_info = []

for cluster_id, nodes in self.best_triclusters.items():

genes = [self.gene_names[node] for node in nodes]

submatrix = [Link][[node for node in nodes], :, :]

msr = self.calculate_msr(submatrix)

volume = len(nodes) * [Link][1] * [Link][2]

tricluster_info.append({
'id': cluster_id,
'genes': genes,
'gene_indices': list(nodes),
'size': len(nodes),
'msr': msr,
'volume': volume,
'tqi': msr / volume if volume > 0 else float('inf')
})

tricluster_info.sort(key=lambda x: x['size'], reverse=True)

return tricluster_info

def calculate_msr(self, submatrix):

"""
Calculate the mean square residue (MSR) for a tricluster.

Parameters:
-----------
submatrix : [Link]
3D submatrix of gene expression values for a tricluster.
 Returns:
--------
float
The MSR value.
"""

if [Link] == 0:
return float('inf')

n_genes, n_conditions, n_time_points = [Link]

mgct = [Link](submatrix)

mgc = [Link](submatrix, axis=0)

mgt = [Link](submatrix, axis=1)

mct = [Link](submatrix, axis=2)

mg = [Link](submatrix, axis=(1, 2))[:, [Link], [Link]]

mc = [Link](submatrix, axis=(0, 2))[[Link], :, [Link]]

mt = [Link](submatrix, axis=(0, 1))[[Link], [Link], :]

r = 0
for i in range(n_genes):
for j in range(n_conditions):
for k in range(n_time_points):
r_ijk = (submatrix[i, j, k] + mgct - mg[i, 0, 0] - mc[0, j, 0]
- mt[0, 0, k])**2
r += r_ijk

msr = r / (n_genes * n_conditions * n_time_points)

return msr

# === VISUALIZATION METHODS ===

def visualize_correlation_matrix(self, cluster_labels=None):

"""
Visualize the correlation matrix with optional cluster ordering.

Parameters:
-----------
cluster_labels : list or dict, optional
Cluster assignments for genes. Can be a list of labels or a dictionary
mapping node indices to cluster IDs.
"""

if self.correlation_matrix is None:
if [Link] is None:
self.construct_gene_co_expression_network()
else:
print("Correlation matrix not found. Recomputing...")
n_genes = [Link][0]
self.correlation_matrix = [Link]((n_genes, n_genes))
for i in range(n_genes):
for j in range(i+1, n_genes):
 gene_i_profile = [Link][i].reshape(-1)
gene_j_profile = [Link][j].reshape(-1)
corr, _ = pearsonr(gene_i_profile, gene_j_profile)
self.correlation_matrix[i, j] = corr
self.correlation_matrix[j, i] = corr

[Link](figsize=(10, 8))

if cluster_labels is not None and isinstance(cluster_labels, dict):

labels = [Link](len([Link]()), dtype=int)
for cluster_id, nodes in cluster_labels.items():
for node in nodes:
labels[node] = cluster_id
cluster_labels = labels

if cluster_labels is not None:

sorted_indices = [Link](cluster_labels)
reordered_matrix = self.correlation_matrix[sorted_indices, :][:,
sorted_indices]

boundaries = []
last_cluster = cluster_labels[sorted_indices[0]]
for i, idx in enumerate(sorted_indices):
if cluster_labels[idx] != last_cluster:
[Link](i - 0.5)
last_cluster = cluster_labels[idx]

[Link](reordered_matrix, cmap='coolwarm', vmin=-1, vmax=1)

for boundary in boundaries:

[Link](y=boundary, color='black', linestyle='-', linewidth=1)
[Link](x=boundary, color='black', linestyle='-', linewidth=1)

[Link]('Gene Correlation Matrix (Ordered by Cluster)')

else:
[Link](self.correlation_matrix, cmap='coolwarm', vmin=-1, vmax=1)
[Link]('Gene Correlation Matrix')

plt.tight_layout()
[Link]()

def visualize_gene_expression_patterns(self, cluster_labels=None,

condition_idx=0):
"""
Visualize gene expression patterns by cluster.

Parameters:
-----------
cluster_labels : list or dict, optional
Cluster assignments for genes. If None, uses best_triclusters.
condition_idx : int
Index of condition to display.
"""

if [Link] is None:
print("No data loaded. Load data first.")
return

if cluster_labels is None:
 if self.best_triclusters is None:
print("No clusters available. Run RNSC first or provide cluster
labels.")
return

labels = [Link]([Link][0], dtype=int)

for cluster_id, nodes in self.best_triclusters.items():
for node in nodes:
labels[node] = cluster_id
cluster_labels = labels

if isinstance(cluster_labels, dict):
labels = [Link]([Link][0], dtype=int)
for cluster_id, nodes in cluster_labels.items():
for node in nodes:
labels[node] = cluster_id
cluster_labels = labels

n_clusters = len([Link](cluster_labels))
time_points = range([Link][2])

[Link](figsize=(15, 4 * n_clusters))

for cluster_id in range(n_clusters):

ax = [Link](n_clusters, 1, cluster_id + 1)

genes_in_cluster = [Link](cluster_labels == cluster_id)[0]

for gene_idx in genes_in_cluster:

gene_expr = [Link][gene_idx, condition_idx, :]
[Link](time_points, gene_expr, alpha=0.3, linewidth=1)

if len(genes_in_cluster) > 0:
mean_expr = [Link]([Link][genes_in_cluster, condition_idx, :],
axis=0)
[Link](time_points, mean_expr, 'k-', linewidth=2.5,
label=f'Mean (n={len(genes_in_cluster)})')

[Link](f'Cluster {cluster_id} (size: {len(genes_in_cluster)})',

fontsize=14)
[Link]('Time Points', fontsize=12)
[Link]('Expression Value', fontsize=12)
[Link](True, alpha=0.3)
[Link](loc='best')

plt.tight_layout()
[Link](f'Gene Expression Patterns by Cluster (Condition:
{self.condition_names[condition_idx]})',
fontsize=16, y=1.02)
[Link]()

def visualize_gcn_2d(self, cluster_labels=None, pos=None, node_size_factor=300,

edge_alpha=0.6):
"""
Visualize the gene co-expression network in 2D.

Parameters:
-----------
cluster_labels : list or dict, optional
 Cluster assignments for genes. If None, uses best_triclusters.
pos : dict, optional
Node positions. If None, spring layout is used.
node_size_factor : float
Scaling factor for node sizes
edge_alpha : float
Transparency for edges
"""

if [Link] is None:
self.construct_gene_co_expression_network()

if len([Link]()) == 0:
print("The graph has no edges to visualize. Try lowering the
correlation threshold.")
return

[Link](figsize=(12, 10))

if pos is None:
pos = nx.spring_layout([Link], seed=42)

edge_weights = nx.get_edge_attributes([Link], 'weight')

edge_widths = [3 * weight for u, v, weight in
[Link](data='weight')]

node_sizes = [node_size_factor * (1 + [Link](n)) for n in

[Link]()]

if cluster_labels is None and self.best_triclusters is not None:

labels = [Link](len([Link]()), dtype=int)
for cluster_id, nodes in self.best_triclusters.items():
for node in nodes:
labels[node] = cluster_id
cluster_labels = labels

if isinstance(cluster_labels, dict):
labels = [Link](len([Link]()), dtype=int)
for cluster_id, nodes in cluster_labels.items():
for node in nodes:
labels[node] = cluster_id
cluster_labels = labels

if cluster_labels is not None:

nx.draw_networkx_nodes([Link], pos, node_color=cluster_labels,
node_size=node_sizes, alpha=0.8,
cmap=[Link].tab10)

unique_clusters = [Link](cluster_labels)
for cluster in unique_clusters:
[Link]([0], [0], 'o', color=[Link].tab10(cluster / max(9,
max(unique_clusters))),
label=f'Cluster {cluster}', markersize=10, alpha=0)
[Link](loc='upper right')
else:
node_degrees = dict([Link]())
node_colors = [node_degrees[n] for n in [Link]()]
nx.draw_networkx_nodes([Link], pos, node_color=node_colors,
node_size=node_sizes, alpha=0.8,
cmap=[Link])

nx.draw_networkx_edges([Link], pos, width=edge_widths,

alpha=edge_alpha, edge_color='gray')

if len([Link]()) <= 20:

nx.draw_networkx_labels([Link], pos, font_size=10)

[Link]('Gene Co-expression Network (2D)', fontsize=16)

[Link]('off')
plt.tight_layout()
[Link]()

def visualize_gcn_3d(self, cluster_labels=None, edge_threshold=None,

use_pca=False):
"""
Visualize the gene co-expression network in 3D.

Parameters:
-----------
cluster_labels : list or dict, optional
Cluster assignments for genes. If None, uses best_triclusters.
edge_threshold : float, optional
Minimum edge weight to display. If None, uses [Link].
use_pca : bool
Whether to use PCA for node positioning.
"""

if [Link] is None:
self.construct_gene_co_expression_network()

if len([Link]()) == 0:
print("The graph has no edges to visualize. Try lowering the
correlation threshold.")
return

if edge_threshold is None:
edge_threshold = [Link]

fig = [Link](figsize=(14, 12))

ax = fig.add_subplot(111, projection='3d')

if use_pca and [Link] is not None:

flattened_data = [Link]([Link][0], -1)

pca = PCA(n_components=3)
positions = pca.fit_transform(flattened_data)

pos = {node: positions[node] for node in [Link]()}

print(f"PCA explained variance: {pca.explained_variance_ratio_}")

else:
pos = nx.spring_layout([Link], dim=3, seed=42)

node_xyz = [Link]([pos[v] for v in [Link]()])

if cluster_labels is None and self.best_triclusters is not None:

labels = [Link](len([Link]()), dtype=int)
for cluster_id, nodes in self.best_triclusters.items():
 for node in nodes:
labels[node] = cluster_id
cluster_labels = labels

if isinstance(cluster_labels, dict):
labels = [Link](len([Link]()), dtype=int)
for cluster_id, nodes in cluster_labels.items():
for node in nodes:
labels[node] = cluster_id
cluster_labels = labels

if cluster_labels is not None:

node_colors = cluster_labels
color_map = [Link].tab10
color_norm = Normalize(vmin=min(cluster_labels),
vmax=max(cluster_labels))
else:
node_colors = [[Link](n) for n in [Link]()]
color_map = [Link]
if max(node_colors) > min(node_colors):
color_norm = Normalize(vmin=min(node_colors),
vmax=max(node_colors))
else:
color_norm = Normalize(vmin=0, vmax=1)

node_sizes = [100 * (1 + [Link](n))**0.5 for n in

[Link]()]

scatter = [Link](
node_xyz[:, 0], node_xyz[:, 1], node_xyz[:, 2],
c=node_colors, cmap=color_map, norm=color_norm,
s=node_sizes, alpha=0.8, edgecolors='w', linewidth=0.5
)

edge_weights = nx.get_edge_attributes([Link], 'weight')

for u, v in [Link]():
weight = edge_weights.get((u, v), 0)
if weight >= edge_threshold:
x = [Link]([pos[u][0], pos[v][0]])
y = [Link]([pos[u][1], pos[v][1]])
z = [Link]([pos[u][2], pos[v][2]])

edge_color = [Link](weight)
line_width = weight * 3

[Link](x, y, z, color=edge_color, alpha=0.6, linewidth=line_width)

cbar = [Link](scatter, ax=ax, shrink=0.7, pad=0.1)

if cluster_labels is not None:

cbar.set_label('Cluster ID')
else:
cbar.set_label('Node Degree')

ax.set_xlabel('X', fontweight='bold')
ax.set_ylabel('Y', fontweight='bold')
ax.set_zlabel('Z', fontweight='bold')
ax.set_title('Gene Co-expression Network (3D)', fontsize=16)
 if cluster_labels is not None:
unique_clusters = [Link](cluster_labels)
legend_elements = []
for cluster in unique_clusters:
color = color_map(color_norm(cluster))
legend_elements.append(
plt.Line2D([0], [0], marker='o', color='w',
markerfacecolor=color,
markersize=10, label=f'Cluster {cluster}')
)
[Link](handles=legend_elements, loc='upper right')

ax.view_init(elev=20, azim=30)

plt.tight_layout()
[Link]()

def perform_go_analysis(self, top_n=10, save_dir="go_results"):

if self.best_triclusters is None:
print("No triclusters found. Run the RNSC algorithm first.")
return

[Link](save_dir, exist_ok=True)

print("\nStarting GO Term Analysis...")

try:
import gseapy as gp
except ImportError:
print("ERROR: gseapy package is not installed. Install it using: pip
install gseapy")
return

successful_clusters = 0

for cluster_id, nodes in self.best_triclusters.items():

if len(nodes) < 3:
print(f"Skipping cluster {cluster_id} - too few genes
({len(nodes)})")
continue # Skip very small clusters

gene_list = [self.gene_names[node] for node in nodes]

print(f"Processing cluster {cluster_id} with {len(gene_list)} genes")

print(f"First few genes in cluster: {gene_list[:5]}")

try:
organisms = ['Human', 'Yeast', 'Mouse']
success = False

for organism in organisms:

try:
print(f"Trying GO enrichment with organism: {organism}")
enr = [Link](
gene_list=gene_list,
gene_sets=['GO_Biological_Process_2021'],
organism=organism,
outdir=None,
 no_plot=True,
cutoff=1.0 # More permissive p-value cutoff
)

if enr.res2d is not None and not [Link]:

success = True
print(f"GO enrichment successful with organism:
{organism}")
break
except Exception as e:
print(f"Failed with organism {organism}: {str(e)}")

if not success:
print(f"Could not perform GO enrichment for cluster
{cluster_id} with any organism")
with
open(f"{save_dir}/cluster_{cluster_id}_go_enrichment_failed.txt", 'w') as f:
[Link](f"GO enrichment failed for cluster {cluster_id}\n")
[Link](f"Genes in cluster: {', '.join(gene_list)}\n")
continue

if enr.res2d is None or [Link]:

print(f"No enriched terms found for cluster {cluster_id}")
with
open(f"{save_dir}/cluster_{cluster_id}_go_enrichment_no_results.txt", 'w') as f:
[Link](f"No enriched GO terms found for cluster
{cluster_id}\n")
[Link](f"Genes in cluster: {', '.join(gene_list)}\n")
continue

# Save results

enr.res2d.to_csv(f"{save_dir}/cluster_{cluster_id}_go_enrichment.csv", index=False)
print(f"Saved GO enrichment results to CSV for cluster
{cluster_id}")

# Plot top enriched GO terms

try:
top_terms = enr.res2d.sort_values('Adjusted P-
value').head(top_n)

if not top_terms.empty:
[Link](figsize=(10, 6))
[Link](top_terms['Term'].[Link](0, 50) + '...',
-np.log10(top_terms['Adjusted P-value']),
color='skyblue')
[Link]('-log10(Adjusted P-value)')
[Link]('GO Terms')
[Link](f'Top GO Terms (Cluster {cluster_id})')
[Link]().invert_yaxis()
plt.tight_layout()
[Link](f"{save_dir}/cluster_{cluster_id}_go_plot.png")
[Link]()
print(f"Created GO plot for cluster {cluster_id}")
successful_clusters += 1
except Exception as e:
print(f"Error creating GO plot for cluster {cluster_id}:
{str(e)}")
 except Exception as e:
print(f"GO Analysis failed for cluster {cluster_id}: {str(e)}")
with open(f"{save_dir}/cluster_{cluster_id}_go_error.log", 'w') as
f:
[Link](f"Error: {str(e)}\n")
[Link](f"Genes: {', '.join(gene_list)}\n")

print(f"GO Term Analysis Completed. Successful analyses for

{successful_clusters} clusters.")

def perform_kegg_analysis(self, top_n=10, save_dir="kegg_results"):

if self.best_triclusters is None:
print("No triclusters found. Run the RNSC algorithm first.")
return

[Link](save_dir, exist_ok=True)

print("\nStarting KEGG Pathway Analysis...")

# Check if gseapy is properly installed

try:
import gseapy as gp
except ImportError:
print("ERROR: gseapy package is not installed. Install it using: pip
install gseapy")
return

successful_clusters = 0

for cluster_id, nodes in self.best_triclusters.items():

if len(nodes) < 3:
print(f"Skipping cluster {cluster_id} - too few genes
({len(nodes)})")
continue

gene_list = [self.gene_names[node] for node in nodes]

print(f"Processing cluster {cluster_id} with {len(gene_list)} genes")

print(f"First few genes in cluster: {gene_list[:5]}")

try:
organisms = ['Human', 'Yeast', 'Mouse']
gene_set_options = ['KEGG_2021_Human', 'KEGG_2019_Human',
'KEGG_2019_Mouse']
success = False

for organism in organisms:

for gene_set in gene_set_options:
try:
print(f"Trying KEGG enrichment with organism:
{organism}, gene set: {gene_set}")
enr = [Link](
gene_list=gene_list,
gene_sets=[gene_set],
organism=organism,
outdir=None,
no_plot=True,
cutoff=1.0 # More permissive p-value cutoff
 )

if enr.res2d is not None and not [Link]:

success = True
print(f"KEGG enrichment successful with organism:
{organism}, gene set: {gene_set}")
break
except Exception as e:
print(f"Failed with organism {organism}, gene set
{gene_set}: {str(e)}")

if success:
break

if not success:
print(f"Could not perform KEGG enrichment for cluster
{cluster_id} with any organism/gene set")
with
open(f"{save_dir}/cluster_{cluster_id}_kegg_enrichment_failed.txt", 'w') as f:
[Link](f"KEGG enrichment failed for cluster {cluster_id}\
n")
[Link](f"Genes in cluster: {', '.join(gene_list)}\n")
continue

if enr.res2d is None or [Link]:

print(f"No enriched pathways found for cluster {cluster_id}")
with
open(f"{save_dir}/cluster_{cluster_id}_kegg_enrichment_no_results.txt", 'w') as f:
[Link](f"No enriched KEGG pathways found for cluster
{cluster_id}\n")
[Link](f"Genes in cluster: {', '.join(gene_list)}\n")
continue

enr.res2d.to_csv(f"{save_dir}/cluster_{cluster_id}_kegg_enrichment.csv",
index=False)
print(f"Saved KEGG enrichment results to CSV for cluster
{cluster_id}")

try:
top_paths = enr.res2d.sort_values('Adjusted P-
value').head(top_n)

if not top_paths.empty:
[Link](figsize=(10, 6))
[Link](top_paths['Term'].[Link](0, 50) + '...',
-np.log10(top_paths['Adjusted P-value']),
color='lightgreen')
[Link]('-log10(Adjusted P-value)')
[Link]('KEGG Pathways')
[Link](f'Top KEGG Pathways (Cluster {cluster_id})')
[Link]().invert_yaxis()
plt.tight_layout()

[Link](f"{save_dir}/cluster_{cluster_id}_kegg_plot.png")
[Link]()
print(f"Created KEGG plot for cluster {cluster_id}")
successful_clusters += 1
except Exception as e:
 print(f"Error creating KEGG plot for cluster {cluster_id}:
{str(e)}")

except Exception as e:
print(f"KEGG Analysis failed for cluster {cluster_id}: {str(e)}")
with open(f"{save_dir}/cluster_{cluster_id}_kegg_error.log", 'w')
as f:
[Link](f"Error: {str(e)}\n")
[Link](f"Genes: {', '.join(gene_list)}\n")

print(f"KEGG Pathway Analysis Completed. Successful analyses for

{successful_clusters} clusters.")

def load_dataset(filename, n_genes=10, n_conditions=5, n_time_points=8,

n_clusters=2):
"""
Load gene expression data from file or generate synthetic data.
Includes comprehensive error checking and debugging.
"""
if filename is not None:
print(f"Attempting to load data from {filename}...")
try:
try:
df = pd.read_csv('dataset_yeast.csv')
print(f"Successfully read CSV file. Columns:
{[Link]()}")
print(f"First few rows: \n{[Link](2)}")
except FileNotFoundError:
print(f"Error: File '{filename}' not found.")
return None, None, None
except Exception as e:
print(f"Error reading CSV file: {str(e)}")
return None, None, None

if 'Gene' not in [Link]:

print("Error: 'Gene' column not found in the CSV file.")
return None, None, None
if 'Expression_values' not in [Link]:
print("Error: 'Expression_values' column not found in the CSV
file.")
return None, None, None

gene_names = df['Gene'].tolist()
print(f"Found {len(gene_names)} genes.")

if not [Link]:
print(f"First expression value (raw):
{df['Expression_values'].iloc[0][:100]}...")

data = []

for idx, row in [Link]():

try:
expression_values = row['Expression_values']

if idx == 0:
print(f"Type of expression_values:
{type(expression_values)}")
 if isinstance(expression_values, str):
try:
expression_values = [Link](expression_values)
if idx == 0:
print("Successfully parsed as JSON")
except [Link]:
try:
expression_values =
ast.literal_eval(expression_values)
if idx == 0:
print("Successfully parsed as Python literal")
except (SyntaxError, ValueError) as e:
print(f"Error parsing expression values for gene
{row['Gene']}: {str(e)}")
print(f"Raw value: {expression_values[:100]}...")
expression_values = []

[Link](expression_values)

if idx == 0:
print(f"Processed first gene. Data structure shape:
{[Link](expression_values).shape if expression_values else 'empty'}")

except Exception as e:
print(f"Error processing row {idx}: {str(e)}")

try:
data = [Link](data, dtype=float)
print(f"Converted data to numpy array with shape: {[Link]}")
except Exception as e:
print(f"Error converting to numpy array: {str(e)}")
print(f"Data structure: {[type(item) for item in data[:3]]}")
return None, None, None

try:
print("Clustering genes...")
true_clusters = cluster_genes(data, n_clusters)
print(f"Generated {len(true_clusters)} cluster labels")
except Exception as e:
print(f"Error during clustering: {str(e)}")
true_clusters = [i % n_clusters for i in range(len(gene_names))]

print(f"Data successfully loaded from {filename}")

return data, gene_names, true_clusters

except Exception as e:
print(f"Unexpected error loading data from {filename}:")
traceback.print_exc()
return None, None, None

def cluster_genes(data, n_clusters):

"""
Cluster genes based on their expression patterns.

Parameters:
-----------
data : [Link]
 Gene expression data with shape (n_genes, n_conditions, n_time_points)
n_clusters : int
Number of clusters to form

Returns:
--------
[Link]
Cluster labels for each gene
"""

from [Link] import KMeans

import numpy as np

n_genes = [Link][0]
flat_data = [Link](n_genes, -1)

kmeans = KMeans(n_clusters=n_clusters, random_state=42)

cluster_labels = kmeans.fit_predict(flat_data)

return cluster_labels

def main():

# Parameters
n_genes = 10
n_conditions = 5
n_time_points = 8
n_clusters = 2
correlation_threshold = 0.7

# Loading dataset
print("Loading input dataset...")
data, gene_names, true_clusters = load_dataset(
filename="dataset_yeast.csv",
n_genes=n_genes,
n_conditions=n_conditions,
n_time_points=n_time_points,
n_clusters=n_clusters
)

print(f"Generated data for {n_genes} genes across {n_conditions} conditions and

{n_time_points} time points")
print(f"True clusters: {[[Link](true_clusters == i) for i in
range(n_clusters)]}")

# Initialize TriRNSC
trirnsc = TriRNSC(
threshold=correlation_threshold,
tabu_length_ratio=0.02,
naive_stopping_tolerance=15,
scaled_stopping_tolerance=15,
diversification_frequency=50,
diversification_length=10,
max_experiments=3
)

trirnsc.load_data(data, gene_names)
 print("\nVisualizing true expression patterns...")
trirnsc.visualize_gene_expression_patterns(true_clusters)

print("\nConstructing and visualizing GCN...")

trirnsc.construct_gene_co_expression_network()

print("\nVisualizing correlation matrix...")

trirnsc.visualize_correlation_matrix()
trirnsc.visualize_correlation_matrix(true_clusters)

print("\nVisualizing GCN with true clusters...")

trirnsc.visualize_gcn_2d(true_clusters)
trirnsc.visualize_gcn_3d(true_clusters)

print("\nRunning TriRNSC algorithm...")

trirnsc.run_rnsc()

tricluster_info = trirnsc.get_tricluster_info()
print("\nTop 3 triclusters:")
for i, tc in enumerate(tricluster_info[:3]):
print(f"Tricluster {tc['id']}: {tc['size']} genes, MSR={tc['msr']:.6f},
TQI={tc['tqi']:.6e}")

print("\nVisualizing GCN with discovered clusters...")

trirnsc.visualize_gcn_2d()
trirnsc.visualize_gcn_3d()

print("\nVisualizing expression patterns of discovered clusters...")

trirnsc.visualize_gene_expression_patterns()

trirnsc.perform_go_analysis(top_n=10)
trirnsc.perform_kegg_analysis(top_n=10)

if __name__ == "__main__":
main()