diff --git a/bin/workflow_glue/report.py b/bin/workflow_glue/report.py index 306a54b..fe5dbef 100644 --- a/bin/workflow_glue/report.py +++ b/bin/workflow_glue/report.py @@ -4,23 +4,36 @@ import json import math from pathlib import Path -from dominate.tags import div, h4, p, pre, strong +from bokeh.resources import INLINE as BOKEH_INLINE +from dominate.tags import div, h3, h4, p, pre, script, strong from dominate.util import raw from ezcharts.components import fastcat +from ezcharts.components.ezchart import EZChart from ezcharts.components.reports import labs +from ezcharts.components.theme import LAB_head_resources +from ezcharts.layout.resource import Resource as EZC_Resource from ezcharts.layout.snippets import Tabs from ezcharts.layout.snippets.table import DataTable import pandas as pd from .util import get_named_logger, wf_parser # noqa: ABS101 +from .volcano import volcano # noqa: ABS101 -def _find_table(directory, pattern): - matches = sorted(Path(directory).glob(pattern)) - return matches[0] if matches else None +def get_bokeh_widgets_js(): + """Return the inline Bokeh widgets JavaScript bundle.""" + widgets_index = BOKEH_INLINE.components_for("js").index("bokeh-widgets") + return raw(BOKEH_INLINE.js_raw[widgets_index]) + + +def get_bokeh_tables_js(): + """Return the inline Bokeh tables JavaScript bundle.""" + tables_index = BOKEH_INLINE.components_for("js").index("bokeh-tables") + return raw(BOKEH_INLINE.js_raw[tables_index]) def _read_table(path, **kwargs): + """Read a TSV file into a DataFrame, returning None if path is absent.""" if path is None or not Path(path).exists(): return None return pd.read_csv(path, sep="\t", **kwargs) @@ -56,6 +69,7 @@ def _format_ratio_value(value): def _cohort_summary(cohort_dir): + """Return cohort-level metrics and transcript class counts DataFrames.""" tx_meta = _read_table(Path(cohort_dir) / "transcript_metadata.tsv") if tx_meta is None: return None, None @@ -79,6 +93,7 @@ def _cohort_summary(cohort_dir): def _sample_summaries(samples_dir): + """Return a dict of per-sample metrics DataFrames keyed by sample name.""" summaries = {} for sample_dir in sorted(Path(samples_dir).iterdir()): if not sample_dir.is_dir(): @@ -104,6 +119,7 @@ def _sample_summaries(samples_dir): def _sqanti_tables(sqanti_dir): + """Return a dict of SQANTI3 classification summary DataFrames keyed by label.""" tables = {} for summary in sorted(Path(sqanti_dir).rglob("classification_summary.tsv")): label = summary.parent.name @@ -111,7 +127,8 @@ def _sqanti_tables(sqanti_dir): return tables -def _top_results(de_dir, filename, n=20): +def _contrast_results(de_dir, filename, n=None): + """Return a dict of per-contrast result DataFrames read from filename.""" tables = {} for contrast_dir in sorted(Path(de_dir).iterdir()): if not contrast_dir.is_dir(): @@ -119,7 +136,10 @@ def _top_results(de_dir, filename, n=20): table = _read_table(contrast_dir / filename) if table is None or table.empty: continue - tables[contrast_dir.name] = table.head(n) + data = table + if n is not None: + data = data.head(n) + tables[contrast_dir.name] = data return tables @@ -266,6 +286,10 @@ def main(args): args.params, args.versions, args.wf_version, + head_resources=[ + *LAB_head_resources, + EZC_Resource(func=get_bokeh_widgets_js, tag=script), + EZC_Resource(func=get_bokeh_tables_js, tag=script)] ) with open(args.metadata, "r") as handle: @@ -869,7 +893,9 @@ def main(args): with report.add_section("Differential gene expression", "DGE"): tabs = Tabs() - for contrast, table in _top_results(args.de_dir, "results_dge.tsv").items(): + for contrast, table in _contrast_results( + args.de_dir, "results_dge.tsv", n=20 + ).items(): with tabs.add_tab(contrast): # Check for contrast-specific warnings if de_qc and contrast in de_qc.get("contrasts", {}): @@ -900,9 +926,35 @@ def main(args): DataTable.from_pandas(table, use_index=False) + h3("Gene expression volcano Plot") + vol, class_table, selected_table = volcano(table) + EZChart(vol, width="100%", height="550") + raw(""" + + """) + with div(_class="volcano-table-grid"): + EZChart(class_table, width="100%", height="auto") + EZChart(selected_table, width="100%", height="auto") + with report.add_section("Differential transcript usage", "DTU"): tabs = Tabs() - dtu_tables = _top_results(args.de_dir, "results_dtu_transcript.tsv") + dtu_tables = _contrast_results( + args.de_dir, "results_dtu_transcript.tsv", n=20) for contrast in sorted(Path(args.de_dir).iterdir()): if not contrast.is_dir(): @@ -940,6 +992,9 @@ def main(args): dtu_tables[contrast_name], use_index=False, ) + h3("Transcript expression volcano Plot") + # EZChart(volcano(dtu_tables[contrast_name])) + else: p("No DTU results available for this contrast.") diff --git a/bin/workflow_glue/tests/common/test_report.py b/bin/workflow_glue/tests/common/test_report.py index 3a9f7e0..bda75f0 100644 --- a/bin/workflow_glue/tests/common/test_report.py +++ b/bin/workflow_glue/tests/common/test_report.py @@ -94,13 +94,16 @@ def _build_report_args(tmp_path, de_qc=None): contrast_dir.mkdir() _write( contrast_dir / "results_dge.tsv", - "GENEID\tlog2FoldChange\tpadj\n" - "gene1\t1.0\t0.01\n", + ( + "GENEID\tnewGeneClass\tgene_name\tbaseMean\tlog2FoldChange\tlfcSE" + "\tstat\tpvlaue\tpadj\n" + "gene1\tannotation\tgene2\t100.0\t1.0\t-0.02\t-8.5\t0.001\t0.05\n" + ) ) _write( contrast_dir / "results_dtu_transcript.tsv", - "featureID\tgroupID\tpadj\n" - "tx1\tgene1\t0.05\n", + "featureID\tgroupID\tlog2FoldChange\tpvalue\tpadj\texonBaseMean\n" + "tx1\tgene1\t1.0\t0.01\t0.05\t20.0\n" ) out_report = tmp_path / "wf-transcriptomes-report.html" diff --git a/bin/workflow_glue/volcano.py b/bin/workflow_glue/volcano.py new file mode 100644 index 0000000..a4115e9 --- /dev/null +++ b/bin/workflow_glue/volcano.py @@ -0,0 +1,1269 @@ +"""Interactive volcano plot visualization for differential expression analysis.""" + +from bokeh.events import Tap +from bokeh.layouts import column as bokeh_column, row as bokeh_row +from bokeh.models import ( + AutocompleteInput, + Button, + ColumnDataSource, + CustomJS, + CustomJSHover, + DataTable, + HoverTool, + HTMLTemplateFormatter, + LabelSet, + NumberFormatter, + Range1d, + Slider, + Span, + TableColumn, + TextInput, + Toggle +) +from bokeh.plotting import figure +from ezcharts.plots import BokehPlot +import numpy as np +import pandas as pd + +SIGNIFICANCE_CLASSES = { + "significant_high_effect": { + "color": "#d62728", + "marker": "circle", + "symbol": "●", + "label": "Significant, high effect", + }, + "significant_low_effect": { + "color": "#1f77b4", + "marker": "square", + "symbol": "■", + "label": "Significant, low effect", + }, + "non_significant_high_effect": { + "color": "#2ca02c", + "marker": "triangle", + "symbol": "▲", + "label": "Not significant, high effect", + }, + "non_significant_low_effect": { + "color": "#8a8a8a", + "marker": "x", + "symbol": "✕", + "label": "Not significant, low effect", + }, +} + + +def _class_attr(attr,): + return {k: v[attr] for k, v in SIGNIFICANCE_CLASSES.items()} + + +def _class_attr_list(attr): + return [v[attr] for k, v in SIGNIFICANCE_CLASSES.items()] + + +def _volcano_source_data(data, fold_threshold=1, p_threshold=0.05): + """Prepare the minimal browser payload for the volcano/MA plot.""" + data = data.copy() + input_columns = set(data.columns) + data["log2FoldChange"] = pd.to_numeric(data["log2FoldChange"], errors="coerce") + data["padj"] = pd.to_numeric(data["padj"], errors="coerce") + is_transcript_plot = "featureID" in data.columns + if is_transcript_plot: + data.rename(columns={ + "featureID": "TXNAME", + "groupID": "GENEID", + "exonBaseMean": "mean_expression" + }, + inplace=True) + identifier_column = "TXNAME" + + else: + identifier_column = "GENEID" + data.rename(columns={ + "baseMean": "mean_expression" + }, inplace=True) + + if identifier_column not in data.columns: + raise ValueError( + "Volcano plot data must contain 'TXNAME' for transcript plots or " + "'GENEID' for gene plots." + ) + + if "mean_expression" not in data.columns: + raise ValueError( + "Volcano/MA plot data must contain 'baseMean' or 'exonBaseMean' column." + ) + + data["mean_expression"] = pd.to_numeric( + data['mean_expression'], errors="coerce" + ) + + has_gene_name = "gene_name" in data.columns + if not has_gene_name: + data["gene_name"] = "" + + data["gene_group"] = data["GENEID"] if "GENEID" in data.columns else "" + if is_transcript_plot: + if "gene_name" in data.columns: + named_group = data["gene_name"].astype(bool) + data.loc[named_group, "gene_group"] = data.loc[named_group, "gene_name"] + else: + named_group = data["GENEID"].astype(bool) + data.loc[named_group, "gene_group"] = data.loc[named_group, "GENEID"] + + data["selected_label"] = "" + + # Transformations and thresholds + data["neg_log10_padj"] = -np.log10(data["padj"]) + y_threshold = -np.log10(p_threshold) + outside_fold_threshold = data["log2FoldChange"].abs() >= fold_threshold + above_p_threshold = data["neg_log10_padj"] >= y_threshold + + # Setup initial significance classes + data["volcano_class"] = "non_significant_low_effect" + data.loc[ + outside_fold_threshold & (data["neg_log10_padj"] < y_threshold), + "volcano_class", + ] = "non_significant_high_effect" + data.loc[ + outside_fold_threshold & above_p_threshold, + "volcano_class", + ] = "significant_high_effect" + data.loc[ + ~outside_fold_threshold & above_p_threshold, + "volcano_class", + ] = "significant_low_effect" + data["color"] = data["volcano_class"].map(_class_attr("color")) + data["marker"] = data["volcano_class"].map(_class_attr("marker")) + + # All the columns we need to send to the browser for interactivity and export + source_columns = [ + identifier_column, + "gene_group", + "gene_name", + "log2FoldChange", + "mean_expression", + "marker", + "neg_log10_padj", + "padj", + "selected_label", + "color", + "volcano_class", + ] + + if is_transcript_plot: + source_columns.insert(1, "GENEID") + + source_data = {} + for column in source_columns: + values = data[column] + if column in ( + "log2FoldChange", + "mean_expression", + "padj", + "neg_log10_padj", + ): + values = values.astype(np.float32) + source_data[column] = values.reset_index(drop=True) + + export_columns = [identifier_column] + if is_transcript_plot: + export_columns.append("GENEID") + if "gene_name" in input_columns: + export_columns.append("gene_name") + export_columns.extend([ + "log2FoldChange", + "mean_expression", + "padj", + "volcano_class", + ]) + + return source_data, export_columns + + +def _tap_selection_callback_code(x_field, y_field): + """Build the shared JS tap-selection callback for volcano-style plots.""" + return f""" + const data = source.data; + const tapX = cb_obj.sx; + const tapY = cb_obj.sy; + const maxPixelDistanceSquared = 20 * 20; + function findPlotView(view, plotId) {{ + if (!view) {{ + return null; + }} + if (view.model && view.model.id === plotId && view.frame) {{ + return view; + }} + const childViews = view.child_views || view._child_views; + if (!childViews) {{ + return null; + }} + const children = childViews.values ? childViews.values() : childViews; + for (const child of children) {{ + const found = findPlotView(child, plotId); + if (found) {{ + return found; + }} + }} + return null; + }} + function selectNearestPointByPixelDistance() {{ + let plotView = Bokeh.index[plot.id]; + if (!plotView || !plotView.frame) {{ + for (const rootView of Object.values(Bokeh.index)) {{ + plotView = findPlotView(rootView, plot.id); + if (plotView) {{ + break; + }} + }} + }} + if ( + !plotView + || !Number.isFinite(tapX) + || !Number.isFinite(tapY) + ) {{ + return; + }} + const xScale = plotView.frame.x_scale; + const yScale = plotView.frame.y_scale; + let bestIndex = -1; + let bestDistance = Infinity; + let i = 0; + while (i !== data.{x_field}.length) {{ + const pointX = xScale.compute(data.{x_field}[i]); + const pointY = yScale.compute(data.{y_field}[i]); + if (!Number.isFinite(pointX) || !Number.isFinite(pointY)) {{ + i += 1; + continue; + }} + const dx = pointX - tapX; + const dy = pointY - tapY; + const distance = dx * dx + dy * dy; + if (Math.sign(bestDistance - distance) === 1) {{ + bestDistance = distance; + bestIndex = i; + }} + i += 1; + }} + if (bestIndex !== -1 && bestDistance <= maxPixelDistanceSquared) {{ + if (gene_select_toggle.active) {{ + const geneGroup = data.gene_group[bestIndex]; + const selected = []; + let j = 0; + while (j !== data.gene_group.length) {{ + if (data.gene_group[j] === geneGroup) {{ + selected.push(j); + }} + j += 1; + }} + selection_state.data.indices[0] = selected; + }} else {{ + const selected = new Set(selection_state.data.indices[0]); + selected.add(bestIndex); + selection_state.data.indices[0] = Array.from(selected).sort( + function(a, b) {{ + return a - b; + }}); + }} + update_callback.execute(source); + }} + }} + selectNearestPointByPixelDistance(); + """ + + +def volcano(data, fold_threshold=1, p_threshold=0.05): + """Build an interactive volcano plot with selection and filtering widgets.""" + source_data, original_columns = _volcano_source_data( + data, fold_threshold=fold_threshold, p_threshold=p_threshold + ) + is_transcript_plot = "TXNAME" in source_data + identifier_col = "TXNAME" if is_transcript_plot else "GENEID" + has_gene_name = "gene_name" in data.columns + counts = pd.Series(source_data["volcano_class"]).value_counts().to_dict() + x_min = float(source_data["log2FoldChange"].min()) + x_max = float(source_data["log2FoldChange"].max()) + y_max = float(source_data["neg_log10_padj"].max()) + # Exclude the top 0.1% from the initial axis scaling to prevent outliers dominating + y_init_max = ( + np.nanpercentile( + source_data["neg_log10_padj"] + .replace([np.inf, -np.inf], np.nan), 99.9) + ) + mean_min = float(source_data["mean_expression"].min()) + mean_max = float(source_data["mean_expression"].max()) + if x_min == x_max: + x_min -= 1 + x_max += 1 + if y_max <= 0: + y_max = 1 + y_init_max = 1 + if mean_min == mean_max: + mean_min *= 0.9 + mean_max *= 1.1 + # Log scale requires strictly positive lower bound; baseMean can be 0 + _pos_means = source_data["mean_expression"][source_data["mean_expression"] > 0] + mean_min_log = float(_pos_means.min()) if len(_pos_means) > 0 else 0.01 + y_threshold = -np.log10(p_threshold) + + source = ColumnDataSource(source_data) + selection_state = ColumnDataSource({"indices": [[]]}) + highlight_source = ColumnDataSource({ + "log2FoldChange": [], + "neg_log10_padj": [], + "mean_expression": [], + "gene_name": [], + }) + selected_source = ColumnDataSource({ + "source_index": [], + "owner_id": [], + "remove": [], + "GENEID": [], + "TXNAME": [], + "gene_name": [], + "log2FoldChange": [], + "padj": [], + "volcano_class": [], + }) + legend_source = ColumnDataSource({ + "symbol": _class_attr_list("symbol"), + "color": _class_attr_list("color"), + "label": _class_attr_list("label"), + "count": [counts.get(k, 0) for k in SIGNIFICANCE_CLASSES], + }) + + hover_content_formatter = CustomJSHover( + args=dict( + source=source, + identifier_col=identifier_col, + is_transcript_plot=is_transcript_plot, + has_gene_name=has_gene_name, + ), + code=""" + const n = special_vars.indices.length; + if (n > 3) { + if (special_vars.index !== special_vars.indices[0]) { return ''; } + return ( + '' + + n + ' points \u2014 zoom in to resolve' + + '' + ); + } + if (!n) { return ''; } + const data = source.data; + + function fmtPadj(v) { + const num = Number(v); + if (!isFinite(num)) { return String(v); } + if (num < 0.0001) { return num.toExponential(4); } + return num.toFixed(4); + } + + const th = 'style="text-align:right;padding-right:4px;color:#666;"'; + function row(label, value, valStyle) { + const td = valStyle + ? '' + value + '' + : '' + value + ''; + return '' + label + '' + td + ''; + } + function renderPoint(i) { + const idLabel = is_transcript_plot ? 'TX' : 'Gene'; + const idValue = is_transcript_plot + ? String(data.TXNAME ? data.TXNAME[i] || '' : '') + : (has_gene_name + ? String(data.gene_name[i] || '') + : String(data.GENEID ? data.GENEID[i] || '' : '')); + let h = ''; + h += row(idLabel, idValue, 'color:#6bb5d6;font-weight:bold;'); + if (is_transcript_plot) { + const gene = has_gene_name + ? String(data.gene_name[i] || '') + : String(data.GENEID ? data.GENEID[i] || '' : ''); + h += row('Gene', gene); + } + h += row('log2FC', Number(data.log2FoldChange[i]).toFixed(3)); + h += row('Mean expr', Number(data.mean_expression[i]).toFixed(3)); + h += row('padj', fmtPadj(data.padj[i])); + h += row('Class', String(data.volcano_class[i] || '')); + h += '
'; + return h; + } + + return renderPoint(special_vars.indices[0]); + """, + ) + hover = HoverTool( + tooltips=f"@{identifier_col}{{custom}}", + formatters={f"@{identifier_col}": hover_content_formatter}, + ) + + vol_fig = figure( + x_axis_label="log2 fold change", + y_axis_label="-log10 adjusted p-value", + height=480, + sizing_mode="stretch_width", + tools="pan,wheel_zoom,box_zoom,reset,save", + x_range=Range1d(x_min, x_max, bounds=(x_min, x_max)), + y_range=Range1d(0, y_init_max, bounds=(0, y_max)), + ) + vol_fig.toolbar.logo = None + + points_renderer = vol_fig.scatter( + x="log2FoldChange", + y="neg_log10_padj", + marker="marker", + source=source, + color="color", + size=7, + fill_alpha=0.25, + line_alpha=0.40, + selection_color="#111111", + selection_alpha=0.9, + nonselection_alpha=0.12, + ) + hover.renderers = [points_renderer] + hover.callback = CustomJS( + args=dict(selected_source=selected_source), + code=""" + const hovered = cb_data.index.indices; + if (!hovered.length || hovered.length > 3) { + selected_source.selected.indices = []; + selected_source.selected.change.emit(); + selected_source.change.emit(); + return; + } + const sourceIndex = hovered[0]; + const selectedRows = []; + let i = 0; + while (i !== selected_source.data.source_index.length) { + if (selected_source.data.source_index[i] === sourceIndex) { + selectedRows.push(i); + break; + } + i += 1; + } + selected_source.selected.indices = selectedRows; + selected_source.selected.change.emit(); + selected_source.change.emit(); + """, + ) + vol_fig.add_tools(hover) + vol_fig.circle( + x="log2FoldChange", + y="neg_log10_padj", + source=highlight_source, + size=15, + fill_alpha=0, + line_alpha=1, + line_color="#111111", + line_width=2, + ) + labels = LabelSet( + x="log2FoldChange", + y="neg_log10_padj", + text="selected_label", + source=source, + x_offset=7, + y_offset=7, + text_font_size="10px", + text_color="#111111", + ) + vol_fig.add_layout(labels) + + ma_hover = HoverTool( + tooltips=f"@{identifier_col}{{custom}}", + formatters={f"@{identifier_col}": hover_content_formatter}, + ) + ma_fig = figure( + x_axis_label="mean expression", + y_axis_label="log2 fold change", + height=480, + sizing_mode="stretch_width", + tools="pan,wheel_zoom,box_zoom,reset,save", + x_axis_type="log", + x_range=Range1d(mean_min_log, mean_max, bounds=(mean_min_log, mean_max)), + y_range=Range1d(x_min, x_max, bounds=(x_min, x_max)), + ) + ma_fig.toolbar.logo = None + ma_fig.visible = False + ma_points_renderer = ma_fig.scatter( + x="mean_expression", + y="log2FoldChange", + marker="marker", + source=source, + color="color", + size=7, + fill_alpha=0.25, + line_alpha=0.40, + selection_color="#111111", + selection_alpha=0.9, + nonselection_alpha=0.12, + ) + ma_hover.renderers = [ma_points_renderer] + ma_hover.callback = hover.callback + ma_fig.add_tools(ma_hover) + ma_fig.circle( + x="mean_expression", + y="log2FoldChange", + source=highlight_source, + size=15, + fill_alpha=0, + line_alpha=1, + line_color="#111111", + line_width=2, + ) + ma_labels = LabelSet( + x="mean_expression", + y="log2FoldChange", + text="selected_label", + source=source, + x_offset=7, + y_offset=7, + text_font_size="10px", + text_color="#111111", + ) + ma_fig.add_layout(ma_labels) + ma_zero_line = Span( + location=0, + dimension="width", + line_dash="dashed", + line_color="#444444", + line_width=1, + ) + ma_fig.add_layout(ma_zero_line) + + # Setup the selection table + table_columns = [ + TableColumn( + field="remove", + title="", + width=28, + sortable=False, + formatter=HTMLTemplateFormatter( + template=( + '' + ), + ), + ), + ] + if is_transcript_plot: + table_columns.extend([ + TableColumn(field="TXNAME", title="TXNAME"), + TableColumn(field="GENEID", title="GENEID"), + ]) + else: + table_columns.append(TableColumn(field="GENEID", title="GENEID")) + if has_gene_name: + table_columns.append(TableColumn(field="gene_name", title="gene_name")) + table_columns.extend([ + TableColumn( + field="log2FoldChange", + title="log2FC", + formatter=NumberFormatter(format="0.000"), + ), + TableColumn( + field="padj", + title="padj", + formatter=NumberFormatter(format="0.000e"), + ), + ]) + selected_table = DataTable( + source=selected_source, + columns=table_columns, + autosize_mode="force_fit", + height=310, + index_position=None, + ) + legend_table = DataTable( + source=legend_source, + columns=[ + TableColumn( + field="label", + title="Class", + width=150, + formatter=HTMLTemplateFormatter( + template=( + '<%= symbol %><%= value %>' + ), + ), + ), + TableColumn(field="count", title="n", width=50), + ], + height=130, + autosize_mode="force_fit", + index_position=None, + editable=False, + sortable=False, + styles={"padding-bottom": "4px"}, + ) + view_toggle = Toggle( + label="Toggle highlights", + button_type="default", + active=True, + width=120, + ) + gene_select_toggle = Toggle( + label="Select by gene", + button_type="default", + active=False, + width=120, + visible=is_transcript_plot, + ) + plot_mode_toggle = Button( + label="Show MA plot", + button_type="default", + width=120, + ) + left_fold_line = Span( + location=-fold_threshold, + dimension="height", + line_dash="dashed", + line_color="#444444", + line_width=1, + ) + right_fold_line = Span( + location=fold_threshold, + dimension="height", + line_dash="dashed", + line_color="#444444", + line_width=1, + ) + p_line = Span( + location=y_threshold, + dimension="width", + line_dash="dashed", + line_color="#444444", + line_width=1, + ) + vol_fig.add_layout(left_fold_line) + vol_fig.add_layout(right_fold_line) + vol_fig.add_layout(p_line) + + fold_slider = Slider( + title="Absolute log2 fold-change threshold", + value=fold_threshold, + start=0, + end=max(5, np.ceil(np.abs(source_data["log2FoldChange"]).max())), + step=0.1, + show_value=False, + sizing_mode="stretch_width", + ) + fold_input = TextInput( + title="", + value=f"{fold_threshold:.1f}", + width=80, + ) + p_slider = Slider( + title="Adjusted p-value threshold", + value=y_threshold, + start=0, + end=max(5, np.ceil(y_max)), + step=0.1, + show_value=False, + sizing_mode="stretch_width", + ) + p_input = TextInput( + title="", + value=f"{p_threshold:.3f}" if p_threshold >= 0.001 else f"{p_threshold:.2e}", + width=80, + ) + # Setup callbacks for interactivity + update_callback = CustomJS( + args=dict( + source=source, + selection_state=selection_state, + highlight_source=highlight_source, + selected_source=selected_source, + legend_source=legend_source, + fold_slider=fold_slider, + fold_input=fold_input, + p_slider=p_slider, + p_input=p_input, + left_fold_line=left_fold_line, + right_fold_line=right_fold_line, + p_line=p_line, + colors=_class_attr("color"), + markers=_class_attr("marker"), + labels=_class_attr("label"), + view_toggle=view_toggle, + identifier_col=identifier_col, + ), + code=""" + const data = source.data; + const fold = fold_slider.value; + const yThreshold = p_slider.value; + const pThreshold = Math.pow(10, -yThreshold); + fold_input.value = fold.toFixed(1); + const pThresholdLabel = pThreshold >= 0.001 + ? pThreshold.toFixed(3) + : pThreshold.toExponential(2); + p_input.value = pThresholdLabel; + const log2fc = data.log2FoldChange; + const negLog10Padj = data.neg_log10_padj; + const color = data.color; + const marker = data.marker; + const volcanoClass = data.volcano_class; + const selectedLabel = data.selected_label; + const counts = { + significant_high_effect: 0, + significant_low_effect: 0, + non_significant_high_effect: 0, + non_significant_low_effect: 0, + }; + + let i = 0; + while (i !== log2fc.length) { + const outsideFold = Math.sign(Math.abs(log2fc[i]) - fold) !== -1; + const aboveP = Math.sign(negLog10Padj[i] - yThreshold) !== -1; + let klass; + if (aboveP) { + if (outsideFold) { + klass = "significant_high_effect"; + } else { + klass = "significant_low_effect"; + } + } else { + if (outsideFold) { + klass = "non_significant_high_effect"; + } else { + klass = "non_significant_low_effect"; + } + } + volcanoClass[i] = klass; + counts[klass] += 1; + color[i] = colors[klass]; + marker[i] = markers[klass]; + i += 1; + } + const selected = Array.from( + new Set(selection_state.data.indices[0]) + ).sort(function(a, b) { + return a - b; + }); + selection_state.data.indices[0] = selected; + source.selected.indices = []; + + const selectedData = { + source_index: [], + owner_id: [], + remove: [], + GENEID: [], + TXNAME: [], + gene_name: [], + log2FoldChange: [], + padj: [], + volcano_class: [], + }; + i = 0; + while (i !== selectedLabel.length) { + selectedLabel[i] = ""; + i += 1; + } + selected.forEach(function(index) { + if (view_toggle.active) { + selectedLabel[index] = \ + data.gene_name[index] || data[identifier_col][index]; + } + selectedData.source_index.push(index); + selectedData.owner_id.push(selected_source.id); + selectedData.remove.push("x"); + selectedData.GENEID.push(data.GENEID ? data.GENEID[index] : ""); + selectedData.TXNAME.push(data.TXNAME ? data.TXNAME[index] : ""); + selectedData.gene_name.push(data.gene_name[index]); + selectedData.log2FoldChange.push(data.log2FoldChange[index]); + selectedData.padj.push(data.padj[index]); + selectedData.volcano_class.push(data.volcano_class[index]); + }); + selected_source.data = selectedData; + if (view_toggle.active) { + const highlightData = { + log2FoldChange: [], + neg_log10_padj: [], + mean_expression: [], + gene_name: [], + }; + selected.forEach(function(index) { + highlightData.log2FoldChange.push(data.log2FoldChange[index]); + highlightData.neg_log10_padj.push(data.neg_log10_padj[index]); + highlightData.mean_expression.push(data.mean_expression[index]); + highlightData.gene_name.push(data.gene_name[index]); + }); + highlight_source.data = highlightData; + } + + left_fold_line.location = -fold; + right_fold_line.location = fold; + p_line.location = yThreshold; + legend_source.data.count = [ + counts.significant_high_effect, + counts.significant_low_effect, + counts.non_significant_high_effect, + counts.non_significant_low_effect, + ]; + legend_source.change.emit(); + source.change.emit(); + selection_state.change.emit(); + highlight_source.change.emit(); + selected_source.change.emit(); + """, + ) + vol_fig.js_on_event(Tap, CustomJS( + args=dict( + plot=vol_fig, + source=source, + selection_state=selection_state, + gene_select_toggle=gene_select_toggle, + update_callback=update_callback, + ), + code=_tap_selection_callback_code("log2FoldChange", "neg_log10_padj"), + )) + ma_fig.js_on_event(Tap, CustomJS( + args=dict( + plot=ma_fig, + source=source, + selection_state=selection_state, + gene_select_toggle=gene_select_toggle, + update_callback=update_callback, + ), + code=_tap_selection_callback_code("mean_expression", "log2FoldChange"), + )) + fold_slider.js_on_change("value", update_callback) + fold_input.js_on_change("value", CustomJS( + args=dict( + fold_input=fold_input, + fold_slider=fold_slider, + update_callback=update_callback, + source=source, + ), + code=""" + const rawValue = fold_input.value.trim(); + if (!rawValue) { + return; + } + const parsed = Number(rawValue); + if (!Number.isFinite(parsed) || parsed < 0) { + return; + } + fold_slider.value = parsed; + update_callback.execute(source); + """, + )) + p_slider.js_on_change("value", update_callback) + p_input.js_on_change("value", CustomJS( + args=dict( + p_input=p_input, + p_slider=p_slider, + update_callback=update_callback, + source=source, + ), + code=""" + const rawValue = p_input.value.trim(); + if (!rawValue) { + return; + } + const parsed = Number(rawValue); + if (!Number.isFinite(parsed) || parsed <= 0 || parsed > 1) { + return; + } + p_slider.value = -Math.log10(parsed); + update_callback.execute(source); + """, + )) + view_toggle.js_on_change("active", CustomJS( + args=dict( + source=source, + selection_state=selection_state, + highlight_source=highlight_source, + view_toggle=view_toggle, + identifier_col=identifier_col, + ), + code=""" + const selectedLabel = source.data.selected_label; + let i = 0; + while (i !== selectedLabel.length) { + selectedLabel[i] = ""; + i += 1; + } + source.selected.indices = []; + if (view_toggle.active) { + const highlightData = { + log2FoldChange: [], + neg_log10_padj: [], + mean_expression: [], + gene_name: [], + }; + selection_state.data.indices[0].forEach(function(index) { + selectedLabel[index] = \ + source.data.gene_name[index] || source.data[identifier_col][index]; + highlightData.log2FoldChange.push(source.data.log2FoldChange[index]); + highlightData.neg_log10_padj.push(source.data.neg_log10_padj[index]); + highlightData.mean_expression.push(source.data.mean_expression[index]); + highlightData.gene_name.push(source.data.gene_name[index]); + }); + highlight_source.data = highlightData; + } else { + highlight_source.data = { + log2FoldChange: [], + neg_log10_padj: [], + mean_expression: [], + gene_name: [], + }; + } + source.change.emit(); + highlight_source.change.emit(); + """, + )) + gene_select_toggle.js_on_change("active", CustomJS( + args=dict( + source=source, + selection_state=selection_state, + gene_select_toggle=gene_select_toggle, + update_callback=update_callback, + ), + code=""" + if (!gene_select_toggle.active) { + return; + } + const selected = selection_state.data.indices[0]; + if (!selected.length) { + return; + } + const geneGroup = source.data.gene_group[selected[0]]; + const groupIndices = []; + let i = 0; + while (i !== source.data.gene_group.length) { + if (source.data.gene_group[i] === geneGroup) { + groupIndices.push(i); + } + i += 1; + } + selection_state.data.indices[0] = groupIndices; + update_callback.execute(source); + """, + )) + plot_mode_toggle.js_on_click(CustomJS( + args=dict( + volcano_fig=vol_fig, + ma_fig=ma_fig, + plot_mode_toggle=plot_mode_toggle, + ), + code=""" + const showingVolcano = volcano_fig.visible; + volcano_fig.visible = !showingVolcano; + ma_fig.visible = showingVolcano; + plot_mode_toggle.label = showingVolcano + ? "Show volcano plot" + : "Show MA plot"; + """, + )) + selected_source.selected.js_on_change("indices", CustomJS( + args=dict( + source=source, + selected_source=selected_source, + highlight_source=highlight_source, + view_toggle=view_toggle, + ), + code=""" + if (view_toggle.active) { + return; + } + const rows = selected_source.selected.indices; + if (rows.length) { + const row = rows[0]; + const sourceIndex = selected_source.data.source_index[row]; + highlight_source.data = { + log2FoldChange: [source.data.log2FoldChange[sourceIndex]], + neg_log10_padj: [source.data.neg_log10_padj[sourceIndex]], + mean_expression: [source.data.mean_expression[sourceIndex]], + gene_name: [source.data.gene_name[sourceIndex]], + }; + } else { + highlight_source.data = { + log2FoldChange: [], + neg_log10_padj: [], + mean_expression: [], + gene_name: [], + }; + } + highlight_source.change.emit(); + """, + )) + selected_source.js_on_change("data", CustomJS( + args=dict( + source=source, + selection_state=selection_state, + selected_source=selected_source, + highlight_source=highlight_source, + gene_select_toggle=gene_select_toggle, + update_callback=update_callback, + ), + code=""" + const listenerKey = "volcano_delete_listener_" + selected_source.id; + if (!window[listenerKey]) { + window[listenerKey] = true; + document.addEventListener("click", function(event) { + const path = event.composedPath ? event.composedPath() : []; + const deleteButton = path.find(function(element) { + return element.classList + && element.classList.contains("volcano-row-delete"); + }); + if (!deleteButton) { + return; + } + if (deleteButton.dataset.selectedSourceId !== selected_source.id) { + return; + } + event.preventDefault(); + event.stopPropagation(); + + const sourceIndex = Number(deleteButton.dataset.sourceIndex); + if (gene_select_toggle.active) { + selection_state.data.indices[0] = []; + } else { + const kept = []; + const selected = selection_state.data.indices[0]; + let i = 0; + while (i !== selected.length) { + const index = selected[i]; + if (index !== sourceIndex) { + kept.push(index); + } + i += 1; + } + selection_state.data.indices[0] = kept; + } + selected_source.selected.indices = []; + highlight_source.data = { + log2FoldChange: [], + neg_log10_padj: [], + mean_expression: [], + gene_name: [], + }; + update_callback.execute(source); + }); + } + """, + )) + + search_input = AutocompleteInput( + placeholder=( + "Search by gene/transcript…" if is_transcript_plot else "Search by gene…"), + completions=[], + min_characters=1, + width=240, + ) + search_input.js_on_change("value_input", CustomJS( + args=dict( + source=source, + search_input=search_input, + ), + code=""" + if (search_input.completions.length) { + return; + } + const data = source.data; + const completions = new Set(); + const fields = [data.gene_name, data.GENEID]; + if (data.TXNAME) { + fields.push(data.TXNAME); + } + fields.forEach(function(field) { + let i = 0; + while (i !== field.length) { + const value = field[i]; + if (value) { + completions.add(String(value)); + } + i += 1; + } + }); + search_input.completions = Array.from(completions).sort(); + """, + )) + search_callback = CustomJS( + args=dict( + source=source, + selection_state=selection_state, + search_input=search_input, + update_callback=update_callback, + ), + code=""" + const query = search_input.value.trim().toLowerCase(); + if (!query) { + return; + } + const data = source.data; + const currentSelected = new Set(selection_state.data.indices[0]); + let i = 0; + while (i !== data.gene_name.length) { + const fields = [data.gene_name, data.GENEID]; + if (data.TXNAME) { + fields.push(data.TXNAME); + } + let matched = false; + let fi = 0; + while (fi !== fields.length) { + const val = fields[fi][i]; + if (val && val.toLowerCase() === query) { + matched = true; + break; + } + fi += 1; + } + if (matched) { + currentSelected.add(i); + } + i += 1; + } + selection_state.data.indices[0] = Array.from(currentSelected).sort( + function(a, b) { return a - b; } + ); + search_input.value = ""; + update_callback.execute(source); + """, + ) + search_input.js_on_change("value", search_callback) + + select_significant_button = Button( + label="Select significant", + button_type="primary", + width=150, + ) + select_significant_button.js_on_click(CustomJS( + args=dict( + source=source, + selection_state=selection_state, + fold_slider=fold_slider, + p_slider=p_slider, + update_callback=update_callback, + ), + code=""" + const data = source.data; + const fold = fold_slider.value; + const yThreshold = p_slider.value; + const passing = []; + let i = 0; + while (i !== data.log2FoldChange.length) { + const outsideFold = Math.abs(data.log2FoldChange[i]) >= fold; + const aboveP = data.neg_log10_padj[i] >= yThreshold; + if (outsideFold && aboveP) { + passing.push(i); + } + i += 1; + } + selection_state.data.indices[0] = passing; + update_callback.execute(source); + """, + )) + + clear_button = Button( + label="Clear selection", + button_type="default", + width=130, + ) + clear_button.js_on_click(CustomJS( + args=dict( + source=source, + selection_state=selection_state, + update_callback=update_callback, + ), + code=""" + selection_state.data.indices[0] = []; + update_callback.execute(source); + """, + )) + + export_button = Button( + label="\u2B07 TSV", + button_type="light", + width=55 + ) + + export_button.js_on_click(CustomJS( + args=dict( + source=source, + selection_state=selection_state, + original_columns=original_columns, + filename=( + "selected_transcripts.tsv" + if is_transcript_plot + else "selected_genes.tsv" + ), + ), + code=""" + const indices = selection_state.data.indices[0]; + if (!indices.length) { + return; + } + const cols = original_columns; + const data = source.data; + const rows = [cols.join('\\t')]; + indices.forEach(function(i) { + rows.push(cols.map(function(col) { + const v = data[col] !== undefined ? data[col][i] : ''; + return v !== null && v !== undefined ? String(v) : ''; + }).join('\\t')); + }); + const blob = new Blob([rows.join('\\n') + '\\n'], { + type: 'text/tab-separated-values', + }); + const url = URL.createObjectURL(blob); + const link = document.createElement('a'); + link.href = url; + link.download = filename; + document.body.appendChild(link); + link.click(); + document.body.removeChild(link); + URL.revokeObjectURL(url); + """, + )) + + volcano_ma_plot = BokehPlot() + controls = bokeh_row( + fold_slider, fold_input, p_slider, p_input, sizing_mode="stretch_width") + toggle_row = bokeh_row( + plot_mode_toggle, view_toggle, + gene_select_toggle, select_significant_button + ) + volcano_ma_plot._fig = bokeh_column( + bokeh_row(controls, sizing_mode="stretch_width"), + bokeh_row( + bokeh_column( + toggle_row, + vol_fig, + ma_fig, + sizing_mode="stretch_width", + ), + sizing_mode="stretch_width", + ), + sizing_mode="stretch_width", + ) + + classes_table = BokehPlot() + classes_table._fig = bokeh_column( + legend_table, + sizing_mode="stretch_width", + styles={ + "padding-top": "4px", + "padding-left": "35px", + "padding-right": "20px" + }, + ) + selected_plot = BokehPlot() + selected_plot._fig = bokeh_column( + bokeh_row( + bokeh_row( + search_input, clear_button, export_button) + ), + selected_table, sizing_mode="stretch_width", + styles={ + "padding-top": "25px", + }) + return volcano_ma_plot, classes_table, selected_plot