From 759cb9d868eab4678c9db09b7adedb9456efc272 Mon Sep 17 00:00:00 2001 From: Kiah McIntosh Date: Thu, 21 May 2026 16:19:18 +0000 Subject: [PATCH] Skip cohort processes when n=1 [CW-7208] --- .gitlab-ci.yml | 13 +- bin/workflow_glue/report.py | 73 +++++--- bin/workflow_glue/tests/common/test_report.py | 165 +++++++++++++++++- main.nf | 4 +- subworkflows/transcriptome.nf | 1 + 5 files changed, 216 insertions(+), 40 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index db34764..c196ad6 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -134,15 +134,14 @@ docker-run: jq -e '.contrasts["condition_treated_vs_control"].dge_status | IN("SUCCESS", "FAILED")' ${CI_PROJECT_NAME}/de_analysis/de_qc_stats.json >/dev/null && jq -e '.contrasts["condition_treated_vs_control"].dtu_status | IN("SUCCESS", "FAILED")' ${CI_PROJECT_NAME}/de_analysis/de_qc_stats.json >/dev/null - # Smoke: quick discover-mode sanity check for core cohort and per-sample outputs. + # Smoke: quick discover-mode sanity check for core cohort outputs. - if: $MATRIX_NAME == "smoke_discover" variables: NF_BEFORE_SCRIPT: ":" NF_WORKFLOW_OPTS: "--fastq test_data/smoke/reads.fastq --sample sampleA --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf" AFTER_NEXTFLOW_CMD: > - test -f ${CI_PROJECT_NAME}/cohort/transcripts.gtf && - test -f ${CI_PROJECT_NAME}/cohort/cohort.transcriptome.fa && - test -f ${CI_PROJECT_NAME}/samples/sampleA/transcripts.gtf + test -f ${CI_PROJECT_NAME}/samples/sampleA/transcripts.gtf && + test -f ${CI_PROJECT_NAME}/samples/sampleA/sampleA.transcriptome.fa # Smoke: fixed-annotation path sanity check for quantification outputs. - if: $MATRIX_NAME == "smoke_fixed" @@ -150,8 +149,8 @@ docker-run: NF_BEFORE_SCRIPT: ":" NF_WORKFLOW_OPTS: "--fastq test_data/smoke/reads.fastq --sample sampleA --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf --transcriptome_mode fixed_annotation" AFTER_NEXTFLOW_CMD: > - test -f ${CI_PROJECT_NAME}/cohort/transcripts.gtf && - test -f ${CI_PROJECT_NAME}/cohort/transcript_counts.tsv + test -f ${CI_PROJECT_NAME}/samples/sampleA/transcripts.gtf && + test -f ${CI_PROJECT_NAME}/samples/sampleA/transcript_counts.tsv # Smoke: direct-RNA alignment profile and downstream SQANTI output presence. - if: $MATRIX_NAME == "smoke_direct_rna" @@ -160,7 +159,7 @@ docker-run: NF_WORKFLOW_OPTS: "--fastq test_data/smoke/reads.fastq --sample sampleA --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf --direct_rna" AFTER_NEXTFLOW_CMD: > test -f ${CI_PROJECT_NAME}/cohort/alignments/sampleA/reads.bam && - test -f ${CI_PROJECT_NAME}/cohort/sqanti_cohort/classification_summary.tsv + test -f ${CI_PROJECT_NAME}/samples/sampleA/sampleA_sqanti/classification_summary.tsv # Smoke: end-to-end DE/DTU wiring and expected contrast output files. - if: $MATRIX_NAME == "smoke_de" diff --git a/bin/workflow_glue/report.py b/bin/workflow_glue/report.py index fe5dbef..a1510c8 100644 --- a/bin/workflow_glue/report.py +++ b/bin/workflow_glue/report.py @@ -68,9 +68,9 @@ def _format_ratio_value(value): return f"{numeric:.2f}x" -def _cohort_summary(cohort_dir): - """Return cohort-level metrics and transcript class counts DataFrames.""" - tx_meta = _read_table(Path(cohort_dir) / "transcript_metadata.tsv") +def _transcriptome_summary(transcriptome_dir): + """Return transcriptome metrics and transcript class counts DataFrames.""" + tx_meta = _read_table(Path(transcriptome_dir) / "transcript_metadata.tsv") if tx_meta is None: return None, None @@ -95,7 +95,10 @@ 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()): + samples_path = Path(samples_dir) + if not samples_path.exists() or not samples_path.is_dir(): + return summaries + for sample_dir in sorted(samples_path.iterdir()): if not sample_dir.is_dir(): continue tx_meta = _read_table(sample_dir / "transcript_metadata.tsv") @@ -143,9 +146,9 @@ def _contrast_results(de_dir, filename, n=None): return tables -def _load_bambu_qc(cohort_dir): +def _load_bambu_qc(bambu_dir): """Load bambu QC statistics JSON.""" - qc_file = Path(cohort_dir) / "bambu_qc_stats.json" + qc_file = Path(bambu_dir) / "bambu_qc_stats.json" if qc_file.exists(): with open(qc_file) as f: return json.load(f) @@ -345,8 +348,6 @@ def main(args): use_index=False, ) - # Load bambu QC statistics - bambu_qc = _load_bambu_qc(args.cohort_dir) annotation_reference_summary = _load_annotation_reference_summary(args.cohort_dir) if annotation_reference_summary: @@ -437,6 +438,16 @@ def main(args): h4("Unstranded Annotation Examples") pre("\n".join(examples)) + # Setup for using cohort or single sample bambu results + is_single_sample = len(metadata) == 1 + primary_label = metadata[0]["alias"] if is_single_sample else "Cohort" + + # Load bambu QC statistics + bambu_dir = ( + Path(args.samples_dir) / metadata[0]["alias"] if is_single_sample + else args.cohort_dir + ) + bambu_qc = _load_bambu_qc(bambu_dir) # Add Bambu QC section with warnings if bambu_qc: with report.add_section("Bambu Quality Control", "Bambu QC"): @@ -556,38 +567,48 @@ def main(args): use_index=False, ) - with report.add_section("Cohort transcriptome", "Cohort"): - cohort_metrics, cohort_classes = _cohort_summary(args.cohort_dir) - if cohort_metrics is not None: + with report.add_section( + f"{primary_label} transcriptome", + f"{primary_label} transcriptome" + ): + transcriptome_metrics, transcriptome_classes = _transcriptome_summary(bambu_dir) + if transcriptome_metrics is not None: DataTable.from_pandas( - cohort_metrics, + transcriptome_metrics, paging=False, searchable=False, use_index=False, ) - if cohort_classes is not None: + if transcriptome_classes is not None: DataTable.from_pandas( - cohort_classes, + transcriptome_classes, paging=False, searchable=False, use_index=False, ) - tx_counts = _read_table(Path(args.cohort_dir) / "transcript_counts.tsv") + tx_counts = _read_table(Path(bambu_dir) / "transcript_counts.tsv") if tx_counts is not None and not tx_counts.empty: - p("Top transcript rows from the cohort abundance table.") + p( + "Top transcript rows from the " + f"{'sample' if is_single_sample else 'cohort'} abundance table." + ) DataTable.from_pandas(tx_counts.head(20), use_index=False) - with report.add_section("Per-sample transcriptomes", "Per sample"): - tabs = Tabs() - for sample, summary_df in _sample_summaries(args.samples_dir).items(): - with tabs.add_tab(sample): - DataTable.from_pandas( - summary_df, - paging=False, - searchable=False, - use_index=False, - ) + if not is_single_sample: + with report.add_section( + "Per-sample transcriptomes", + "Per-sample transcriptomes" + ): + tabs = Tabs() + for sample, summary_df in _sample_summaries(args.samples_dir).items(): + with tabs.add_tab(sample): + DataTable.from_pandas( + summary_df, + paging=False, + searchable=False, + use_index=False, + ) if args.alignment_stats_dir and Path(args.alignment_stats_dir).exists(): with report.add_section("Alignment statistics", "Alignments"): diff --git a/bin/workflow_glue/tests/common/test_report.py b/bin/workflow_glue/tests/common/test_report.py index bda75f0..d7a51ef 100644 --- a/bin/workflow_glue/tests/common/test_report.py +++ b/bin/workflow_glue/tests/common/test_report.py @@ -277,8 +277,12 @@ def test_report_main_handles_degenerate_bambu_qc_and_read_summary( } ), ) + samples = tmp_path / "samples" + samples.mkdir() + sample_a = samples / "sampleA" + sample_a.mkdir() _write( - cohort / "bambu_qc_stats.json", + sample_a / "bambu_qc_stats.json", json.dumps( { "samples": 1, @@ -298,10 +302,6 @@ def test_report_main_handles_degenerate_bambu_qc_and_read_summary( ), ) - samples = tmp_path / "samples" - samples.mkdir() - (samples / "OPTIONAL_FILE").touch() - sqanti = tmp_path / "sqanti" sqanti.mkdir() (sqanti / "OPTIONAL_FILE").touch() @@ -346,6 +346,161 @@ def test_report_main_handles_degenerate_bambu_qc_and_read_summary( ) +def test_report_main_uses_cohort_bambu_qc_for_multi_sample_inputs( + monkeypatch, + tmp_path, +): + """Multi-sample runs should render bambu QC from cohort-level outputs.""" + tables = [] + + monkeypatch.setattr(report.labs, "LabsReport", _FakeReport) + monkeypatch.setattr(report, "Tabs", _FakeTabs) + monkeypatch.setattr(report, "p", lambda *args, **kwargs: None) + monkeypatch.setattr(report, "pre", lambda *args, **kwargs: None) + monkeypatch.setattr(report, "_create_warning_banner", lambda *args, **kwargs: None) + monkeypatch.setattr(report.fastcat, "SeqSummary", lambda *args, **kwargs: None) + monkeypatch.setattr( + report.DataTable, + "from_pandas", + staticmethod(lambda table, *args, **kwargs: tables.append(table.copy())), + ) + + metadata = _write( + tmp_path / "metadata.json", + json.dumps( + [ + {"alias": "sampleA", "has_stats": False}, + {"alias": "sampleB", "has_stats": False}, + ] + ), + ) + params = _write(tmp_path / "params.json", "{}") + versions = tmp_path / "versions" + versions.mkdir() + _write(versions / "versions.txt", "tool,1.0\n") + + cohort = tmp_path / "cohort" + cohort.mkdir() + _write( + cohort / "bambu_qc_stats.json", + json.dumps( + { + "samples": 2, + "library_sizes": {"sampleA": 1200, "sampleB": 900}, + "min_library_size": 900, + "max_library_size": 1200, + "median_library_size": 1050, + "library_size_ratio": 1.3333, + "total_transcripts_before_filter": 100, + "total_transcripts_after_filter": 80, + "transcripts_filtered": 20, + "median_transcripts_detected": 70, + "total_genes_after_filter": 60, + "transcriptome_mode": "discover", + "ndr_used": 0.1, + } + ), + ) + + samples = tmp_path / "samples" + samples.mkdir() + sample_a = samples / "sampleA" + sample_a.mkdir() + _write( + sample_a / "bambu_qc_stats.json", + json.dumps({"samples": 1, "library_sizes": {"sampleA": 5}}), + ) + _write( + sample_a / "transcript_metadata.tsv", + "TXNAME\tGENEID\n" + "tx1\tgene1\n" + "tx2\tgene1\n" + "tx3\tgene2\n", + ) + sample_b = samples / "sampleB" + sample_b.mkdir() + _write( + sample_b / "bambu_qc_stats.json", + json.dumps({"samples": 1, "library_sizes": {"sampleB": 7}}), + ) + _write( + sample_b / "transcript_metadata.tsv", + "TXNAME\tGENEID\n" + "txA\tgeneA\n" + "txB\tgeneB\n", + ) + + sqanti = tmp_path / "sqanti" + sqanti.mkdir() + (sqanti / "OPTIONAL_FILE").touch() + + alignment_stats = tmp_path / "alignment_stats" + alignment_stats.mkdir() + (alignment_stats / "OPTIONAL_FILE").touch() + + out_report = tmp_path / "wf-transcriptomes-report.html" + args = report.argparser().parse_args( + [ + str(out_report), + "--metadata", + str(metadata), + "--alignment_stats_dir", + str(alignment_stats), + "--cohort_dir", + str(cohort), + "--samples_dir", + str(samples), + "--sqanti_dir", + str(sqanti), + "--versions", + str(versions), + "--params", + str(params), + ] + ) + + report.main(args) + + assert out_report.exists() + assert any( + "Samples analyzed" in table.to_string() and "2" in table.to_string() + for table in tables + ) + assert any( + "Library size ratio (max/min)" in table.to_string() + and "1.33x" in table.to_string() + for table in tables + ) + assert any( + "Sample" in table.columns + and "Library Size" in table.columns + and {"sampleA", "sampleB"}.issubset(set(table["Sample"].tolist())) + and "1,200" in table.to_string() + and "900" in table.to_string() + for table in tables + ) + per_sample_metric_tables = [ + table + for table in tables + if list(table.columns) == ["Metric", "Value"] + and set(table["Metric"].tolist()) == {"Transcripts", "Genes"} + ] + assert any( + set(zip(table["Metric"], table["Value"])) == { + ("Transcripts", 3), + ("Genes", 2), + } + for table in per_sample_metric_tables + ) + assert any( + set(zip(table["Metric"], table["Value"])) == { + ("Transcripts", 2), + ("Genes", 2), + } + for table in per_sample_metric_tables + ) + + def test_report_main_renders_statistical_methods_and_warnings( monkeypatch, tmp_path, diff --git a/main.nf b/main.nf index 769e9c8..66d288c 100644 --- a/main.nf +++ b/main.nf @@ -39,7 +39,7 @@ process makeReport { tuple val(metadata), path(stats, stageAs: "stats_*") path "versions/*" path "params.json" - path cohort_dir, stageAs: "cohort" + path cohort_dir, stageAs: "cohort/*" path sample_dirs, stageAs: "samples/*" path sqanti_dirs, stageAs: "sqanti/*" path de_files @@ -198,7 +198,7 @@ workflow pipeline { report_input, software_versions, workflow_params, - transcriptome.joint_dir, + transcriptome.joint_dir.ifEmpty(OPTIONAL_FILE), sample_dirs_for_report, sqanti_dirs_for_report, de_dir, diff --git a/subworkflows/transcriptome.nf b/subworkflows/transcriptome.nf index 221b2fd..d139918 100644 --- a/subworkflows/transcriptome.nf +++ b/subworkflows/transcriptome.nf @@ -222,6 +222,7 @@ workflow transcriptome_analysis { joint_discover = runJointBambuDiscover( alignments .toSortedList { a, b -> a[0].alias <=> b[0].alias } + .filter { rows -> rows.size() > 1 } .map { rows -> // transform [meta, bam, bai] rows to // [meta, [alias1...aliasN], [bam1...bamN], [bai1...baiN], sample_sheet]