diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 0a22d18..9801cf5 100644 --- a/.gitlab-ci.yml +++ b/.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" diff --git a/README.md b/README.md index 58dc786..052d99c 100644 --- a/README.md +++ b/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//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` -directories are published under `samples//_sqanti/`. +results live under `cohort/sqanti/`, while per-sample `SQANTI3` +directories are published under `samples//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 | diff --git a/bin/workflow_glue/report.py b/bin/workflow_glue/report.py index 21e647c..0e86613 100644 --- a/bin/workflow_glue/report.py +++ b/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): - """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 - tables[label] = _read_table(summary) - return tables +def _sqanti_table(sqanti_dir): + """Return a SQANTI3 classification summary DataFrame.""" + rows = [] + summaries = [] + for root, _, files in os.walk(sqanti_dir, followlinks=True): + if "classification_summary.tsv" in files: + 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) - if sqanti_tables: + sqanti_table = _sqanti_table(args.sqanti_dir) + if sqanti_table is not None and not sqanti_table.empty: with report.add_section("SQANTI3 classification", "SQANTI3"): - tabs = Tabs() - for label, table in sqanti_tables.items(): - with tabs.add_tab(label): - DataTable.from_pandas(table, use_index=False) + p("Summary of structural classification of isoforms using SQANTI3.") + DataTable.from_pandas(sqanti_table, use_index=False) + with p(): + for category, description in classification_categories.items(): + small(strong(f"{category}: ")) + small(raw(f"{description}
")) if args.de_dir and Path(args.de_dir).exists(): # Load DE QC statistics diff --git a/docs/06_pipeline_overview.md b/docs/06_pipeline_overview.md index a7f763c..0e73039 100644 --- a/docs/06_pipeline_overview.md +++ b/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//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` -directories are published under `samples//_sqanti/`. +results live under `cohort/sqanti/`, while per-sample `SQANTI3` +directories are published under `samples//sqanti/`. ### 7. Optional DE and DTU analysis diff --git a/docs/08_outputs.md b/docs/08_outputs.md index bdfddd2..d7cd4a5 100644 --- a/docs/08_outputs.md +++ b/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 | diff --git a/main.nf b/main.nf index 5682b04..88053f9 100644 --- a/main.nf +++ b/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] }) diff --git a/output_definition.json b/output_definition.json index 13dd96b..171da6c 100644 --- a/output_definition.json +++ b/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", diff --git a/subworkflows/transcriptome.nf b/subworkflows/transcriptome.nf index c41d45b..7500635 100644 --- a/subworkflows/transcriptome.nf +++ b/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 "sqanti_cohort/classification_summary.tsv", emit: summary + path "cohort", emit: dir + 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 \ - --output sqanti_cohort/classification_summary.tsv + workflow-glue summarise_sqanti --sqanti_dir cohort \ + --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}_sqanti/classification_summary.tsv"), emit: summary + tuple val(meta), path("${meta.alias}"), emit: dir + 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" \ - --output "${meta.alias}_sqanti/classification_summary.tsv" + workflow-glue summarise_sqanti --sqanti_dir "${meta.alias}" \ + --output "${meta.alias}/classification_summary.tsv" """ }