diff --git a/CHANGELOG.md b/CHANGELOG.md
index ceabb35..883f32f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,9 +13,9 @@ Users of wf-transcriptomes v2.0.0 who have encountered issues during discovery a
### Fixed
- "Error in full_join" encountered during `runPerSampleBambuQuant` when all read classes have no compatible transcript assignment. An empty quant table is correctly emitted instead.
- "unable to find an inherited method for function 'rowData'" encountered during `runJointBambuDiscover` when providing many samples. The workflow now correctly handles data spilled to disk by bambu discover.
+- Volcano plot class counts incorrect when `log2FoldChange` or `padj` columns contained NA values.
- Adjusted p-values below 0.001 in the volcano selection table are now shown in scientific notation instead of being rounded to 0.000.
-
## [v2.0.0]
This release refreshes `wf-transcriptomes` around a new reference-guided transcriptomics workflow built on `bambu`, with `SQANTI3` transcript classification and QC, `DESeq2` for differential gene expression, `DEXSeq` for differential transcript usage, and per-sample modified base summarisation with `modkit` when modification tags are present in aligned BAMs.
diff --git a/bin/workflow_glue/report.py b/bin/workflow_glue/report.py
index ddc7ea4..e8299fa 100644
--- a/bin/workflow_glue/report.py
+++ b/bin/workflow_glue/report.py
@@ -548,7 +548,7 @@ def _heatmap_style():
grid-template-columns: 1fr;
}
}
- .clustering-info {
+ .info-text {
font-size: 11px;
}"""
@@ -1384,7 +1384,7 @@ def main(args):
EZChart(hierarchical_result.heatmap, width="100%")
EZChart(hierarchical_result.pca, width="100%")
EZChart(hierarchical_result.distance, width="100%")
- with div(cls="clustering-info"):
+ with div(cls="info-text"):
br()
clustering_info('gene')
tabs = Tabs()
@@ -1422,19 +1422,22 @@ def main(args):
_round_de_table(table.head(args.de_table_size)),
use_index=False
)
- with div(cls="clustering-info"):
+ with div(cls="info-text"):
raw(
f"Table showing the top {args.de_table_size} genes sorted "
"by adjusted p-value.
"
)
h3("Gene expression volcano Plot")
- gn_vol, gn_class_table, gn_selected_table = volcano(table)
+ gn_vol, gn_class_table, gn_selected_table, vol_text = volcano(table)
with div(_class="volcano-plot-grid"):
EZChart(gn_vol, width="100%", height="550")
with div(_class="volcano-side-panel"):
EZChart(gn_class_table, width="100%", height="auto")
EZChart(gn_selected_table, width="100%", height="auto")
+ if vol_text is not None:
+ with div(cls="info-text"):
+ raw(vol_text)
with report.add_section("Differential transcript usage", "DTU"):
p(
@@ -1467,7 +1470,7 @@ def main(args):
EZChart(hierarchical_result.heatmap, width="100%")
EZChart(hierarchical_result.pca, width="100%")
EZChart(hierarchical_result.distance, width="100%")
- with div(cls="clustering-info"):
+ with div(cls="info-text"):
br()
clustering_info('transcript')
@@ -1511,19 +1514,24 @@ def main(args):
_round_de_table(dtu_table.head(args.de_table_size)),
use_index=False
)
- with div(cls="clustering-info"):
+ with div(cls="info-text"):
raw(
f"Table showing the top {args.de_table_size} "
"transcripts sorted by adjusted p-value.
"
)
h3("Transcript expression volcano Plot")
- tr_vol, tr_class_table, tr_selected_table = volcano(dtu_table)
+ (
+ tr_vol, tr_class_table, tr_selected_table, vol_text
+ ) = volcano(dtu_table)
with div(_class="volcano-plot-grid"):
EZChart(tr_vol, width="100%", height="550")
with div(_class="volcano-side-panel"):
EZChart(tr_class_table, width="100%", height="auto")
EZChart(tr_selected_table, width="100%", height="auto")
+ if vol_text is not None:
+ with div(cls="info-text"):
+ raw(vol_text)
else:
p("No DTU results available for this contrast.")
diff --git a/bin/workflow_glue/volcano.py b/bin/workflow_glue/volcano.py
index 9597b94..016c29a 100644
--- a/bin/workflow_glue/volcano.py
+++ b/bin/workflow_glue/volcano.py
@@ -67,6 +67,10 @@ def _volcano_source_data(data, fold_threshold=1, p_threshold=0.05):
input_columns = set(data.columns)
data["log2FoldChange"] = pd.to_numeric(data["log2FoldChange"], errors="coerce")
data["padj"] = pd.to_numeric(data["padj"], errors="coerce")
+
+ data.replace([np.inf, -np.inf], np.nan, inplace=True)
+ data.dropna(subset=["log2FoldChange", "padj"], inplace=True)
+
is_transcript_plot = "featureID" in data.columns
if is_transcript_plot:
data.rename(columns={
@@ -267,6 +271,9 @@ def _tap_selection_callback_code(x_field, y_field):
def volcano(data, fold_threshold=1, p_threshold=0.05):
"""Build an interactive volcano plot with selection and filtering widgets."""
+ input_size = len(data)
+ source_data = data.dropna(subset=["log2FoldChange", "padj"])
+ n_filtered = input_size - len(source_data)
source_data, original_columns = _volcano_source_data(
data, fold_threshold=fold_threshold, p_threshold=p_threshold
)
@@ -1321,4 +1328,15 @@ def volcano(data, fold_threshold=1, p_threshold=0.05):
}
""",
))
- return volcano_ma_plot, classes_table, selected_plot
+ message = None
+ if n_filtered > 0:
+ tx_file = (
+ "results_dtu_transcript.tsv" if is_transcript_plot else "results_dge.tsv")
+ message = (
+ f"{n_filtered} features were omitted from the volcano plot because "
+ "they do not have plottable log2 fold-change or adjusted p-values. "
+ "
This is expected for low-information features filtered by the "
+ "differential analysis method. Full results are available in "
+ f"{tx_file}."
+ )
+ return volcano_ma_plot, classes_table, selected_plot, message