NeDis: Network-based disruption analysis of biological coordination

random_state = 43

from nedis.data.synthetic import load_example
from nedis.visualization import visualize_data

# load data
X, y, entities, labels = load_example(random_state=random_state)

# visualize
fig, axes, correlation_matrices, coordinates = visualize_data(
    X, y, entities, mode="network", random_state=random_state);

from nedis.cordis.default \
    import DefaultCorrelationDisruptionFeatureTransformer \
    as DefaultCorrelationDisruptionTransformer

# configure and fir correlation disruption transformer
disruption_transformer = DefaultCorrelationDisruptionTransformer(
    default_optimization_separation_score="spearman",
    default_derive_features_aggregation="mean")
disruption_transformer.fit(X, y, groups=entities, subset_masks="y");

# transform samples
disruption_values = disruption_transformer.transform(X)

# we get one aggregated disruption value per cluster for each sample
disruption_values.shape

import numpy as np
import scipy.stats

import matplotlib.pyplot as plt
import seaborn as sns

from nedis.visualization import plot_cordis_cluster as plot_cluster

y_unique = np.unique(y)

for i_cluster, cluster in enumerate(disruption_transformer.selected_clusters_):
    
    values = disruption_values[:, i_cluster]
    r, p = scipy.stats.spearmanr(values, y)
    
    fig, axes = plt.subplots(
        1, len(y_unique) + 1, 
        figsize=(4 * 1 * (len(y_unique) + 1), 4 * 1), 
        dpi=300)
 
    # correlation disruption plot
    ax = axes[0]
    x_rank = scipy.stats.rankdata(y, method="dense")
    sns.lineplot(x=x_rank, y=values, hue=entities, color="blue", alpha=0.1, ax=ax)
    sns.lineplot(x=x_rank, y=values, ax=ax)
    ax.set(
        xlabel="timepoints",
        xticks=np.unique(x_rank),
        xticklabels= y_unique,
        ylabel="disruption",
        title=f"r={r:.02f}"
    )
    ax.get_legend().remove()
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    
    # cluster visualization
    for i, yy in enumerate(y_unique):
        ax = axes[i + 1]
        ax.axis("off")
        plot_cluster(
            cluster, 
            coordinates, 
            correlation_matrices[yy], 
            correlation_threshold=0, 
            verbose=0,
            ax=ax)
    fig.suptitle(
        f"Cluster: Reference data={cluster['reference_label']}, Id={cluster['id']}")
    
    plt.show()

from nedis.cluster.leidenalg import WeightedLeidenClustering
from nedis.cordis.clustering import ReferenceCorrelationMatrixClusteringStep
from nedis.cordis.optimization import GreedyRefinementOptimizationStep
from nedis.cordis.disruption import CorrelationDisruption

# NOTE: you can use any cluster algorithm that assigns a cluster label to each feature
clustering_algorithm = WeightedLeidenClustering(random_state=random_state)

clustering_step = ReferenceCorrelationMatrixClusteringStep(
    clustering_algorithm=clustering_algorithm,
#     clustering_absolute_correlation=True,
#     correlation_function="spearman",
#     feature_filters=None
)

optimization_step = GreedyRefinementOptimizationStep(
    separation_score="spearman",
#     separation_score_comparison='all',
#     refinement_mode="rows-and-columns",
#     correlation_function="spearman",
#     disruption_metric="direction",
#     disruption_robustness='loo',
#     disruption_aggregation='mean',
#     max_runs=-1
)

codi = CorrelationDisruption(
    clustering_step=clustering_step,
    cluster_optimization_step=optimization_step,
    # filtering
    filter_coverage_threshold=0.5, 
    separation_score_threshold=("auto", 1)
)
codi.fit(X, y, subset_masks="y")

from nedis.cordis.clustering import (
    ListClusteringStep, ReferenceFeatureLabelClusteringStep, init_cluster)
from nedis.cordis.optimization import (
    ReferenceScoreOptimizationStep)

# use to define clusters by feature labels
# `labels` is defined in the data creation above
clustering_step = ReferenceFeatureLabelClusteringStep(labels)

# # ALTERNATIVE: completely custom clusters
# clustering_step = ListClusteringStep([
#     init_cluster(reference_label=0, reference_shape=X.shape[1], rows=np.arange(0,5)),
#     init_cluster(reference_label=4, reference_shape=X.shape[1], rows=np.arange(5,15)),
#     init_cluster(reference_label=0, reference_shape=X.shape[1], rows=np.arange(15,25)),
# ])

# use this to prevent optimization 
# ALTERNATIVE: use the GreedyRefinementOptimizationStep from above
optimization_step = ReferenceScoreOptimizationStep(
    separation_score="spearman",
)

codi = CorrelationDisruption(
    clustering_step=clustering_step,
    cluster_optimization_step=optimization_step,
    # filtering
    filter_coverage_threshold=0.5, 
    separation_score_threshold=("auto", 1)
)
codi.fit(X, y, subset_masks="y")