Merge pull request #655 from FZJ-INM1-BDA/maint_paper_nbs

Maint paper nbs

Author: AhmetNSimsek

CITATION.cff

@@ -15,9 +15,9 @@ license: Apache-2.0
 
 repository-code: https://github.com/FZJ-INM1-BDA/siibra-python
 
-version: v1.0.1-alpha.2
+version: v1.0.1-alpha.4
 
-date-released: 2025-03-12
+date-released: 2025-04-10
 
 doi: 10.5281/zenodo.7885728
 

examples/tutorials/2025-paper-fig3.py

@@ -27,7 +27,6 @@
 """
 
 # %%
-import re
 import numpy as np
 import pandas as pd
 import matplotlib.pyplot as plt
@@ -39,6 +38,7 @@
 
 sns.set_style("dark")
 
+
 # %%
 # Input: Activation map or other feature distribution as image volume in MNI space
 # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -72,7 +72,7 @@
 
 # %%
 # Split input volume into cluster components
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 #
 # There are many ways to get components out of a feature map. Here we use siibra to
 # - draw random points from the distribution encoded by the input volume, then
@@ -105,11 +105,11 @@
 )
 
 # %%
-# Assign peaks and clusters to cytoarchitectonic regions
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+# Assign clusters to cytoarchitectonic regions
+# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 #
 # To assign the clusters to brain regions, we build feature maps from each cluster
-# and assign them to the Julich-Brain probabilistic maps. The assignemint is one-to-many
+# and assign them to the Julich-Brain probabilistic maps. The assignment is one-to-many
 # since the structures in the image and parcellation are continuous. Assignments
 # report correlation, intersection over union, and some other measures which we can use
 # to filter and sort them. The result is an assignment table from cluster components
@@ -120,24 +120,8 @@
 pmaps = siibra.get_map(
     parcellation="julich 3.0.3", space="mni152", maptype="statistical"
 )
-assignments = []
-
-# assign peaks to regions
-peaks = input_volume.find_peaks(mindist=5, sigma_mm=0)
-with siibra.QUIET:
-    df = pmaps.assign(peaks)
-df.query(f"`map value` >= {min_map_value}", engine="python", inplace=True)
-df["type"] = "peak"
-df["id"] = df["input structure"]
-assignments.append(df[["type", "id", "region", "map value"]])
-
-view = plotting.plot_glass_brain(
-    input_volume.fetch(), alpha=1, cmap="RdBu", symmetric_cbar=True
-)
-view.add_markers(peaks.as_list(), marker_size=30)
 
-# %%
-# assign clusters to regions
+assignments = []
 for cl in clusterlabels:
     clustermap = siibra.volumes.from_pointcloud(samples, label=cl, target=input_volume)
     plotting.plot_glass_brain(
@@ -150,12 +134,11 @@
     with siibra.QUIET:
         df = pmaps.assign(clustermap)
     df.query(f"correlation >= {min_correlation}", engine="python", inplace=True)
-    df["type"] = "cluster"
-    df["id"] = cl
-    assignments.append(df[["type", "id", "region", "correlation", "map value"]])
+    df['id'] = cl
+    assignments.append(df[["id", "region", "correlation", "map value"]])
 
 all_assignments = pd.concat(assignments).sort_values(by="correlation", ascending=False)
-all_assignments
+all_assignments[:9].round(2)
 
 # %%
 # plot the three primary assigned probability maps
@@ -174,64 +157,58 @@
 # Find features
 # ^^^^^^^^^^^^^
 #
-# To demonstrate multimodal feature profiling, we only choose the first connected region.
-
-
-def shortname(region):
-    return (
-        re.sub("\s*\(.*\)\s*|\s*Areaa\s*", " ", region.name)
-        .replace("left", "L")
-        .replace("right", "R")
-        .strip()
-    )
-
-
-def filter_features(feats, region):
-    return [f for f in feats if any(r.assign(region) for r in f.anchor.regions)]
+# To demonstrate multimodal feature profiling, we choose one of the connected regions.
+selected_region = siibra.get_region("julich 3.0.3", "Area hOc1 (V1, 17, CalcS) left")
 
+# %%
+# Let us first query receptor density fingerprint for this region
+receptor_fingerprints = siibra.features.get(selected_region, "receptor density fingerprint")
+fp = receptor_fingerprints[0]
+print("\n".join(fp.urls))
+fig, axs = plt.subplots(1, 1, figsize=(4, 3.2))
+fp.plot(ax=axs, capsize=4)
 
-def plot_receptors(region, ax):
-    fts = filter_features(siibra.features.get(region, "ReceptorDensityFingerprint"), region)
-    fts[0].plot(ax=ax)
 
+# %%
+# Now, query for gene expressions for the same region
+genes = ["gabarapl1", "gabarapl2", "maoa", "tac1"]
+gene_expressions = siibra.features.get(selected_region, "gene expressions", gene=genes)
+print("\n".join(gene_expressions[0].urls))
+fig, axs = plt.subplots(1, 1, figsize=(4, 3.5))
+gene_expressions[0].plot(ax=axs)
 
-def plot_celldensities(region, ax):
-    fts = filter_features(siibra.features.get(region, "LayerwiseCellDensity"), region)
-    fts[0].plot(ax=ax)
 
+# %%
+# We can next obtain cell body density values for this region
+cell_densities = siibra.features.get(selected_region, "LayerwiseCellDensity")
+print("\n".join(cell_densities[0].urls))
+fig, axs = plt.subplots(1, 1, figsize=(4, 2.8))
+cell_densities[0].plot(ax=axs, textwrap=30)
+# axs.set_title(axs.get_title().replace("in", "\n"))
 
-def plot_gene_expressions(region, ax, genes):
-    fts = siibra.features.get(region, "GeneExpressions", gene=genes)
-    fts[0].plot(ax=ax)
+# %%
+# Lastly, we can obtain the regional profile of streamline count type
+# parcellation-based connectivity matrices
+connectivity_matrices = siibra.features.get(selected_region, "StreamlineCounts")
+conn = connectivity_matrices[0]  # select the first set of matrices
+print("\n".join(conn.urls))
 
 
-def plot_connectivity(region, ax):
-    fts = siibra.features.get(region, "StreamlineCounts")
-    conndata = fts[0][0].get_profile(region).data
-    conndata.rename(index={r: shortname(r) for r in conndata.index}, inplace=True)
-    conndata[:15].plot(kind="bar", ax=ax)
-    plt.xticks(ha="right")
-    plt.tight_layout()
-    plt.grid(True)
-    plt.title(f"Connectivity for {region.name}")
+def shorten_name(region):
+    # to simplify readability
+    return (
+        region.replace("Area ", "")
+        .replace(" (GapMap)", "")
+        .replace("left", "L")
+        .replace("right", "R")
+    )
 
 
-# %%
-selected_region = siibra.get_region("julich 3.0.3", "Area hOc1 (V1, 17, CalcS) left")
-plot_funcs = [
-    lambda r, a: plot_receptors(r, a),
-    lambda r, a: plot_celldensities(r, a),
-    lambda r, a: plot_connectivity(r, a),
-    lambda r, a: plot_gene_expressions(
-        r, a, ["gabarapl1", "gabarapl2", "maoa", "tac1"]
-    ),
-]
-
-fig, axs = plt.subplots(1, len(plot_funcs), figsize=(3.5 * len(plot_funcs), 4))
-for ax, func in zip(axs.ravel(), plot_funcs):
-    func(r=selected_region, a=ax)
-    ax.grid(True)
-    xtl = ax.get_xticklabels()
-    ax.set_xticklabels(xtl, rotation=55, ha="right")
-plt.suptitle("")
-plt.tight_layout()
+fig, axs = plt.subplots(1, 1, figsize=(4, 4.5))
+conn.plot(selected_region, ax=axs, max_rows=15, kind="bar", rot=50, width=0.7)
+axs.set_ylabel(
+    f"Mean of {conn.modality} \u00b1 std \n in {len(conn.elements)} {conn.indexing_attributes[0]}s"
+)
+axs.xaxis.set_ticklabels([shorten_name(t.get_text()) for t in axs.xaxis.get_majorticklabels()])
+plt.grid(True, 'minor')
+plt.title(f"Connectivity for {selected_region.name}")

examples/tutorials/2025-paper-fig4.py

@@ -23,8 +23,6 @@
 # %%
 import matplotlib.pyplot as plt
 from nilearn import plotting
-import pandas as pd
-import re
 import seaborn as sns
 import siibra
 
@@ -38,6 +36,7 @@
 
 jubrain = siibra.parcellations.JULICH_BRAIN_CYTOARCHITECTONIC_ATLAS_V3_0_3
 pmaps = jubrain.get_map(space="mni152", maptype="statistical")
+print(jubrain.publications[0]['citation'])
 
 # %%
 # Obtain definitions and probabilistic maps in MNI asymmetric space of area IFG 44
@@ -53,73 +52,59 @@
         colorbar=False,
         annotate=False,
         symmetric_cbar=True,
+        title=r.name,
     )
-    plt.savefig(f"{r.key}.png")
 
 # %%
-# For each of the two brain areas, retrieve average densities of a selection of monogenetic
-# neurotransmitter receptors, layer-specific distributions of cell bodies, as well as expressions
-# of a selection of genes coding for these receptors.
+# For each of the two brain areas, query for layer-specific distributions of cell bodies.
+fig, axs = plt.subplots(1, len(regions), sharey=True, figsize=(8, 2.7))
+for i, region in enumerate(regions):
+    layerwise_cellbody_densities = siibra.features.get(region, "layerwise cell density")
+    layerwise_cellbody_densities[0].plot(ax=axs[i], textwrap=25)
+    print(layerwise_cellbody_densities[0].urls)
+    axs[i].set_ylim(25, 115)
 
-receptors = ["M1", "M2", "M3", "D1", "5-HT1A", "5-HT2"]
-genes = ["CHRM1", "CHRM2", "CHRM3", "HTR1A", "HTR2A", "DRD1"]
-modalities = [
-    ("receptor density fingerprint", {}, {"receptors": receptors, "rot": 90}),
-    ("layerwise cell density", {}, {"rot": 0}),
-    ("gene expressions", {"gene": genes}, {"rot": 90}),
-]
-fig, axs = plt.subplots(
-    len(modalities) + 1, len(regions), figsize=(4 * len(regions), 11), sharey="row"
-)
-ymax = [1000, 150, None]
+# %%
+# Next, retrieve average densities of a selection of monogenetic
+# neurotransmitter receptors.
+receptors = ["M1", "M2", "M3", "5-HT1A", "5-HT2", "D1"]
+fig, axs = plt.subplots(1, len(regions), sharey=True, figsize=(8, 3.5))
+for i, region in enumerate(regions):
+    receptor_fingerprints = siibra.features.get(region, "receptor density fingerprint")
+    receptor_fingerprints[0].plot(receptors=receptors, ax=axs[i])
+    print(receptor_fingerprints[0].urls)
+    axs[i].set_ylim(0, 1000)
+    axs[i].title.set_size(12)
+    transmitters = [
+        siibra.vocabularies.RECEPTOR_SYMBOLS[r]['neurotransmitter']['name']
+        for r in receptors
+    ]
+    axs[i].set_xticklabels([f"{r}\n({n})" for r, n in zip(receptors, transmitters)])
 
+# %%
+# Now, for further insight, query for expressions of a selection of genes coding
+# for these receptors.
+genes = ["CHRM1", "CHRM2", "CHRM3", "HTR1A", "HTR2A", "DRD1"]
+fig, axs = plt.subplots(1, len(regions), sharey=True, figsize=(7, 3))
 for i, region in enumerate(regions):
-    axs[0, i].imshow(plt.imread(f"{region.key}.png"))
-    axs[0, i].set_title(f'{region.name.replace("Area ", "")}')
-    axs[0, i].axis("off")
-    for j, (modality, kwargs, plotargs) in enumerate(modalities):
-        fts = siibra.features.get(region, modality, **kwargs)
-        assert len(fts) > 0
-        if len(fts) > 1:
-            print(f"More than one feature found for {modality}, {region.name}")
-        f = fts[0]
-        f.plot(ax=axs[j + 1, i], **plotargs)
-        if modality == "receptor density fingerprint":
-            # add neurotransmitter names to  receptor names in xtick labels
-            transmitters = [re.sub(r"(^.*\()|(\))", "", n) for n in f.neurotransmitters]
-            axs[j + 1, i].set_xticklabels(
-                [f"{r}\n({n})" for r, n in zip(receptors, transmitters)]
-            )
-        if ymax[j] is not None:
-            axs[j + 1, i].set_ylim(0, ymax[j])
-        if "std" in axs[j + 1, i].yaxis.get_label_text():
-            axs[j + 1, i].set_ylabel(
-                axs[j + 1, i].yaxis.get_label_text().replace("std", "std\n")
-            )
-        axs[j + 1, i].set_title(f"{fts[0].modality}")
-fig.suptitle("")
-fig.tight_layout()
+    gene_expressions = siibra.features.get(region, "gene expressions", gene=genes)
+    assert len(gene_expressions) == 1
+    gene_expressions[0].plot(ax=axs[i])
+    axs[i].title.set_size(11)
+    print(gene_expressions[0].urls)
 
 # %%
 # For each of the two brain areas, collect functional connectivity profiles referring to
 # temporal correlation of fMRI timeseries of several hundred subjects from the Human Connectome
 # Project. We show the strongest connections per brain area for the average connectivity patterns
-fts = siibra.features.get(regions[0], "functional connectivity")
-conn = fts[1]
-
-# aggregate connectivity profiles for first region across subjects
-D1 = (
-    pd.concat([c.get_profile(regions[0]).data for c in conn], axis=1)
-    .agg(["mean", "std"], axis=1)
-    .sort_values(by="mean", ascending=False)
-)
-
-# aggregate connectivity profiles for second region across subjects
-D2 = (
-    pd.concat([c.get_profile(regions[1]).data for c in conn], axis=1)
-    .agg(["mean", "std"], axis=1)
-    .sort_values(by="mean", ascending=False)
-)
+fts = siibra.features.get(jubrain, "functional connectivity")
+for cf in fts:
+    if cf.cohort != "HCP":
+        continue
+    if cf.paradigm == "Resting state (EmpCorrFC concatenated)":
+        conn = cf
+        break
+print(conn.urls)
 
 
 # plot both average connectivity profiles to target regions
@@ -132,35 +117,32 @@ def shorten_name(n):
     )
 
 
-fig, (a1, a2) = plt.subplots(1, 2, sharey=True, figsize=(3.6 * len(regions), 4.1))
-kwargs = {"kind": "bar", "width": 0.85, "logy": True}
-D1.iloc[:17]["mean"].plot(
-    **kwargs, yerr=D1.iloc[:17]["std"], ax=a1, ylabel=shorten_name(regions[0].name)
-)
-D2.iloc[:17]["mean"].plot(
-    **kwargs, yerr=D2.iloc[:17]["std"], ax=a2, ylabel=shorten_name(regions[1].name)
-)
-a1.set_ylabel("temporal correlation")
-a1.xaxis.set_ticklabels(
-    [shorten_name(t.get_text()) for t in a1.xaxis.get_majorticklabels()]
-)
-a2.xaxis.set_ticklabels(
-    [shorten_name(t.get_text()) for t in a2.xaxis.get_majorticklabels()]
-)
-a1.grid(True)
-a2.grid(True)
-a1.set_title(regions[0].name.replace("Area ", ""))
-a2.set_title(regions[1].name.replace("Area ", ""))
+plotkwargs = {
+    "kind": "bar",
+    "width": 0.85,
+    "logscale": True,
+    "xlabel": "",
+    "ylabel": "temporal correlation",
+    "rot": 90,
+    "ylim": (0.3, 1.1),
+    "capsize": 2,
+}
+fig, axs = plt.subplots(1, len(regions), sharey=True, figsize=(7, 4.5))
+for i, region in enumerate(regions):
+    plotkwargs["ax"] = axs[i]
+    conn.plot(region, max_rows=17, **plotkwargs)
+    axs[i].xaxis.set_ticklabels([shorten_name(t.get_text()) for t in axs[i].xaxis.get_majorticklabels()])
+    axs[i].set_title(region.name.replace("Area ", ""), fontsize=8)
+    axs[i].grid(True, 'minor')
 plt.suptitle("Functional Connectivity")
-plt.tight_layout()
 
 # %%
 # For both brain areas, sample representative cortical image patches at 1µm
 # resolution that were automatically extracted from scans of BigBrain sections.
 # The image patches display clearly the differences in laminar structure of the two regions
 
 selected_sections = [4950, 1345]
-fig, axs = plt.subplots(1, len(regions))
+fig, axs = plt.subplots(1, len(regions), sharey=True)
 for r, sn, ax in zip(regions, selected_sections, axs):
     # find 1 micron sections intersecting the region
     pmap = pmaps.get_volume(r)
@@ -176,6 +158,4 @@ def shorten_name(n):
     ax.axis("off")
     ax.set_title(r.name)
 
-plt.tight_layout()
-
 # sphinx_gallery_thumbnail_number = 2