BraiAn: a demonstration

This notebook assumes you have done the following steps:

• alignment of brain slices in ABBA, exported to a QuPath project as annotations.
• used qupath-extension-braian to run the image analysis comprising of cell detection and. if needed, cell classification and cell overlaps
• tagged the regions to exclude from the analysis due to tissue, imaging or alignment problems.
  This is done for each image slice of the QuPath projects, as explained here)
• Exported detections countings and the list of regions to exclude using qupath-extension-braian into a results and a regions_to_exclude.
  The first should be exported to .csv files (one per slice). The latter to .txt files (one per slice).

All these steps should be done for each brain of the group to analyse. For this demostration, we will be using this dataset, placed in the same folder as the notebook.

Before we start ...

We set some parameters for the analysis

import braian
import braian.plot as bap
import braian.stats as bas
import os
import plotly.io as pio

# This ensures BraiAn's figures works in multiple places:
pio.renderers.default = "plotly_mimetype+notebook"

# Where do we read data and save the results of the computation
DATA_INPUT_PATH  = "/mnt/tenibre/projects/data/IEGs/cFos_Arc1/QuPath_output"
DATA_OUTPUT_PATH = "./BraiAn_output"
if not(os.path.exists(DATA_OUTPUT_PATH)):
    os.makedirs(DATA_OUTPUT_PATH, exist_ok=True)

# Here we define which exported folders, in DATA_INPUT_PATH, are corresponding to which animal group
group1 = dict(
    name="HC",  # Home Cage
    dirs=["287HC", "342HC", "343HC", "346HC", "371HC"]
)
group2 = dict(
    name="CTX", # Context
    dirs=["329CTX", "331CTX", "355CTX", "400CTX", "401CTX", "402CTX"]
)

# Parameters for the Allen Brain atlas
ALLEN_VERSION = "CCFv3"                                             # this version depends on the annotation used in ABBA when aligning the brain slices
EXCLUDED_BRANCHES = ["fiber tracts", "VS", "grv", "retina", "CB"]   # some regions can be excluded from the start for every slice since we are not interested in them

import braian import braian.plot as bap import braian.stats as bas import os import plotly.io as pio # This ensures BraiAn's figures works in multiple places: pio.renderers.default = "plotly_mimetype+notebook" # Where do we read data and save the results of the computation DATA_INPUT_PATH = "/mnt/tenibre/projects/data/IEGs/cFos_Arc1/QuPath_output" DATA_OUTPUT_PATH = "./BraiAn_output" if not(os.path.exists(DATA_OUTPUT_PATH)): os.makedirs(DATA_OUTPUT_PATH, exist_ok=True) # Here we define which exported folders, in DATA_INPUT_PATH, are corresponding to which animal group group1 = dict( name="HC", # Home Cage dirs=["287HC", "342HC", "343HC", "346HC", "371HC"] ) group2 = dict( name="CTX", # Context dirs=["329CTX", "331CTX", "355CTX", "400CTX", "401CTX", "402CTX"] ) # Parameters for the Allen Brain atlas ALLEN_VERSION = "CCFv3" # this version depends on the annotation used in ABBA when aligning the brain slices EXCLUDED_BRANCHES = ["fiber tracts", "VS", "grv", "retina", "CB"] # some regions can be excluded from the start for every slice since we are not interested in them

Allen's brain region ontology

Since we work with mice, we place our data in Allen Institute's Common Coordinate Framework

path_to_allen_json = "AllenMouseBrainOntology.json"
braian.cache(path_to_allen_json, "http://api.brain-map.org/api/v2/structure_graph_download/1.json")
brain_ontology = braian.AllenBrainOntology(path_to_allen_json, EXCLUDED_BRANCHES, version=ALLEN_VERSION)

path_to_allen_json = "AllenMouseBrainOntology.json" braian.cache(path_to_allen_json, "http://api.brain-map.org/api/v2/structure_graph_download/1.json") brain_ontology = braian.AllenBrainOntology(path_to_allen_json, EXCLUDED_BRANCHES, version=ALLEN_VERSION)

We can also plot the ontology to make sure that we excluded the we excluded the correct branches.
It can also be handy in future to choose on which non-overlapping brain regions do the analysis!

bap.hierarchy(brain_ontology)

bap.hierarchy(brain_ontology)

self.groups = [BraiAnConfig.GroupDirectory(
                    id=int(group[len("group"):]),
                    name=self.config["experiment"][group]["name"],
                    animal_directories=self.config["experiment"][group]["dirs"]
                ) for group in self.config["experiment"]
                  if group.startswith("group") and group[len("group"):].isdigit()
              ]

self.groups = [BraiAnConfig.GroupDirectory( id=int(group[len("group"):]), name=self.config["experiment"][group]["name"], animal_directories=self.config["experiment"][group]["dirs"] ) for group in self.config["experiment"] if group.startswith("group") and group[len("group"):].isdigit() ]

config = braian.BraiAnConfig(os.path.join(rootdir_path, "data"), config_file)
selected_regions = config.brain_ontology.get_regions(config.config["comparison"]["regions-to-plot"])
print(f"You selected {len(selected_regions)} regions to plot.")

#parent_region = config.brain_ontology.parent_region
#direct_subregions = config.brain_ontology.direct_subregions
#full_name = config.brain_ontology.full_name
#regions = config.brain_ontology.list_all_subregions("root", mode="depth")

config = braian.BraiAnConfig(os.path.join(rootdir_path, "data"), config_file) selected_regions = config.brain_ontology.get_regions(config.config["comparison"]["regions-to-plot"]) print(f"You selected {len(selected_regions)} regions to plot.") #parent_region = config.brain_ontology.parent_region #direct_subregions = config.brain_ontology.direct_subregions #full_name = config.brain_ontology.full_name #regions = config.brain_ontology.list_all_subregions("root", mode="depth")

REMOVE_HIGH_CV_REGIONS = False
CVAR_THRESHOLD = 1

# ###################################### PLOT OPTIONS ######################################
PLOT_ALLENBRAIN_HIERARCHY = True
PLOT_ANIMALS_ROOTS = False
PLOT_COEFFICIENT_OF_VARIATION = True
PLOT_COEFFICIENT_OF_VARIATION_THRESHOLD = 1

SAVE_ANIMALS = True
SAVE_GROUPS = True

REMOVE_HIGH_CV_REGIONS = False CVAR_THRESHOLD = 1 # ###################################### PLOT OPTIONS ###################################### PLOT_ALLENBRAIN_HIERARCHY = True PLOT_ANIMALS_ROOTS = False PLOT_COEFFICIENT_OF_VARIATION = True PLOT_COEFFICIENT_OF_VARIATION_THRESHOLD = 1 SAVE_ANIMALS = True SAVE_GROUPS = True

Script's code

The Allen Brain Atlas

We start by importing the mouse Allen Brain Atlas, in which we find information about all brain regions (their parent region and children regions in the brain hierarchy, for example).

# Plot brain region hierarchy
if PLOT_ALLENBRAIN_HIERARCHY:
    fig = bap.hierarchy(config.brain_ontology)
    fig.show()

# Plot brain region hierarchy if PLOT_ALLENBRAIN_HIERARCHY: fig = bap.hierarchy(config.brain_ontology) fig.show()

Load data

Now, we're ready to read the .csv files with the cell counts, and also the exclusion files (if there were regions to exclude).
Below, you have to specify:

• animals_root: Absolute path to the folder that contains the animal folders.
• group_1_dirs: A list of names of the folders corresponding to animals in Group 1 (e.g., Control group). Indeed, it is necessary to store the results in individual folders for each animal.
• group_2_dirs: A list of names of the folders corresponding to animals in Group 2 (e.g., Stress group).
• group_1_name: A meaningful string for Group 1.
• group_2_name: A meaningful string for Group 2.
• area_key: A string of the column in the .csv files that refers to the size of a brain areatra
• tracer_key: A string of the column in the .csv files that refers to the tracer number used to highlight the marker
• marker: A string of the marker we would like to highlight (e.g. CFos)

Provare a modificar per ottenere densita in mm^2 (da micron)

Now, we load the Control and Stress results seperately in two pandas dataframes, and save the results.

Note: regions to exclude are automatically excluded.

config.read_groups(DATA_INPUT_PATH)
groups_slices = [[sliced_brain for sliced_brain in group.sliced_brains] for group in config.groups]

config.read_groups(DATA_INPUT_PATH) groups_slices = [[sliced_brain for sliced_brain in group.sliced_brains] for group in config.groups]

if PLOT_ANIMALS_ROOTS:
    region_name = "root"
    root_plot = bap.plot_region_density(region_name, *groups_slices, width=1000, height=500)
    root_plot.show()

if PLOT_ANIMALS_ROOTS: region_name = "root" root_plot = bap.plot_region_density(region_name, *groups_slices, width=1000, height=500) root_plot.show()

# print("N regions above threshold:", sum([(brain.data > cv_threshold).sum() for brain in cvar_brains]))
# print("N regions below threshold:", sum([(brain.data <= cv_threshold).sum() for brain in cvar_brains]))
if PLOT_COEFFICIENT_OF_VARIATION:
    cvar_plot = bap.plot_cv_above_threshold(config.brain_ontology, *groups_slices, cv_threshold=PLOT_COEFFICIENT_OF_VARIATION_THRESHOLD, width=1000, height=500)
    cvar_plot.show()

# print("N regions above threshold:", sum([(brain.data > cv_threshold).sum() for brain in cvar_brains])) # print("N regions below threshold:", sum([(brain.data <= cv_threshold).sum() for brain in cvar_brains])) if PLOT_COEFFICIENT_OF_VARIATION: cvar_plot = bap.plot_cv_above_threshold(config.brain_ontology, *groups_slices, cv_threshold=PLOT_COEFFICIENT_OF_VARIATION_THRESHOLD, width=1000, height=500) cvar_plot.show()

animal_name = "224.3CM"
region_acronym = "ENTl6a"
# for a in config.groups[-1].brains:
config.check_animal_region(animal_name, region_acronym, marker=None) #, marker="cFos")
config.check_animal_region_slices(animal_name, region_acronym)

animal_name = "224.3CM" region_acronym = "ENTl6a" # for a in config.groups[-1].brains: config.check_animal_region(animal_name, region_acronym, marker=None) #, marker="cFos") config.check_animal_region_slices(animal_name, region_acronym)

# NOTE: brains are being written WITH Left/Right discrimination
# If you desire to save them without, call AnimalBrain with hemisphere_distinction=False

config.reduce_slices()
if SAVE_ANIMALS:
    for group in config.groups:
        for animal in group.brains:
            print(f"{animal} saved to {animal.to_csv(DATA_OUTPUT_PATH, sep="\t")}")

# NOTE: brains are being written WITH Left/Right discrimination # If you desire to save them without, call AnimalBrain with hemisphere_distinction=False config.reduce_slices() if SAVE_ANIMALS: for group in config.groups: for animal in group.brains: print(f"{animal} saved to {animal.to_csv(DATA_OUTPUT_PATH, sep="\t")}")

groups_density: list[braian.AnimalGroup] = config.to_groups("density")
if SAVE_GROUPS:
    for dgroup in groups_density:
        dgroup.to_csv(DATA_OUTPUT_PATH, sep="\t", overwrite=True)

groups_density: list[braian.AnimalGroup] = config.to_groups("density") if SAVE_GROUPS: for dgroup in groups_density: dgroup.to_csv(DATA_OUTPUT_PATH, sep="\t", overwrite=True)

# compute each comparison 
for comparison in config.comparisons:
    print("comparison."+comparison.id+":", end=" ")
    comparison.apply()
    print("COMPUTED!")

# compute each comparison for comparison in config.comparisons: print("comparison."+comparison.id+":", end=" ") comparison.apply() print("COMPUTED!")

for comparison in config.comparisons:
    comparison.plot_heatmaps(PLOTS_ROOT, n=20, orientation="frontal", show_text=False, cmap="magma_r") #, title=f"{right_m} & {left_m} {str(metric)}")

for comparison in config.comparisons: comparison.plot_heatmaps(PLOTS_ROOT, n=20, orientation="frontal", show_text=False, cmap="magma_r") #, title=f"{right_m} & {left_m} {str(metric)}")

# SIMILARITY
comparison = next(comp for comp in config.comparisons if comp.metric == "Overlapping")
comparison.to_braindata()["CTX"][0].data_name, comparison.to_braindata()["CTX"][1].data_name

# SIMILARITY comparison = next(comp for comp in config.comparisons if comp.metric == "Overlapping") comparison.to_braindata()["CTX"][0].data_name, comparison.to_braindata()["CTX"][1].data_name

# SIMILARITY
comparison = next(comp for comp in config.comparisons if comp.metric == "Similarity")
comparison.braindata["cFos+Arc1"][0].data["CA1"], comparison.braindata["cFos+Arc1"][1].data["CA1"]
comparison.result[0].mean["cFos+Arc1"].data["CA1"], comparison.result[1].mean["cFos+Arc1"].data["CA1"]

# SIMILARITY comparison = next(comp for comp in config.comparisons if comp.metric == "Similarity") comparison.braindata["cFos+Arc1"][0].data["CA1"], comparison.braindata["cFos+Arc1"][1].data["CA1"] comparison.result[0].mean["cFos+Arc1"].data["CA1"], comparison.result[1].mean["cFos+Arc1"].data["CA1"]

# PLS
comparison = next(comp for comp in config.comparisons if comp.id == "6")
pls_ = comparison.to_braindata()
pls_ = {m: d[0] for m,d in pls_.items()}
marker1, marker2 = comparison.markers
braian.plot_gridgroups(comparison.result,
                       comparison.regions_to_plot,
                       marker1, marker2,
                       config.brain_onthology,
                       markers_salience_scores=pls_,
                       plot_scatter=True,
                       height=5000, width=None,
                       barplot_width=0.7, space_between_markers=0.02,
                       groups_marker1_colours=["LightCoral", "SandyBrown"],
                       groups_marker2_colours=["IndianRed", "Orange"],
                       color_heatmap="deep_r")

# PLS comparison = next(comp for comp in config.comparisons if comp.id == "6") pls_ = comparison.to_braindata() pls_ = {m: d[0] for m,d in pls_.items()} marker1, marker2 = comparison.markers braian.plot_gridgroups(comparison.result, comparison.regions_to_plot, marker1, marker2, config.brain_onthology, markers_salience_scores=pls_, plot_scatter=True, height=5000, width=None, barplot_width=0.7, space_between_markers=0.02, groups_marker1_colours=["LightCoral", "SandyBrown"], groups_marker2_colours=["IndianRed", "Orange"], color_heatmap="deep_r")

# if REMOVE_HIGH_CV_REGIONS:
#     for i, group_slices in enumerate(groups_slices):
#         for animal_brain, slices in zip(groups_sum_brains[i], group_slices):
#             cvars = braian.AnimalBrain.from_slices(slices, mode="cvar", hemisphere_distinction=animal_brain.is_split, min_slices=0)
#             # TODO: currently there is no differentiation between real markers and overlapping markers.
#             # This bad workaround excludes all those markers having a '+' in the name.
#             real_markers = [m for m in cvars.markers if "+" not in m]
#             cvars_data = cvars.to_pandas()
#             disperse_regions = cvars_data.index[(cvars_data > CVAR_THRESHOLD)[real_markers].any(axis=1)]
#             print(f"removing {len(disperse_regions)}/{len(cvars_data)} dispersive regions from '{slices.name}'")
#             animal_brain.remove_region(*disperse_regions)

# if REMOVE_HIGH_CV_REGIONS: # for i, group_slices in enumerate(groups_slices): # for animal_brain, slices in zip(groups_sum_brains[i], group_slices): # cvars = braian.AnimalBrain.from_slices(slices, mode="cvar", hemisphere_distinction=animal_brain.is_split, min_slices=0) # # TODO: currently there is no differentiation between real markers and overlapping markers. # # This bad workaround excludes all those markers having a '+' in the name. # real_markers = [m for m in cvars.markers if "+" not in m] # cvars_data = cvars.to_pandas() # disperse_regions = cvars_data.index[(cvars_data > CVAR_THRESHOLD)[real_markers].any(axis=1)] # print(f"removing {len(disperse_regions)}/{len(cvars_data)} dispersive regions from '{slices.name}'") # animal_brain.remove_region(*disperse_regions)

Experimenal code

import plotly.graph_objects as go
from plotly.subplots import make_subplots

groups_perc = braian.AnimalGroup("+".join((g.name for g in groups_sum)),
                                 [a for g in groups_sum for a in g.animals], 
                                 metric="percentage", onthology=brain_onthology, merge_hemispheres=True)
fig = make_subplots(rows=1, cols=3, specs=[[{'type':'domain'}, {'type':'domain'}, {'type':'domain'}]],
                    subplot_titles=groups_perc.markers)
for i,marker in enumerate(groups_perc.markers):
    data = groups_perc.mean[marker].data[braian.MAJOR_DIVISIONS].copy()
    data["other"] = 1-data.sum()
    # groups_perc[0].mean["Arc1"].data[braian.MAJOR_DIVISIONS].sum()
    allen_colours = brain_onthology.get_region_colors()
    allen_colours["other"] = "lightgrey"
    fig.add_trace(
        go.Pie(
            labels=[brain_onthology.full_name[r] if r in braian.MAJOR_DIVISIONS else r for r,p in data.items() if p != 0],
            values=[p for p in data if p != 0],
            name=marker,
            marker=dict(
                colors=[allen_colours[r] for r,p in data.items() if p != 0],
                line=dict(color="#000000", width=2)
            ),
            sort=False,
            textfont=dict(size=12),
            hole=0.3,
            textposition="outside", textinfo="percent+label",
            showlegend=False
        ),
        1, i+1
    )
    
fig.update_layout(
    title_text=f"Positive percentage by areas - {groups_perc.name}",
)
fig.show()

import plotly.graph_objects as go from plotly.subplots import make_subplots groups_perc = braian.AnimalGroup("+".join((g.name for g in groups_sum)), [a for g in groups_sum for a in g.animals], metric="percentage", onthology=brain_onthology, merge_hemispheres=True) fig = make_subplots(rows=1, cols=3, specs=[[{'type':'domain'}, {'type':'domain'}, {'type':'domain'}]], subplot_titles=groups_perc.markers) for i,marker in enumerate(groups_perc.markers): data = groups_perc.mean[marker].data[braian.MAJOR_DIVISIONS].copy() data["other"] = 1-data.sum() # groups_perc[0].mean["Arc1"].data[braian.MAJOR_DIVISIONS].sum() allen_colours = brain_onthology.get_region_colors() allen_colours["other"] = "lightgrey" fig.add_trace( go.Pie( labels=[brain_onthology.full_name[r] if r in braian.MAJOR_DIVISIONS else r for r,p in data.items() if p != 0], values=[p for p in data if p != 0], name=marker, marker=dict( colors=[allen_colours[r] for r,p in data.items() if p != 0], line=dict(color="#000000", width=2) ), sort=False, textfont=dict(size=12), hole=0.3, textposition="outside", textinfo="percent+label", showlegend=False ), 1, i+1 ) fig.update_layout( title_text=f"Positive percentage by areas - {groups_perc.name}", ) fig.show()

if "comparison" in config:
    comparisons = dict()
    for comp_name, comp in config["comparison"].items():
        if not isinstance(comp, dict):
            continue
        comparisons[comp_name] = []
        metric = braian.BrainMetrics(comp["metric"]) if comp["metric"].lower() != "correlation" else "correlation"
        kwargs = dict(brain_onthology=brain_onthology, merge_hemispheres=True)
        for i, group_sum in enumerate(groups_sum): 
            if "markers" in comp: #len(comp["markers"]) == 2
                markers = (group_sum.markers[i-1] for i in comp["markers"])
            else:
                markers = group_sum.markers
            match metric:
                case "correlation":
                    data = bas.markers_correlation(*markers,
                                                   group=braian.AnimalGroup(group_sum.name, group_sum.animals, "density", **kwargs))
                case braian.BrainMetrics.OVERLAPPING | braian.BrainMetrics.SIMILARITY_INDEX | braian.BrainMetrics.DENSITY_DIFFERENCE:
                    m1, m2 = markers
                    data = braian.AnimalGroup(group_sum.name, group_sum.animals, metric, marker1=m1, marker2=m2, **kwargs)
                case _:
                    data = braian.AnimalGroup(group_sum.name, group_sum.animals, metric, **kwargs)
            comparisons[comp_name].append(data)

if "comparison" in config: comparisons = dict() for comp_name, comp in config["comparison"].items(): if not isinstance(comp, dict): continue comparisons[comp_name] = [] metric = braian.BrainMetrics(comp["metric"]) if comp["metric"].lower() != "correlation" else "correlation" kwargs = dict(brain_onthology=brain_onthology, merge_hemispheres=True) for i, group_sum in enumerate(groups_sum): if "markers" in comp: #len(comp["markers"]) == 2 markers = (group_sum.markers[i-1] for i in comp["markers"]) else: markers = group_sum.markers match metric: case "correlation": data = bas.markers_correlation(*markers, group=braian.AnimalGroup(group_sum.name, group_sum.animals, "density", **kwargs)) case braian.BrainMetrics.OVERLAPPING | braian.BrainMetrics.SIMILARITY_INDEX | braian.BrainMetrics.DENSITY_DIFFERENCE: m1, m2 = markers data = braian.AnimalGroup(group_sum.name, group_sum.animals, metric, marker1=m1, marker2=m2, **kwargs) case _: data = braian.AnimalGroup(group_sum.name, group_sum.animals, metric, **kwargs) comparisons[comp_name].append(data)

def make_filename(*ss: str):
    return "_".join((s.replace(' ', '_') for s in ss if s != ""))

if "comparison" in config:
    for comp_name, comp in config["comparison"].items():
        if not isinstance(comp, dict):
            continue
        plots_output_dir = os.path.join(PLOTS_ROOT, comp["dir"])
        os.makedirs(plots_output_dir, exist_ok=True)
        metric_data = []
        if "groups" in comp:
            if len(comp["groups"]) != 2:
                print(f"WARNING: you can compare only two groups at the time!")
                continue
            group1, group2 = (comparisons[comp_name][i-1] for i in comp["groups"])
            if isinstance(group1, braian.AnimalGroup) and isinstance(group2, braian.AnimalGroup):
                if not group1.is_comparable(group2):
                    print(f"WARNING: '{group1}' is not comparable with '{group2}'!")
                    continue
                for marker in group1.markers:
                    right_data = group1.mean[marker]
                    left_data = group2.mean[marker]
                    metric_data.append((right_data, left_data, marker, f"{group1.name}+{group2.name}"))
            elif isinstance(group1, braian.BrainData) and isinstance(group2, braian.BrainData):
                metric_data.append((group1, group2, "", f"{right_data.data_name}+{left_data.data_name}"))
            else:
                print(f"WARNING: '{group1}' and '{group2}' are of an unexpected type!")
                continue
        elif "markers" in comp:
            groups = comparisons[comp_name]
            for group in groups:
                if isinstance(group, braian.BrainData):
                    right_data = group
                    left_data = None
                    metric_data.append((right_data, left_data, right_data.data_name, ""))
                    continue
                try:
                    marker1, marker2 = (group.markers[i-1] for i in comp["markers"])
                except IndexError:
                    print(f"WARNING: be sure that metric='{str(group.metric)}' is asymmetric.", end=" ")
                    print(", ".join((str(braian.BrainMetrics.SIMILARITY_INDEX), str(braian.BrainMetrics.DENSITY_DIFFERENCE), "correlation"))+" are not!")
                    right_data = group.mean[group.markers[0]]
                    left_data = None
                    metric_data.append((right_data, left_data, group.name, ""))
                    continue
                right_data = group.mean[marker1]
                right_data.data_name = marker1
                left_data = group.mean[marker2]
                left_data.data_name = marker2
                metric_data.append((right_data, left_data, group.name, f"{marker1}+{marker2}"))
        print("comparison", comp_name, "---", comp["metric"])
        for right_data, left_data, common_str, comparison_str in metric_data:
            # print("\tRIGHT:" if left_data is not None else "\tBOTH", right_data.data_name, right_data.metric, right_data.units)
            # if left_data is not None:
            #     print("\tLEFT:", left_data.data_name, left_data.metric, left_data.units)
            metric = right_data.metric.lower().split(" ")[0]
            filename = make_filename(metric, common_str, comparison_str)
            if metric.endswith("-corr") or metric.startswith(str(braian.BrainMetrics.DENSITY_DIFFERENCE)):
                centered_cmap = True
            else:
                centered_cmap = False
            if metric.endswith("-corr"):
                cmin, cmax = -1, 1
            elif metric.startswith(str(braian.BrainMetrics.SIMILARITY_INDEX)) or \
                 metric.startswith(str(braian.BrainMetrics.OVERLAPPING)):
                cmin, cmax = 0, 1
            else:
                cmin, cmax = None, None
            print(f"\t{filename}: ", end="")
            bap.heatmap(right_data, selected_regions,
                        plots_output_dir, filename, other=left_data, n=20,
                        cmin=cmin, cmax=cmax, cmap="magma_r", centered_cmap=centered_cmap,
                        orientation="frontal", show_text=False) #, title=f"{right_m} & {left_m} {str(metric)}")

def make_filename(*ss: str): return "_".join((s.replace(' ', '_') for s in ss if s != "")) if "comparison" in config: for comp_name, comp in config["comparison"].items(): if not isinstance(comp, dict): continue plots_output_dir = os.path.join(PLOTS_ROOT, comp["dir"]) os.makedirs(plots_output_dir, exist_ok=True) metric_data = [] if "groups" in comp: if len(comp["groups"]) != 2: print(f"WARNING: you can compare only two groups at the time!") continue group1, group2 = (comparisons[comp_name][i-1] for i in comp["groups"]) if isinstance(group1, braian.AnimalGroup) and isinstance(group2, braian.AnimalGroup): if not group1.is_comparable(group2): print(f"WARNING: '{group1}' is not comparable with '{group2}'!") continue for marker in group1.markers: right_data = group1.mean[marker] left_data = group2.mean[marker] metric_data.append((right_data, left_data, marker, f"{group1.name}+{group2.name}")) elif isinstance(group1, braian.BrainData) and isinstance(group2, braian.BrainData): metric_data.append((group1, group2, "", f"{right_data.data_name}+{left_data.data_name}")) else: print(f"WARNING: '{group1}' and '{group2}' are of an unexpected type!") continue elif "markers" in comp: groups = comparisons[comp_name] for group in groups: if isinstance(group, braian.BrainData): right_data = group left_data = None metric_data.append((right_data, left_data, right_data.data_name, "")) continue try: marker1, marker2 = (group.markers[i-1] for i in comp["markers"]) except IndexError: print(f"WARNING: be sure that metric='{str(group.metric)}' is asymmetric.", end=" ") print(", ".join((str(braian.BrainMetrics.SIMILARITY_INDEX), str(braian.BrainMetrics.DENSITY_DIFFERENCE), "correlation"))+" are not!") right_data = group.mean[group.markers[0]] left_data = None metric_data.append((right_data, left_data, group.name, "")) continue right_data = group.mean[marker1] right_data.data_name = marker1 left_data = group.mean[marker2] left_data.data_name = marker2 metric_data.append((right_data, left_data, group.name, f"{marker1}+{marker2}")) print("comparison", comp_name, "---", comp["metric"]) for right_data, left_data, common_str, comparison_str in metric_data: # print("\tRIGHT:" if left_data is not None else "\tBOTH", right_data.data_name, right_data.metric, right_data.units) # if left_data is not None: # print("\tLEFT:", left_data.data_name, left_data.metric, left_data.units) metric = right_data.metric.lower().split(" ")[0] filename = make_filename(metric, common_str, comparison_str) if metric.endswith("-corr") or metric.startswith(str(braian.BrainMetrics.DENSITY_DIFFERENCE)): centered_cmap = True else: centered_cmap = False if metric.endswith("-corr"): cmin, cmax = -1, 1 elif metric.startswith(str(braian.BrainMetrics.SIMILARITY_INDEX)) or \ metric.startswith(str(braian.BrainMetrics.OVERLAPPING)): cmin, cmax = 0, 1 else: cmin, cmax = None, None print(f"\t{filename}: ", end="") bap.heatmap(right_data, selected_regions, plots_output_dir, filename, other=left_data, n=20, cmin=cmin, cmax=cmax, cmap="magma_r", centered_cmap=centered_cmap, orientation="frontal", show_text=False) #, title=f"{right_m} & {left_m} {str(metric)}")

group_sum = config.groups[0] #groups_sum[0]
group_diff = braian.AnimalGroup(name=group_sum.name, animals=group_sum.brains, metric="ddiff", marker1="Arc1", marker2="cFos", merge_hemispheres=True, onthology=config.brain_onthology)
data = group_diff.mean["Arc1+cFos"]

for i in range(5,7):
    config.brain_onthology.select_at_structural_level(i)
    depth_i = config.brain_onthology.get_selected_regions()
    try:
        f = bap.heatmap(data, depth_i, None, None, depth=6682, centered_cmap=True, orientation="frontal", show_text=True, title=f"structural_level={i}")
    except ValueError:
        print(f"ERROR: structural_level={i}")

group_sum = config.groups[0] #groups_sum[0] group_diff = braian.AnimalGroup(name=group_sum.name, animals=group_sum.brains, metric="ddiff", marker1="Arc1", marker2="cFos", merge_hemispheres=True, onthology=config.brain_onthology) data = group_diff.mean["Arc1+cFos"] for i in range(5,7): config.brain_onthology.select_at_structural_level(i) depth_i = config.brain_onthology.get_selected_regions() try: f = bap.heatmap(data, depth_i, None, None, depth=6682, centered_cmap=True, orientation="frontal", show_text=True, title=f"structural_level={i}") except ValueError: print(f"ERROR: structural_level={i}")

for i in range(4,8):
    config.brain_onthology.select_at_depth(i)
    depth_i = config.brain_onthology.get_selected_regions()
    try:
        f = bap.heatmap(data, depth_i, None, None, depth=6682, centered_cmap=True, orientation="frontal", show_text=True, title=f"depth={i}")
    except ValueError:
        print(f"ERROR: depth={i}")

for i in range(4,8): config.brain_onthology.select_at_depth(i) depth_i = config.brain_onthology.get_selected_regions() try: f = bap.heatmap(data, depth_i, None, None, depth=6682, centered_cmap=True, orientation="frontal", show_text=True, title=f"depth={i}") except ValueError: print(f"ERROR: depth={i}")

marker = groups_density[0].markers[0]
pls = bas.PLS(selected_regions, groups_density[0], groups_density[1], marker=marker)
pls.random_permutation
pls.bootstrap_salience_scores(n=5000)
pls.X

marker = groups_density[0].markers[0] pls = bas.PLS(selected_regions, groups_density[0], groups_density[1], marker=marker) pls.random_permutation pls.bootstrap_salience_scores(n=5000) pls.X

def make_filename(*ss: str):
    return "_".join((s.replace(' ', '_') for s in ss if s != ""))

filename = make_filename("pls_salient_scores", "CTX-FC", "cFos+Arc1")
pls["cFos"].data_name = "cFos"
pls["Arc1"].data_name = "Arc1"
bap.heatmap(pls["cFos"], selected_regions, os.path.join(PLOTS_ROOT, "CTX-FC"), filename, other=pls["Arc1"], n=20, centered_cmap=True, orientation="frontal", show_text=True)

def make_filename(*ss: str): return "_".join((s.replace(' ', '_') for s in ss if s != "")) filename = make_filename("pls_salient_scores", "CTX-FC", "cFos+Arc1") pls["cFos"].data_name = "cFos" pls["Arc1"].data_name = "Arc1" bap.heatmap(pls["cFos"], selected_regions, os.path.join(PLOTS_ROOT, "CTX-FC"), filename, other=pls["Arc1"], n=20, centered_cmap=True, orientation="frontal", show_text=True)

import scipy.stats
import matplotlib.colors
import numpy as np

threshold = bas.PLS.to_zscore01, two_tailed=True)
print(f"THRESHOLD: ~{threshold:.2f}")
marker1 = "cFos"
marker2 = "Arc1"

marker1_salients = pls[marker1].data.abs() >= threshold
marker2_salients = pls[marker2].data.abs() >= threshold
assert (marker1_salients.index == marker2_salients.index).all(), "PLS results must be on the same brain regions!"
print(f"#salient regions ({marker1}):", sum(marker1_salients))
print(f"#salient regions ({marker2}):", sum(marker2_salients))
marker1_only_salients = marker1_salients & (~marker2_salients)
marker2_only_salients = marker2_salients & (~marker1_salients)
both_markers_salients = marker1_salients & marker2_salients

salient_regions = pd.Series(0, index=marker1_salients.index)    # 0 -> not salient for neither of the markers
salient_regions[marker1_only_salients] = 1                      # 1 -> salient only for marker1
salient_regions[marker2_only_salients] = 2                      # 2 -> salient only for marker2
salient_regions[both_markers_salients] = 3                      # 3 -> salient for both marker1 and marker2
salient_regions[np.isnan(pls[marker1].data)] = np.nan
salient_regions = braian.BrainData(salient_regions, pls[marker1].data_name, f"salient regions (|z-score|≥{threshold:.2f})", "class", fill_nan=False)

def make_filename(*ss: str):
    return "_".join((s.replace(' ', '_') for s in ss if s != ""))

bounds = np.array([-.5,.5,1.5,2.5,3.5])
#norm = matplotlib.colors.BoundaryNorm(boundaries=(bounds-bounds.min())/(bounds.max()-bounds.min()), ncolors=256) #, extend='both')
#cmap = braian.brain_data.NormalizedColormap(cmap=matplotlib.colormaps['viridis'], norm=norm)
cmap = matplotlib.colors.ListedColormap(["#ffffcc", "#a1dab4", "#ff8686", "#2c7fb8"], name="my_cmap")
# bap.heatmap(salient_regions, selected_regions, None, None, depth=6500, cmap=cmap, centered_cmap=False, orientation="frontal", show_text=True)
filename = make_filename("pls_salient_regions", "CTX-FC", "cFos+Arc1")
bap.heatmap(salient_regions, selected_regions, os.path.join(PLOTS_ROOT, "CTX-FC"), filename, n=20, cmap=cmap, centered_cmap=False, orientation="frontal", show_text=False)

import scipy.stats import matplotlib.colors import numpy as np threshold = bas.PLS.to_zscore01, two_tailed=True) print(f"THRESHOLD: ~{threshold:.2f}") marker1 = "cFos" marker2 = "Arc1" marker1_salients = pls[marker1].data.abs() >= threshold marker2_salients = pls[marker2].data.abs() >= threshold assert (marker1_salients.index == marker2_salients.index).all(), "PLS results must be on the same brain regions!" print(f"#salient regions ({marker1}):", sum(marker1_salients)) print(f"#salient regions ({marker2}):", sum(marker2_salients)) marker1_only_salients = marker1_salients & (~marker2_salients) marker2_only_salients = marker2_salients & (~marker1_salients) both_markers_salients = marker1_salients & marker2_salients salient_regions = pd.Series(0, index=marker1_salients.index) # 0 -> not salient for neither of the markers salient_regions[marker1_only_salients] = 1 # 1 -> salient only for marker1 salient_regions[marker2_only_salients] = 2 # 2 -> salient only for marker2 salient_regions[both_markers_salients] = 3 # 3 -> salient for both marker1 and marker2 salient_regions[np.isnan(pls[marker1].data)] = np.nan salient_regions = braian.BrainData(salient_regions, pls[marker1].data_name, f"salient regions (|z-score|≥{threshold:.2f})", "class", fill_nan=False) def make_filename(*ss: str): return "_".join((s.replace(' ', '_') for s in ss if s != "")) bounds = np.array([-.5,.5,1.5,2.5,3.5]) #norm = matplotlib.colors.BoundaryNorm(boundaries=(bounds-bounds.min())/(bounds.max()-bounds.min()), ncolors=256) #, extend='both') #cmap = braian.brain_data.NormalizedColormap(cmap=matplotlib.colormaps['viridis'], norm=norm) cmap = matplotlib.colors.ListedColormap(["#ffffcc", "#a1dab4", "#ff8686", "#2c7fb8"], name="my_cmap") # bap.heatmap(salient_regions, selected_regions, None, None, depth=6500, cmap=cmap, centered_cmap=False, orientation="frontal", show_text=True) filename = make_filename("pls_salient_regions", "CTX-FC", "cFos+Arc1") bap.heatmap(salient_regions, selected_regions, os.path.join(PLOTS_ROOT, "CTX-FC"), filename, n=20, cmap=cmap, centered_cmap=False, orientation="frontal", show_text=False)

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats

np.random.seed(42)  # Set the random seed for reproducibility
data = np.random.normal(0, 1, 100_000)
mean = np.mean(data)
std = np.std(data)
print("(μ+2σ) percentile:", 1-scipy.stats.norm.pdf((mean + 2*std)))

# Plot normal
x = np.linspace(scipy.stats.norm.ppf(0.01), scipy.stats.norm.ppf(0.99), 100)
plt.plot(x, scipy.stats.norm.pdf(x), 'r-', lw=5, alpha=0.6, label='norm pdf')

# Plot histogram
plt.hist(data, bins=30, density=True, alpha=0.7)

# Plot the mean and standard deviations
plt.axvline(mean, color='r', linestyle='dashed', linewidth=1, label='Mean')
plt.axvline(mean - std, color='g', linestyle='dashed', linewidth=1, label='1 STD')
plt.axvline(mean + std, color='g', linestyle='dashed', linewidth=1)
plt.axvline(mean - 2*std, color='b', linestyle='dashed', linewidth=1, label='2 STD')
plt.axvline(mean + 2*std, color='b', linestyle='dashed', linewidth=1)
plt.axvline(mean - 3*std, color='m', linestyle='dashed', linewidth=1, label='3 STD')
plt.axvline(mean + 3*std, color='m', linestyle='dashed', linewidth=1)

plt.legend()
plt.xlabel('Value')
plt.ylabel('Density')
plt.title('Histogram of the Dataset')
plt.show()

import numpy as np import matplotlib.pyplot as plt import scipy.stats np.random.seed(42) # Set the random seed for reproducibility data = np.random.normal(0, 1, 100_000) mean = np.mean(data) std = np.std(data) print("(μ+2σ) percentile:", 1-scipy.stats.norm.pdf((mean + 2*std))) # Plot normal x = np.linspace(scipy.stats.norm.ppf(0.01), scipy.stats.norm.ppf(0.99), 100) plt.plot(x, scipy.stats.norm.pdf(x), 'r-', lw=5, alpha=0.6, label='norm pdf') # Plot histogram plt.hist(data, bins=30, density=True, alpha=0.7) # Plot the mean and standard deviations plt.axvline(mean, color='r', linestyle='dashed', linewidth=1, label='Mean') plt.axvline(mean - std, color='g', linestyle='dashed', linewidth=1, label='1 STD') plt.axvline(mean + std, color='g', linestyle='dashed', linewidth=1) plt.axvline(mean - 2*std, color='b', linestyle='dashed', linewidth=1, label='2 STD') plt.axvline(mean + 2*std, color='b', linestyle='dashed', linewidth=1) plt.axvline(mean - 3*std, color='m', linestyle='dashed', linewidth=1, label='3 STD') plt.axvline(mean + 3*std, color='m', linestyle='dashed', linewidth=1) plt.legend() plt.xlabel('Value') plt.ylabel('Density') plt.title('Histogram of the Dataset') plt.show()

f = bap.heatmap(bas.pls_regions_salience(groups_density[0], groups_density[1], selected_regions, n_bootstrap=5000, marker="cFos"),
                selected_regions, None, None, depth=6500, centered_cmap=True, orientation="frontal", show_text=True)

f = bap.heatmap(bas.pls_regions_salience(groups_density[0], groups_density[1], selected_regions, n_bootstrap=5000, marker="cFos"), selected_regions, None, None, depth=6500, centered_cmap=True, orientation="frontal", show_text=True)

import bgheatmaps as bgh
import vedo as vd
# import vtk
vd.embedWindow(None)
# vd.embedWindow("k3d")
from braian.brain_data import CenteredColormap

group_sum = groups_sum[0]
group_diff = braian.AnimalGroup(name=group_sum.name, animals=group_sum.animals, metric="diff", marker1="Arc1", marker2="cFos", merge_hemispheres=True, brain_onthology=brain_onthology)
data = group_diff.mean["Arc1+cFos"]
# data = data.data[~data.data.isna()]
data = data[braian.MAJOR_DIVISIONS+selected_regions]

bgh.Heatmap(
    data.to_dict(),
    position=6658,
    thickness=1000,
    orientation="frontal",
    title="difference",
    cmap=CenteredColormap("RdBu", data.min(), data.max()),
    vmin=data.min(),
    vmax=data.max(),
    format="3D",
    hemisphere="right"
).show()

import bgheatmaps as bgh import vedo as vd # import vtk vd.embedWindow(None) # vd.embedWindow("k3d") from braian.brain_data import CenteredColormap group_sum = groups_sum[0] group_diff = braian.AnimalGroup(name=group_sum.name, animals=group_sum.animals, metric="diff", marker1="Arc1", marker2="cFos", merge_hemispheres=True, brain_onthology=brain_onthology) data = group_diff.mean["Arc1+cFos"] # data = data.data[~data.data.isna()] data = data[braian.MAJOR_DIVISIONS+selected_regions] bgh.Heatmap( data.to_dict(), position=6658, thickness=1000, orientation="frontal", title="difference", cmap=CenteredColormap("RdBu", data.min(), data.max()), vmin=data.min(), vmax=data.max(), format="3D", hemisphere="right" ).show()

import bgheatmaps as bgh
import math
import vedo as vd
# vd.embedWindow("k3d")
vd.embedWindow(None)

brain_data = groups_density[0].mean[groups_density[0].markers[0]].data
brain_data = brain_data[brain_data.index.isin(selected_regions)]

heatmap = bgh.Heatmap(
            brain_data[~brain_data.isna() | (brain_data.index == "RSPd4")].to_dict(),
            position=7545,
            orientation="frontal",
            thickness=1000,
            title=f"{groups_density[0].markers[0]} density",
            cmap="magma_r",
            vmin=math.floor(brain_data[brain_data != np.inf].min(axis=None)),
            vmax=math.ceil(brain_data[brain_data != np.inf].max(axis=None)),
            format="2D",
            hemisphere="both"
        )
heatmap.show()

import bgheatmaps as bgh import math import vedo as vd # vd.embedWindow("k3d") vd.embedWindow(None) brain_data = groups_density[0].mean[groups_density[0].markers[0]].data brain_data = brain_data[brain_data.index.isin(selected_regions)] heatmap = bgh.Heatmap( brain_data[~brain_data.isna() | (brain_data.index == "RSPd4")].to_dict(), position=7545, orientation="frontal", thickness=1000, title=f"{groups_density[0].markers[0]} density", cmap="magma_r", vmin=math.floor(brain_data[brain_data != np.inf].min(axis=None)), vmax=math.ceil(brain_data[brain_data != np.inf].max(axis=None)), format="2D", hemisphere="both" ) heatmap.show()

import importlib
import sys
__imported_modules = sys.modules.copy()
for module_name, module in __imported_modules.items():
    if not module_name.startswith("braian"): # and not module_name.startswith("bgheatmaps"):
        continue
    try:
        # print("reaload:", module_name)
        importlib.reload(module)
    except ModuleNotFoundError:
        continue

import importlib import sys __imported_modules = sys.modules.copy() for module_name, module in __imported_modules.items(): if not module_name.startswith("braian"): # and not module_name.startswith("bgheatmaps"): continue try: # print("reaload:", module_name) importlib.reload(module) except ModuleNotFoundError: continue