hanhou/patchseq
0
1"""2Scatter plot component for the visualization app.3"""4 5import logging6from typing import Any, Dict, List, Tuple7 8import numpy as np9import matplotlib.pyplot as plt10import seaborn as sns11import pandas as pd12import panel as pn13from bokeh.layouts import gridplot14from bokeh.models import BoxZoomTool, ColumnDataSource, DatetimeTickFormatter, HoverTool, Legend15from bokeh.plotting import figure16from scipy import stats17from scipy.stats import mannwhitneyu18from sklearn.metrics import silhouette_score19from sklearn.mixture import GaussianMixture20from itertools import combinations21 22from LCNE_patchseq_analysis.pipeline_util.s3 import get_public_url_cell_summary23from LCNE_patchseq_analysis.pipeline_util.s3 import load_mesh_from_s324from LCNE_patchseq_analysis.data_util.mesh import trimesh_to_bokeh_data25 26from components.color_mapping import ColorMapping27from components.size_mapping import SizeMapping28 29# Set seaborn style30sns.set_context("paper")31 32logger = logging.getLogger(__name__)33 34# Define available color palettes35COLOR_PALETTES = [36 "Viridis256",37 "Plasma256",38 "Magma256",39 "Inferno256",40 "Cividis256",41 "Turbo256",42 "Set3",43 "Category10",44 "Category20",45 "Category20b",46 "Category20c",47]48 49 50class ScatterPlot:51 """Handles scatter plot creation and updates."""52 53 def __init__(self, df_meta: pd.DataFrame, data_holder: Any):54 """Initialize with metadata dataframe."""55 self.df_meta = df_meta56 self.color_mapping = ColorMapping(df_meta)57 self.size_mapping = SizeMapping(df_meta)58 self.data_holder = data_holder59 # Add cell summary URLs to dataframe60 self._add_cell_summary_urls()61 self.controls = self.create_plot_controls(width=300)62 self._latest_figures = {}63 64 def _add_cell_summary_urls(self):65 """Add cell summary URLs to the dataframe."""66 # Create a new column for cell summary URLs67 self.df_meta["cell_summary_url"] = None68 69 # Get URLs for each ephys_roi_id70 for idx, row in self.df_meta.iterrows():71 ephys_roi_id = str(int(row["ephys_roi_id"]))72 try:73 url = get_public_url_cell_summary(ephys_roi_id, if_check_exists=False)74 self.df_meta.at[idx, "cell_summary_url"] = url75 except Exception as e:76 logger.warning(f"Could not get URL for ephys_roi_id {ephys_roi_id}: {e}")77 self.df_meta.at[idx, "cell_summary_url"] = None78 79 def create_plot_controls(self, width: int = 180) -> Dict[str, Any]:80 """Create the control widgets for the scatter plot."""81 # Get numeric and categorical columns82 numeric_cols = self.df_meta.select_dtypes(include=["number"]).columns.tolist()83 categorical_cols = self.df_meta.select_dtypes(include=["object"]).columns.tolist()84 available_cols = sorted(numeric_cols + categorical_cols)85 86 # Append [valid N] to the available_cols for display purposes87 available_cols = [f"{col} [valid {self.df_meta[col].count()}]" for col in available_cols]88 all_cols = ["None"] + available_cols89 90 controls = {91 "x_axis_select": pn.widgets.Select(92 name="X Axis",93 options=all_cols,94 value=[col for col in all_cols if "Y" in col][0],95 sizing_mode="stretch_width",96 ),97 "y_axis_select": pn.widgets.Select(98 name="Y Axis",99 options=all_cols,100 value=[101 col102 for col in all_cols103 if "ipfx_tau" in col104 ][0],105 sizing_mode="stretch_width",106 ),107 "color_col_select": pn.widgets.Select(108 name="Color By",109 options=all_cols,110 value=[col for col in all_cols if "injection region" in col][0],111 sizing_mode="stretch_width",112 ),113 "color_palette_select": pn.widgets.Select(114 name="Color Palette",115 options=COLOR_PALETTES,116 value="Viridis256",117 sizing_mode="stretch_width",118 ),119 "size_col_select": pn.widgets.Select(120 name="Size By",121 options=all_cols,122 value="None",123 sizing_mode="stretch_width",124 ),125 "size_range_slider": pn.widgets.RangeSlider(126 name="Size Range",127 start=5,128 end=40,129 value=(10, 30),130 step=1,131 sizing_mode="stretch_width",132 ),133 "size_gamma_slider": pn.widgets.FloatSlider(134 name="Size Gamma",135 start=0.1,136 end=5,137 value=1,138 step=0.1,139 sizing_mode="stretch_width",140 ),141 "alpha_slider": pn.widgets.FloatSlider(142 name="Alpha",143 start=0.1,144 end=1,145 value=0.7,146 step=0.1,147 sizing_mode="stretch_width",148 ),149 "width_slider": pn.widgets.IntSlider(150 name="Width",151 start=400,152 end=1200,153 value=800,154 step=50,155 sizing_mode="stretch_width",156 ),157 "height_slider": pn.widgets.IntSlider(158 name="Height",159 start=400,160 end=1200,161 value=600,162 step=50,163 sizing_mode="stretch_width",164 ),165 "bins_slider": pn.widgets.IntSlider(166 name="Histogram bins",167 start=10,168 end=100,169 value=50,170 step=1,171 sizing_mode="stretch_width",172 ),173 "show_gmm": pn.widgets.Checkbox(174 name="Show Gaussian Mixture Model",175 value=True,176 sizing_mode="stretch_width",177 ),178 "show_linear_fit": pn.widgets.Checkbox(179 name="Show Linear Fit",180 value=True,181 sizing_mode="stretch_width",182 ),183 "n_components_x": pn.widgets.IntSlider(184 name="Number of components (X)",185 start=1,186 end=5,187 value=2,188 step=1,189 disabled=False,190 sizing_mode="stretch_width",191 ),192 "n_components_y": pn.widgets.IntSlider(193 name="Number of components (Y)",194 start=1,195 end=5,196 value=1,197 step=1,198 disabled=False,199 sizing_mode="stretch_width",200 ),201 "hist_height_slider": pn.widgets.IntSlider(202 name="Distribution plot height",203 start=50,204 end=300,205 value=150,206 step=10,207 sizing_mode="stretch_width",208 ),209 "font_size_slider": pn.widgets.IntSlider(210 name="Font Size",211 start=10,212 end=30,213 value=15,214 sizing_mode="stretch_width",215 ),216 }217 218 # Link the GMM checkbox to enable/disable the component sliders219 def toggle_gmm_components(event):220 controls["n_components_x"].disabled = not event.new221 controls["n_components_y"].disabled = not event.new222 223 controls["show_gmm"].param.watch(toggle_gmm_components, "value")224 225 # Initialize the disabled state based on the initial checkbox value226 controls["n_components_x"].disabled = not controls["show_gmm"].value227 controls["n_components_y"].disabled = not controls["show_gmm"].value228 229 return controls230 231 def sync_controls_to_url(self):232 """Sync scatter plot controls to URL query parameters."""233 234 location = pn.state.location235 mapping = {236 "x_axis_select": ("value", "scatter_x"),237 "y_axis_select": ("value", "scatter_y"),238 "color_col_select": ("value", "scatter_color"),239 "color_palette_select": ("value", "scatter_palette"),240 "size_col_select": ("value", "scatter_size"),241 "size_range_slider": ("value", "scatter_size_range"),242 "size_gamma_slider": ("value", "scatter_gamma"),243 "alpha_slider": ("value", "scatter_alpha"),244 "width_slider": ("value", "scatter_width"),245 "height_slider": ("value", "scatter_height"),246 "bins_slider": ("value", "scatter_bins"),247 "show_gmm": ("value", "scatter_gmm"),248 "show_linear_fit": ("value", "scatter_linear_fit"),249 "n_components_x": ("value", "scatter_components_x"),250 "n_components_y": ("value", "scatter_components_y"),251 "hist_height_slider": ("value", "scatter_hist_height"),252 "font_size_slider": ("value", "scatter_font_size"),253 }254 for control_name, (param_name, url_param) in mapping.items():255 location.sync(self.controls[control_name], {param_name: url_param})256 257 def create_tooltips(258 self, x_col: str, y_col: str, color_col: str, size_col: str259 ) -> List[Tuple[str, str]]:260 """Create tooltips for the hover tool."""261 262 tooltips = f"""263 <div style="text-align: left; flex: auto; white-space: nowrap; margin: 0 10px;264 border: 2px solid black; padding: 10px;">265 <span style="font-size: 17px;">266 <b>@Date_str, @{{injection region}}, @{{ephys_roi_id}},267 @{{jem-id_cell_specimen}}</b><br>268 <b>X = @{{{x_col}}}</b> [{x_col}]<br>269 <b>Y = @{{{y_col}}}</b> [{y_col}]<br>270 <b> Color = @{{{color_col}}}</b> [{color_col}]<br>271 <b> Size = @{{{size_col}}}</b> [{size_col}]<br>272 </span>273 <img src="@cell_summary_url{{safe}}" alt="Cell Summary"274 style="width: 800px; height: auto;">275 </div>276 """277 278 return tooltips279 280 def create_marginal_histogram(281 self,282 data: pd.Series,283 orientation: str,284 width: int,285 height: int,286 alpha: float,287 bins: int,288 show_gmm: bool = False,289 n_components: int = 1,290 ) -> figure:291 """Create a histogram for marginal distribution with optional GMM overlay."""292 # Remove NaN values and convert to numeric293 clean_data = pd.to_numeric(data, errors="coerce").dropna()294 295 # If no valid data, create an empty plot296 if clean_data.empty:297 p = figure(298 height=height,299 width=width,300 tools="",301 toolbar_location=None,302 x_range=(0, 1),303 y_range=(0, 1),304 )305 p.text(306 x=0.5,307 y=0.5,308 text=["No valid data"],309 text_align="center",310 text_baseline="middle",311 )312 return p313 314 # Calculate histogram data (independent of orientation)315 hist, edges = np.histogram(clean_data, bins=bins, density=True)316 317 # Set axis ranges and quad parameters based on orientation318 if orientation == "x":319 x_range = (edges[0], edges[-1])320 y_range = (0, hist.max() * 1.1)321 else: # "y" orientation322 x_range = (0, hist.max() * 1.1)323 y_range = (edges[0], edges[-1])324 325 # Create the figure326 p = figure(327 height=height,328 width=width,329 tools="",330 toolbar_location=None,331 x_range=x_range,332 y_range=y_range,333 )334 335 # Plot the histogram using vbar/hbar (Bokeh's bar plot) instead of quad336 if orientation == "x":337 # Use vbar for x-orientation338 p.vbar(339 x=[(edges[i] + edges[i + 1]) / 2 for i in range(len(edges) - 1)],340 top=hist,341 width=(edges[1] - edges[0]) * 0.9, # Slightly narrower than bin width342 fill_color="gray",343 line_color="white",344 alpha=0.9,345 )346 else: # "y" orientation347 # Use hbar for y-orientation348 p.hbar(349 y=[(edges[i] + edges[i + 1]) / 2 for i in range(len(edges) - 1)],350 right=hist,351 height=(edges[1] - edges[0]) * 0.9, # Slightly narrower than bin width352 fill_color="gray",353 line_color="white",354 alpha=0.9,355 )356 357 # Optional: Plot Gaussian Mixture Model overlay358 if show_gmm:359 360 gmm = GaussianMixture(n_components=n_components, random_state=42)361 gmm.fit(clean_data.values.reshape(-1, 1))362 domain = np.linspace(edges[0], edges[-1], 1000)363 density = np.exp(gmm.score_samples(domain.reshape(-1, 1)))364 365 # Calculate evaluation metrics366 if n_components > 1:367 labels = gmm.predict(clean_data.values.reshape(-1, 1))368 silhouette = silhouette_score(clean_data.values.reshape(-1, 1), labels)369 370 # Calculate BIC and AIC371 bic = gmm.bic(clean_data.values.reshape(-1, 1))372 else:373 bic = np.nan374 silhouette = np.nan375 376 p.line(377 *((domain, density) if orientation == "x" else (density, domain)),378 line_color="black",379 line_width=4,380 alpha=0.9,381 )382 383 # Plot individual components384 for i in range(n_components):385 mean = gmm.means_[i][0]386 std = np.sqrt(gmm.covariances_[i][0][0])387 weight = gmm.weights_[i]388 comp_density = (389 weight390 * np.exp(-0.5 * ((domain - mean) / std) ** 2)391 / (std * np.sqrt(2 * np.pi))392 )393 p.line(394 *((domain, comp_density) if orientation == "x" else (comp_density, domain)),395 line_color="black",396 line_width=2,397 alpha=0.9,398 line_dash="dashed",399 )400 401 # Add metrics to the plot title402 axis_to_show = p.xaxis if orientation == "x" else p.yaxis403 404 axis_to_show.axis_label = f"Silhouette: {silhouette:.3f}, " f"BIC: {bic:.3f}"405 406 # Font size407 axis_to_show.axis_label_text_font_size = "10pt"408 axis_to_show.major_label_text_font_size = "0pt"409 410 # Hide axes and grid411 axis_to_hide = p.xaxis if orientation == "y" else p.yaxis412 axis_to_hide.visible = False413 p.grid.visible = False414 return p415 416 def add_lc_mesh_overlay(self, p: figure, x_col: str, y_col: str) -> None:417 """Add LC mesh overlay to the plot based on axis column names.418 419 Args:420 p: Bokeh figure to add the mesh to421 x_col: X-axis column name422 y_col: Y-axis column name423 """424 # Determine mesh direction based on column names425 direction = None426 if x_col.startswith("X ") and y_col.startswith("Y "):427 direction = "sagittal"428 elif x_col.startswith("Z "):429 direction = "coronal"430 431 if direction is not None:432 try:433 # Load and add LC mesh overlay434 mesh = load_mesh_from_s3()435 lc_mesh_bokeh = trimesh_to_bokeh_data(mesh, direction=direction)436 mesh_source = ColumnDataSource(lc_mesh_bokeh)437 p.patches(438 source=mesh_source,439 xs="xs",440 ys="ys",441 fill_alpha=0.3,442 line_color=None,443 fill_color="lightgray",444 level="underlay",445 nonselection_fill_alpha=0.3,446 nonselection_line_alpha=0,447 selection_fill_alpha=0.3,448 selection_line_alpha=0,449 muted_alpha=0.3,450 )451 except Exception as e:452 logger.warning(f"Could not add LC mesh overlay: {e}")453 454 def update_scatter_plot( # noqa: C901455 self,456 x_col: str,457 y_col: str,458 color_col: str,459 color_palette: str,460 size_col: str,461 size_range: tuple,462 size_gamma: float,463 alpha: float,464 width: int,465 height: int,466 font_size: int = 14,467 bins: int = 30,468 hist_height_slider: int = 100,469 show_gmm: bool = False,470 n_components_x: int = 2,471 n_components_y: int = 1,472 show_linear_fit: bool = True,473 df_meta: pd.DataFrame = None,474 ) -> gridplot:475 """Update the scatter plot with new parameters."""476 # Use provided dataframe if supplied, otherwise use instance df_meta477 df_to_use = df_meta if df_meta is not None else self.df_meta478 479 # Strip off [valid N] from the column name480 x_col = x_col.split(" [valid ")[0]481 y_col = y_col.split(" [valid ")[0]482 color_col = color_col.split(" [valid ")[0]483 size_col = size_col.split(" [valid ")[0]484 485 # Create a new figure for the main scatter plot486 p = figure(487 x_axis_label=x_col,488 y_axis_label=y_col,489 tools="pan,wheel_zoom,box_zoom,reset,tap",490 height=height,491 width=width,492 )493 494 # Create ColumnDataSource from the dataframe495 source = ColumnDataSource(df_to_use)496 497 # If any column is Date, convert it to datetime498 if x_col == "Date":499 source.data[x_col] = pd.to_datetime(pd.Series(source.data[x_col]), errors="coerce")500 p.xaxis.formatter = DatetimeTickFormatter(501 years="%Y",502 months="%Y-%m",503 days="%Y-%m-%d",504 )505 506 # Create temporary color mapping for this specific dataframe507 temp_color_mapping = ColorMapping(df_to_use)508 # Determine color mapping509 color = temp_color_mapping.determine_color_mapping(510 color_col, color_palette, p, font_size=font_size511 )512 513 # Create temporary size mapping for this specific dataframe514 temp_size_mapping = SizeMapping(df_to_use)515 # Determine size mapping516 size = temp_size_mapping.determine_size_mapping(517 size_col, source, min_size=size_range[0], max_size=size_range[1], gamma=size_gamma518 )519 520 # Add scatter glyph using the data source521 scatter_glyph = p.scatter(x=x_col, y=y_col, source=source, size=size, color=color, alpha=alpha)522 523 # Add linear regression if requested and both columns are numeric524 if show_linear_fit and x_col != "Date" and x_col != "None" and y_col != "None":525 # Get clean numeric data526 # Convert to numeric and drop rows where either x or y is NA527 df_clean = df_to_use[[x_col, y_col]].apply(pd.to_numeric, errors="coerce").dropna()528 x_data = df_clean[x_col]529 y_data = df_clean[y_col]530 531 # Only proceed if we have valid data532 if not x_data.empty and not y_data.empty:533 # Perform linear regression534 slope, intercept, r_value, p_value, std_err = stats.linregress(x_data, y_data)535 536 # Calculate fitted line points over dense grid for smooth CI curve537 x_min, x_max = x_data.min(), x_data.max()538 x_vals = np.linspace(x_min, x_max, 200)539 y_fit = slope * x_vals + intercept540 541 # Add confidence band for the regression line542 n_points = len(x_data)543 if n_points > 2:544 residuals = y_data - (slope * x_data + intercept)545 mse = np.sum(residuals ** 2) / (n_points - 2)546 x_mean = x_data.mean()547 Sxx = np.sum((x_data - x_mean) ** 2)548 if Sxx > 0:549 t_val = stats.t.ppf(0.975, n_points - 2)550 se_fit = np.sqrt(551 mse * (1 / n_points + (x_vals - x_mean) ** 2 / Sxx)552 )553 ci_upper = y_fit + t_val * se_fit554 ci_lower = y_fit - t_val * se_fit555 band_source = ColumnDataSource(556 {557 "x": np.concatenate([x_vals, x_vals[::-1]]),558 "y": np.concatenate([ci_upper, ci_lower[::-1]]),559 }560 )561 p.patch(562 x="x",563 y="y",564 source=band_source,565 fill_color="lightgray",566 fill_alpha=0.3,567 line_alpha=0,568 level="underlay",569 )570 571 # Add fitted line572 setting = (573 {"line_width": 3, "line_dash": "solid"}574 if p_value < 0.05575 else {"line_width": 2, "line_dash": "dashed"}576 )577 line = p.line(x_vals, y_fit, line_color="black", **setting)578 579 # Add legend with R² and p-value580 legend_items = [581 (f"Linear Fit (p = {p_value:.3e}, R² = {r_value**2:.3f})", [line]),582 ]583 legend = Legend(584 items=legend_items, 585 location="top_left", 586 label_text_font_size=f"{font_size-2}pt"587 )588 legend.click_policy = "hide"589 p.add_layout(legend, 'above')590 591 # Flip the y-axis if y_col is depth592 if y_col == "Y (D --> V)":593 p.y_range.flipped = True594 595 # Add HoverTool with tooltips596 tooltips = self.create_tooltips(x_col, y_col, color_col, size_col)597 hovertool = HoverTool(598 tooltips=tooltips,599 attachment="right", # Fix tooltip to the right of the plot600 formatters={"@Date": "datetime"},601 renderers=[scatter_glyph], # Only apply hover to scatter points, not mesh patches602 )603 604 p.add_tools(hovertool)605 606 # Define callback to update ephys_roi_id on point tap607 def update_ephys_roi_id(attr, old, new):608 if new:609 selected_index = new[0]610 ephys_roi_id = str(int(df_to_use.iloc[selected_index]["ephys_roi_id"]))611 logger.info(f"Selected ephys_roi_id: {ephys_roi_id}")612 # Update the data holder's ephys_roi_id613 if hasattr(self, "data_holder"):614 self.data_holder.ephys_roi_id_selected = ephys_roi_id615 616 # Attach the callback to the selection changes617 source.selected.on_change("indices", update_ephys_roi_id)618 619 # Set the default tool activated on drag to be box zoom620 p.toolbar.active_drag = p.select_one(BoxZoomTool)621 622 # Set axis label font sizes623 p.xaxis.axis_label_text_font_size = f"{font_size}pt"624 p.yaxis.axis_label_text_font_size = f"{font_size}pt"625 626 # Set major tick label font sizes627 p.xaxis.major_label_text_font_size = f"{font_size*0.9}pt"628 p.yaxis.major_label_text_font_size = f"{font_size*0.9}pt"629 630 # Add LC mesh overlay if appropriate columns are selected631 self.add_lc_mesh_overlay(p, x_col, y_col)632 633 # Create marginal histograms634 x_hist = None635 try:636 if x_col != "Date" and x_col != "None": # Skip histogram for Date column637 x_hist = self.create_marginal_histogram(638 df_to_use[x_col],639 "x",640 width=width,641 height=hist_height_slider,642 alpha=alpha,643 bins=bins,644 show_gmm=show_gmm,645 n_components=n_components_x,646 )647 x_hist.x_range = p.x_range # Link x ranges648 except Exception as e:649 logger.warning(f"Could not create x histogram: {e}")650 x_hist = None651 652 y_hist = None653 try:654 if y_col != "Date" and y_col != "None": # Skip histogram for Date column655 y_hist = self.create_marginal_histogram(656 df_to_use[y_col],657 "y",658 width=hist_height_slider,659 height=height,660 alpha=alpha,661 bins=bins,662 show_gmm=show_gmm,663 n_components=n_components_y,664 )665 y_hist.y_range = p.y_range # Link y ranges666 except Exception as e:667 logger.warning(f"Could not create y histogram: {e}")668 y_hist = None669 670 # Count non-NaN values grouped by "injection region"671 count_non_nan = df_to_use.groupby("injection region")[[x_col, y_col]].count().T672 count_non_nan.insert(0, "Total", count_non_nan.sum(axis=1))673 count_non_nan.index = pd.Index(["X", "Y"], name="Valid N")674 675 # Count NaN values (missing data) grouped by "injection region"676 count_nan = df_to_use.groupby("injection region")[[x_col, y_col]].apply(lambda x: x.isna().sum()).T677 count_nan.insert(0, "Total", count_nan.sum(axis=1))678 count_nan.index = pd.Index(["X", "Y"], name="Missing Data")679 680 # If the color column is categorical, generate violin plot and pairwise statistical tests681 if color_col in self.df_meta.select_dtypes(include=["object"]).columns:682 # --- Create marginalized histograms to compare across colors ---683 # marginalized_histograms, _ = self.create_marginalized_histograms(684 # df_to_use, y_col, color_col, color_palette, temp_color_mapping, p, font_size685 # )686 687 # --- Create violin plot to compare across injection regions ---688 violin_plot = self.create_violin_plot(689 df_to_use, y_col, color_col, color_palette, temp_color_mapping, p, font_size690 )691 692 # Perform pairwise statistical tests693 # Drop NA for y_col and color_col for statistical tests694 plot_df_for_stats = df_to_use[[y_col, color_col]].dropna() if (y_col != "Date" and y_col != "None" and color_col != "None") else pd.DataFrame()695 pvalues_table = self.perform_pairwise_statistical_tests(plot_df_for_stats, y_col, color_col) if not plot_df_for_stats.empty else pn.pane.Markdown("**No statistical tests available**")696 else:697 violin_plot = pn.pane.Markdown("**Violin plot only available for categorical color columns**")698 pvalues_table = pn.pane.Markdown("")699 700 # Create grid layout701 layout = pn.Row(702 pn.Column(703 gridplot(704 [[y_hist, p], [None, x_hist]],705 toolbar_location="right",706 merge_tools=False,707 toolbar_options={"logo": None},708 ),709 pn.pane.Markdown(count_non_nan.to_markdown()),710 pn.pane.Markdown(count_nan.to_markdown()),711 ),712 pn.Column(713 violin_plot,714 pvalues_table,715 pn.Spacer(height=20),716 # marginalized_histograms,717 ),718 )719 720 # Store figures for export721 self._latest_figures = {722 "scatter_plot": p,723 }724 if x_hist is not None:725 self._latest_figures["x_histogram"] = x_hist726 if y_hist is not None:727 self._latest_figures["y_histogram"] = y_hist728 # Store violin plot if it's a matplotlib figure (not just a markdown pane)729 if hasattr(violin_plot, 'object') and hasattr(violin_plot.object, 'savefig'):730 self._latest_figures["violin_plot"] = violin_plot731 732 return layout733 734 def create_marginalized_histograms(735 self,736 df_to_use: pd.DataFrame,737 y_col: str,738 color_col: str,739 color_palette: str,740 temp_color_mapping: ColorMapping,741 p: figure,742 font_size: int743 ) -> Tuple[Any, Any]:744 """Create marginalized histograms to compare data across color groups.745 746 Args:747 df_to_use: DataFrame containing the data748 y_col: Column name for y-axis variable749 color_col: Column name for color grouping750 color_palette: Color palette name751 temp_color_mapping: ColorMapping instance for this data752 p: Bokeh figure (used for color mapping extraction)753 font_size: Font size for labels754 755 Returns:756 Tuple of (marginalized_histograms, None) - Creates KDE plots with mean±SEM757 """758 # Prepare marginalized histogram using seaborn's histplot (KDE) for y_col by color_col759 marginalized_histograms = pn.pane.Markdown("No marginalized histogram available.")760 761 try:762 if y_col != "Date" and y_col != "None" and color_col != "None":763 fig, ax = plt.subplots(figsize=(4, 3.5), dpi=300)764 # Drop NA for y_col and color_col765 plot_df = df_to_use[[y_col, color_col]].dropna()766 if not plot_df.empty:767 # Extract color mapping from the scatter plot768 color_palette_dict = None769 color_mapping_result = temp_color_mapping.determine_color_mapping(770 color_col, color_palette, p, font_size=font_size, if_add_color_bar=False771 )772 if isinstance(color_mapping_result, dict) and 'transform' in color_mapping_result:773 color_mapper = color_mapping_result['transform']774 if hasattr(color_mapper, 'factors') and hasattr(color_mapper, 'palette'):775 color_palette_dict = dict(zip(color_mapper.factors, color_mapper.palette))776 777 # Count number of samples per group (valid data)778 group_counts = plot_df[color_col].value_counts().to_dict()779 780 # Count missing data (NaN) per group from original dataframe781 group_nan_counts = {}782 for group in group_counts.keys():783 # Get all rows for this group from original dataframe784 group_mask = df_to_use[color_col] == group785 # Count NaN values in y_col for this group786 nan_count = df_to_use.loc[group_mask, y_col].isna().sum()787 group_nan_counts[group] = nan_count788 789 # Create a mapping from original group name to "group (n = valid, nan = missing)"790 group_labels = {791 group: f"{group} (n = {count}, missing {group_nan_counts.get(group, 0)})" 792 for group, count in group_counts.items()793 }794 # Add a new column for legend labels795 plot_df["_legend_label"] = plot_df[color_col].map(group_labels)796 797 sns.kdeplot(798 data=plot_df,799 x=y_col,800 hue="_legend_label",801 common_norm=False,802 fill=False,803 ax=ax,804 palette=color_palette_dict if color_palette_dict is None else {805 group_labels[group]: color_palette_dict[group] for group in group_labels if group in color_palette_dict806 },807 )808 sns.despine(trim=True)809 ax.set_xlabel(y_col)810 811 # Compute mean ± SEM for each group and add as dot + errorbar812 y_positions = [] # Track y positions for staggering813 for i, (group, data_subset) in enumerate(plot_df.groupby(color_col)):814 values = data_subset[y_col].dropna()815 if len(values) > 0:816 # Ensure values are numeric and convert to float817 try:818 numeric_values = pd.to_numeric(values, errors='coerce').dropna()819 if len(numeric_values) > 0:820 mean_val = float(np.mean(numeric_values))821 # Use numpy's std with ddof=1 to calculate SEM manually822 if len(numeric_values) > 1:823 sem_val = float(np.std(numeric_values, ddof=1) / np.sqrt(len(numeric_values)))824 else:825 sem_val = 0.0826 827 # Get color for this group828 group_color = color_palette_dict.get(group, 'black') if color_palette_dict else 'black'829 830 # Stagger y position slightly for each group831 # Compute y position as 10% + i * 10% of the current ylim range832 ylim = ax.get_ylim()833 y_pos = ylim[0] + 0.05 * (ylim[1] - ylim[0]) + i * 0.05 * (ylim[1] - ylim[0])834 y_positions.append(y_pos)835 836 # Add dot for mean837 ax.plot(mean_val, y_pos, 'o', color=group_color, markersize=4, 838 markeredgewidth=1)839 840 # Add error bar for SEM841 if sem_val > 0.0:842 ax.errorbar(mean_val, y_pos, xerr=sem_val, color=group_color, 843 capsize=3, capthick=1.5, elinewidth=1.5, zorder=9)844 except (ValueError, TypeError):845 # Skip non-numeric data846 continue847 848 # Adjust y-axis limits to accommodate the error bars849 current_ylim = ax.get_ylim()850 if y_positions:851 max_y_pos = max(y_positions)852 ax.set_ylim(current_ylim[0], max(current_ylim[1], max_y_pos + 0.02))853 854 # Move legend to top of the plot855 y_lim = ax.get_ylim()856 sns.move_legend(857 ax,858 loc="center left",859 bbox_to_anchor=(1.01, 0.5),860 ncol=1,861 frameon=False,862 fontsize="small",863 title=color_col,864 )865 # Use Panel's matplotlib pane instead of manual base64 conversion866 marginalized_histograms = pn.pane.Matplotlib(fig, dpi=300, tight=True, width=400)867 868 # Close the figure to prevent memory leaks869 plt.close(fig)870 871 except Exception as e:872 logger.warning(f"Could not create marginalized KDE histogram: {e}")873 marginalized_histograms = pn.pane.Markdown("Marginalized histogram error.")874 875 return marginalized_histograms, None876 877 def perform_pairwise_statistical_tests(878 self,879 plot_df: pd.DataFrame,880 y_col: str,881 color_col: str882 ) -> Any:883 """Perform pairwise Mann-Whitney U tests between groups.884 885 Args:886 plot_df: DataFrame containing the cleaned data887 y_col: Column name for y-axis variable888 color_col: Column name for color grouping889 890 Returns:891 Panel markdown object with statistical test results892 """893 pvalues_table = pn.pane.Markdown("**No statistical tests available**")894 895 try:896 # Perform pairwise Mann-Whitney U tests897 pairwise_tests = {}898 groups = list(plot_df[color_col].unique())899 900 for group1, group2 in combinations(groups, 2):901 data1 = pd.to_numeric(plot_df[plot_df[color_col] == group1][y_col], errors='coerce').dropna()902 data2 = pd.to_numeric(plot_df[plot_df[color_col] == group2][y_col], errors='coerce').dropna()903 904 if len(data1) > 0 and len(data2) > 0:905 try:906 statistic, p_value = mannwhitneyu(data1, data2, alternative='two-sided')907 pairwise_tests[f"{group1} vs {group2}"] = p_value908 except Exception as e:909 logger.warning(f"Could not perform Mann-Whitney U test for {group1} vs {group2}: {e}")910 pairwise_tests[f"{group1} vs {group2}"] = np.nan911 912 # Create a table of p-values913 if pairwise_tests:914 pvalues_df = pd.DataFrame(list(pairwise_tests.items()), columns=['Comparison', 'p-value'])915 pvalues_df['p-value'] = pvalues_df['p-value'].apply(lambda x: f"{x:.3e}" if not pd.isna(x) else "NaN")916 pvalues_table = pn.pane.Markdown(917 f"**Mann-Whitney U Test (pairwise comparisons)**\n\n{pvalues_df.to_markdown(index=False)}"918 )919 else:920 pvalues_table = pn.pane.Markdown("**No pairwise comparisons available**")921 922 except Exception as e:923 logger.warning(f"Could not perform pairwise statistical tests: {e}")924 pvalues_table = pn.pane.Markdown("**Statistical test error**")925 926 return pvalues_table927 928 def create_violin_plot(929 self,930 df_to_use: pd.DataFrame,931 y_col: str,932 color_col: str,933 color_palette: str,934 temp_color_mapping: ColorMapping,935 p: figure,936 font_size: int937 ) -> Any:938 """Create violin plot to compare data distributions across injection regions.939 940 Args:941 df_to_use: DataFrame containing the data942 y_col: Column name for y-axis variable943 color_col: Column name for color grouping (injection regions)944 color_palette: Color palette name945 temp_color_mapping: ColorMapping instance for this data946 p: Bokeh figure (used for color mapping extraction)947 font_size: Font size for labels948 949 Returns:950 Panel matplotlib object with violin plot951 """952 violin_plot = pn.pane.Markdown("No violin plot available.")953 954 try:955 if y_col != "Date" and y_col != "None" and color_col != "None":956 fig, ax = plt.subplots(figsize=(5, 4), dpi=300)957 # Drop NA for y_col and color_col958 plot_df = df_to_use[[y_col, color_col]].dropna()959 if not plot_df.empty:960 # Extract color mapping from the scatter plot961 color_palette_dict = None962 color_mapping_result = temp_color_mapping.determine_color_mapping(963 color_col, color_palette, p, font_size=font_size, if_add_color_bar=False964 )965 if isinstance(color_mapping_result, dict) and 'transform' in color_mapping_result:966 color_mapper = color_mapping_result['transform']967 if hasattr(color_mapper, 'factors') and hasattr(color_mapper, 'palette'):968 color_palette_dict = dict(zip(color_mapper.factors, color_mapper.palette))969 970 # Count number of samples per group (valid data)971 group_counts = plot_df[color_col].value_counts().to_dict()972 973 # Count missing data (NaN) per group from original dataframe974 group_nan_counts = {}975 for group in group_counts.keys():976 # Get all rows for this group from original dataframe977 group_mask = df_to_use[color_col] == group978 # Count NaN values in y_col for this group979 nan_count = df_to_use.loc[group_mask, y_col].isna().sum()980 group_nan_counts[group] = nan_count981 982 # Convert y_col to numeric983 plot_df[y_col] = pd.to_numeric(plot_df[y_col], errors='coerce')984 plot_df = plot_df.dropna(subset=[y_col])985 986 if not plot_df.empty:987 # Get the order of groups to ensure consistency between violin plot and overlays988 # Use sorted order to match seaborn's default behavior989 groups_order = sorted(plot_df[color_col].unique())990 991 # Create violin plot using seaborn with explicit order992 sns.violinplot(993 data=plot_df,994 x=color_col,995 y=y_col,996 hue=color_col,997 ax=ax,998 palette=color_palette_dict,999 inner="quart", # Show quartiles as inner elements1000 alpha=0.6,1001 cut=0, # No extension beyond the data range1002 order=groups_order, # Explicitly set the order1003 width=0.5,1004 legend=False1005 )1006 1007 # Overlay raw data points with strip plot using the same order1008 sns.stripplot(1009 data=plot_df,1010 x=color_col,1011 y=y_col,1012 ax=ax,1013 color='black',1014 size=2,1015 alpha=0.5,1016 jitter=True,1017 order=groups_order # Use the same order1018 )1019 1020 # Calculate and plot mean ± SEM for each group using the same order1021 groups = groups_order # Use the same order as seaborn plots1022 x_positions = np.arange(len(groups))1023 1024 for i, group in enumerate(groups):1025 group_data = pd.to_numeric(plot_df[plot_df[color_col] == group][y_col], errors='coerce').dropna()1026 if len(group_data) > 0:1027 mean_val = float(np.mean(group_data))1028 if len(group_data) > 1:1029 sem_val = float(np.std(group_data, ddof=1) / np.sqrt(len(group_data)))1030 else:1031 sem_val = 0.01032 1033 # Plot mean as a larger point1034 group_color = color_palette_dict.get(group, 'black') if color_palette_dict else 'black'1035 ax.plot(i+0.45, mean_val, 'o', color=group_color, markersize=5, 1036 markeredgecolor='black', markeredgewidth=1, zorder=10)1037 1038 # Add error bar for SEM1039 if sem_val > 0.0:1040 ax.errorbar(i+0.45, mean_val, yerr=sem_val, color='black', 1041 capsize=5, capthick=1, elinewidth=1, zorder=9,1042 fmt='none')1043 1044 # Set x-axis labels with sample counts1045 group_labels_with_counts = [1046 f"{group}\n(n={group_counts.get(group, 0)}, missing {group_nan_counts.get(group, 0)})" 1047 for group in groups1048 ]1049 ax.set_xticks(x_positions)1050 ax.set_xticklabels(group_labels_with_counts, rotation=30, ha='right')1051 ax.set_ylabel(y_col)1052 ax.set_xlabel(color_col)1053 1054 # Add grid and styling1055 sns.despine(trim=True)1056 1057 # Adjust layout to prevent label cutoff1058 plt.tight_layout()1059 1060 # Use Panel's matplotlib pane1061 violin_plot = pn.pane.Matplotlib(fig, dpi=300, tight=True, width=400)1062 1063 # Close the figure to prevent memory leaks1064 plt.close(fig)1065 1066 except Exception as e:1067 logger.warning(f"Could not create violin plot: {e}")1068 violin_plot = pn.pane.Markdown("Violin plot error.")1069 1070 return violin_plot1071