Merge branch 'cw-7238' into 'dev'

SQANTI3 results in report [CW-7238]

See merge request epi2melabs/workflows/wf-transcriptomes!271

.gitlab-ci.yml

@@ -160,7 +160,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}/samples/sampleA/alignment/reads.bam &&
-                  test -f ${CI_PROJECT_NAME}/samples/sampleA/sampleA_sqanti/classification_summary.tsv
+                  test -f ${CI_PROJECT_NAME}/samples/sampleA/sqanti/classification_summary.tsv
 
         # Smoke: end-to-end DE/DTU wiring and expected contrast output files.
         - if: $MATRIX_NAME == "smoke_de"

README.md

@@ -220,7 +220,7 @@ Each sample is aligned to the supplied reference genome with
 [`minimap2`](https://github.com/lh3/minimap2), then sorted and indexed with
 [`samtools`](https://www.htslib.org/). The aligned BAMs under
 `samples/<alias>/alignment/` are the main alignment files used for transcriptome
-analysis, optional `SQANTI3` QC, and optional IGV viewing.
+analysis, optional [`SQANTI3`](https://github.com/conesalab/SQANTI3) QC, and optional IGV viewing.
 
 ### 4. Cohort transcriptome construction
 
@@ -245,8 +245,8 @@ for DE/DTU.
 Transcript FASTA files are derived from GTF plus genome using `gffread`.
 When `--skip_sqanti` is not set, `SQANTI3` classifies the cohort and per-sample
 transcriptomes and produces structural QC summaries. The cohort `SQANTI3`
-results live under `cohort/sqanti_cohort/`, while per-sample `SQANTI3`
+results live under `cohort/sqanti/`, while per-sample `SQANTI3`
-directories are published under `samples/<alias>/<alias>_sqanti/`.
+directories are published under `samples/<alias>/sqanti/`.
 
 ### 7. Optional DE and DTU analysis
 
@@ -368,13 +368,13 @@ Output files may be aggregated including information for all samples or provided
 | Cohort transcript counts | cohort/transcript_counts.tsv | Transcript-level count matrix produced by bambu. | aggregated |
 | Cohort gene counts | cohort/gene_counts.tsv | Gene-level count matrix derived from bambu output. | aggregated |
 | Cohort transcript metadata | cohort/transcript_metadata.tsv | Transcript annotations and bambu transcript classes for the cohort model. | aggregated |
-| Cohort SQANTI3 summary | cohort/sqanti_cohort/classification_summary.tsv | SQANTI3 classification summary for the cohort transcriptome when SQANTI3 QC is enabled. | aggregated |
+| Cohort SQANTI3 summary | cohort/sqanti/classification_summary.tsv | SQANTI3 classification summary for the cohort transcriptome when SQANTI3 QC is enabled. | aggregated |
 | Per-sample transcriptome GTF | samples/{{ alias }}/transcripts.gtf | Independent bambu transcript model for an individual sample. | per-sample |
 | Per-sample transcriptome FASTA | samples/{{ alias }}/{{ alias }}.transcriptome.fa | Transcript sequences derived from the per-sample GTF. | per-sample |
 | Per-sample transcript counts | samples/{{ alias }}/transcript_counts.tsv | Transcript-level abundance estimates for the per-sample bambu model. | per-sample |
 | Per-sample gene counts | samples/{{ alias }}/gene_counts.tsv | Gene-level abundance estimates for the per-sample bambu model. | per-sample |
 | Per-sample transcript metadata | samples/{{ alias }}/transcript_metadata.tsv | Transcript annotations and bambu transcript classes for the per-sample model. | per-sample |
-| Per-sample SQANTI3 summary | samples/{{ alias }}/{{ alias }}_sqanti/classification_summary.tsv | SQANTI3 classification summary for the per-sample transcriptome when SQANTI3 QC is enabled. | per-sample |
+| Per-sample SQANTI3 summary | samples/{{ alias }}/sqanti/classification_summary.tsv | SQANTI3 classification summary for the per-sample transcriptome when SQANTI3 QC is enabled. | per-sample |
 | Differential gene expression results | de_analysis/{{ contrast }}/results_dge.tsv | DESeq2 gene-level differential expression results for one contrast. | aggregated |
 | Differential gene expression plots | de_analysis/{{ contrast }}/results_dge.pdf | PDF plots generated during DESeq2 analysis for one contrast. | aggregated |
 | Differential transcript usage results | de_analysis/{{ contrast }}/results_dtu_transcript.tsv | Transcript-level DTU results for one contrast. | aggregated |

bin/workflow_glue/report.py

@@ -2,10 +2,13 @@
 
 import json
 import math
+import os
 from pathlib import Path
 
 from bokeh.resources import INLINE as BOKEH_INLINE
-from dominate.tags import br, div, h3, h4, p, pre, script, strong, style as dom_style
+from dominate.tags import (
+    br, div, h3, h4, p, pre, script, small, strong, style as dom_style
+)
 from dominate.util import raw
 from ezcharts.components import fastcat
 from ezcharts.components.ezchart import EZChart
@@ -21,6 +24,35 @@ from .util import get_named_logger, wf_parser  # noqa: ABS101
 from .volcano import volcano  # noqa: ABS101
 
 
+classification_categories = {
+    "Full splice match": (
+        "Reference and query isoforms have the same number of exons and "
+        "all internal junctions agree."
+    ),
+    "Incomplete splice match": (
+        "Query isoform has fewer 5′ exons than the reference, with "
+        "matching internal junctions."
+    ),
+    "Novel in catalog": (
+        "No full or incomplete splice match, but uses a "
+        "combination of known donor/acceptor splice sites."
+    ),
+    "Novel not in catalog": (
+        "No full or incomplete splice match, with at least "
+        "one unannotated donor or acceptor splice site."
+    ),
+    "Antisense": (
+        "No same-strand reference overlap, but antisense to an "
+        "annotated gene."
+    ),
+    "Genic intron": (
+        "Query isoform is fully contained within an annotated intron."
+    ),
+    "Genic": "Query isoform overlaps introns and exons.",
+    "Intergenic": "Query isoform lies in an intergenic region.",
+}
+
+
 def get_bokeh_widgets_js():
     """Return the inline Bokeh widgets JavaScript bundle."""
     widgets_index = BOKEH_INLINE.components_for("js").index("bokeh-widgets")
@@ -69,6 +101,11 @@ def _format_ratio_value(value):
     return f"{numeric:.2f}x"
 
 
+def _format_classification_label(name):
+    """Return canonical report label for a summary classification value."""
+    return str(name).strip().replace("-", "_").replace("_", " ").capitalize()
+
+
 def _transcriptome_summary(transcriptome_dir):
     """Return transcriptome metrics and transcript class counts DataFrames."""
     tx_meta = _read_table(Path(transcriptome_dir) / "transcript_metadata.tsv")
@@ -122,13 +159,40 @@ def _sample_summaries(samples_dir):
     return summaries
 
 
-def _sqanti_tables(sqanti_dir):
+def _sqanti_table(sqanti_dir):
-    """Return a dict of SQANTI3 classification summary DataFrames keyed by label."""
+    """Return a SQANTI3 classification summary DataFrame."""
-    tables = {}
+    rows = []
-    for summary in sorted(Path(sqanti_dir).rglob("classification_summary.tsv")):
+    summaries = []
-        label = summary.parent.name
+    for root, _, files in os.walk(sqanti_dir, followlinks=True):
-        tables[label] = _read_table(summary)
+        if "classification_summary.tsv" in files:
-    return tables
+            summaries.append(Path(root) / "classification_summary.tsv")
+    for summary in sorted(summaries):
+        table = _read_table(summary)
+        if table is None or table.empty:
+            continue
+
+        sample = summary.parent.name
+
+        sample_counts = {"Sample": sample}
+        for _, row in table.iterrows():
+            feature = _format_classification_label(row["structural_category"])
+            count = _coerce_float(row["count"])
+            sample_counts[feature] = int(round(count)) if count is not None else 0
+        rows.append(sample_counts)
+
+    if not rows:
+        return None
+
+    sqanti_df = pd.DataFrame(rows).fillna(0)
+    feature_cols = list(classification_categories.keys())
+    for col in feature_cols:
+        if col not in sqanti_df.columns:
+            sqanti_df[col] = 0
+        sqanti_df[col] = sqanti_df[col].astype(int)
+    sqanti_df = sqanti_df[["Sample"] + feature_cols]
+    is_cohort = sqanti_df["Sample"].str.lower().eq("cohort")
+    sqanti_df = sqanti_df.assign(_is_cohort=is_cohort)
+    return sqanti_df.sort_values(["_is_cohort", "Sample"]).drop(columns="_is_cohort")
 
 
 def _contrast_results(de_dir, filename, n=None):
@@ -686,13 +750,15 @@ def main(args):
                 with tabs.add_tab(stats_file.stem.replace(".flagstat", "")):
                     pre(stats_file.read_text())
 
-    sqanti_tables = _sqanti_tables(args.sqanti_dir)
+    sqanti_table = _sqanti_table(args.sqanti_dir)
-    if sqanti_tables:
+    if sqanti_table is not None and not sqanti_table.empty:
         with report.add_section("SQANTI3 classification", "SQANTI3"):
-            tabs = Tabs()
+            p("Summary of structural classification of isoforms using SQANTI3.")
-            for label, table in sqanti_tables.items():
+            DataTable.from_pandas(sqanti_table, use_index=False)
-                with tabs.add_tab(label):
+            with p():
-                    DataTable.from_pandas(table, use_index=False)
+                for category, description in classification_categories.items():
+                    small(strong(f"{category}: "))
+                    small(raw(f"{description}<br>"))
 
     if args.de_dir and Path(args.de_dir).exists():
         # Load DE QC statistics

docs/06_pipeline_overview.md

@@ -65,7 +65,7 @@ Each sample is aligned to the supplied reference genome with
 [`minimap2`](https://github.com/lh3/minimap2), then sorted and indexed with
 [`samtools`](https://www.htslib.org/). The aligned BAMs under
 `samples/<alias>/alignment/` are the main alignment files used for transcriptome
-analysis, optional `SQANTI3` QC, and optional IGV viewing.
+analysis, optional [`SQANTI3`](https://github.com/conesalab/SQANTI3) QC, and optional IGV viewing.
 
 ### 4. Cohort transcriptome construction
 
@@ -90,8 +90,8 @@ for DE/DTU.
 Transcript FASTA files are derived from GTF plus genome using `gffread`.
 When `--skip_sqanti` is not set, `SQANTI3` classifies the cohort and per-sample
 transcriptomes and produces structural QC summaries. The cohort `SQANTI3`
-results live under `cohort/sqanti_cohort/`, while per-sample `SQANTI3`
+results live under `cohort/sqanti/`, while per-sample `SQANTI3`
-directories are published under `samples/<alias>/<alias>_sqanti/`.
+directories are published under `samples/<alias>/sqanti/`.
 
 ### 7. Optional DE and DTU analysis
 

docs/08_outputs.md

@@ -13,13 +13,13 @@ Output files may be aggregated including information for all samples or provided
 | Cohort transcript counts | cohort/transcript_counts.tsv | Transcript-level count matrix produced by bambu. | aggregated |
 | Cohort gene counts | cohort/gene_counts.tsv | Gene-level count matrix derived from bambu output. | aggregated |
 | Cohort transcript metadata | cohort/transcript_metadata.tsv | Transcript annotations and bambu transcript classes for the cohort model. | aggregated |
-| Cohort SQANTI3 summary | cohort/sqanti_cohort/classification_summary.tsv | SQANTI3 classification summary for the cohort transcriptome when SQANTI3 QC is enabled. | aggregated |
+| Cohort SQANTI3 summary | cohort/sqanti/classification_summary.tsv | SQANTI3 classification summary for the cohort transcriptome when SQANTI3 QC is enabled. | aggregated |
 | Per-sample transcriptome GTF | samples/{{ alias }}/transcripts.gtf | Independent bambu transcript model for an individual sample. | per-sample |
 | Per-sample transcriptome FASTA | samples/{{ alias }}/{{ alias }}.transcriptome.fa | Transcript sequences derived from the per-sample GTF. | per-sample |
 | Per-sample transcript counts | samples/{{ alias }}/transcript_counts.tsv | Transcript-level abundance estimates for the per-sample bambu model. | per-sample |
 | Per-sample gene counts | samples/{{ alias }}/gene_counts.tsv | Gene-level abundance estimates for the per-sample bambu model. | per-sample |
 | Per-sample transcript metadata | samples/{{ alias }}/transcript_metadata.tsv | Transcript annotations and bambu transcript classes for the per-sample model. | per-sample |
-| Per-sample SQANTI3 summary | samples/{{ alias }}/{{ alias }}_sqanti/classification_summary.tsv | SQANTI3 classification summary for the per-sample transcriptome when SQANTI3 QC is enabled. | per-sample |
+| Per-sample SQANTI3 summary | samples/{{ alias }}/sqanti/classification_summary.tsv | SQANTI3 classification summary for the per-sample transcriptome when SQANTI3 QC is enabled. | per-sample |
 | Differential gene expression results | de_analysis/{{ contrast }}/results_dge.tsv | DESeq2 gene-level differential expression results for one contrast. | aggregated |
 | Differential gene expression plots | de_analysis/{{ contrast }}/results_dge.pdf | PDF plots generated during DESeq2 analysis for one contrast. | aggregated |
 | Differential transcript usage results | de_analysis/{{ contrast }}/results_dtu_transcript.tsv | Transcript-level DTU results for one contrast. | aggregated |

main.nf

@@ -225,8 +225,6 @@ workflow pipeline {
             .concat(transcriptome.sample_gene_rds.map { meta, rds -> [rds, "samples/${meta.alias}"] })
             .concat(transcriptome.sample_metadata.map { meta, metadata -> [metadata, "samples/${meta.alias}"] })
             .concat(generated_alignment_outputs)
-            .concat(transcriptome.joint_sqanti_dir.map { [it, "cohort"] })
-            .concat(transcriptome.sample_sqanti_dirs.map { meta, sqanti_dir -> [sqanti_dir, "samples/${meta.alias}"] })
 
         if (params.de_analysis) {
             results = results.concat(de_dir.map { [it, null] })

output_definition.json

@@ -89,7 +89,7 @@
       "type": "aggregated"
     },
     "cohort-sqanti-summary": {
-      "filepath": "cohort/sqanti_cohort/classification_summary.tsv",
+      "filepath": "cohort/sqanti/classification_summary.tsv",
       "title": "Cohort SQANTI3 summary",
       "description": "SQANTI3 classification summary for the cohort transcriptome when SQANTI3 QC is enabled.",
       "mime-type": "text/tab-separated-values",
@@ -137,7 +137,7 @@
       "type": "per-sample"
     },
     "sample-sqanti-summary": {
-      "filepath": "samples/{{ alias }}/{{ alias }}_sqanti/classification_summary.tsv",
+      "filepath": "samples/{{ alias }}/sqanti/classification_summary.tsv",
       "title": "Per-sample SQANTI3 summary",
       "description": "SQANTI3 classification summary for the per-sample transcriptome when SQANTI3 QC is enabled.",
       "mime-type": "text/tab-separated-values",

subworkflows/transcriptome.nf

@@ -140,18 +140,21 @@ process runJointSqanti {
     label "wf_transcriptomes_sqanti"
     cpus { params.threads ?: 4 }
     memory "24 GB"
+    publishDir "${params.out_dir}/cohort",
+        mode: "copy",
+        saveAs: { "sqanti" }
     input:
         path gtf
         path annotation, stageAs: "annotation/*"
         tuple path(reference, stageAs: "reference/reference.fa"), path(ref_fai, stageAs: "reference/reference.fai")
     output:
-        path "sqanti_cohort", emit: dir
+        path "cohort", emit: dir
-        path "sqanti_cohort/classification_summary.tsv", emit: summary
+        path "cohort/classification_summary.tsv", emit: summary
     script:
         String extra = params.sqanti_extra_args ?: ""
         String skip_orf = params.sqanti_skip_orf ? "--skipORF" : ""
     """
-    mkdir sqanti_cohort
+    mkdir cohort
     sqanti3_qc.py \
         --isoforms "${gtf}" \
         --refGTF "${annotation}" \
@@ -160,11 +163,11 @@ process runJointSqanti {
         --force_id_ignore \
         --report skip \
         -t ${task.cpus} \
-        -d sqanti_cohort \
+        -d cohort \
         -o cohort \
         ${extra}
-    workflow-glue summarise_sqanti --sqanti_dir sqanti_cohort \
+    workflow-glue summarise_sqanti --sqanti_dir cohort \
-        --output sqanti_cohort/classification_summary.tsv
+        --output cohort/classification_summary.tsv
     """
 }
 
@@ -173,18 +176,21 @@ process runPerSampleSqanti {
     label "wf_transcriptomes_sqanti"
     cpus { params.threads ?: 4 }
     memory "24 GB"
+    publishDir "${params.out_dir}/samples/${meta.alias}",
+        mode: "copy",
+        saveAs: { "sqanti"}
     input:
         tuple val(meta), path(gtf)
         path annotation, stageAs: "annotation/*"
         tuple path(reference, stageAs: "reference/reference.fa"), path(ref_fai, stageAs: "reference/reference.fai")
     output:
-        tuple val(meta), path("${meta.alias}_sqanti"), emit: dir
+        tuple val(meta), path("${meta.alias}"), emit: dir
-        tuple val(meta), path("${meta.alias}_sqanti/classification_summary.tsv"), emit: summary
+        tuple val(meta), path("${meta.alias}/classification_summary.tsv"), emit: summary
     script:
         String extra = params.sqanti_extra_args ?: ""
         String skip_orf = params.sqanti_skip_orf ? "--skipORF" : ""
     """
-    mkdir "${meta.alias}_sqanti"
+    mkdir "${meta.alias}"
     sqanti3_qc.py \
         --isoforms "${gtf}" \
         --refGTF "${annotation}" \
@@ -193,11 +199,11 @@ process runPerSampleSqanti {
         --force_id_ignore \
         --report skip \
         -t ${task.cpus} \
-        -d "${meta.alias}_sqanti" \
+        -d "${meta.alias}" \
         -o "${meta.alias}" \
         ${extra}
-    workflow-glue summarise_sqanti --sqanti_dir "${meta.alias}_sqanti" \
+    workflow-glue summarise_sqanti --sqanti_dir "${meta.alias}" \
-        --output "${meta.alias}_sqanti/classification_summary.tsv"
+        --output "${meta.alias}/classification_summary.tsv"
     """
 }