WIP: QC dashboard

Author: CBroz1

element_array_ephys/ephys_no_curation.py

@@ -34,7 +34,6 @@ def activate(
         ephys_schema_name (str): A string containing the name of the ephys schema.
         probe_schema_name (str): A string containing the name of the probe scehma.
         create_schema (bool): If True, schema will be created in the database.
-        create_tables (bool): If True, tables related to the schema will be created in the database.
         linking_module (str): A string containing the module name or module containing the required dependencies to activate the schema.
 
     Dependencies:

element_array_ephys/ephys_report.py

@@ -2,7 +2,8 @@
 import datetime
 import datajoint as dj
 import typing as T
-import json
+
+from element_interface.utils import insert1_skip_full_duplicates
 
 schema = dj.schema()
 
@@ -136,6 +137,89 @@ def make(self, key):
         )
 
 
+@schema
+class QualityMetricCutoffs(dj.Lookup):
+    definition = """
+    cutoffs_id                    : smallint
+    ---
+    amplitude_cutoff_maximum=null : float # Defualt null, no cutoff applied
+    presence_ratio_minimum=null   : float # Defualt null, no cutoff applied
+    isi_violations_maximum=null   : float # Defualt null, no cutoff applied
+    """
+
+    contents = [(0, None, None, None), (1, 0.1, 0.9, 0.5)]
+
+    def insert_new_cutoffs(
+        cls,
+        cutoffs_id: int,
+        amplitude_cutoff_maximum: float,
+        presence_ratio_minimum: float,
+        isi_violations_maximum: float,
+    ):
+        insert1_skip_full_duplicates(  # depends on element-interface/pull/43
+            cls,
+            dict(
+                cutoffs_id=cutoffs_id,
+                amplitude_cutoff_maximum=amplitude_cutoff_maximum,
+                presence_ratio_minimum=presence_ratio_minimum,
+                isi_violations_maximum=isi_violations_maximum,
+            ),
+        )
+
+
+@schema
+class QualityMetricSet(dj.Manual):
+    definition = """
+    -> ephys.QualityMetrics
+    -> QualityMetricCutoffs
+    """
+
+
+@schema
+class QualityMetricReport(dj.Computed):
+    definition = """
+    -> QualityMetricSet
+    """
+
+    class Cluster(dj.Part):
+        definition = """
+        -> master
+        ---
+        firing_rate_plot    : longblob
+        presence_ratio_plot : longblob
+        amp_cutoff_plot     : longblob
+        isi_violation_plot  : longblob
+        snr_plot            : longblob
+        iso_dist_plot       : longblob
+        d_prime_plot        : longblob
+        nn_hit_plot         : longblob
+        """
+
+    def make(self, key):
+        from .plotting.qc import QualityMetricFigs
+
+        cutoffs = (QualityMetricCutoffs & key).fetch1()
+        qc_key = ephys.QualityMetrics & key
+        qc_figs = QualityMetricFigs(qc_key, **cutoffs)
+
+        # ADD SAVEFIG?
+
+        self.insert1(key)
+        self.Cluster.insert1(
+            dict(
+                **key,
+                firing_rate_plot=qc_figs.firing_rate_plot().to_plotly_json(),
+                presence_ratio_plot=qc_figs.presence_ratio_plot().to_plotly_json(),
+                amp_cutoff_plot=qc_figs.amp_cutoff_plot().to_plotly_json(),
+                isi_violation_plot=qc_figs.isi_violation_plot().to_plotly_json(),
+                snr_plot=qc_figs.snr_plot().to_plotly_json(),
+                iso_dist_plot=qc_figs.iso_dist_plot().to_plotly_json(),
+                d_prime_plot=qc_figs.d_prime_plot().to_plotly_json(),
+                nn_hit_plot=qc_figs.nn_hit_plot().to_plotly_json(),
+            )
+        )
+
+
 def _make_save_dir(root_dir: pathlib.Path = None) -> pathlib.Path:
     if root_dir is None:
         root_dir = pathlib.Path().absolute()

element_array_ephys/plotting/qc.py

@@ -0,0 +1,131 @@
+import numpy as np
+import pandas as pd
+import plotly.graph_objs as go
+from scipy.ndimage import gaussian_filter1d
+from .. import ephys_no_curation as ephys
+
+
+class QualityMetricFigs(object):
+    def __init__(self, key=None, **kwargs) -> None:
+        self._key = key
+        self._amplitude_cutoff_max = kwargs.get("amplitude_cutoff_maximum", None)
+        self._presence_ratio_min = kwargs.get("presence_ratio_minimum", None)
+        self._isi_violations_max = kwargs.get("isi_violations_maximum", None)
+        self._units = pd.DataFrame()
+
+    @property
+    def units(self):
+        assert self._key, "Must use key when retrieving units for QC figures"
+        if self._units.empty:
+            restrictions = ["TRUE"]
+            if self._amplitude_cutoff_max:
+                restrictions.append(f"amplitude_cutoff < {self._amplitude_cutoff_max}")
+            if self._presence_ratio_min:
+                restrictions.append(f"presence_ratio > {self._presence_ratio_min}")
+            if self._isi_violations_max:
+                restrictions.append(f"isi_violation < {self._isi_violations_max}")
+            " AND ".join(restrictions)
+            return (
+                ephys.QualityMetrics
+                * ephys.QualityMetrics.Cluster
+                * ephys.QualityMetrics.Waveform
+                & self._key
+                & restrictions
+            ).fetch(format="frame")
+        return self._units
+
+    def _plot_metric(
+        self,
+        data,
+        bins,
+        x_axis_label,
+        scale=1,
+        vline=None,
+    ):
+        fig = go.Figure()
+        fig.update_layout(
+            xaxis_title=x_axis_label,
+            template="plotly_dark",  # "simple_white",
+            width=350 * scale,
+            height=350 * scale,
+            margin=dict(l=20 * scale, r=20 * scale, t=20 * scale, b=20 * scale),
+            xaxis=dict(showgrid=False, zeroline=False, linewidth=2, ticks="outside"),
+            yaxis=dict(showgrid=False, linewidth=0, zeroline=True, visible=False),
+        )
+        if data.isnull().all():
+            return fig.add_annotation(text="No data available", showarrow=False)
+
+        histogram, histogram_bins = np.histogram(data, bins=bins, density=True)
+
+        fig.add_trace(
+            go.Scatter(
+                x=histogram_bins[:-1],
+                y=gaussian_filter1d(histogram, 1),  # TODO: remove smoothing
+                mode="lines",
+                line=dict(color="rgb(0, 160, 223)", width=2 * scale),  # DataJoint Blue
+                hovertemplate="%{x:.2f}<br>%{y:.2f}<extra></extra>",
+            )
+        )
+
+        if vline:
+            fig.add_vline(x=vline, line_width=2 * scale, line_dash="dash")
+
+        return fig
+
+    def empty_fig(self):  # TODO: Remove before submission?
+        return self._plot_metric(
+            pd.Series(["nan"]), np.linspace(0, 0, 0), "This fig left blank"
+        )
+
+    def firing_rate_plot(self):
+        return self._plot_metric(
+            np.log10(self.units["firing_rate"]),
+            np.linspace(-3, 2, 100),  # If linear, use np.linspace(0, 50, 100)
+            "log<sub>10</sub> firing rate (Hz)",
+        )
+
+    def presence_ratio_plot(self):
+        return self._plot_metric(
+            self.units["presence_ratio"],
+            np.linspace(0, 1, 100),
+            "Presence ratio",
+            vline=0.9,
+        )
+
+    def amp_cutoff_plot(self):
+        return self._plot_metric(
+            self.units["amplitude_cutoff"],
+            np.linspace(0, 0.5, 200),
+            "Amplitude cutoff",
+            vline=0.1,
+        )
+
+    def isi_violation_plot(self):
+        return self._plot_metric(
+            np.log10(self.units["isi_violation"] + 1e-5),  # Offset b/c log(0)
+            np.linspace(-6, 2.5, 100),  # If linear np.linspace(0, 10, 200)
+            "log<sub>10</sub> ISI violations",
+            vline=np.log10(0.5),
+        )
+
+    def snr_plot(self):
+        return self._plot_metric(self.units["snr"], np.linspace(0, 10, 100), "SNR")
+
+    def iso_dist_plot(self):
+        return self._plot_metric(
+            self.units["isolation_distance"],
+            np.linspace(0, 170, 50),
+            "Isolation distance",
+        )
+
+    def d_prime_plot(self):
+        return self._plot_metric(
+            self.units["d_prime"], np.linspace(0, 15, 50), "d-prime"
+        )
+
+    def nn_hit_plot(self):
+        return self._plot_metric(
+            self.units["nn_hit_rate"],
+            np.linspace(0, 1, 100),
+            "Nearest-neighbors hit rate",
+        )

element_array_ephys/plotting/unit_level.py

@@ -26,7 +26,7 @@ def plot_waveform(waveform: np.ndarray, sampling_rate: float) -> go.Figure:
             x=waveform_df["timestamp"],
             y=waveform_df["waveform"],
             mode="lines",
-            line=dict(color="rgb(51, 76.5, 204)", width=2),
+            line=dict(color="rgb(0, 160, 223)", width=2),  # DataJoint Blue
             hovertemplate="%{y:.2f} μV<br>" + "%{x:.2f} ms<extra></extra>",
         )
     )

element_array_ephys/plotting/widget.py

@@ -1,11 +1,15 @@
+import pathlib
+from skimage import io
 from modulefinder import Module
 import ipywidgets as widgets
-import pathlib
+from ipywidgets import widgets as wg
 from IPython.display import display
-from .. import ephys_report
+from plotly.io import from_json
+from plotly.subplots import make_subplots
 import plotly.graph_objs as go
 import plotly.express as px
-from skimage import io
+
+from .. import ephys_report
 
 
 def main(ephys: Module) -> widgets:
@@ -123,3 +127,125 @@ def plot_unit_widget(unit):
     unit_widget = widgets.interactive(plot_unit_widget, unit=unit_dropdown_wg)
 
     return widgets.VBox([probe_widget, unit_widget])
+
+
+def qc_widget(ephys: Module) -> widgets:
+    from .qc import QualityMetricFigs
+
+    title_button = wg.Button(
+        description="Ephys Quality Control Metrics",
+        button_style="info",
+        layout=wg.Layout(
+            height="auto", width="auto", grid_area="title_button", border="solid"
+        ),
+        style=wg.ButtonStyle(button_color="#00a0df"),
+        disabled=True,
+    )
+
+    cluster_dropdown = wg.Dropdown(
+        options=ephys.QualityMetrics.fetch("KEY"),
+        description="Clusters:",
+        description_tooltip='Press "Load" to visualize the clusters identified.',
+        disabled=False,
+        layout=wg.Layout(
+            width="95%",
+            display="flex",
+            flex_flow="row",
+            justify_content="space-between",
+            grid_area="cluster_dropdown",
+        ),
+        style={"description_width": "80px"},
+    )
+
+    cutoff_dropdown = wg.Dropdown(
+        options=ephys_report.QualityMetricCutoffs.fetch("KEY"),
+        description="Cutoffs:",
+        description_tooltip='Press "Load" to visualize the clusters identified.',
+        disabled=False,
+        layout=wg.Layout(
+            width="95%",
+            display="flex",
+            flex_flow="row",
+            justify_content="space-between",
+            grid_area="cutoff_dropdown",
+        ),
+        style={"description_width": "80px"},
+    )
+
+    fig = make_subplots(
+        rows=1,
+        cols=2,
+        shared_yaxes=False,
+        horizontal_spacing=0.01,
+        vertical_spacing=0,
+        column_titles=["Firing", "Title2"],
+    )
+    fwg = go.FigureWidget(fig)
+
+    # figure_output = wg.VBox(
+    #     [QualityMetricFigs.empty_fig()],
+    #     layout=wg.Layout(width="95%", grid_area="figure_output"),
+    # )
+    # figure_output.add_class("box_style")
+
+    load_button = wg.Button(
+        description="Load",
+        tooltip="Load figures.",
+        layout=wg.Layout(width="auto", grid_area="load_button"),
+    )
+
+    def response(change, usedb=False):  # TODO: Accept cutoff vals?
+        global firing_rate_plot
+        if usedb:
+            if cluster_dropdown.value not in ephys_report.QualityMetricReport():
+                ephys_report.QualityMetricReport.populate(cluster_dropdown.value)
+
+            firing_rate_plot = from_json(
+                (ephys_report.QualityMetricReport & cluster_dropdown.value).fetch1(
+                    "firing_rate_plot"
+                )
+            )
+
+            presence_ratio_plot = from_json(
+                (ephys_report.QualityMetricReport & cluster_dropdown.value).fetch1(
+                    "presence_ratio_plot"
+                )
+            )
+
+        else:
+            qc_figs = QualityMetricFigs(cluster_dropdown)
+            firing_rate_plot = qc_figs.empty_fig()
+            presence_ratio_plot = qc_figs.empty_fig()
+
+        with fwg.batch_update():
+            fwg.data[0] = firing_rate_plot
+            fwg.data[1] = presence_ratio_plot
+
+    figure_output = wg.VBox(
+        [fwg], layout=wg.Layout(width="95%", grid_area="figure_output")
+    )
+    figure_output.add_class("box_style")
+
+    load_button.on_click(response)
+
+    main_container = wg.GridBox(
+        children=[
+            title_button,
+            cluster_dropdown,
+            cutoff_dropdown,
+            load_button,
+            figure_output,
+        ],
+        layout=wg.Layout(
+            grid_template_areas="""
+            "title_button title_button title_button"
+            "cluster_dropdown . load_button"
+            "cutoff_dropdown . load_button"
+            "figure_output figure_output figure_output"
+            """
+        ),
+        grid_template_rows="auto auto auto auto",
+        grid_template_columns="auto auto auto",
+    )
+
+    return main_container