Microns example
Note: while these packages are not required for the use of MorphLink, they are used here to demonstrate the functionality of the package with real data and integration with multiple other packages. This tutorial requires:
CAVEclient,skeletor,pcg_skel,pyvista,
all of which are available on PyPI.
Initialization¶
Start a CAVEclient
from caveclient import CAVEclient
client = CAVEclient("minnie65_public")
Create a MorphLink instance. This is the object we'll use for keeping track of
different features and representations of morphology.
from morphlink import MorphLink
morphology = MorphLink()
Adding a mesh¶
We'll start by getting the morphology of a single neuron as a mesh.
root_id = 864691135361314119
cv = client.info.segmentation_cloudvolume(progress=False)
mesh = cv.mesh.get(
root_id, deduplicate_chunk_boundaries=False, remove_duplicate_vertices=False
)[root_id]
We can add this to our morphology representation using the add_mesh method.
morphology.add_mesh(mesh, "mesh")
morphology.layers
Adding point annotations¶
Next, we can add the nucleus location for this neuron. This is a new type of layer, in this case just a single point.
nuc_info = client.materialize.query_table(
"nucleus_detection_v0",
filter_equal_dict={"pt_root_id": root_id},
split_positions=True,
desired_resolution=[1, 1, 1],
)
nuc_loc = nuc_info[["pt_position_x", "pt_position_y", "pt_position_z"]].values.squeeze()
nuc_loc
morphology.add_points(nuc_loc, "nucleus")
morphology.layers
polydatas = morphology.to_pyvista()
import pyvista as pv
def set_camera(plotter, zoom=10):
plotter.camera.focal_point = nuc_loc
plotter.camera_position = "zy"
plotter.camera.up = (0, -1, 0)
plotter.camera.focal_point = nuc_loc
plotter.camera.zoom(zoom)
WINDOW_SIZE = [2400, 1500]
pv.set_jupyter_backend("static")
plotter = pv.Plotter(window_size=WINDOW_SIZE)
plotter.add_mesh(polydatas["mesh"], color="lightgrey", opacity=0.5)
plotter.add_points(
polydatas["nucleus"], color="red", point_size=60, render_points_as_spheres=True
)
set_camera(plotter, zoom=5)
plotter.show()
Now, let's link these objects together. The simplest mapping is to
annotate the nucleus location on the mesh as the closest point.
Under the hood, the add_link method will find the closest point on the mesh and save
that mapping.
morphology.add_link("nucleus", "mesh", mapping="closest")
morphology.links
We can retreive the mapping between the nucleus and the mesh using the
get_link method. This returns a DataFrame with the mapping as its two columns.
morphology.get_link("nucleus", "mesh")
We can also ask for the specific mapping for a point in our nucleus layer. Since we only have one point in this layer, we just get one item back from the mapping, denoting the closest point on the mesh.
mesh_nuc_index = morphology.get_mapping("nucleus", "mesh")
mesh_nuc_index
Similarly, we can add another layer that consists of many points - in this case, synapses onto this neuron.
post_synapses = client.materialize.query_table(
"synapses_pni_2",
filter_equal_dict={"post_pt_root_id": root_id},
split_positions=True,
desired_resolution=[1, 1, 1],
)
post_synapses.set_index("id", inplace=True)
spatial_cols = ["ctr_pt_position_x", "ctr_pt_position_y", "ctr_pt_position_z"]
post_synapse_locs = post_synapses[spatial_cols]
post_synapse_info = post_synapses.drop(columns=spatial_cols)
Specifically, we'll specify the synapse locations using a points layer, which allows future
methods to do space-aware operations on this layer.
morphology.add_points(
post_synapse_locs,
"post_synapses",
)
Additionally, we'll add the rest of the synapse info using the add_table method, which
will add a space-independent table of data.
morphology.add_table(
post_synapse_info,
"post_synapse_info",
)
Finally, to make sure that the above information about the synapses gets mapped to
the synapse locations, we can use the add_link method again.
morphology.add_link(
"post_synapse_info", "post_synapses", mapping="index", reciprocal=True
)
plotter = pv.Plotter(window_size=WINDOW_SIZE)
plotter.add_mesh(polydatas["mesh"], color="lightgrey")
plotter.add_points(
morphology.post_synapses.to_pyvista(),
color="blue",
point_size=10,
render_points_as_spheres=True,
)
set_camera(plotter, zoom=5)
plotter.show()
Let's add another link, here again finding the closest point on the mesh for each synapse.
morphology.add_link("post_synapses", "mesh")
morphology.links
It can be helpful to visualize these mappings. Conceptually, it is as if we have created
edges from our points in the synapse layer to the corresponding closest points in the
mesh layer. We include a get_link_as_layer method to visualize these edges (note that
this does not actually add a new layer to the MorphLink object).
plotter = pv.Plotter(window_size=WINDOW_SIZE)
plotter.add_mesh(polydatas["mesh"], color="lightgrey", opacity=0.1)
plotter.add_points(
morphology.post_synapses.to_pyvista(),
color="blue",
point_size=10,
render_points_as_spheres=True,
)
plotter.add_mesh(
morphology.get_link_as_layer("post_synapses", "mesh").to_pyvista(),
color="black",
show_vertices=True,
line_width=3,
point_size=5,
)
set_camera(plotter, zoom=80)
plotter.show()
Adding a skeleton¶
Now, we can do something a bit more interesting. Let's skeletonize the mesh, using the nucleus as the source points for the skeletonization.
Note: requires skeletor package to be installed.
import time
from skeletor.skeletonize import by_wavefront
currtime = time.time()
out = by_wavefront(mesh, origins=mesh_nuc_index.to_list(), progress=True)
print(f"{time.time() - currtime:.3f} seconds elapsed.")
This skeletonization process stores the mapping between mesh vertices and the new, collapsed vertices from the skeletonization. First, let's add the skeleton to our morphology representations.
morphology.add_graph(out, "skeleton")
plotter = pv.Plotter(window_size=WINDOW_SIZE)
plotter.add_mesh(morphology.mesh.to_pyvista(), color="lightgrey", opacity=0.3)
plotter.add_mesh(
morphology.post_synapses.to_pyvista(),
color="blue",
point_size=10,
render_points_as_spheres=True,
)
plotter.add_mesh(morphology.skeleton.to_pyvista(), color="red", line_width=2)
set_camera(plotter, zoom=5)
plotter.show()
Then, we can add the mapping between the mesh and the skeleton. This mapping is
stored in the mesh_map attribute of the skeletonization output.
morphology.add_link("mesh", "skeleton", mapping=out.mesh_map)
morphology.links
More complex mappings¶
Now, we might also be interested in where synapses are located along the skeleton. Even though the skeletonization map doesn't have a direct mapping between the mesh and the synapses, we can first map synapses to their points on the mesh, and then map those mesh points to their points on the skeleton.
Fortunately, this kind of transitive mapping is handled automatically by MorphLink
under the hood. Internally, there is a graph that denotes relationships between
different layers, and the link_path method can be used to find the path from a
source layer to a target layer (if one exists).
Many times, this will just be a direct link.
morphology.get_link_path("post_synapses", "mesh")
But as in the case of the synapses and the skeleton, it will find the path that involves mapping synapses to the mesh, and then the mesh to the skeleton.
morphology.get_link_path("post_synapses", "skeleton")
synapse_skeleton_ids = morphology.get_mapping("post_synapses", "skeleton")
synapse_skeleton_ids
skeleton_synapse_points = morphology.skeleton.nodes.iloc[synapse_skeleton_ids]
plotter = pv.Plotter(window_size=WINDOW_SIZE)
plotter.add_mesh(morphology.skeleton.to_pyvista(), color="red", line_width=2)
plotter.add_points(
skeleton_synapse_points.values,
color="blue",
point_size=5,
render_points_as_spheres=True,
)
set_camera(plotter, zoom=5)
plotter.show()
Working with a level2 graph¶
import numpy as np
edgelist = client.chunkedgraph.level2_chunk_graph(root_id)
edgelist = np.array(edgelist)
level2_ids = np.unique(edgelist)
l2_data = client.l2cache.get_l2data(
level2_ids,
attributes=["area_nm2", "max_dt_nm", "mean_dt_nm", "size_nm3", "rep_coord_nm"],
)
import pandas as pd
l2_nodes = pd.DataFrame(l2_data).T
l2_nodes.index = l2_nodes.index.astype(int)
l2_nodes["x"] = l2_nodes["rep_coord_nm"].apply(lambda x: x[0])
l2_nodes["y"] = l2_nodes["rep_coord_nm"].apply(lambda x: x[1])
l2_nodes["z"] = l2_nodes["rep_coord_nm"].apply(lambda x: x[2])
l2_node_locs = l2_nodes[["x", "y", "z"]]
l2_node_info = l2_nodes.drop(columns=["x", "y", "z", "rep_coord_nm"])
morphology.add_graph((l2_node_locs, edgelist), "l2_graph")
morphology.add_table(l2_node_info, "l2_info")
morphology.add_link("l2_info", "l2_graph", mapping="index", reciprocal=True)
morphology.layers
Skeletonizing the level2 graph can also provide a quick way to get a skeleton, albeit at the loss of some spatial resolution.
from pcg_skel import pcg_skeleton_direct
Note that the edges were supplied as named index pairs, not on positional node indices.
morphology.l2_graph.edges
But edges_positional will return the edges with positional node indices.
morphology.l2_graph.edges_positional
pcg_skeleton = pcg_skeleton_direct(
morphology.l2_graph.vertices,
morphology.l2_graph.edges_positional,
collapse_soma=True,
root_point=np.squeeze(morphology.nucleus.vertices),
)
morphology.add_graph(pcg_skeleton, "pcg_skeleton")
plotter = pv.Plotter(window_size=WINDOW_SIZE)
# pv.set_jupyter_backend('client')
# plotter = pv.Plotter()
plotter.add_mesh(morphology.mesh.to_pyvista(), color="lightgrey", opacity=0.3)
plotter.add_mesh(morphology.l2_graph.to_pyvista(), color="lime", line_width=3)
plotter.add_mesh(morphology.pcg_skeleton.to_pyvista(), color="darkred", line_width=3)
set_camera(plotter, zoom=15)
plotter.show()
Again, let's add some links, including the one created by the skeletonization process. The rest will just be again based on nearest point mappings, for now.
morphology.add_link("l2_graph", "pcg_skeleton", mapping=pcg_skeleton.mesh_to_skel_map)
morphology.add_link("nucleus", "pcg_skeleton", mapping="closest")
# these ones in particular are not perfect and could be refined with CAVEclient calls,
# but let's see how it does
morphology.add_link("mesh", "pcg_skeleton", mapping="closest")
morphology.add_link("l2_graph", "mesh", mapping="closest")
morphology.links
Masking and filtering¶
A common operation one might want to do on a neuron is to mask out the axon, or extract
just the axon, etc. morphlink aims to be agnostic to how you want to compute something
like "where is the axon on this neuron", and is only here to help you keep track of how
that labeling implicitly applies to other parts of a morphology.
As a concrete example, computing this kind of split of axon/dendrite can be done using the skeleton and its synapses.
pre_synapses = client.materialize.query_table(
"synapses_pni_2",
filter_equal_dict={"pre_pt_root_id": root_id},
split_positions=True,
desired_resolution=[1, 1, 1],
)
pre_synapses = pre_synapses.query("pre_pt_root_id != post_pt_root_id")
pre_synapses.set_index("id", inplace=True)
morphology.add_points(
pre_synapses[spatial_cols],
"pre_synapses",
)
morphology.add_table(pre_synapses.drop(columns=spatial_cols), "pre_synapse_info")
morphology.add_link("pre_synapse_info", "pre_synapses", mapping="index")
morphology.add_link("pre_synapses", "mesh")
morphology.links
morphology.get_mapping("post_synapses", "pcg_skeleton")
morphology.get_mapping("pre_synapses", "pcg_skeleton")
from scipy.sparse.linalg import eigsh
from scipy.sparse import diags_array
adj = morphology.pcg_skeleton.to_adjacency()
adj = adj + adj.T
degrees = adj.sum(axis=0)
degree_matrix = diags_array([degrees], offsets=[0], shape=adj.shape)
laplacian = degree_matrix - adj
eigvals, eigvecs = eigsh(laplacian, k=5, sigma=-1e-12)
post_inds = morphology.get_mapping("post_synapses", "pcg_skeleton")
pre_inds = morphology.get_mapping("pre_synapses", "pcg_skeleton")
X_post = eigvecs[post_inds]
X_pre = eigvecs[pre_inds]
y_post = np.zeros(len(post_inds))
y_pre = np.ones(len(pre_inds))
X = np.concatenate([X_post, X_pre])
y = np.concatenate([y_post, y_pre])
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
cols = np.arange(1, eigvecs.shape[1])
lda = LinearDiscriminantAnalysis()
lda.fit(X[:, cols], y)
pred_labels = lda.predict(eigvecs[:, cols])
pcg_skeleton_info = pd.DataFrame(
index=morphology.pcg_skeleton.nodes.index, data=pred_labels, columns=["pred_labels"]
)
morphology.add_table(pcg_skeleton_info, "pcg_skeleton_info")
morphology.add_link(
"pcg_skeleton_info", "pcg_skeleton", mapping="index", reciprocal=True
)
plotter = pv.Plotter(window_size=WINDOW_SIZE)
plotter.add_mesh(morphology.mesh.to_pyvista(), color="lightgrey", opacity=0.3)
plotter.add_mesh(
morphology.pcg_skeleton.to_pyvista(),
color="darkred",
scalars=morphology.pcg_skeleton_info.nodes["pred_labels"].values,
cmap="coolwarm",
line_width=3,
show_scalar_bar=False,
)
plotter.add_points(
morphology.post_synapses.to_pyvista(),
color="blue",
point_size=5,
render_points_as_spheres=True,
)
plotter.add_points(
morphology.pre_synapses.to_pyvista(),
color="red",
point_size=5,
render_points_as_spheres=True,
)
set_camera(plotter, zoom=4)
plotter.show()
Now, the more interesting part: independent of whatever algorithm we chose for doing this split for us, let's see how to mask the rest of the neuron and all of its layers based on this split. Even though the split will be defined in terms of the pcg_skeleton_info layer, the masking will be applied to all layers, and will leverage the direct and transitive mappings we've set up so far.
To do this, we'll use the query_nodes method, and specify the layer we want to mask.
morphology.drop_layer("skeleton") # dropping because no path to pcg_skeleton
masked_morphology = morphology.query_nodes("pred_labels == 0", "pcg_skeleton_info")
masked_morphology
plotter = pv.Plotter(window_size=WINDOW_SIZE, shape=(1, 2))
plotter.subplot(0, 0)
plotter.add_mesh(morphology.mesh.to_pyvista(), color="lightgrey", opacity=0.3)
plotter.add_mesh(
morphology.pcg_skeleton.to_pyvista(),
scalars=morphology.pcg_skeleton_info.nodes["pred_labels"],
cmap="coolwarm",
line_width=2,
show_scalar_bar=False,
)
plotter.add_points(
morphology.pre_synapses.to_pyvista(),
color="red",
point_size=5,
render_points_as_spheres=True,
)
plotter.add_points(
morphology.post_synapses.to_pyvista(),
color="blue",
point_size=5,
render_points_as_spheres=True,
)
plotter.subplot(0, 1)
plotter.add_mesh(masked_morphology.mesh.to_pyvista(), color="lightgrey", opacity=0.3)
plotter.add_mesh(
masked_morphology.pcg_skeleton.to_pyvista(),
scalars=masked_morphology.pcg_skeleton_info.nodes["pred_labels"],
cmap="coolwarm",
line_width=2,
show_scalar_bar=False,
)
plotter.add_points(
masked_morphology.pre_synapses.to_pyvista(),
color="red",
point_size=5,
render_points_as_spheres=True,
)
plotter.add_points(
masked_morphology.post_synapses.to_pyvista(),
color="blue",
point_size=5,
render_points_as_spheres=True,
)
plotter.link_views()
set_camera(plotter, zoom=2)
plotter.show()
As a sanity check, note that now very few presynaptic sites are left in our morphology, as expected since we masked out the axon.
masked_morphology.pre_synapses.nodes.shape[0]
We can also see how the volume of our morphology changed after masking.
pre_volume = morphology.l2_info.nodes["size_nm3"].sum()
post_volume = masked_morphology.l2_info.nodes["size_nm3"].sum()
print(pre_volume)
print(post_volume)
print(post_volume / pre_volume)