<jupyter_start><jupyter_text>Methods for Multiscale Comparative ConnectomicsThis demo shows you how to use methods in `graspologic` to analyze patternsin brain connectivity in connectomics datasets. We specifically demonstratemethods for identifying differences in edges and vertices across subjects.<jupyter_code>import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import warnings
pd.options.mode.chained_assignment = None
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Load the Duke mouse brain datasetDataset of 32 mouse connectomes derived from whole-brain diffusionmagnetic resonance imaging of four distinct mouse genotypes:BTBR T+ Itpr3tf/J (BTBR), C57BL/6J(B6), CAST/EiJ (CAST), and DBA/2J (DBA2).For each strain, connectomes were generated from eight age-matched mice(N = 8 per strain), with a sex distribution of four males and four females.Each connectome was parcellated using asymmetric Waxholm Space, yielding avertex set with a total of 332 regions of interest (ROIs) symmetricallydistributed across the left and right hemispheres. Within a givenhemisphere, there are seven superstructures consisting up multiple ROIs,resulting in a total of 14 distinct communities in each connectome.<jupyter_code>from graspologic.datasets import load_mice
# Load the full mouse dataset
mice = load_mice()

# Stack all adjacency matrices in a 3D numpy array
graphs = np.array(mice.graphs)

# Get sample parameters
n_subjects = mice.meta["n_subjects"]
n_vertices = mice.meta["n_vertices"]
# Split the set of graphs by genotype
btbr = graphs[mice.labels == "BTBR"]
b6 = graphs[mice.labels == "B6"]
cast = graphs[mice.labels == "CAST"]
dba2 = graphs[mice.labels == "DBA2"]

connectomes = [btbr, b6, cast, dba2]<jupyter_output><empty_output><jupyter_text>Identifying Signal EdgesThe simplest approach for comparing connectomes is to treat them as a _bag of edges_ without considering interactions between the edges.Serially performing univariate statistical analyses at each edge enables the discovery of _signal edges_ whose neurological connectivity differs across categorical or dimensional phenotypes.Here, we demonstrate the possibility of using Distance Correlation (`Dcorr`), a nonparametric universally consistent test, to detect changes in edges.In this model, we assume that each edge in the connectome is independently andidentically sampled from some distribution $F_i$, where $i$ represents the groupto which the given connectome belongs. In this setting the groups are the mousegenotypes. `Dcorr` allows us to test the following hypothesis:\begin{align*}H_0:\ &F_1 = F_2 = \ldots F_k \\H_A:\ &\exists \ j \neq j' \text{ s.t. } F_j \neq F_{j'}\end{align*}<jupyter_code>from hyppo.ksample import KSample<jupyter_output><empty_output><jupyter_text>Since the connectomes in this dataset are undirected, we only need to do edgecomparisons on the upper triangle of the adjacency matrices.<jupyter_code># Make iterator for traversing the upper triangle of the connectome
indices = zip(*np.triu_indices(n_vertices, 1))
edge_pvals = []

for roi_i, roi_j in indices:
    
    # Get the (i,j)-th edge for each connectome
    samples = [genotype[:, roi_i, roi_j] for genotype in connectomes]
    
    # Calculate the p-value for the (i,j)-th edge
    try:
        statistic, pvalue = KSample("Dcorr").test(*samples, reps=10000, workers=-1)
    except ValueError:
        # A ValueError is thrown when any of the samples have equal edge 
        # weights (i.e. one of the inputs has 0 variance)
        statistic = np.nan
        pvalue = 1

    edge_pvals.append([roi_i, roi_j, statistic, pvalue])<jupyter_output><empty_output><jupyter_text>Connectomes are a high-dimensional dataype.Thus, statistical tests on components of the connectome (e.g. edge and vertices)results in multiple comparisons.We recommend correcting for multiple comparisons using the Holm-Bonferroni (HB)corection.The algorithm is described below:1. Sort the p-values from lowest-to-highest $P_1, P_2, \dots, P_n$, where $n$ is the number of tests2. Correct the p-value as $P_1(n), P_2(n-1), \dots, P_n(1)$3. If any corrected p-value is $>1$, replace with $1$4. If the corrected p-value is less than a significance level $\alpha$, reject<jupyter_code># Convert the nested list to a dataframe
edge_pvals = pd.DataFrame(edge_pvals, columns=["ROI_1", "ROI_2", "stat", "pval"])

# Correct p-values using the Holm-Bonferroni correction
edge_pvals.sort_values(by="pval", inplace=True, ignore_index=True)
pval_rank = edge_pvals["pval"].rank(ascending=False, method="max")
edge_pvals["holm_pval"] = edge_pvals["pval"].multiply(pval_rank)
edge_pvals["holm_pval"] = edge_pvals["holm_pval"].apply(
    lambda pval: 1 if pval > 1 else pval
)

# Test for significance using alpha=0.05
alpha = 0.05
edge_pvals["significant"] = (edge_pvals["holm_pval"] < alpha)

# Get the top 10 strongest signal edges
edge_pvals_top = edge_pvals.head(10)

# Replace ROI indices with actual names
def lookup_roi_name(index):
    hemisphere = "R" if index // 166 else "L"
    index = index % 166
    roi_name = mice.atlas.query(f"ROI == {index+1}")["Structure"].item()
    return f"{roi_name} ({hemisphere})"

edge_pvals_top["ROI_1"] = edge_pvals_top["ROI_1"].apply(lookup_roi_name)
edge_pvals_top["ROI_2"] = edge_pvals_top["ROI_2"].apply(lookup_roi_name)

edge_pvals_top.head()<jupyter_output><empty_output><jupyter_text>Note that none of the edges achieve significance at $\alpha=0.05$ following Holm-Bonferroni correction.We might expect this, given that we are correcting for $N=54,946$ tests.Instead of the magnitude, the **ranking of the p-values** can be used to determine signal edges.  Identifying Signal VerticesA sample of connectomes can be jointly embedded in a low-dimensional Euclidean space using the omnibus embedding (`omni`).A host of machine learning tasks can be accomplished with this joint embedded representation of the connectome, such as clustering or classification of vertices.Here, we use the embedding to formulate a statistical test that can be used to identify vertices that are strongly associated with given phenotypes.According to a Central Limit Theorem for `omni`, these latent positions are universally consistent and asymptotically normal.<jupyter_code>from graspologic.embed import OmnibusEmbed
from graspologic.plot import pairplot
# Jointly embed graphs using omnibus embedding
embedder = OmnibusEmbed()
omni_embedding = embedder.fit_transform(graphs)
print(omni_embedding.shape)<jupyter_output><empty_output><jupyter_text>For each of the 32 mice, `omni` embeds each vertex in the connectome to alatent position vector $x_i \in \mathbb{R}^5$.We test for differences in the distribution of vertex latent positions using the `MGC` independence test implemented in `hyppo`.This test essentially acts as a nonparametric MANOVA test.<jupyter_code>with warnings.catch_warnings():
    warnings.filterwarnings("ignore", module="hyppo")
    warnings.filterwarnings("ignore", category=RuntimeWarning)
    vertex_pvals = []

    for vertex_i in range(n_vertices):
    
    # Get the embedding of the i-th vertex for 
        samples = [
            omni_embedding[mice.labels==genotype, vertex_i, :] 
            for genotype in np.unique(mice.labels)
        ]
    
    # Calculate the p-value for the i-th vertex
        statistic, pvalue, _ = KSample("MGC").test(
            *samples, reps=250, workers=-1
        )
        vertex_pvals.append([vertex_i, statistic, pvalue])<jupyter_output><empty_output><jupyter_text>**Important:** When running on data, make sure `reps` $\geq \frac{|V|}{\alpha}$ where $|V|$ is the cardinality of the vertex set and $\alpha$ is the significance level for your test.Again, this has to do with p-value correction:because `MGC` is a permutation test, you need many replications to ensure accurate computation of the p-value.Because this notebook is primarily for demonstration purposes, we use a small value of `reps=250`. For this dataset, a more appropriate number of permutations would be around $10^6$.<jupyter_code># Convert the nested list to a dataframe
vertex_pvals = pd.DataFrame(vertex_pvals, columns=["ROI", "stat", "pval"])

# Correct p-values using the Holm-Bonferroni correction
vertex_pvals.sort_values(by="pval", inplace=True, ignore_index=True)
pval_rank = vertex_pvals["pval"].rank(ascending=False, method="max")
vertex_pvals["holm_pval"] = vertex_pvals["pval"].multiply(pval_rank)
vertex_pvals["holm_pval"] = vertex_pvals["holm_pval"].apply(
    lambda pval: 1 if pval > 1 else pval
)

# Test for significance using alpha=0.05
alpha = 0.05
vertex_pvals["significant"] = (vertex_pvals["holm_pval"] < alpha)

# Get the top 10 strongest signal edges
vertex_pvals_top = vertex_pvals.head(10)

# Replace ROI indices with actual names
vertex_pvals_top["ROI"] = vertex_pvals_top["ROI"].apply(lookup_roi_name)

vertex_pvals_top.head()<jupyter_output><empty_output><jupyter_text>We can use pairplots to visualize the embeddings of specific vertices.Below are pairsplots of the corpus callosum from the left and right hemispheres.<jupyter_code>left = pairplot(
    omni_embedding[:, 121 - 1, :3], 
    mice.labels, 
    palette=["#e7298a", "#1b9e77", "#d95f02", "#7570b3"],
    title="Left Corpus Callosum"
)
right = pairplot(
    omni_embedding[:, 121 + 166 - 1, :3], 
    mice.labels, 
    palette=["#e7298a", "#1b9e77", "#d95f02", "#7570b3"],
    title="Right Corpus Callosum",
)<jupyter_output><empty_output><jupyter_text>Whole-brain ComparisonsWe can use the results of the omnibus embedding to perform whole-brain comparisonsacross subjects from different phenotypes.<jupyter_code>from graspologic.embed import ClassicalMDS
from graspologic.plot import heatmap<jupyter_output><empty_output><jupyter_text>Two-dimensional representations of each connectome were obtained by using Classical Multidimensional Scaling (`cmds`) to reduce the dimensionality of the embeddings obtained by `omni`.<jupyter_code># Further reduce embedding dimensionality using cMDS
cmds = ClassicalMDS(2)
cmds_embedding = cmds.fit_transform(omni_embedding)
cmds_embedding = pd.DataFrame(cmds_embedding, columns=["Dimension 1", "Dimension 2"])
cmds_embedding["Genotype"] = mice.labels

# Embedding with BTBR
sns.scatterplot(
    x="Dimension 1",
    y="Dimension 2",
    hue="Genotype",
    data=cmds_embedding,
    palette=["#e7298a", "#1b9e77", "#d95f02", "#7570b3"],
)
plt.show()
# Find the dissimilarity between subjects' connectomes using the cMDS embedding
dis = cmds.dissimilarity_matrix_
scaled_dissimilarity = dis / np.max(dis)

heatmap(scaled_dissimilarity,
        context="paper",
        inner_hier_labels=mice.labels,)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09290178155022852, 'fraction_numerical': 0.016136554425893106, 'mean_word_length': 5.382142857142857, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 48, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '89205', 'license': 'MIT', 'n_tokens_mistral': 3485, 'n_tokens_neox': 3242, 'n_words': 1103, 'path': 'docs/tutorials/connectomics/mcc.ipynb', 'repo_name': 'aaitbr/graspologic', 'size': '10721'}