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