From c50f98547769eb87fc4c7d229774c7f33408f679 Mon Sep 17 00:00:00 2001 From: Syed Shujaat Ali Zaidi Date: Mon, 31 Aug 2026 11:29:00 -0500 Subject: [PATCH 1/7] Single-cell: interactive master summary, TCRdist3 + TCRi fixes, GIANA clonotype dedup MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Single-cell modality work. Both routes verified end to end on 8 samples / 2 patients: VDJ-only 123 tasks, full SC 126 tasks, 0 failures, 12/12 modules contributing. ## Fixes **TCRdist3 was silently dropped from every full-SC run.** `enrich_seurat.R` required `rhdf5`, which is absent from the single-cell container (only `hdf5r` is present), so every sample fell through to a `.csv` branch that never exists in sparse mode. Because `tcrdist_export` is `optional: true`, nothing failed — the method just vanished. Now reads via whichever HDF5 package is available: **0 → 22,496 cells annotated (27%)**. **TCRi was dead code.** `modules/scratch/TCRI/` and `subworkflows/scratch/tcri.nf` existed but nothing included them, while `run_tcri = true` in the config implied otherwise. Wired `TCRI_SW` into the SC workflow, added the missing `tables`/`figures` emits, and guarded an unbound `LD_LIBRARY_PATH` that aborted the task under `set -u` before any R ran. **17,843 / 18,610 cells scored (95.9%).** **Master Summary rebuilt.** 21 of the 30 file params the `.qmd` declared were never passed, so 7 sections rendered blank; `stageAs` was also stripping file extensions, so no figures embedded at all. Replaced the ~30 hand-wired params with per-module table staging (`intables//*`), so new upstream tables no longer need process-input changes. Now 17 interactive plotly figures and 12 sortable tables. Figures use native `plot_ly()` rather than `ggplotly()`, which is broken against ggplot2 4.x in this container. **New `CLUSTER_ROLLUP` step.** GIANA/GLIPH2/TCRdist3 produced raw per-patient output but no rollup, so those modules always reported absent. Adds rollups plus method coverage, and radius-based TCRdist clustering for the VDJ-only route. **`VDJ_QC`**: added `pairing_bar_by_sample`, which the Master Summary requested by name but no module produced. ## Please review: one change reaches the bulk engine `modules/local/compare/giana.nf` — collapses to one row per clonotype before GIANA. `PATIENT_CONCATENATE` pools a patient's samples by stacking rows, so GIANA was clustering cross-sample duplicates rather than similar sequences: **76 of 77 clusters on the single-cell test set contained a single distinct CDR3b**, and single-sample patients produced no output at all. This affects bulk too, so it is **gated behind `giana_dedup_clonotypes`, default `false`** — bulk results are unchanged unless explicitly enabled. `params_singlecell.yml` sets it true. Happy to drop the change entirely if you would rather handle it separately. No other bulk-side file is modified, and nothing is deleted. ## Not covered - Bulk mode has not been executed against this branch (`-preview` only). The gate makes it a no-op by default, but it is unverified. - No regression tests added. The `rhdf5` issue in particular is invisible when it breaks. - Tested on one dataset only. --- .gitignore | 11 + bin/cluster_rollup.py | 390 ++++++ bin/enrich_seurat.R | 30 +- modules/bridges/cluster_rollup.nf | 75 + modules/local/compare/giana.nf | 33 +- .../MASTER_SUMMARY/Master_Summary_Report.qmd | 1226 +++++++++++------ modules/scratch/MASTER_SUMMARY/main.nf | 57 +- modules/scratch/TCRI/main.nf | 5 +- modules/scratch/VDJ_QC/VDJ_QC_analysis.qmd | 26 + nextflow.config | 8 + params_singlecell.yml | 4 + subworkflows/scratch/master_summary.nf | 44 +- subworkflows/scratch/tcri.nf | 4 + workflows/tcrtoolkit_sc.nf | 93 +- 14 files changed, 1504 insertions(+), 502 deletions(-) create mode 100755 bin/cluster_rollup.py create mode 100644 modules/bridges/cluster_rollup.nf diff --git a/.gitignore b/.gitignore index e9ff6c1..f1f5157 100644 --- a/.gitignore +++ b/.gitignore @@ -38,3 +38,14 @@ tmp .vscode/* .e2e_test_tmp/ + +## Run artifacts — existing rules only match bare `work` and `results*`, +## so custom -w / --outdir names slip through. +work*/ +test_results*/ +test_v2port*/ +logs/ + +## Let the report notebooks through: the `notebooks/*` rule above would +## otherwise silently skip every template_*.qmd, which the bulk route needs. +!notebooks/template_*.qmd diff --git a/bin/cluster_rollup.py b/bin/cluster_rollup.py new file mode 100755 index 0000000..76ad1d6 --- /dev/null +++ b/bin/cluster_rollup.py @@ -0,0 +1,390 @@ +#!/usr/bin/env python3 +""" +CLUSTER_ROLLUP + +Summarises the clonotype-clustering modules (GIANA / GLIPH2 / TCRdist3) into the +rollup tables the Master Summary expects. Those three modules write raw per-patient +/ per-sample output but no rollup, so MASTER_SUMMARY's giana_summary_file / +gliph2_summary_file / tcrdist3_summary_file had nothing to read. + +Outputs (written into --outdir, matching the .qmd's expected filenames): + + giana_summary_rollup.tsv metric / value + gliph2_summary_rollup.tsv metric / value + tcrdist3_summary_rollup.tsv metric / value + + method_presence_summary.tsv method / frac_cells (fallback: see below) + method_cluster_counts.tsv method / n_clusters (fallback: see below) + + annotation_giana_summary.tsv annot / frac_giana (only when informative) + annotation_gliph2_summary.tsv annot / frac_gliph + annotation_tcrdist3_summary.tsv annot / frac_tcrdist + +The two method_* tables are also produced by CONSENSUS_CLUSTERING, which sees CoNGA +and consensus labels too and is therefore the better source. The workflow prefers the +CONSENSUS copies when consensus ran and falls back to these otherwise (VDJ-only). + +A file is simply not written when its inputs are absent, so the Master Summary's +0-byte/NO_FILE handling treats it as missing rather than empty. +""" + +import argparse +import os +import sys + +import pandas as pd + +# CDR3 amino-acid column, in resolution order, across the various export flavours. +CDR3_COLS = ["cdr3b", "CDR3b", "junction_aa", "cdr3_b_aa", "cdr3"] +ANNOT_COLS = ["annot", "Annotation", "celltype", "CellType", "seurat_clusters"] + + +def usable(path): + """A staged NO_FILE placeholder is 0 bytes; treat it (and absent paths) as missing.""" + return bool(path) and os.path.isfile(path) and os.path.getsize(path) > 0 + + +def tag_from_name(path, suffix): + """Patient/sample identity is carried in the filename prefix, e.g. HRS371754_giana.txt.""" + base = os.path.basename(path) + return base[: -len(suffix)] if base.endswith(suffix) else os.path.splitext(base)[0] + + +def fmt(v): + """Keep the metric/value column as text so integer counts don't render as '263.0'.""" + if isinstance(v, float) and v.is_integer(): + return str(int(v)) + if isinstance(v, float): + return f"{v:g}" + return str(v) + + +def write_rollup(outdir, name, rows): + df = pd.DataFrame([(m, fmt(v)) for m, v in rows], columns=["metric", "value"]) + df.to_csv(os.path.join(outdir, name), sep="\t", index=False) + print(f"[cluster_rollup] wrote {name} ({len(df)} metrics)") + + +def read_giana(paths): + """GIANA writes two '##' comment lines before the real header.""" + frames = [] + for p in paths: + if not usable(p): + continue + try: + df = pd.read_csv(p, sep="\t", comment="#", low_memory=False) + except Exception as exc: # a truncated/empty per-patient file must not sink the run + print(f"[cluster_rollup] WARN: could not read {p}: {exc}", file=sys.stderr) + continue + if df.empty: + continue + df["__group"] = tag_from_name(p, "_giana.txt") + frames.append(df) + return pd.concat(frames, ignore_index=True) if frames else None + + +def read_gliph2(paths): + frames = [] + for p in paths: + if not usable(p): + continue + try: + df = pd.read_csv(p, sep="\t", low_memory=False) + except Exception as exc: + print(f"[cluster_rollup] WARN: could not read {p}: {exc}", file=sys.stderr) + continue + if df.empty: + continue + df["__group"] = tag_from_name(p, "_cluster_member_details.txt") + frames.append(df) + return pd.concat(frames, ignore_index=True) if frames else None + + +def read_tcrdist(paths): + # SAMPLE emits clone_df.csv alongside the .hdf5/.csv distance matrices; only the + # clone tables are parseable here, so ignore everything else rather than warning on it. + frames = [] + for p in paths: + if not usable(p) or not p.endswith("_clone_df.csv"): + continue + try: + df = pd.read_csv(p, low_memory=False) + except Exception as exc: + print(f"[cluster_rollup] WARN: could not read {p}: {exc}", file=sys.stderr) + continue + if df.empty: + continue + df["__group"] = tag_from_name(p, "_clone_df.csv") + frames.append(df) + return pd.concat(frames, ignore_index=True) if frames else None + + +def tcrdist_clusters(matrix_paths, radius): + """ + Threshold each sample's pairwise distance matrix at `radius` and count connected + components of size >= 2 - the same definition enrich_seurat.R uses. + + Only the Full-SC route ran that clustering (it lives in enrich_seurat.R), so on the + VDJ-only route TCRdist3 produced distance matrices that were never thresholded and the + method reported 0 clusters. Doing it here makes TCRdist3 comparable to GIANA/GLIPH2 on + both routes. + + Returns (n_clusters_total, {sample: n_clusters}, n_clustered_clones). + """ + try: + import numpy as np + import scipy.sparse as sp + from scipy.sparse.csgraph import connected_components + except ImportError as exc: + print(f"[cluster_rollup] WARN: scipy/numpy unavailable, skipping tcrdist clustering: {exc}", + file=sys.stderr) + return 0, {}, 0 + + per_sample = {} + total_clustered = 0 + + for p in matrix_paths: + if not usable(p): + continue + name = os.path.basename(p) + sample = name.split("_distance_matrix")[0] + ext = os.path.splitext(p)[1].lower() + try: + if ext == ".csv": + full = np.loadtxt(p, delimiter=",") + adj = sp.csr_matrix((full <= radius) & (full > 0)) + elif ext == ".hdf5": + import h5py + with h5py.File(p, "r") as f: + m = sp.csr_matrix( + (f["data"][:], f["indices"][:], f["indptr"][:]), + shape=tuple(f["shape"][:]), + ) + # Stored sparse: an absent entry means "not within the sparsity cutoff", + # so only explicit non-zero entries can be within radius. + coo = m.tocoo() + keep = (coo.data <= radius) & (coo.data > 0) + adj = sp.coo_matrix( + (np.ones(keep.sum()), (coo.row[keep], coo.col[keep])), shape=m.shape + ).tocsr() + else: + continue + except Exception as exc: + print(f"[cluster_rollup] WARN: could not read {p}: {exc}", file=sys.stderr) + continue + + n_comp, labels = connected_components(adj, directed=False) + sizes = pd.Series(labels).value_counts() + n_real = int((sizes >= 2).sum()) + per_sample[sample] = n_real + total_clustered += int(sizes[sizes >= 2].sum()) + + return sum(per_sample.values()), per_sample, total_clustered + + +def first_col(df, candidates): + for c in candidates: + if c in df.columns: + return c + return None + + +def resolve_cdr3(df): + """Return a Series of CDR3b strings, falling back to the '_' clone_id form.""" + col = first_col(df, CDR3_COLS) + if col is not None: + return df[col].astype("string").str.strip() + for key in ("clone_id", "CTaa", "clonotype"): + if key in df.columns: + return df[key].astype("string").str.split("_", n=1).str[0].str.strip() + return None + + +def cluster_key(df, cluster_col): + """ + GIANA and GLIPH2 both number/label clusters *within* a patient, so the same id + reappears across patients for entirely unrelated clusters (45 of 77 collide in the + 8-sample test set). Counting on the raw id silently merges them - scope it by the + patient the file came from. + """ + if cluster_col not in df.columns: + return None + return df["__group"].astype("string") + "::" + df[cluster_col].astype("string") + + +def cluster_series(df, cdr3, cluster_col): + """Map CDR3 -> set of CDR3s that the method actually assigned to a cluster.""" + if cluster_col not in df.columns or cdr3 is None: + return set() + keep = df[cluster_col].notna() & (df[cluster_col].astype("string").str.strip() != "") + return set(cdr3[keep].dropna().unique()) + + +def annotation_table(export, annot_col, cdr3_export, clustered, frac_name): + """Per-annotation fraction of cells whose clonotype was clustered by this method.""" + tmp = pd.DataFrame({"annot": export[annot_col].astype("string"), "cdr3": cdr3_export}) + tmp = tmp.dropna(subset=["annot"]) + if tmp.empty: + return None + tmp["hit"] = tmp["cdr3"].isin(clustered) + out = ( + tmp.groupby("annot", dropna=True)["hit"] + .mean() + .reset_index() + .rename(columns={"hit": frac_name}) + .sort_values(frac_name, ascending=False) + ) + return out + + +def main(): + ap = argparse.ArgumentParser(description="Roll up GIANA / GLIPH2 / TCRdist3 outputs") + ap.add_argument("--outdir", default=".") + ap.add_argument("--giana", nargs="*", default=[]) + ap.add_argument("--gliph2", nargs="*", default=[]) + ap.add_argument("--tcrdist", nargs="*", default=[]) + ap.add_argument("--tcrdist-matrix", nargs="*", default=[], + help="per-sample distance matrices (.hdf5/.csv) to threshold at --tcrdist-radius") + ap.add_argument("--tcrdist-radius", type=float, default=24.0) + ap.add_argument("--export-cells", default=None) + args = ap.parse_args() + + os.makedirs(args.outdir, exist_ok=True) + + giana = read_giana(args.giana) + gliph = read_gliph2(args.gliph2) + tcrd = read_tcrdist(args.tcrdist) + + export = None + cdr3_export = None + annot_col = None + if usable(args.export_cells): + try: + export = pd.read_csv(args.export_cells, sep="\t", low_memory=False) + cdr3_export = resolve_cdr3(export) + annot_col = first_col(export, ANNOT_COLS) + except Exception as exc: + print(f"[cluster_rollup] WARN: could not read export_cells: {exc}", file=sys.stderr) + export = None + + clustered = {} + counts = {} + + # ── GIANA ──────────────────────────────────────────────────────────────── + if giana is not None: + cdr3 = resolve_cdr3(giana) + gkey = cluster_key(giana, "cluster") + n_clusters = gkey.nunique() if gkey is not None else 0 + sizes = gkey.value_counts() if gkey is not None else pd.Series(dtype=int) + write_rollup(args.outdir, "giana_summary_rollup.tsv", [ + ("Patients with clusters", giana["__group"].nunique()), + ("Clonotypes clustered", int(len(giana))), + ("Unique CDR3b clustered", int(cdr3.nunique()) if cdr3 is not None else 0), + ("Clusters detected", int(n_clusters)), + ("Largest cluster size", int(sizes.max()) if len(sizes) else 0), + ("Median cluster size", float(sizes.median()) if len(sizes) else 0.0), + ("Samples represented", int(giana["sample"].nunique()) if "sample" in giana.columns else 0), + ]) + clustered["GIANA"] = cluster_series(giana, cdr3, "cluster") + counts["GIANA"] = int(n_clusters) + + # ── GLIPH2 ─────────────────────────────────────────────────────────────── + if gliph is not None: + cdr3 = resolve_cdr3(gliph) + pkey = cluster_key(gliph, "tag") + n_clusters = pkey.nunique() if pkey is not None else 0 + sizes = pkey.value_counts() if pkey is not None else pd.Series(dtype=int) + write_rollup(args.outdir, "gliph2_summary_rollup.tsv", [ + ("Patients with clusters", gliph["__group"].nunique()), + ("Clonotype-cluster memberships", int(len(gliph))), + ("Unique CDR3b clustered", int(cdr3.nunique()) if cdr3 is not None else 0), + ("Motif clusters detected", int(n_clusters)), + ("Largest cluster size", int(sizes.max()) if len(sizes) else 0), + ("Median cluster size", float(sizes.median()) if len(sizes) else 0.0), + ("Samples represented", int(gliph["sample"].nunique()) if "sample" in gliph.columns else 0), + ]) + clustered["GLIPH2"] = cluster_series(gliph, cdr3, "tag") + counts["GLIPH2"] = int(n_clusters) + + # ── TCRdist3 ───────────────────────────────────────────────────────────── + # clone_df carries no cluster labels (those come from thresholding the distance + # matrix in enrich_seurat.R), so cluster counts are taken from the per-cell export + # when CLUSTER_TO_SC has already written a tcrdist_cluster column. + if tcrd is not None: + cdr3 = resolve_cdr3(tcrd) + + # Prefer the per-cell tcrdist_cluster column when CLUSTER_TO_SC already wrote it + # (Full SC); otherwise threshold the distance matrices here so the VDJ-only route + # reports real cluster counts instead of 0. + n_tcrdist_clusters = 0 + per_sample = {} + n_clustered_clones = 0 + if export is not None and "tcrdist_cluster" in export.columns: + n_tcrdist_clusters = int(export["tcrdist_cluster"].nunique(dropna=True)) + elif args.tcrdist_matrix: + n_tcrdist_clusters, per_sample, n_clustered_clones = tcrdist_clusters( + args.tcrdist_matrix, args.tcrdist_radius + ) + + rows = [ + ("Samples analysed", tcrd["__group"].nunique()), + ("Clones in distance matrices", int(len(tcrd))), + ("Unique CDR3b", int(cdr3.nunique()) if cdr3 is not None else 0), + ("Median clones per sample", float(tcrd.groupby("__group").size().median())), + ("Max clones per sample", int(tcrd.groupby("__group").size().max())), + (f"Clusters detected (radius {args.tcrdist_radius:g})", n_tcrdist_clusters), + ] + if n_clustered_clones: + rows.append(("Clones in a cluster", n_clustered_clones)) + if per_sample: + rows.append(("Median clusters per sample", + float(pd.Series(list(per_sample.values())).median()))) + write_rollup(args.outdir, "tcrdist3_summary_rollup.tsv", rows) + + if cdr3 is not None: + clustered["TCRdist3"] = set(cdr3.dropna().unique()) + if n_tcrdist_clusters: + counts["TCRdist3"] = n_tcrdist_clusters + + # ── method_* fallbacks (CONSENSUS supersedes these when it runs) ────────── + if counts: + pd.DataFrame( + sorted(counts.items()), columns=["method", "n_clusters"] + ).to_csv(os.path.join(args.outdir, "method_cluster_counts.tsv"), sep="\t", index=False) + print(f"[cluster_rollup] wrote method_cluster_counts.tsv ({len(counts)} methods)") + + if clustered and export is not None and cdr3_export is not None: + n_cells = len(export) + rows = [ + (m, float(cdr3_export.isin(s).sum()) / n_cells if n_cells else 0.0) + for m, s in sorted(clustered.items()) + ] + pd.DataFrame(rows, columns=["method", "frac_cells"]).to_csv( + os.path.join(args.outdir, "method_presence_summary.tsv"), sep="\t", index=False + ) + print(f"[cluster_rollup] wrote method_presence_summary.tsv ({len(rows)} methods)") + + # ── per-annotation coverage ────────────────────────────────────────────── + # Only meaningful when the export actually carries >1 distinct annotation; + # in VDJ-only mode every cell is "Unannotated" and the panel would be a single bar. + if export is not None and annot_col is not None and cdr3_export is not None: + if export[annot_col].nunique(dropna=True) > 1: + for method, frac_name, fname in ( + ("GIANA", "frac_giana", "annotation_giana_summary.tsv"), + ("GLIPH2", "frac_gliph", "annotation_gliph2_summary.tsv"), + ("TCRdist3", "frac_tcrdist", "annotation_tcrdist3_summary.tsv"), + ): + if method not in clustered: + continue + tbl = annotation_table(export, annot_col, cdr3_export, clustered[method], frac_name) + if tbl is not None and not tbl.empty: + tbl.to_csv(os.path.join(args.outdir, fname), sep="\t", index=False) + print(f"[cluster_rollup] wrote {fname} ({len(tbl)} annotations)") + else: + print("[cluster_rollup] single annotation value in export; skipping annotation panels") + + print("[cluster_rollup] done") + + +if __name__ == "__main__": + main() diff --git a/bin/enrich_seurat.R b/bin/enrich_seurat.R index faeb9be..77196cb 100755 --- a/bin/enrich_seurat.R +++ b/bin/enrich_seurat.R @@ -164,11 +164,31 @@ if (length(clone_df_files) > 0) { mat_hdf5 <- paste0(sample_name, "_distance_matrix.hdf5") mat_csv <- paste0(sample_name, "_distance_matrix.csv") - if (file.exists(mat_hdf5) && requireNamespace("rhdf5", quietly = TRUE)) { - data <- rhdf5::h5read(mat_hdf5, "data") - indices <- rhdf5::h5read(mat_hdf5, "indices") - indptr <- rhdf5::h5read(mat_hdf5, "indptr") - shape <- as.integer(rhdf5::h5read(mat_hdf5, "shape")) + # The SC container ships hdf5r, not rhdf5, and TCRDIST3_MATRIX only writes a + # .csv matrix in dense mode - so requiring rhdf5 meant every sample fell through + # to the missing-CSV branch and TCRdist3 was silently dropped from the whole run + # ("all samples failed or had no distance matrix"). Read with whichever HDF5 + # package is available. + h5_reader <- if (requireNamespace("rhdf5", quietly = TRUE)) "rhdf5" + else if (requireNamespace("hdf5r", quietly = TRUE)) "hdf5r" + else NA_character_ + + read_h5 <- function(path, keys) { + if (identical(h5_reader, "rhdf5")) { + stats::setNames(lapply(keys, function(k) rhdf5::h5read(path, k)), keys) + } else { + fh <- hdf5r::H5File$new(path, mode = "r") + on.exit(fh$close_all(), add = TRUE) + stats::setNames(lapply(keys, function(k) fh[[k]]$read()), keys) + } + } + + if (file.exists(mat_hdf5) && !is.na(h5_reader)) { + h5 <- read_h5(mat_hdf5, c("data", "indices", "indptr", "shape")) + data <- h5$data + indices <- h5$indices + indptr <- h5$indptr + shape <- as.integer(h5$shape) mat <- Matrix::sparseMatrix( i = as.integer(indices) + 1L, p = as.integer(indptr), diff --git a/modules/bridges/cluster_rollup.nf b/modules/bridges/cluster_rollup.nf new file mode 100644 index 0000000..4edd0be --- /dev/null +++ b/modules/bridges/cluster_rollup.nf @@ -0,0 +1,75 @@ +/* + * CLUSTER_ROLLUP + * + * GIANA, GLIPH2 and TCRdist3 write raw per-patient / per-sample output but no rollup + * table, so MASTER_SUMMARY's giana_summary_file / gliph2_summary_file / + * tcrdist3_summary_file had nothing to read and those modules always reported as + * absent. This process summarises them into the filenames the .qmd expects. + * + * Every output is `optional: true`: a route where a method did not run simply yields + * no file, and the workflow substitutes the NO_FILE placeholder that the .qmd already + * treats as missing. + */ +process CLUSTER_ROLLUP { + tag "${project_name}" + label 'process_single' + container "${params.container}" + + publishDir "${params.outdir}/bridge/cluster_rollup", mode: 'copy', overwrite: true + + input: + path giana_files, stageAs: 'giana/*' + path gliph2_files, stageAs: 'gliph2/*' + path tcrdist_files, stageAs: 'tcrdist/*' + path tcrdist_mats, stageAs: 'tcrdistmat/*' + path export_cells, stageAs: 'in_export_cells.tsv' + // Staged rather than invoked as ${projectDir}/bin/... so its contents take part in + // the task hash: otherwise editing the script leaves -resume silently reusing the + // previous output. + path rollup_script, stageAs: 'cluster_rollup.py' + val project_name + + output: + path "giana_summary_rollup.tsv", emit: giana_summary, optional: true + path "gliph2_summary_rollup.tsv", emit: gliph2_summary, optional: true + path "tcrdist3_summary_rollup.tsv", emit: tcrdist3_summary, optional: true + path "method_presence_summary.tsv", emit: method_presence, optional: true + path "method_cluster_counts.tsv", emit: method_cluster_counts, optional: true + path "annotation_giana_summary.tsv", emit: annotation_giana, optional: true + path "annotation_gliph2_summary.tsv", emit: annotation_gliph2, optional: true + path "annotation_tcrdist3_summary.tsv", emit: annotation_tcrdist3, optional: true + + script: + // nullglob so an empty staging directory expands to nothing rather than a literal + // "giana/*". Guards are written as if-blocks, not `test && assign`: Nextflow runs the + // script under `set -e`, where a failing `&&` chain would abort the whole task. + """ + shopt -s nullglob + + giana_args=() + gliph_args=() + tcrdist_args=() + matrix_args=() + + g=( giana/* ) + if [ \${#g[@]} -gt 0 ]; then giana_args=( --giana "\${g[@]}" ); fi + + p=( gliph2/* ) + if [ \${#p[@]} -gt 0 ]; then gliph_args=( --gliph2 "\${p[@]}" ); fi + + t=( tcrdist/* ) + if [ \${#t[@]} -gt 0 ]; then tcrdist_args=( --tcrdist "\${t[@]}" ); fi + + m=( tcrdistmat/* ) + if [ \${#m[@]} -gt 0 ]; then matrix_args=( --tcrdist-matrix "\${m[@]}" ); fi + + python3 cluster_rollup.py \\ + --outdir . \\ + --export-cells in_export_cells.tsv \\ + --tcrdist-radius ${params.tcrdist_radius} \\ + \${matrix_args[@]+"\${matrix_args[@]}"} \\ + \${giana_args[@]+"\${giana_args[@]}"} \\ + \${gliph_args[@]+"\${gliph_args[@]}"} \\ + \${tcrdist_args[@]+"\${tcrdist_args[@]}"} + """ +} diff --git a/modules/local/compare/giana.nf b/modules/local/compare/giana.nf index fdaea2c..ec9f52d 100644 --- a/modules/local/compare/giana.nf +++ b/modules/local/compare/giana.nf @@ -13,12 +13,43 @@ process GIANA_CALC { path "${patient}_giana.txt", emit: 'giana_output' // path "giana_EncodingMatrix.txt" - script: + script: + def dedup_clonotypes = (params.giana_dedup_clonotypes ?: false) ? 'True' : 'False' + """ python3 - < {len(df)} unique clonotypes", flush=True) + df.to_csv("concat_cdr3_renamed.tsv", sep="\t", index=False) EOF diff --git a/modules/scratch/MASTER_SUMMARY/Master_Summary_Report.qmd b/modules/scratch/MASTER_SUMMARY/Master_Summary_Report.qmd index ef08f82..c7ecbc5 100644 --- a/modules/scratch/MASTER_SUMMARY/Master_Summary_Report.qmd +++ b/modules/scratch/MASTER_SUMMARY/Master_Summary_Report.qmd @@ -1,64 +1,35 @@ --- -title: "SCRATCH-TCR: Master Summary Report" -author: "Syed Shujaat Ali Zaidi" +title: "TCR-Toolkit-SCRATCH: Master Summary" format: html: toc: true - toc-depth: 3 + toc-depth: 2 + toc-location: left number-sections: true code-fold: true code-summary: "Show code" embed-resources: true - theme: defualt - df-print: paged + theme: cosmo + fig-width: 10 + fig-height: 6 execute: echo: false warning: false message: false params: - project_name: "SCRATCH-TCR Project" + project_name: "TCR-Toolkit-SCRATCH" outdir: "Master_Summary_Report" - vdj_qc_per_sample_compact_file: "" - vdj_qc_before_after_summary_file: "" - vdj_qc_sample_sheet_resolved_file: "" - vdj_qc_clone_rank_abundance_file: "" - - vdj_qc_before_after_retention_fig: "" - vdj_qc_pairing_bar_fig: "" - vdj_qc_clone_rank_abundance_fig: "" - vdj_qc_multiple_chains_fig: "" - - # core inputs - seurat_rds: "data/seurat_with_consensus_clonotype_clusters.rds" - export_cells_file: "data/consensus_export_cells.tsv" - - # optional module rollups - vdj_qc_summary_file: "data/vdj_qc_per_sample_compact.tsv" - tcell_integration_summary_file: "data/tcell_summary_rollup.tsv" - tcri_summary_file: "data/tcri_summary_rollup.tsv" - conga_summary_file: "data/conga_summary_rollup.tsv" - gliph2_summary_file: "data/gliph2_summary_rollup.tsv" - tcrdist3_summary_file: "data/tcrdist3_summary_rollup.tsv" - giana_summary_file: "data/giana_summary_rollup.tsv" - consensus_summary_file: "data/consensus_summary_rollup.tsv" - repertoire_summary_file: "data/repertoire_summary_rollup.tsv" - - # optional richer inputs - diversity_by_sample_file: "data/diversity_by_sample.tsv" - clone_burden_file: "data/clone_burden_cells_by_sample.tsv" - sample_overlap_matrix_file: "data/sample_overlap_matrix.tsv" - method_presence_file: "data/method_presence_summary.tsv" - method_cluster_counts_file: "data/method_cluster_counts.tsv" - consensus_cluster_summary_file: "data/consensus_cluster_summary.tsv" - annotation_tcri_summary_file: "data/annotation_tcri_summary.tsv" - annotation_conga_summary_file: "data/annotation_conga_summary.tsv" - annotation_gliph2_summary_file: "data/annotation_gliph2_summary.tsv" - annotation_tcrdist3_summary_file: "data/annotation_tcrdist3_summary.tsv" - annotation_giana_summary_file: "data/annotation_giana_summary.tsv" - annotation_consensus_summary_file: "data/annotation_consensus_summary.tsv" - - # mapping + # Root of the staged per-module table directories (tables_dir//). + # Replaces the previous ~30 individually-wired file params: modules publish whole + # tables/ directories, so new upstream tables become available here without another + # round of process-input plumbing. + tables_dir: "intables" + + seurat_rds: "" + export_cells_file: "" + + # Column mapping / resolution label_col: "" sample_col: "sample" patient_col: "patient" @@ -76,538 +47,913 @@ params: clone_id_candidates: "clone_id;CTaa;clonotype;clone" clone_size_candidates: "clone_size;CloneSize;clone_n" consensus_cluster_candidates: "consensus_cluster;consensus;consensus_id" - - - figure_format: "png" - figure_width: 10 - figure_height: 7 - figure_dpi: 300 - base_size: 12 + tcrdist_radius: 24 top_n_annotations: 20 - top_n_consensus_clusters: 20 + top_n_clusters: 20 + top_n_genes: 25 + base_size: 12 --- -# setup ```{r} -#| label: setup-chunk +#| label: setup suppressPackageStartupMessages({ - library(Seurat) - library(SeuratObject) library(data.table) library(dplyr) library(tidyr) library(ggplot2) library(forcats) library(scales) - library(knitr) - library(kableExtra) - library(ComplexHeatmap) - library(circlize) - library(patchwork) + library(plotly) + library(reactable) + library(htmltools) }) -dir.create(params$outdir, recursive = TRUE, showWarnings = FALSE) -dir.create(file.path(params$outdir, "tables"), recursive = TRUE, showWarnings = FALSE) +dir.create(file.path(params$outdir, "tables"), recursive = TRUE, showWarnings = FALSE) dir.create(file.path(params$outdir, "figures"), recursive = TRUE, showWarnings = FALSE) `%||%` <- function(a, b) if (!is.null(a) && length(a) > 0 && !all(is.na(a))) a else b -sanitize_param_string <- function(x) { - if (is.null(x) || length(x) == 0) return(NULL) - trimws(gsub("\\u00A0", " ", as.character(x))) +# ── Input access ──────────────────────────────────────────────────────────── +# Every upstream module stages its tables under //. A module that +# did not run on this route simply has no directory, so T() returns NULL and the +# corresponding section renders an explanatory note rather than silently nothing. +T <- function(module, name) { + p <- file.path(params$tables_dir, module, name) + if (!file.exists(p) || file.size(p) == 0) return(NULL) + out <- tryCatch({ + if (tolower(tools::file_ext(name)) == "csv") fread(p) else fread(p, sep = "\t") + }, error = function(e) NULL) + if (is.null(out) || nrow(out) == 0) return(NULL) + as.data.frame(out) } -normalize_colnames <- function(df) { - colnames(df) <- trimws(gsub("\\u00A0", " ", colnames(df))) - df +# `-P label_col=` (an empty params.label_col) reaches R as NULL, and nzchar(NULL) is +# logical(0), which makes if() error rather than return FALSE. Normalise every incoming +# string param to a length-1 character before testing it. +as_str1 <- function(x) { + if (is.null(x) || length(x) == 0) return("") + x <- as.character(x)[1] + if (is.na(x)) "" else trimws(x) } -split_candidate_string <- function(x) { - if (is.null(x) || length(x) == 0 || is.na(x) || x == "") return(character()) - x <- gsub('"', "", as.character(x)) - x <- trimws(unlist(strsplit(x, ";", fixed = TRUE))) - x[nzchar(x)] +# First match wins across several modules. CONSENSUS computes method_presence / +# method_cluster_counts over CoNGA and consensus labels as well, so its copies are +# strictly richer than CLUSTER_ROLLUP's GIANA/GLIPH2/TCRdist-only fallback. +# Count distinct values, ignoring NA and empty strings. VDJ_QC leaves patient_id blank +# when no GEX object is supplied, so a naive n_distinct() on that column returns 1 for +# "" rather than falling through to a source that actually carries the value. +n_distinct_real <- function(x) { + if (is.null(x)) return(NA_integer_) + x <- trimws(as.character(x)) + x <- x[!is.na(x) & nzchar(x)] + if (!length(x)) return(NA_integer_) + length(unique(x)) } -safe_read_table <- function(path) { - # Treat empty placeholder files (staged NO_FILE, 0 bytes) as absent. - if (is.null(path) || is.na(path) || path == "" || !file.exists(path)) return(NULL) - if (isTRUE(file.size(path) == 0)) return(NULL) - ext <- tolower(tools::file_ext(path)) - if (ext %in% c("tsv", "tab", "txt")) fread(path, sep = "\t") else fread(path) +# First source with a real value wins. +coalesce_count <- function(...) { + for (v in list(...)) if (!is.na(v) && v > 0) return(as.character(v)) + NA_character_ } -first_existing_col <- function(df, preferred, candidates = character()) { - df <- normalize_colnames(df) - preferred <- sanitize_param_string(preferred) - - if (length(candidates) == 1 && is.character(candidates)) { - candidates <- split_candidate_string(candidates) +Tfirst <- function(modules, name) { + for (m in modules) { + d <- T(m, name) + if (!is.null(d)) return(d) } - candidates <- trimws(gsub("\\u00A0", " ", as.character(candidates))) - - if (!is.null(preferred) && preferred != "" && preferred %in% colnames(df)) return(preferred) - - hits <- intersect(candidates, colnames(df)) - if (length(hits) > 0) hits[[1]] else NULL -} - -save_plot_safe <- function(plot_obj, filename, - width = params$figure_width, - height = params$figure_height, - dpi = params$figure_dpi) { - ggsave( - filename = file.path(params$outdir, "figures", filename), - plot = plot_obj, - width = width, - height = height, - dpi = dpi, - bg = "white", - limitsize = FALSE - ) + NULL } -save_table_safe <- function(df, filename) { - fwrite(df, file.path(params$outdir, "tables", filename), sep = "\t") +have <- function(x) !is.null(x) && nrow(x) > 0 + +note <- function(txt) knitr::asis_output(paste0("\n::: {.callout-note appearance=\"simple\"}\n", + txt, "\n:::\n\n")) + +save_table <- function(df, filename) { + if (have(df)) fwrite(df, file.path(params$outdir, "tables", filename), sep = "\t") + invisible(df) } -theme_scratch_pub <- function(base_size = 12) { +theme_ms <- function(base_size = params$base_size) { theme_bw(base_size = base_size) + theme( - plot.title = element_text(face = "bold", size = base_size + 2, hjust = 0), - plot.subtitle = element_text(size = base_size, hjust = 0), - axis.title = element_text(face = "bold"), - axis.text = element_text(color = "black"), + plot.title = element_text(face = "bold", size = base_size + 1, hjust = 0), + axis.title = element_text(face = "bold"), panel.grid.minor = element_blank(), panel.grid.major = element_line(linewidth = 0.2, color = "grey90"), strip.background = element_rect(fill = "grey95", color = "grey80"), - strip.text = element_text(face = "bold"), - legend.title = element_text(face = "bold") + strip.text = element_text(face = "bold"), + legend.title = element_text(face = "bold") ) } -``` +# ── Interactive output helpers ────────────────────────────────────────────── +# Figures are built with plot_ly() directly rather than ggplotly(). The container pairs +# ggplot2 4.0.2 with plotly 4.10.4, and that plotly predates ggplot2's 4.x internals, so +# ggplotly() dies with "subscript out of bounds" on most plots. Native plot_ly() calls do +# not touch ggplot internals and are version-robust. +# +# Every figure supports: hover for exact values, drag to zoom, double-click to reset, +# click legend entries to toggle series. +pal <- c("#2c7fb8", "#d95f02", "#7570b3", "#1b9e77", + "#e7298a", "#66a61e", "#e6ab02", "#a6761d") + +pconf <- function(p) { + plotly::config(p, displaylogo = FALSE, + modeBarButtonsToRemove = c("select2d", "lasso2d"), + toImageButtonOptions = list(format = "png", scale = 2)) +} -# input loading -```{r} -seu <- if (file.exists(params$seurat_rds) && isTRUE(file.size(params$seurat_rds) > 0)) readRDS(params$seurat_rds) else NULL -export_cells <- safe_read_table(params$export_cells_file) +ttl <- function(t) list(text = t, x = 0, xanchor = "left", + font = list(size = params$base_size + 4)) + +# Horizontal bar. df: cat (chr), val (num), hover (chr) +hbar <- function(df, title, xlab, colour = "#2c7fb8", height = 440, tickformat = NULL) { + df <- df[order(df$val), , drop = FALSE] + df$cat <- factor(df$cat, levels = unique(df$cat)) + pconf(plotly::plot_ly(df, x = ~val, y = ~cat, type = "bar", orientation = "h", + marker = list(color = colour), + text = ~hover, hovertemplate = "%{text}") %>% + plotly::layout(title = ttl(title), height = height, + margin = list(l = 140, t = 60), + xaxis = list(title = xlab, tickformat = tickformat), + yaxis = list(title = ""))) +} -vdj_qc_per_sample_compact <- safe_read_table(params$vdj_qc_per_sample_compact_file) -vdj_qc_before_after_summary <- safe_read_table(params$vdj_qc_before_after_summary_file) -vdj_qc_sample_sheet_resolved <- safe_read_table(params$vdj_qc_sample_sheet_resolved_file) -vdj_qc_clone_rank_abundance <- safe_read_table(params$vdj_qc_clone_rank_abundance_file) -vdj_qc_summary <- safe_read_table(params$vdj_qc_summary_file) -tcell_integration_summary <- safe_read_table(params$tcell_integration_summary_file) -tcri_summary <- safe_read_table(params$tcri_summary_file) -conga_summary <- safe_read_table(params$conga_summary_file) -gliph2_summary <- safe_read_table(params$gliph2_summary_file) -tcrdist3_summary <- safe_read_table(params$tcrdist3_summary_file) -giana_summary <- safe_read_table(params$giana_summary_file) -consensus_summary <- safe_read_table(params$consensus_summary_file) -repertoire_summary <- safe_read_table(params$repertoire_summary_file) +# Stacked horizontal bar. df: cat, val, grp, hover +sbar <- function(df, title, xlab, height = 470, tickformat = ".0%") { + pconf(plotly::plot_ly(df, x = ~val, y = ~cat, color = ~grp, colors = pal, + type = "bar", orientation = "h", + text = ~hover, hovertemplate = "%{text}") %>% + plotly::layout(barmode = "stack", title = ttl(title), height = height, + margin = list(l = 140, t = 60), + xaxis = list(title = xlab, tickformat = tickformat), + yaxis = list(title = ""), + legend = list(orientation = "h", y = -0.15))) +} -diversity_by_sample <- safe_read_table(params$diversity_by_sample_file) -clone_burden <- safe_read_table(params$clone_burden_file) -sample_overlap_matrix <- safe_read_table(params$sample_overlap_matrix_file) -method_presence <- safe_read_table(params$method_presence_file) -method_cluster_counts <- safe_read_table(params$method_cluster_counts_file) -consensus_cluster_summary <- safe_read_table(params$consensus_cluster_summary_file) +# Grouped vertical bar. df: cat, val, grp, hover +gbar <- function(df, title, xlab, ylab, height = 470) { + pconf(plotly::plot_ly(df, x = ~cat, y = ~val, color = ~grp, colors = pal, + type = "bar", + text = ~hover, hovertemplate = "%{text}") %>% + plotly::layout(barmode = "group", title = ttl(title), height = height, + margin = list(t = 60), + xaxis = list(title = xlab), yaxis = list(title = ylab), + legend = list(orientation = "h", y = -0.2))) +} -annotation_tcri_summary <- safe_read_table(params$annotation_tcri_summary_file) -annotation_conga_summary <- safe_read_table(params$annotation_conga_summary_file) -annotation_gliph2_summary <- safe_read_table(params$annotation_gliph2_summary_file) -annotation_tcrdist3_summary <- safe_read_table(params$annotation_tcrdist3_summary_file) -annotation_giana_summary <- safe_read_table(params$annotation_giana_summary_file) -annotation_consensus_summary <- safe_read_table(params$annotation_consensus_summary_file) +# Box plot with points. df: grp, val, hover +boxby <- function(df, title, ylab, height = 470, logy = FALSE, angle = 0) { + pconf(plotly::plot_ly(df, x = ~grp, y = ~val, color = ~grp, colors = pal, + type = "box", boxpoints = "all", jitter = 0.4, pointpos = 0, + marker = list(size = 4, opacity = 0.55), + text = ~hover, hoverinfo = "text") %>% + plotly::layout(title = ttl(title), height = height, showlegend = FALSE, + margin = list(t = 60, b = 110), + xaxis = list(title = "", tickangle = angle), + yaxis = list(title = ylab, + type = if (logy) "log" else "linear"))) +} + +# Scatter. df: x, y, grp, hover (grp may be a constant) +scat <- function(df, title, xlab, ylab, height = 470, size = 11) { + pconf(plotly::plot_ly(df, x = ~x, y = ~y, color = ~grp, colors = pal, + type = "scatter", mode = "markers", + marker = list(size = size, opacity = 0.85), + text = ~hover, hovertemplate = "%{text}") %>% + plotly::layout(title = ttl(title), height = height, margin = list(t = 60), + xaxis = list(title = xlab), yaxis = list(title = ylab))) +} + +# Line/step. df: x, y, grp, hover +lplot <- function(df, title, xlab, ylab, logx = TRUE, logy = TRUE, height = 470) { + pconf(plotly::plot_ly(df, x = ~x, y = ~y, color = ~grp, colors = pal, + type = "scatter", mode = "lines", + text = ~hover, hovertemplate = "%{text}") %>% + plotly::layout(title = ttl(title), height = height, margin = list(t = 60), + xaxis = list(title = xlab, type = if (logx) "log" else "linear"), + yaxis = list(title = ylab, type = if (logy) "log" else "linear"))) +} + +itable <- function(df, caption = NULL, pageSize = 10, ...) { + if (!have(df)) return(note("No data available.")) + reactable( + df, + searchable = TRUE, sortable = TRUE, filterable = FALSE, + highlight = TRUE, compact = TRUE, striped = TRUE, + defaultPageSize = pageSize, showPageSizeOptions = TRUE, + pageSizeOptions = c(5, 10, 25, 50, 100), + defaultColDef = colDef( + format = colFormat(digits = 3), + minWidth = 90, + headerStyle = list(background = "#f7f7f8", fontWeight = 600) + ), + ... + ) +} ``` -# Resolving Columns ```{r} -if (is.null(export_cells) && !is.null(seu)) { - export_cells <- seu@meta.data %>% tibble::rownames_to_column("cell_id") +#| label: load-all +# ── VDJ QC ── +vdj_compact <- T("vdj_qc", "vdj_qc_per_sample_compact.tsv") +vdj_retention <- T("vdj_qc", "qc_contigs_before_after_summary.tsv") +vdj_rank <- T("vdj_qc", "clone_rank_abundance.tsv") +vdj_pairing <- T("vdj_qc", "pairing_status_all.tsv") +vdj_cdr3len <- T("vdj_qc", "cdr3_lengths_all.tsv") +vdj_vusage <- T("vdj_qc", "v_gene_usage.tsv") +vdj_jusage <- T("vdj_qc", "j_gene_usage.tsv") +vdj_samplesheet<- T("vdj_qc", "sample_sheet_resolved.tsv") + +# ── Pseudobulk QC gate ── +pb_qc <- T("pseudobulk", "pseudobulk_qc_summary.csv") +pb_vfam <- T("pseudobulk", "pseudobulk_v_family.csv") + +# ── Shared bulk engine ── +sample_stats <- T("sample", "sample_stats.csv") +sample_vfam <- T("sample", "v_family.csv") +sample_jfam <- T("sample", "j_family.csv") +sharing <- T("compare", "cdr3_sharing.tsv") + +# ── Clustering rollups ── +giana_roll <- T("rollup", "giana_summary_rollup.tsv") +gliph2_roll <- T("rollup", "gliph2_summary_rollup.tsv") +tcrdist_roll <- T("rollup", "tcrdist3_summary_rollup.tsv") +method_presence<- Tfirst(c("consensus", "rollup"), "method_presence_summary.tsv") +method_counts <- Tfirst(c("consensus", "rollup"), "method_cluster_counts.tsv") +anno_giana <- T("rollup", "annotation_giana_summary.tsv") +anno_gliph2 <- T("rollup", "annotation_gliph2_summary.tsv") +anno_tcrdist <- T("rollup", "annotation_tcrdist3_summary.tsv") + +# ── Repertoire ── +rep_rollup <- T("repertoire", "repertoire_summary_rollup.tsv") +diversity <- T("repertoire", "diversity_by_sample.tsv") +div_patient <- T("repertoire", "diversity_by_patient.tsv") +clone_burden <- T("repertoire", "clone_burden_cells_by_sample.tsv") +overlap_mat <- T("repertoire", "sample_overlap_matrix.tsv") +top_shared <- T("repertoire", "top_shared_clones.tsv") +shareability <- T("repertoire", "clone_shareability.tsv") +pheno_flux <- T("repertoire", "phenotypic_flux.tsv") + +# ── GEX-dependent ── +tcell_roll <- T("tcell", "tcell_summary_rollup.tsv") +tcell_anno <- T("tcell", "tcell_per_annotation_summary.tsv") +tcell_topclone <- T("tcell", "tcell_top_clones.tsv") +conga_roll <- T("conga", "conga_summary_rollup.tsv") +conga_clusters <- T("conga", "conga_cluster_summary.tsv") +anno_conga <- T("conga", "annotation_conga_summary.tsv") +cons_roll <- T("consensus", "consensus_summary_rollup.tsv") +tcri_roll <- T("tcri", "tcri_summary_rollup.tsv") +tcri_scores <- T("tcri", "high_tcri_clone_summary.tsv") +anno_tcri <- T("tcri", "annotation_tcri_summary.tsv") +cons_clusters <- T("consensus", "consensus_cluster_summary.tsv") +anno_cons <- T("consensus", "annotation_consensus_summary.tsv") + +# ── Per-cell export ── +export_cells <- NULL +.ec <- as_str1(params$export_cells_file) +if (nzchar(.ec) && file.exists(.ec) && file.size(.ec) > 0) { + export_cells <- tryCatch(as.data.frame(fread(.ec, sep = "\t")), + error = function(e) NULL) } -label_col <- if (!is.null(export_cells)) first_existing_col(export_cells, params$label_col, params$label_candidates) else NULL -sample_col <- if (!is.null(export_cells)) first_existing_col(export_cells, params$sample_col, params$sample_candidates) else NULL -patient_col <- if (!is.null(export_cells)) first_existing_col(export_cells, params$patient_col, params$patient_candidates) else NULL -condition_col <- if (!is.null(export_cells)) first_existing_col(export_cells, params$condition_col, params$condition_candidates) else NULL -timepoint_col <- if (!is.null(export_cells)) first_existing_col(export_cells, params$timepoint_col, params$timepoint_candidates) else NULL -clone_id_col <- if (!is.null(export_cells)) first_existing_col(export_cells, params$clone_id_col, params$clone_id_candidates) else NULL -clone_size_col <- if (!is.null(export_cells)) first_existing_col(export_cells, params$clone_size_col, params$clone_size_candidates) else NULL -consensus_cluster_col <- if (!is.null(export_cells)) first_existing_col(export_cells, params$consensus_cluster_col, params$consensus_cluster_candidates) else NULL +pick_col <- function(df, preferred, candidates) { + if (is.null(df) || !ncol(df)) return(NULL) + preferred <- as_str1(preferred) + cand <- trimws(unlist(strsplit(gsub('"', "", as_str1(candidates)), ";", fixed = TRUE))) + cand <- cand[nzchar(cand)] + if (nzchar(preferred) && preferred %in% names(df)) return(preferred) + hit <- intersect(cand, names(df)) + if (length(hit)) hit[[1]] else NULL +} +sample_col <- pick_col(export_cells, params$sample_col, params$sample_candidates) +patient_col <- pick_col(export_cells, params$patient_col, params$patient_candidates) +label_col <- pick_col(export_cells, params$label_col, params$label_candidates) +clone_col <- pick_col(export_cells, params$clone_id_col, params$clone_id_candidates) +cons_col <- pick_col(export_cells, params$consensus_cluster_col, params$consensus_cluster_candidates) + +has_gex <- !is.null(tcell_roll) || !is.null(cons_roll) || !is.null(conga_roll) || !is.null(tcri_roll) +route <- if (has_gex) "Full single-cell (GEX + VDJ)" else "VDJ-only (no GEX object)" ``` +::: {.callout-tip appearance="simple"} +**All figures in this report are interactive.** Hover a point or bar for exact values, drag +to zoom, double-click to reset, and click legend entries to show or hide series. Tables are +sortable and searchable, with adjustable page size. +::: + # Overview + +What was analysed, and how the run was routed. The route is detected from whether a GEX +annotated object was supplied: without one, the T-cell integration, CoNGA and consensus +steps are skipped and the sections that depend on them say so explicitly. + ```{r} #| label: overview -#| results: asis -overview_tbl <- tibble::tibble( - Metric = c( - "Project name", - "Cells in export table", - "Samples", - "Patients", - "Conditions", - "Timepoints", - "Unique clonotypes", - "Unique consensus clusters" - ), +ov <- data.frame( + Metric = c("Project", "Analysis route", "Samples", "Patients", + "Cells / clonotype rows", "Unique clonotypes", "TCRdist radius"), Value = c( params$project_name, - if (!is.null(export_cells)) nrow(export_cells) else NA, - if (!is.null(export_cells) && !is.null(sample_col)) length(unique(na.omit(export_cells[[sample_col]]))) else NA, - if (!is.null(export_cells) && !is.null(patient_col)) length(unique(na.omit(export_cells[[patient_col]]))) else NA, - if (!is.null(export_cells) && !is.null(condition_col)) length(unique(na.omit(export_cells[[condition_col]]))) else NA, - if (!is.null(export_cells) && !is.null(timepoint_col)) length(unique(na.omit(export_cells[[timepoint_col]]))) else NA, - if (!is.null(export_cells) && !is.null(clone_id_col)) length(unique(na.omit(export_cells[[clone_id_col]]))) else NA, - if (!is.null(export_cells) && !is.null(consensus_cluster_col)) length(unique(na.omit(export_cells[[consensus_cluster_col]]))) else NA - ) + route, + coalesce_count( + if (have(vdj_compact)) nrow(vdj_compact) else NA_integer_, + if (!is.null(export_cells) && !is.null(sample_col)) n_distinct_real(export_cells[[sample_col]]) else NA_integer_), + coalesce_count( + if (have(vdj_compact) && "patient_id" %in% names(vdj_compact)) n_distinct_real(vdj_compact$patient_id) else NA_integer_, + if (!is.null(export_cells) && !is.null(patient_col)) n_distinct_real(export_cells[[patient_col]]) else NA_integer_), + if (!is.null(export_cells)) format(nrow(export_cells), big.mark = ",") else NA, + if (!is.null(export_cells) && !is.null(clone_col)) format(dplyr::n_distinct(export_cells[[clone_col]]), big.mark = ",") else NA, + as.character(params$tcrdist_radius) + ), + stringsAsFactors = FALSE ) +save_table(ov, "master_overview.tsv") +itable(ov, pageSize = 10) +``` -save_table_safe(overview_tbl, "master_overview.tsv") +## Module coverage -print( - kableExtra::kbl(overview_tbl, caption = "Master project overview.") %>% - kableExtra::kable_styling(full_width = FALSE, bootstrap_options = c("striped","hover","condensed")) +Which analysis modules contributed output to this report. `FALSE` for the GEX-dependent +modules is expected on the VDJ-only route and does not indicate a failure. + +```{r} +#| label: module-availability +mods <- list( + "VDJ QC" = vdj_compact, + "Pseudobulk QC" = pb_qc, + "Sample statistics" = sample_stats, + "Clonotype sharing" = sharing, + "GIANA" = giana_roll, + "GLIPH2" = gliph2_roll, + "TCRdist3" = tcrdist_roll, + "Repertoire" = rep_rollup, + "T-cell integration" = tcell_roll, + "CoNGA" = conga_roll, + "Consensus" = cons_roll, + "TCRi" = tcri_roll ) +avail <- data.frame(module = names(mods), + available = vapply(mods, have, logical(1)), + row.names = NULL) +save_table(avail, "module_output_availability.tsv") + +d <- avail %>% + mutate(status = ifelse(available, "present", "absent"), + cat = module, val = as.integer(available), + hover = paste0(module, ": ", status)) +hbar(d, "Module output availability", "", colour = "#2c7fb8", height = 460) +``` + +# Sample QC and the pseudobulk gate +Before any clustering, every sample passes a QC gate on clonotype and cell counts. Samples +falling below `pseudobulk_qc_min_clones` / `pseudobulk_qc_min_cells` are dropped from all +downstream analysis, so this is the first place to look if a sample is missing later on. + +```{r} +#| label: pseudobulk-gate +if (!have(pb_qc)) { + note("No pseudobulk QC summary was staged for this run.") +} else { + save_table(pb_qc, "master_pseudobulk_qc.tsv") + d <- pb_qc %>% + transmute(x = n_clones, y = n_cells, grp = as.character(pass), + hover = paste0(sample, "
clones: ", n_clones, + "
cells: ", n_cells, "
", pass)) + p <- scat(d, "Pseudobulk QC gate (dashed lines = thresholds)", "Clonotypes", "Cells") + p %>% plotly::layout(shapes = list( + list(type = "line", x0 = pb_qc$min_clones[1], x1 = pb_qc$min_clones[1], + y0 = 0, y1 = max(pb_qc$n_cells) * 1.05, + line = list(dash = "dash", color = "grey")), + list(type = "line", x0 = 0, x1 = max(pb_qc$n_clones) * 1.05, + y0 = pb_qc$min_cells[1], y1 = pb_qc$min_cells[1], + line = list(dash = "dash", color = "grey")))) +} ``` -# Module Rollups ```{r} -#| label: module-rollups -#| results: asis -rollup_list <- list( - VDJ_QC = vdj_qc_summary, - TCell_Integration = tcell_integration_summary, - TCRi = tcri_summary, - CoNGA = conga_summary, - GLIPH2 = gliph2_summary, - TCRdist3 = tcrdist3_summary, - GIANA = giana_summary, - Consensus = consensus_summary, - Repertoire = repertoire_summary -) +#| label: pseudobulk-gate-table +itable(pb_qc, pageSize = 10) +``` -rollup_presence <- tibble::tibble( - module = names(rollup_list), - available = vapply(rollup_list, function(x) !is.null(x), logical(1)) -) +# VDJ QC -save_table_safe(rollup_presence, "module_output_availability.tsv") +Contig-level quality control from Cell Ranger output: how many contigs survived filtering, +how many cells carry a properly paired receptor, and the shape of the clone size +distribution. Low pairing or heavy multi-chain burden here propagates into every downstream +clustering result. -print( - kableExtra::kbl(rollup_presence, caption = "Module output availability.") %>% - kableExtra::kable_styling(full_width = FALSE, bootstrap_options = c("striped","hover","condensed")) -) +## Per-sample summary + +```{r} +#| label: vdj-compact +if (!have(vdj_compact)) note("VDJ QC did not contribute a per-sample summary.") else { + save_table(vdj_compact, "master_vdj_qc_per_sample.tsv") + itable(vdj_compact, pageSize = 10) +} ``` -# VDJ QC +## Contig retention + +Fraction of contigs surviving QC filtering, per sample. A sample retaining markedly less +than its peers usually indicates a library or sequencing-depth problem rather than biology. + ```{r} -#| label: vdj-qc-section -#| results: asis +#| label: vdj-retention +if (!have(vdj_retention)) note("No contig retention table available.") else { + d <- vdj_retention + ycol <- intersect(c("retained_pct", "retained_prop", "pct_retained"), names(d)) + if (!length(ycol)) note("Retention table present but has no recognisable retention column.") else { + v <- as.numeric(d[[ycol[1]]]); if (max(v, na.rm = TRUE) <= 1) v <- 100 * v + dd <- data.frame(cat = d$sample, val = v, + hover = paste0(d$sample, "
retained: ", round(v, 1), "%")) + hbar(dd, "Contig retention after QC", "Retained (%)") + } +} +``` -if (!is.null(vdj_qc_per_sample_compact) && nrow(vdj_qc_per_sample_compact) > 0) { - vdj_qc_per_sample_compact_master <- vdj_qc_per_sample_compact %>% - dplyr::select( - dplyr::any_of(c( - "sample", - "patient_id", - "timepoint", - "n_cells_total", - "n_cells_paired", - "pct_cells_paired", - "n_unique_clones" - )) - ) +## Chain pairing - save_table_safe(vdj_qc_per_sample_compact, "master_vdj_qc_per_sample_compact.tsv") - - print( - kableExtra::kbl( - vdj_qc_per_sample_compact_master, - caption = "VDJ QC per-sample compact summary." - ) %>% - kableExtra::kable_styling( - full_width = FALSE, - bootstrap_options = c("striped","hover","condensed"), - font_size = 11 - ) %>% - kableExtra::scroll_box(width = "100%", height = "420px") - ) +Share of cells by chain-pairing status. Only cells with a paired receptor contribute a +usable clonotype, so this sets the effective ceiling on everything downstream. + +```{r} +#| label: vdj-pairing +if (!have(vdj_pairing) || !all(c("sample", "pairing") %in% names(vdj_pairing))) { + note("No pairing status table available.") +} else { + d <- vdj_pairing %>% + count(sample, pairing, name = "cells") %>% + group_by(sample) %>% mutate(frac = cells / sum(cells)) %>% ungroup() + save_table(d, "master_pairing_by_sample.tsv") + dd <- d %>% transmute(cat = sample, val = frac, grp = pairing, + hover = paste0(sample, "
", pairing, "
", + cells, " cells (", percent(frac, accuracy = 0.1), ")")) + sbar(dd, "Chain pairing status by sample", "Fraction of cells") } +``` -if (!is.null(vdj_qc_before_after_summary) && nrow(vdj_qc_before_after_summary) > 0) { - save_table_safe(vdj_qc_before_after_summary, "master_vdj_qc_before_after_summary.tsv") - - print( - kableExtra::kbl( - vdj_qc_before_after_summary, - caption = "VDJ QC contigs before/after filtering summary." - ) %>% - kableExtra::kable_styling( - full_width = FALSE, - bootstrap_options = c("striped","hover","condensed"), - font_size = 11 - ) %>% - kableExtra::scroll_box(width = "100%", height = "320px") - ) +## CDR3 length distribution + +CDR3β length distribution across all productive contigs. A healthy repertoire is centred +near 14–15 amino acids; sharp spikes can indicate contamination or a dominant clone. + +```{r} +#| label: cdr3-length +if (!have(vdj_cdr3len)) note("No CDR3 length table available.") else { + d <- vdj_cdr3len + lc <- intersect(c("cdr3_length", "length"), names(d)) + if (!length(lc)) note("CDR3 length table has no length column.") else { + d$.len <- suppressWarnings(as.integer(d[[lc[1]]])) + if ("chain" %in% names(d)) d <- d[d$chain %in% c("TRB", "TRA", "TRD", "TRG"), , drop = FALSE] + d <- d[!is.na(d$.len) & d$.len > 0 & d$.len < 40, , drop = FALSE] + if (!nrow(d)) note("No usable CDR3 lengths after filtering.") else { + dd <- if ("chain" %in% names(d)) count(d, .len, chain, name = "n") else + cbind(count(d, .len, name = "n"), chain = "all") + dd <- dd %>% transmute(cat = .len, val = n, grp = chain, + hover = paste0(chain, "
length ", .len, "
", + format(n, big.mark = ","), " contigs")) + gbar(dd, "CDR3 length distribution", "CDR3 length (aa)", "Contigs") + } + } } +``` -if (!is.null(vdj_qc_sample_sheet_resolved) && nrow(vdj_qc_sample_sheet_resolved) > 0) { - save_table_safe(vdj_qc_sample_sheet_resolved, "master_vdj_qc_sample_sheet_resolved.tsv") +## V and J gene usage + +Germline gene usage per sample. Systematic shifts between samples can reflect real +biology, but can equally flag a batch or primer effect worth ruling out. + +```{r} +#| label: gene-usage +usage_box <- function(d, gcol, title) { + if (!have(d)) return(NULL) + fc <- intersect(c("freq", "frequency", "prop", "n"), names(d)) + if (!length(fc) || !(gcol %in% names(d))) return(NULL) + d$.f <- as.numeric(d[[fc[1]]]) + top <- d %>% group_by(.data[[gcol]]) %>% + summarise(m = mean(.f, na.rm = TRUE), .groups = "drop") %>% + slice_max(m, n = params$top_n_genes, with_ties = FALSE) %>% pull(.data[[gcol]]) + dd <- d[d[[gcol]] %in% top, , drop = FALSE] + ord <- dd %>% group_by(.data[[gcol]]) %>% summarise(m = median(.f), .groups = "drop") %>% + arrange(desc(m)) %>% pull(.data[[gcol]]) + out <- data.frame(grp = factor(dd[[gcol]], levels = ord), val = dd$.f, + hover = paste0(dd[[gcol]], "
", dd$sample, "
", signif(dd$.f, 3))) + boxby(out, title, "Usage frequency", height = 520, angle = -45) } +pv <- usage_box(vdj_vusage, "v_gene", paste0("Top ", params$top_n_genes, " V genes")) +if (is.null(pv)) note("No V gene usage table available.") else pv +``` -if (!is.null(vdj_qc_clone_rank_abundance) && nrow(vdj_qc_clone_rank_abundance) > 0) { - top_clone_rank <- vdj_qc_clone_rank_abundance %>% - dplyr::group_by(sample) %>% - dplyr::slice_head(n = 10) %>% - dplyr::ungroup() +```{r} +#| label: j-usage +pj <- usage_box(vdj_jusage, "j_gene", paste0("Top ", params$top_n_genes, " J genes")) +if (is.null(pj)) note("No J gene usage table available.") else pj +``` - save_table_safe(top_clone_rank, "master_vdj_qc_clone_rank_abundance_top.tsv") +# Repertoire diversity and clonal structure - print( - kableExtra::kbl( - top_clone_rank, - caption = "Top ranked clonotypes by sample." - ) %>% - kableExtra::kable_styling( - full_width = FALSE, - bootstrap_options = c("striped","hover","condensed"), - font_size = 11 - ) %>% - kableExtra::scroll_box(width = "100%", height = "320px") - ) +Diversity summarises how evenly the repertoire is distributed across clonotypes. Low +Shannon entropy with high clonality indicates a repertoire dominated by a few expanded +clones — typically the signal of interest in an antigen-driven response. + +## Diversity by sample + +```{r} +#| label: diversity +if (!have(diversity)) { + note("No diversity table available. This is produced by the Repertoire module.") +} else { + save_table(diversity, "master_diversity_by_sample.tsv") + metrics <- intersect(c("shannon", "simpson", "inv_simpson", "richness"), names(diversity)) + d <- diversity %>% + pivot_longer(all_of(metrics), names_to = "metric", values_to = "value") %>% + transmute(grp = metric, val = value, + hover = paste0(sample, "
", metric, ": ", signif(value, 4))) + boxby(d, "Repertoire diversity across samples", "Metric value") } +``` -fig_map <- list( - "QC retention before/after filtering" = params$vdj_qc_before_after_retention_fig, - "Pairing status by sample" = params$vdj_qc_pairing_bar_fig, - "Clone rank abundance" = params$vdj_qc_clone_rank_abundance_fig, - "Multiple chains by sample" = params$vdj_qc_multiple_chains_fig -) +```{r} +#| label: diversity-table +itable(diversity, pageSize = 10) +``` -valid_figs <- fig_map[ - vapply(fig_map, function(x) !is.na(x) && nzchar(x) && file.exists(x) && isTRUE(file.size(x) > 0), logical(1)) -] +## Sample statistics -if (length(valid_figs) > 0) { - cat("### Selected VDJ QC figures\n\n") +Clonality, productive fraction and convergence per sample, from the shared bulk engine. +Convergent recombination — distinct nucleotide sequences encoding the same CDR3 amino acid +sequence — is a useful independent marker of antigen-driven selection. - for (nm in names(valid_figs)) { - fp <- valid_figs[[nm]] - out_name <- basename(fp) - file.copy(fp, file.path(params$outdir, "figures", out_name), overwrite = TRUE) +```{r} +#| label: sample-stats +if (!have(sample_stats)) note("No sample statistics table available.") else { + save_table(sample_stats, "master_sample_stats.tsv") + itable(sample_stats, pageSize = 10) +} +``` - cat("#### ", nm, "\n\n", sep = "") - print(knitr::include_graphics(fp)) - cat("\n\n") +```{r} +#| label: clonality-plot +if (!have(sample_stats) || !("clonality" %in% names(sample_stats))) { + note("No clonality column in the sample statistics table.") +} else { + d <- sample_stats %>% + transmute(cat = sample, val = clonality, + hover = paste0(sample, "
clonality: ", signif(clonality, 3), + "
clonotypes: ", num_clones)) + hbar(d, "Clonality by sample", "Clonality", colour = "#7570b3") +} +``` + +## Clonal burden + +How much of each sample's cell mass sits in small versus expanded clones. A sample shifting +mass into the large-clone bins is undergoing clonal expansion. + +```{r} +#| label: clone-burden +if (!have(clone_burden)) { + note("No clone burden table available. This is produced by the Repertoire module.") +} else { + save_table(clone_burden, "master_clone_burden.tsv") + fc <- intersect(c("frac_cells", "frac_clones"), names(clone_burden)) + if (!("clone_bin" %in% names(clone_burden)) || !length(fc)) { + note("Clone burden table lacks the expected clone_bin / fraction columns.") + } else { + d <- clone_burden + d$.f <- as.numeric(d[[fc[1]]]) + dd <- d %>% transmute(cat = sample, val = .f, grp = clone_bin, + hover = paste0(sample, "
", clone_bin, "
", + percent(.f, accuracy = 0.1))) + sbar(dd, "Clonal burden across samples", "Fraction of cells") } } ``` +## Clone rank abundance + +Clone size against rank. A steep head indicates a few dominant clones; a long flat tail +indicates a diverse, largely unexpanded repertoire. -# Diversity ```{r} -if (!is.null(diversity_by_sample) && nrow(diversity_by_sample) > 0) { - div_long <- diversity_by_sample %>% - pivot_longer(cols = intersect(c("shannon","simpson","inv_simpson","richness"), colnames(diversity_by_sample)), - names_to = "metric", values_to = "value") +#| label: rank-abundance +ra <- if (have(vdj_rank)) vdj_rank else T("repertoire", "clone_rank_abundance.tsv") +if (!have(ra)) note("No clone rank abundance table available.") else { + rc <- intersect(c("rank", "clone_rank"), names(ra)) + nc <- intersect(c("n", "cells", "count", "clone_size", "size"), names(ra)) + if (!length(rc) || !length(nc)) note("Rank abundance table lacks rank/size columns.") else { + d <- ra + d$.r <- as.numeric(d[[rc[1]]]); d$.n <- as.numeric(d[[nc[1]]]) + d <- d[!is.na(d$.r) & !is.na(d$.n) & d$.n > 0 & d$.r > 0, , drop = FALSE] + grp <- if ("sample" %in% names(d)) as.character(d$sample) else "all" + dd <- data.frame(x = d$.r, y = d$.n, grp = grp, + hover = paste0(grp, "
rank ", d$.r, "
", d$.n, " cells")) + lplot(dd, "Clone rank abundance (log-log)", "Clone rank", "Clone size") + } +} +``` - p_div <- ggplot(div_long, aes(x = metric, y = value, fill = metric)) + - geom_boxplot(alpha = 0.85, outlier.size = 0.4) + - geom_jitter(width = 0.12, size = 1.2, alpha = 0.7) + - guides(fill = "none") + - labs( - title = "Repertoire diversity summary", - x = NULL, - y = "Metric value" - ) + - theme_scratch_pub(params$base_size) +# Repertoire overlap and clonotype sharing - p_div - save_plot_safe(p_div, "master_diversity_boxplots.png") +Clonotypes seen in more than one sample. Within a patient this reflects genuine clonal +persistence across timepoints or tissues; across patients it usually reflects public +clonotypes with high generation probability rather than a shared response. + +## Sample overlap + +```{r} +#| label: overlap-heatmap +if (!have(overlap_mat)) { + note("No sample overlap matrix available. This is produced by the Repertoire module.") +} else { + save_table(overlap_mat, "master_sample_overlap.tsv") + m <- overlap_mat + rn <- as.character(m[[1]]); m[[1]] <- NULL + mat <- as.matrix(m); rownames(mat) <- rn + plot_ly( + x = colnames(mat), y = rownames(mat), z = mat, + type = "heatmap", colors = colorRamp(c("white", "#fee391", "#cb181d")), + hovertemplate = "%{y} vs %{x}
overlap: %{z:.3f}" + ) %>% + layout(title = "Sample repertoire overlap", + xaxis = list(title = "", tickangle = -45), yaxis = list(title = "")) %>% + config(displaylogo = FALSE) } ``` -# Clone Birden +## Most shared clonotypes + ```{r} -if (!is.null(clone_burden) && nrow(clone_burden) > 0) { - sample_col_burden <- if ("sample" %in% colnames(clone_burden)) "sample" else colnames(clone_burden)[1] - frac_col <- if ("frac_cells" %in% colnames(clone_burden)) "frac_cells" else if ("frac_clones" %in% colnames(clone_burden)) "frac_clones" else NULL - bin_col <- if ("clone_bin" %in% colnames(clone_burden)) "clone_bin" else NULL +#| label: top-shared +if (!have(top_shared)) note("No top shared clones table available.") else { + save_table(top_shared, "master_top_shared_clones.tsv") + itable(top_shared, pageSize = 10) +} +``` + +## Sharing versus generation probability - if (!is.null(frac_col) && !is.null(bin_col)) { - p_burden <- ggplot(clone_burden, aes_string(x = sample_col_burden, y = frac_col, fill = bin_col)) + - geom_col(position = "fill") + - coord_flip() + - scale_y_continuous(labels = percent_format(accuracy = 1)) + - labs( - title = "Clonal burden across samples", - x = NULL, - y = "Fraction", - fill = "Clone bin" - ) + - theme_scratch_pub(params$base_size) +Each clonotype's sharing breadth against its OLGA generation probability. Clonotypes shared +widely *and* easy to generate are likely public/convergent rather than response-specific — +this plot is the quickest way to tell those apart. - p_burden - save_plot_safe(p_burden, "master_clone_burden.png") +```{r} +#| label: sharing-pgen +if (!have(sharing)) note("No clonotype sharing table available (produced by COMPARE).") else { + d <- sharing + pg <- intersect(c("log10_pgen", "pgen"), names(d)) + ns <- intersect(c("total_samples", "n_samples"), names(d)) + if (!length(pg) || !length(ns)) note("Sharing table lacks pgen / sample-count columns.") else { + d$.p <- as.numeric(d[[pg[1]]]); d$.n <- as.integer(d[[ns[1]]]) + if (pg[1] == "pgen") d$.p <- log10(pmax(d$.p, 1e-300)) + d <- d[is.finite(d$.p) & !is.na(d$.n), , drop = FALSE] + dd <- data.frame(grp = factor(d$.n, levels = sort(unique(d$.n))), val = d$.p, + hover = paste0(d$junction_aa, "
samples: ", d$.n, + "
log10 pgen: ", signif(d$.p, 4))) + boxby(dd, "Generation probability by sharing breadth", "log10 pgen") } } ``` -# Method coverage and cluster counts -```{r} -plot_list <- list() +# Sequence clustering -if (!is.null(method_presence) && nrow(method_presence) > 0) { - mp <- method_presence %>% - mutate(method = fct_reorder(method, frac_cells)) - p1 <- ggplot(mp, aes(x = method, y = frac_cells)) + - geom_col(fill = "steelblue") + - coord_flip() + - scale_y_continuous(labels = percent_format(accuracy = 1)) + - labs(title = "Method coverage", x = NULL, y = "Fraction of cells annotated") + - theme_scratch_pub(params$base_size) - plot_list$coverage <- p1 -} +Three complementary approaches to grouping clonotypes by receptor similarity: **GIANA** +(isometric distance on CDR3β + V gene), **GLIPH2** (shared local motif enrichment) and +**TCRdist3** (pairwise TCR distance, thresholded at radius `r params$tcrdist_radius`). +GIANA and GLIPH2 run per patient; TCRdist3 runs per sample. + +## Method coverage + +Fraction of cells whose clonotype was assigned to a cluster by each method, and how many +clusters each produced. Large disagreements between methods are expected — they optimise +for different notions of similarity — but a method at or near zero is worth investigating. -if (!is.null(method_cluster_counts) && nrow(method_cluster_counts) > 0) { - mc <- method_cluster_counts %>% - mutate(method = fct_reorder(method, n_clusters)) - p2 <- ggplot(mc, aes(x = method, y = n_clusters)) + - geom_col(fill = "darkorange") + - coord_flip() + - scale_y_continuous(labels = comma) + - labs(title = "Clusters detected by method", x = NULL, y = "Number of clusters") + - theme_scratch_pub(params$base_size) - plot_list$counts <- p2 +```{r} +#| label: method-coverage +if (!have(method_presence)) { + note("No method coverage table available.") +} else { + save_table(method_presence, "master_method_presence.tsv") + d <- method_presence %>% + transmute(cat = method, val = frac_cells, + hover = paste0(method, "
", percent(frac_cells, accuracy = 0.1), " of cells")) + hbar(d, "Cell coverage by method", "Fraction of cells", height = 380, tickformat = ".0%") } +``` -if (length(plot_list) == 2) { - wrap_plots(plot_list$coverage, plot_list$counts, ncol = 2) - save_plot_safe(wrap_plots(plot_list$coverage, plot_list$counts, ncol = 2), "master_method_panels.png", width = 14, height = 6) -} else if (length(plot_list) == 1) { - print(plot_list[[1]]) +```{r} +#| label: method-cluster-counts +if (!have(method_counts)) { + note("No cluster count table available.") +} else { + save_table(method_counts, "master_method_cluster_counts.tsv") + d <- method_counts %>% + transmute(cat = method, val = n_clusters, + hover = paste0(method, "
", n_clusters, " clusters")) + hbar(d, "Clusters detected by method", "Clusters", colour = "#d95f02", height = 380) } ``` -# Sample Overlap +## Per-method summaries + ```{r} -if (!is.null(sample_overlap_matrix) && nrow(sample_overlap_matrix) > 0) { - som <- as.data.frame(sample_overlap_matrix) - rownames(som) <- som[[1]] - som[[1]] <- NULL - som <- as.matrix(som) +#| label: cluster-rollups +rolls <- list(GIANA = giana_roll, GLIPH2 = gliph2_roll, TCRdist3 = tcrdist_roll) +rolls <- rolls[vapply(rolls, have, logical(1))] +if (!length(rolls)) note("No clustering rollups available.") else { + combined <- bind_rows(lapply(names(rolls), function(n) + data.frame(method = n, rolls[[n]], stringsAsFactors = FALSE))) + save_table(combined, "master_cluster_rollups.tsv") + itable(combined, pageSize = 15, + groupBy = "method", + columns = list(method = colDef(name = "Method"))) +} +``` - ht <- Heatmap( - som, - name = "Overlap", - col = colorRamp2(c(0, 0.5, 1), c("white", "gold", "firebrick")), - cluster_rows = TRUE, - cluster_columns = TRUE, - row_names_side = "left", - column_title = "Sample repertoire overlap", - heatmap_legend_param = list(title = "Overlap") - ) - draw(ht) +::: {.callout-note appearance="simple"} +GIANA receives one row per unique clonotype. Before this was enforced, the same clonotype +seen in several samples of a patient arrived as duplicate rows and GIANA reported those +duplicates as clusters — inflating cluster counts and producing no output at all for +patients with a single sample. +::: + +## Per-annotation clustering coverage + +Fraction of cells clustered by each method, broken down by cell-type annotation. Requires +cell-type labels from a GEX object. + +```{r} +#| label: annotation-panels +panels <- list(GIANA = anno_giana, GLIPH2 = anno_gliph2, TCRdist3 = anno_tcrdist, + CoNGA = anno_conga, Consensus = anno_cons, TCRi = anno_tcri) +rows <- list() +for (nm in names(panels)) { + d <- panels[[nm]] + if (!have(d)) next + ac <- if ("annot" %in% names(d)) "annot" else names(d)[1] + fc <- intersect(c("frac_giana", "frac_gliph", "frac_tcrdist", "frac_clustered", + "frac_high", "frac_consensus", "frac_tcri", "frac_high_tcri"), names(d)) + if (!length(fc)) next + rows[[nm]] <- data.frame(method = nm, annot = as.character(d[[ac]]), + frac = as.numeric(d[[fc[1]]]), stringsAsFactors = FALSE) +} +if (!length(rows)) { + note(paste("No per-annotation summaries available. These require cell-type labels from a", + "GEX object, so this section is empty on the VDJ-only route.")) +} else { + d <- bind_rows(rows) %>% + group_by(method) %>% slice_max(frac, n = params$top_n_annotations, with_ties = FALSE) %>% + ungroup() + save_table(d, "master_annotation_coverage.tsv") + dd <- d %>% transmute(cat = annot, val = frac, grp = method, + hover = paste0(method, "
", annot, "
", + percent(frac, accuracy = 0.1))) + gbar(dd, "Clustering coverage by cell-type annotation", "", "Fraction of cells", height = 560) } ``` -# Annotation-centered comparison +# GEX-dependent analyses + +These sections require an annotated Seurat object (`--input_annotated_object`) and are +skipped entirely on the VDJ-only route. + +## T-cell integration + +Fusion of VDJ contigs with the GEX object, matching receptors to transcriptome-defined cell +states. + ```{r} -annotation_panels <- list( - TCRi = annotation_tcri_summary, - CoNGA = annotation_conga_summary, - GLIPH2 = annotation_gliph2_summary, - TCRdist3 = annotation_tcrdist3_summary, - GIANA = annotation_giana_summary, - Consensus = annotation_consensus_summary -) +#| label: tcell +if (!have(tcell_roll)) { + note("T-cell integration did not run on this route (no GEX object supplied).") +} else { + save_table(tcell_roll, "master_tcell_rollup.tsv") + itable(tcell_roll, pageSize = 15) +} +``` -anno_comp <- lapply(names(annotation_panels), function(nm) { - df <- annotation_panels[[nm]] - if (is.null(df) || nrow(df) == 0) return(NULL) +```{r} +#| label: tcell-annot +if (!have(tcell_anno)) { + note("No per-annotation T-cell integration summary available.") +} else { + d <- tcell_anno + nc <- intersect(c("n_cells", "cells", "n"), names(d)) + ac <- if ("annot" %in% names(d)) "annot" else names(d)[1] + if (!length(nc)) note("T-cell annotation summary lacks a cell-count column.") else { + dd <- data.frame(cat = as.character(d[[ac]]), val = as.numeric(d[[nc[1]]])) + dd$hover <- paste0(dd$cat, "
", format(dd$val, big.mark = ","), " cells") + hbar(dd, "Cells per annotation", "Cells", height = 560) + } +} +``` - annot_col <- if ("annot" %in% colnames(df)) "annot" else colnames(df)[1] - frac_col <- intersect(c("frac_high","frac_clustered","frac_gliph","frac_tcrdist","frac_giana","frac_consensus"), colnames(df)) - if (length(frac_col) == 0) return(NULL) +## CoNGA - df %>% - transmute( - method = nm, - annot = .data[[annot_col]], - frac = .data[[frac_col[1]]] - ) -}) +Graph-level correspondence between TCR similarity and transcriptional similarity. Clusters +scoring highly here contain cells that are neighbours in *both* spaces — the strongest +available evidence of antigen-driven clonal programmes. -anno_comp <- dplyr::bind_rows(anno_comp) +```{r} +#| label: conga +if (!have(conga_roll) && !have(conga_clusters)) { + note("CoNGA did not run on this route (requires a GEX object).") +} else { + if (have(conga_roll)) save_table(conga_roll, "master_conga_rollup.tsv") + if (have(conga_clusters)) save_table(conga_clusters, "master_conga_clusters.tsv") + itable(if (have(conga_roll)) conga_roll else conga_clusters, pageSize = 15) +} +``` -if (nrow(anno_comp) > 0) { - anno_comp2 <- anno_comp %>% - group_by(method) %>% - slice_max(order_by = frac, n = params$top_n_annotations, with_ties = FALSE) %>% - ungroup() +## Consensus clustering - p_anno <- ggplot(anno_comp2, aes(x = reorder(annot, frac), y = frac, fill = method)) + - geom_col() + - coord_flip() + - facet_wrap(~method, scales = "free_y") + - scale_y_continuous(labels = percent_format(accuracy = 1)) + - labs( - title = "Top annotation enrichments across modules", - x = NULL, - y = "Fraction / enrichment summary" - ) + - theme_scratch_pub(params$base_size) +Reconciliation of the GIANA, GLIPH2, TCRdist3 and CoNGA labels into a single consensus +cluster per cell. Consensus clusters supported by several independent methods are the most +trustworthy candidates for follow-up. - p_anno - save_plot_safe(p_anno, "master_annotation_comparison.png", width = 14, height = 10) +```{r} +#| label: consensus +if (!have(cons_clusters)) { + note("Consensus clustering did not run on this route (requires a GEX object).") +} else { + save_table(cons_clusters, "master_consensus_clusters.tsv") + d <- cons_clusters + nc <- intersect(c("n_cells", "cells", "n"), names(d)) + cc <- intersect(c("consensus_cluster", "cluster"), names(d)) + if (!length(nc) || !length(cc)) itable(d, pageSize = 15) else { + dd <- data.frame(cat = as.character(d[[cc[1]]]), val = as.numeric(d[[nc[1]]])) + dd <- dd[order(-dd$val), , drop = FALSE] + dd <- head(dd, params$top_n_clusters) + dd$hover <- paste0("cluster ", dd$cat, "
", dd$val, " cells") + hbar(dd, paste0("Top ", params$top_n_clusters, " consensus clusters"), "Cells", height = 560) + } } ``` -# Concensus Cluster Burden ```{r} -if (!is.null(consensus_cluster_summary) && nrow(consensus_cluster_summary) > 0) { - ccs <- consensus_cluster_summary %>% - slice_head(n = params$top_n_consensus_clusters) %>% - mutate(consensus_cluster = fct_reorder(consensus_cluster, n_cells)) +#| label: consensus-table +itable(cons_clusters, pageSize = 10) +``` - p_cons <- ggplot(ccs, aes(x = consensus_cluster, y = n_cells)) + - geom_col(fill = "steelblue") + - coord_flip() + - scale_y_continuous(labels = comma) + - labs( - title = "Top consensus clonotype clusters", - x = NULL, - y = "Cells" - ) + - theme_scratch_pub(params$base_size) +## TCRi - p_cons - save_plot_safe(p_cons, "master_consensus_clusters.png") +Immunogenicity scoring: how strongly each clone's receptor resembles known +antigen-reactive receptors. High-scoring clones are the candidates worth following up +experimentally. Requires the GEX object, so this section is empty on the VDJ-only route. + +```{r} +#| label: tcri +#| results: asis +if (!have(tcri_roll)) { + note("TCRi did not run on this route (requires a GEX object).") +} else { + save_table(tcri_roll, "master_tcri_rollup.tsv") + itable(tcri_roll, pageSize = 15) } ``` -# key tables ```{r} -if (!is.null(consensus_cluster_summary) && nrow(consensus_cluster_summary) > 0) { - kable(consensus_cluster_summary %>% slice_head(n = params$top_n_consensus_clusters), - caption = "Top consensus clusters.") %>% - kable_styling(full_width = TRUE, bootstrap_options = c("striped","hover","condensed","responsive")) +#| label: tcri-clones +if (!have(tcri_scores)) { + note("No high-TCRi clone summary available.") +} else { + save_table(tcri_scores, "master_tcri_high_clones.tsv") + d <- tcri_scores + nc <- intersect(c("n_cells", "cells", "clone_size", "n"), names(d)) + cc <- intersect(c("clone_id", "CTaa", "clonotype"), names(d)) + if (!length(nc) || !length(cc)) { + itable(d, pageSize = 10) + } else { + dd <- data.frame(cat = as.character(d[[cc[1]]]), val = as.numeric(d[[nc[1]]])) + dd <- head(dd[order(-dd$val), , drop = FALSE], params$top_n_clusters) + dd$hover <- paste0(dd$cat, "
", dd$val, " cells") + hbar(dd, paste0("Top ", params$top_n_clusters, " high-TCRi clones"), "Cells", height = 560) + } } +``` -if (!is.null(diversity_by_sample) && nrow(diversity_by_sample) > 0) { - kable(diversity_by_sample, - caption = "Diversity metrics by sample.") %>% - kable_styling(full_width = TRUE, bootstrap_options = c("striped","hover","condensed","responsive")) +# Reference tables + +Every table rendered above is also written to `Master_Summary_Report/tables/` as TSV for +downstream use. + +```{r} +#| label: written-tables +files <- list.files(file.path(params$outdir, "tables"), pattern = "\\.tsv$") +if (!length(files)) note("No tables were written.") else { + itable(data.frame(table = sort(files), stringsAsFactors = FALSE), pageSize = 25) } ``` -# session info +# Session info + ```{r} -writeLines(capture.output(sessionInfo()), file.path(params$outdir, "sessionInfo.master_summary.txt")) +#| label: session +writeLines(capture.output(sessionInfo()), + file.path(params$outdir, "sessionInfo.master_summary.txt")) sessionInfo() ``` diff --git a/modules/scratch/MASTER_SUMMARY/main.nf b/modules/scratch/MASTER_SUMMARY/main.nf index 282ad04..7309e3b 100644 --- a/modules/scratch/MASTER_SUMMARY/main.nf +++ b/modules/scratch/MASTER_SUMMARY/main.nf @@ -6,22 +6,26 @@ process MASTER_SUMMARY { publishDir "${params.outdir}/Master_Summary", mode: 'copy', overwrite: true input: - // stageAs unique names: several of these are optional and default to the shared NO_FILE - // placeholder in VDJ-only mode. Without distinct staged names, Nextflow errors on - // "input file name collision" when two inputs resolve to the same NO_FILE. The .qmd - // treats any 0-byte staged file as absent. - path seurat_rds, stageAs: 'in_seurat_rds' - path export_cells, stageAs: 'in_export_cells' - - path vdj_qc_per_sample_compact, stageAs: 'in_vdj_qc_per_sample_compact' - path vdj_qc_before_after_summary, stageAs: 'in_vdj_qc_before_after_summary' - path vdj_qc_sample_sheet_resolved, stageAs: 'in_vdj_qc_sample_sheet_resolved' - path vdj_qc_clone_rank_abundance, stageAs: 'in_vdj_qc_clone_rank_abundance' - - path vdj_qc_before_after_retention_fig, stageAs: 'in_vdj_qc_before_after_retention_fig' - path vdj_qc_pairing_bar_fig, stageAs: 'in_vdj_qc_pairing_bar_fig' - path vdj_qc_clone_rank_abundance_fig, stageAs: 'in_vdj_qc_clone_rank_abundance_fig' - path vdj_qc_multiple_chains_fig, stageAs: 'in_vdj_qc_multiple_chains_fig' + // Core per-cell inputs. stageAs with a real extension matters: the .qmd picks its + // reader from the extension, and knitr/quarto infer image type the same way. + path seurat_rds, stageAs: 'in_seurat_rds.rds' + path export_cells, stageAs: 'in_export_cells.tsv' + + // Each upstream module contributes its whole tables/ directory under + // intables//. This replaces the previous ~30 individually-wired file params: + // a module that did not run on this route simply stages nothing, the .qmd finds no + // files, and the affected sections render an explanatory note. Adding a new upstream + // table no longer requires another process input. + path vdj_qc_tables, stageAs: 'intables/vdj_qc/*' + path pseudobulk_tables, stageAs: 'intables/pseudobulk/*' + path sample_tables, stageAs: 'intables/sample/*' + path compare_tables, stageAs: 'intables/compare/*' + path rollup_tables, stageAs: 'intables/rollup/*' + path repertoire_tables, stageAs: 'intables/repertoire/*' + path tcell_tables, stageAs: 'intables/tcell/*' + path conga_tables, stageAs: 'intables/conga/*' + path consensus_tables, stageAs: 'intables/consensus/*' + path tcri_tables, stageAs: 'intables/tcri/*' path qmd val barrier_done @@ -29,33 +33,24 @@ process MASTER_SUMMARY { output: path "Master_Summary_Report.html", emit: report_html - path "Master_Summary_Report/tables/*", emit: tables, optional: true + path "Master_Summary_Report/tables/*", emit: tables, optional: true path "Master_Summary_Report/figures/*", emit: figures, optional: true script: """ - mkdir -p Master_Summary_Report - mkdir -p Master_Summary_Report/data - mkdir -p Master_Summary_Report/tables - mkdir -p Master_Summary_Report/figures + mkdir -p Master_Summary_Report/tables Master_Summary_Report/figures quarto render ${qmd} \\ -P project_name="${project_name}" \\ + -P tables_dir="intables" \\ -P seurat_rds="${seurat_rds}" \\ -P export_cells_file="${export_cells}" \\ - -P vdj_qc_per_sample_compact_file="${vdj_qc_per_sample_compact}" \\ - -P vdj_qc_before_after_summary_file="${vdj_qc_before_after_summary}" \\ - -P vdj_qc_sample_sheet_resolved_file="${vdj_qc_sample_sheet_resolved}" \\ - -P vdj_qc_clone_rank_abundance_file="${vdj_qc_clone_rank_abundance}" \\ - -P vdj_qc_before_after_retention_fig="${vdj_qc_before_after_retention_fig}" \\ - -P vdj_qc_pairing_bar_fig="${vdj_qc_pairing_bar_fig}" \\ - -P vdj_qc_clone_rank_abundance_fig="${vdj_qc_clone_rank_abundance_fig}" \\ - -P vdj_qc_multiple_chains_fig="${vdj_qc_multiple_chains_fig}" \\ + -P outdir="Master_Summary_Report" \\ -P label_col="${params.label_col}" \\ -P sample_col="${params.sample_col}" \\ -P patient_col="${params.patient_col}" \\ -P condition_col="${params.condition_col}" \\ -P timepoint_col="${params.timepoint_col}" \\ - -P outdir="Master_Summary_Report" + -P tcrdist_radius=${params.tcrdist_radius} """ -} \ No newline at end of file +} diff --git a/modules/scratch/TCRI/main.nf b/modules/scratch/TCRI/main.nf index 700a720..850b33b 100644 --- a/modules/scratch/TCRI/main.nf +++ b/modules/scratch/TCRI/main.nf @@ -32,7 +32,10 @@ process TCRI { export QUARTO_PRINT_STACK=true export HOME="\$PWD" - export LD_LIBRARY_PATH="/opt/conda/envs/tcrenv/lib:\$LD_LIBRARY_PATH" + # :- default is required: Nextflow runs task scripts under `set -u`, and + # LD_LIBRARY_PATH is unset in this image, so a bare \$LD_LIBRARY_PATH aborts the + # task with "unbound variable" before anything runs. + export LD_LIBRARY_PATH="/opt/conda/envs/tcrenv/lib:\${LD_LIBRARY_PATH:-}" export RETICULATE_PYTHON="/opt/conda/envs/tcrenv/bin/python" echo "Testing Python environment natively..." diff --git a/modules/scratch/VDJ_QC/VDJ_QC_analysis.qmd b/modules/scratch/VDJ_QC/VDJ_QC_analysis.qmd index bd9cbd4..27c000a 100644 --- a/modules/scratch/VDJ_QC/VDJ_QC_analysis.qmd +++ b/modules/scratch/VDJ_QC/VDJ_QC_analysis.qmd @@ -1144,6 +1144,32 @@ patch_before_after save_plot_safe(patch_before_after, glue("qc_before_after_retention.{params$figure_format}"), width = 16, height = 5) ``` +```{r} +#| label: pairing-bar-by-sample +# Chain-pairing status per sample. This module already computes pairing_all and writes +# pairing_status_all.tsv but never plotted it - and the Master Summary requests exactly +# this figure by name (vdj_qc_pairing_bar_fig), so that panel was always empty. +if (exists("pairing_all") && !is.null(pairing_all) && nrow(pairing_all) > 0) { + pairing_bar <- pairing_all %>% + count(sample, pairing, name = "cells") %>% + group_by(sample) %>% + mutate(frac = cells / sum(cells)) %>% + ungroup() + + p_pairing <- ggplot(pairing_bar, aes(x = sample, y = frac, fill = pairing)) + + geom_col() + + coord_flip() + + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + + labs(title = "Chain pairing status by sample", + x = NULL, y = "Fraction of cells", fill = "Pairing") + + theme_scratch_pub(params$base_size) + + print(p_pairing) + save_plot_safe(p_pairing, glue("pairing_bar_by_sample.{params$figure_format}"), + width = 12, height = 7) +} +``` + ```{r} #| label: read-umi-plots if (isTRUE(params$show_reads_umis) && any(!is.na(contigs_all_post$reads))) { diff --git a/nextflow.config b/nextflow.config index f4a3370..62e3973 100644 --- a/nextflow.config +++ b/nextflow.config @@ -49,6 +49,12 @@ params { // GLIPH2 parameters use_gliph2 = true + // Collapse duplicate clonotype rows before GIANA. PATIENT_CONCATENATE stacks a + // patient's per-sample rows, so GIANA otherwise clusters duplicates rather than + // similar sequences. Default false keeps existing bulk results unchanged; the + // single-cell route sets it true. + giana_dedup_clonotypes = false + local_min_pvalue = "0.001" simulation_depth = "1000" kmer_min_depth = "3" @@ -118,6 +124,8 @@ params { run_repertoire = true run_master_summary = true run_tcri = true + // Optional pre-computed TCRi scores; when unset the report computes them itself. + tcri_scores_file = null // Global plotting / embedding metadata_file = "${projectDir}/assets/NO_FILE" diff --git a/params_singlecell.yml b/params_singlecell.yml index 1f0e581..1fe6c10 100644 --- a/params_singlecell.yml +++ b/params_singlecell.yml @@ -30,6 +30,10 @@ pseudobulk_qc_mode: "drop" # drop | hard_stop workflow_level: "sample,patient,compare" use_gliph2: true +# Collapse duplicate clonotype rows before GIANA. PATIENT_CONCATENATE stacks a patient's +# per-sample rows, so without this GIANA clusters duplicates rather than similar sequences. +giana_dedup_clonotypes: true + # ── Single-cell report toggles (full-SC route) ─────────────────────────────── run_conga: true run_consensus: true diff --git a/subworkflows/scratch/master_summary.nf b/subworkflows/scratch/master_summary.nf index 7318bc4..c54bde0 100644 --- a/subworkflows/scratch/master_summary.nf +++ b/subworkflows/scratch/master_summary.nf @@ -3,20 +3,28 @@ nextflow.enable.dsl = 2 include { MASTER_SUMMARY } from '../../modules/scratch/MASTER_SUMMARY/main.nf' +/* + * MASTER_SUMMARY_SW + * + * Runs on both single-cell routes. Each `*_tables` input is the collected tables/ + * directory of one upstream module, or an empty list on a route where that module did + * not run - the report then explains the absence instead of rendering a blank section. + */ workflow MASTER_SUMMARY_SW { take: seurat_rds export_cells - vdj_qc_per_sample_compact - vdj_qc_before_after_summary - vdj_qc_sample_sheet_resolved - vdj_qc_clone_rank_abundance - - vdj_qc_before_after_retention_fig - vdj_qc_pairing_bar_fig - vdj_qc_clone_rank_abundance_fig - vdj_qc_multiple_chains_fig + vdj_qc_tables + pseudobulk_tables + sample_tables + compare_tables + rollup_tables + repertoire_tables + tcell_tables + conga_tables + consensus_tables + tcri_tables barrier_done project_name @@ -30,14 +38,16 @@ workflow MASTER_SUMMARY_SW { MASTER_SUMMARY( seurat_rds, export_cells, - vdj_qc_per_sample_compact, - vdj_qc_before_after_summary, - vdj_qc_sample_sheet_resolved, - vdj_qc_clone_rank_abundance, - vdj_qc_before_after_retention_fig, - vdj_qc_pairing_bar_fig, - vdj_qc_clone_rank_abundance_fig, - vdj_qc_multiple_chains_fig, + vdj_qc_tables, + pseudobulk_tables, + sample_tables, + compare_tables, + rollup_tables, + repertoire_tables, + tcell_tables, + conga_tables, + consensus_tables, + tcri_tables, ch_notebook, barrier_done, project_name diff --git a/subworkflows/scratch/tcri.nf b/subworkflows/scratch/tcri.nf index 35be62f..3394722 100644 --- a/subworkflows/scratch/tcri.nf +++ b/subworkflows/scratch/tcri.nf @@ -21,4 +21,8 @@ workflow TCRI_SW { report_html = ch_tcri.report_html seurat_with_tcri = ch_tcri.seurat_with_tcri export_cells = ch_tcri.export_cells + // Re-export the module's tables/figures so the Master Summary can pick up + // tcri_summary_rollup.tsv and annotation_tcri_summary.tsv. + tables = ch_tcri.tables + figures = ch_tcri.figures } diff --git a/workflows/tcrtoolkit_sc.nf b/workflows/tcrtoolkit_sc.nf index 163d3f5..da681e1 100644 --- a/workflows/tcrtoolkit_sc.nf +++ b/workflows/tcrtoolkit_sc.nf @@ -27,6 +27,7 @@ include { TCELL_INTEGRATION_SW } from '../subworkflows/scratch/tcell_integration include { CONGA_SW } from '../subworkflows/scratch/conga.nf' include { CONSENSUS_SW } from '../subworkflows/scratch/consensus_clustering.nf' include { REPERTOIRE_SW } from '../subworkflows/scratch/repertoire.nf' +include { TCRI_SW } from '../subworkflows/scratch/tcri.nf' include { MASTER_SUMMARY_SW } from '../subworkflows/scratch/master_summary.nf' // ── Bridges (SC ↔ bulk-engine schema conversion) ────────────────────────── @@ -35,6 +36,7 @@ include { VDJ_TO_BULK_SW } from '../subworkflows/bridges/vdj_to_bulk.nf' include { CLUSTER_TO_SC_SW } from '../subworkflows/bridges/cluster_to_sc.nf' include { SC_SAMPLE_STATS } from '../modules/bridges/sc_sample_stats.nf' include { BULK_TO_EXPORT } from '../modules/bridges/bulk_to_export.nf' +include { CLUSTER_ROLLUP } from '../modules/bridges/cluster_rollup.nf' // ── Shared bulk-TCR analysis engine (main/local — unchanged behavior) ───── include { PSEUDOBULK_QC_SW } from '../subworkflows/local/pseudobulk_qc.nf' @@ -45,6 +47,19 @@ include { BULKTCR_ANALYSIS } from '../subworkflows/local/bulktcr_analysis.nf' // ("`enabled` is not defined") - a plain top-level function is. def enabled(x) { x == null || x == true } +// Pull a named file out of a module's collected tables/ channel, falling back to NO_FILE +// so a route where that module did not run still supplies something the report can read +// as "absent" (0 bytes). +def pickFile(ch, fname, nofile) { + ch.flatten().filter { it.name == fname }.first().ifEmpty(nofile) +} + +// Same with a second source. concat() preserves order, so .first() deterministically +// prefers the primary file and falls back otherwise - unlike mix(), which would race. +def pickFileOr(ch, fname, alt, nofile) { + ch.flatten().filter { it.name == fname }.concat(alt).first().ifEmpty(nofile) +} + workflow TCRTOOLKIT_SC { def nofile = file("${projectDir}/assets/NO_FILE") @@ -123,6 +138,12 @@ workflow TCRTOOLKIT_SC { // clonotype-level export is synthesized instead (below). conga_report = channel.empty() consensus_report = channel.empty() + // Collected tables/ per module for the Master Summary; stay empty where skipped. + tcell_tables = channel.empty() + conga_tables = channel.empty() + consensus_tables = channel.empty() + tcri_tables = channel.empty() + tcri_report = channel.empty() if (!vdjOnly) { // Reuse tcrdist3 outputs computed inside SAMPLE (no second run). @@ -135,10 +156,20 @@ workflow TCRTOOLKIT_SC { BULKTCR_ANALYSIS.out.tcrdist_output.map { _meta, f -> f } ) enriched_seurat = CLUSTER_TO_SC_SW.out.enriched_seurat + tcell_tables = tcell_out.tables if (enabled(params.run_conga)) { conga_out = CONGA_SW( enriched_seurat, tcell_out.export_cells, ch_project_name ) conga_report = conga_out.report_html + conga_tables = conga_out.tables + } + + // TCRi: immunogenicity scoring on the enriched Seurat. GEX-gated - it needs the + // transcriptome object, so it cannot run on the VDJ-only route. + if (enabled(params.run_tcri)) { + tcri_out = TCRI_SW( enriched_seurat, tcell_out.export_cells, ch_project_name ) + tcri_tables = tcri_out.tables + tcri_report = tcri_out.report_html } if (enabled(params.run_consensus)) { @@ -150,6 +181,7 @@ workflow TCRTOOLKIT_SC { gliph2_export, tcrdist_export, giana_export, ch_project_name ) consensus_report = CONSENSUS_SW.out.report_html + consensus_tables = CONSENSUS_SW.out.tables } rep_seurat = enabled(params.run_consensus) ? CONSENSUS_SW.out.seurat_with_consensus : enriched_seurat @@ -167,24 +199,71 @@ workflow TCRTOOLKIT_SC { // Cell-level + full with a GEX object; clonotype-level repertoire and a CoNGA-excluded // summary without one (only CoNGA and the cell-cluster mapping are truly GEX-gated). repertoire_report = channel.empty() + repertoire_tables = channel.empty() if (enabled(params.run_repertoire)) { REPERTOIRE_SW( rep_seurat, rep_export, ch_project_name ) repertoire_report = REPERTOIRE_SW.out.report_html + repertoire_tables = REPERTOIRE_SW.out.tables } + // ── Step 7b: rollups for GIANA / GLIPH2 / TCRdist3 ──────────────────── + // These three write raw per-patient / per-sample output but no rollup table, so the + // Master Summary had nothing to read for them and they always reported as absent. + CLUSTER_ROLLUP( + BULKTCR_ANALYSIS.out.giana_clusters.collect().ifEmpty([]), + BULKTCR_ANALYSIS.out.gliph2_cluster_details.collect().ifEmpty([]), + BULKTCR_ANALYSIS.out.tcrdist_clone_df.collect().ifEmpty([]), + BULKTCR_ANALYSIS.out.tcrdist_output.map { _meta, f -> f }.collect().ifEmpty([]), + rep_export, + channel.fromPath("${projectDir}/bin/cluster_rollup.py", checkIfExists: true), + ch_project_name + ) + if (enabled(params.run_master_summary)) { master_barrier = channel.empty() - .mix(conga_report, consensus_report, repertoire_report) + .mix(conga_report, consensus_report, repertoire_report, tcri_report) .collect() .ifEmpty([nofile]) + // Aggregated per-sample stats, built here rather than emitted from the shared + // engine so no bulk-side file needs changing. + def sample_stats_agg = BULKTCR_ANALYSIS.out.sample_csv + .collectFile(name: "sample_stats.csv", keepHeader: true, skip: 1, sort: true) + + def rollup_tables = channel.empty() + .mix(CLUSTER_ROLLUP.out.giana_summary, + CLUSTER_ROLLUP.out.gliph2_summary, + CLUSTER_ROLLUP.out.tcrdist3_summary, + CLUSTER_ROLLUP.out.method_presence, + CLUSTER_ROLLUP.out.method_cluster_counts, + CLUSTER_ROLLUP.out.annotation_giana, + CLUSTER_ROLLUP.out.annotation_gliph2, + CLUSTER_ROLLUP.out.annotation_tcrdist3) + .collect().ifEmpty([]) + + def pseudobulk_tables = channel.empty() + .mix(PSEUDOBULK_QC_SW.out.qc_summary, PSEUDOBULK_QC_SW.out.v_family) + .collect().ifEmpty([]) + + def sample_tables = channel.empty() + .mix(sample_stats_agg, BULKTCR_ANALYSIS.out.v_family, BULKTCR_ANALYSIS.out.j_family) + .collect().ifEmpty([]) + MASTER_SUMMARY_SW( - rep_seurat, rep_export, - vdj_qc_per_sample_compact, vdj_qc_before_after_summary, - vdj_qc_sample_sheet_resolved, vdj_qc_clone_rank_abundance, - vdj_qc_before_after_retention_fig, vdj_qc_pairing_bar_fig, - vdj_qc_clone_rank_abundance_fig, vdj_qc_multiple_chains_fig, - master_barrier, ch_project_name + rep_seurat, + rep_export, + vdj_qc_out.qc_tables.flatten().collect().ifEmpty([]), + pseudobulk_tables, + sample_tables, + BULKTCR_ANALYSIS.out.shared_cdr3.flatten().collect().ifEmpty([]), + rollup_tables, + repertoire_tables.flatten().collect().ifEmpty([]), + tcell_tables.flatten().collect().ifEmpty([]), + conga_tables.flatten().collect().ifEmpty([]), + consensus_tables.flatten().collect().ifEmpty([]), + tcri_tables.flatten().collect().ifEmpty([]), + master_barrier, + ch_project_name ) } } From 65908f8ddfa8f3bb5becf80807322b2d0de41d38 Mon Sep 17 00:00:00 2001 From: Syed Shujaat Ali Zaidi Date: Mon, 31 Aug 2026 20:14:18 -0500 Subject: [PATCH 2/7] Cirro: expose single-cell params and warn on a missing patient column - process-form.json / process-input.json: add patient_col, giana_dedup_clonotypes and tcrdist_radius to both single-cell profiles; run_tcri on full_gex only (it needs a GEX object). - preprocess.py: warn when the sample sheet has no patient column. PATIENT pools samples by it for GIANA/GLIPH2, so without it every sample becomes its own patient and the clustering is silently wrong. Also logs the resolved sample/patient counts. The check is defensive and cannot fail preprocessing. --- .cirro/singlecell_full_gex/preprocess.py | 27 +++++++++++++++++++ .cirro/singlecell_full_gex/process-form.json | 24 +++++++++++++++++ .cirro/singlecell_full_gex/process-input.json | 4 +++ .cirro/singlecell_vdjonly/preprocess.py | 27 +++++++++++++++++++ .cirro/singlecell_vdjonly/process-form.json | 18 +++++++++++++ .cirro/singlecell_vdjonly/process-input.json | 3 +++ 6 files changed, 103 insertions(+) diff --git a/.cirro/singlecell_full_gex/preprocess.py b/.cirro/singlecell_full_gex/preprocess.py index a4581c1..4897fab 100644 --- a/.cirro/singlecell_full_gex/preprocess.py +++ b/.cirro/singlecell_full_gex/preprocess.py @@ -105,6 +105,33 @@ def prepare_sample_sheet(ds): ds.logger.warning(f"sample_sheet is missing required column '{colname}'. Populating with NaN.") sample_sheet[colname] = np.nan + try: + # The PATIENT step (GIANA + GLIPH2) pools samples by params.patient_col. If that column + # is absent every sample becomes its own patient and the clustering is silently wrong - + # no error, just useless results - so say so loudly here. samplesheet_from_params()'s + # fallback frame in particular carries only sample/path. + # ds.params may be a plain dict or a params object depending on cirro version, and a + # warning must never be the thing that breaks preprocessing - so degrade to the default. + try: + patient_col = dict(ds.params).get("patient_col") or "patient_id" + except Exception: + patient_col = "patient_id" + if patient_col not in sample_sheet.columns: + ds.logger.warning( + f"sample_sheet has no '{patient_col}' column. GIANA and GLIPH2 pool samples per " + "patient, so each sample will be treated as its own patient and cross-sample " + "clustering within a patient will be lost. Add the column to the dataset " + "samplesheet, or set 'Patient column' to one that exists." + ) + elif sample_sheet[patient_col].isna().any(): + missing = sample_sheet.loc[sample_sheet[patient_col].isna(), "sample"].tolist() + ds.logger.warning(f"Samples with no {patient_col}: {missing}. These will not pool with any patient.") + else: + n_pat = sample_sheet[patient_col].nunique() + ds.logger.info(f"{len(sample_sheet)} samples across {n_pat} patients (by '{patient_col}')") + except Exception as e: + ds.logger.warning(f"patient-column check skipped: {e}") + sample_sheet.to_csv('sample_sheet.csv', index=None) ds.add_param('sample_sheet', 'sample_sheet.csv') ds.logger.info(sample_sheet.to_dict()) diff --git a/.cirro/singlecell_full_gex/process-form.json b/.cirro/singlecell_full_gex/process-form.json index 84a8411..a529be6 100644 --- a/.cirro/singlecell_full_gex/process-form.json +++ b/.cirro/singlecell_full_gex/process-form.json @@ -66,6 +66,30 @@ "description": "Render the master summary report.", "title": "Run master summary report", "type": "boolean" + }, + "patient_col": { + "type": "string", + "title": "Patient column", + "description": "Sample-sheet column holding the patient/subject ID. GIANA and GLIPH2 pool samples per patient using this; if it is wrong, every sample is treated as its own patient.", + "default": "patient_id" + }, + "run_tcri": { + "type": "boolean", + "title": "Run TCRi", + "description": "Immunogenicity scoring on the enriched Seurat object. Requires a GEX object. Adds roughly an hour of runtime.", + "default": true + }, + "giana_dedup_clonotypes": { + "type": "boolean", + "title": "Deduplicate clonotypes before GIANA", + "description": "Collapse to one row per clonotype before GIANA. Patient pooling stacks per-sample rows, so without this GIANA clusters cross-sample duplicates rather than similar sequences.", + "default": true + }, + "tcrdist_radius": { + "type": "integer", + "title": "TCRdist radius", + "description": "Distance threshold for TCRdist3 connected-component clustering.", + "default": 24 } } }, diff --git a/.cirro/singlecell_full_gex/process-input.json b/.cirro/singlecell_full_gex/process-input.json index 514359f..496afdf 100644 --- a/.cirro/singlecell_full_gex/process-input.json +++ b/.cirro/singlecell_full_gex/process-input.json @@ -15,5 +15,9 @@ "run_consensus": "$.params.dataset.paramJson.run_consensus", "run_repertoire": "$.params.dataset.paramJson.run_repertoire", "run_master_summary": "$.params.dataset.paramJson.run_master_summary", + "patient_col": "$.params.dataset.paramJson.patient_col", + "run_tcri": "$.params.dataset.paramJson.run_tcri", + "giana_dedup_clonotypes": "$.params.dataset.paramJson.giana_dedup_clonotypes", + "tcrdist_radius": "$.params.dataset.paramJson.tcrdist_radius", "outdir": "$.params.dataset.s3|/data/" } diff --git a/.cirro/singlecell_vdjonly/preprocess.py b/.cirro/singlecell_vdjonly/preprocess.py index a4581c1..4897fab 100644 --- a/.cirro/singlecell_vdjonly/preprocess.py +++ b/.cirro/singlecell_vdjonly/preprocess.py @@ -105,6 +105,33 @@ def prepare_sample_sheet(ds): ds.logger.warning(f"sample_sheet is missing required column '{colname}'. Populating with NaN.") sample_sheet[colname] = np.nan + try: + # The PATIENT step (GIANA + GLIPH2) pools samples by params.patient_col. If that column + # is absent every sample becomes its own patient and the clustering is silently wrong - + # no error, just useless results - so say so loudly here. samplesheet_from_params()'s + # fallback frame in particular carries only sample/path. + # ds.params may be a plain dict or a params object depending on cirro version, and a + # warning must never be the thing that breaks preprocessing - so degrade to the default. + try: + patient_col = dict(ds.params).get("patient_col") or "patient_id" + except Exception: + patient_col = "patient_id" + if patient_col not in sample_sheet.columns: + ds.logger.warning( + f"sample_sheet has no '{patient_col}' column. GIANA and GLIPH2 pool samples per " + "patient, so each sample will be treated as its own patient and cross-sample " + "clustering within a patient will be lost. Add the column to the dataset " + "samplesheet, or set 'Patient column' to one that exists." + ) + elif sample_sheet[patient_col].isna().any(): + missing = sample_sheet.loc[sample_sheet[patient_col].isna(), "sample"].tolist() + ds.logger.warning(f"Samples with no {patient_col}: {missing}. These will not pool with any patient.") + else: + n_pat = sample_sheet[patient_col].nunique() + ds.logger.info(f"{len(sample_sheet)} samples across {n_pat} patients (by '{patient_col}')") + except Exception as e: + ds.logger.warning(f"patient-column check skipped: {e}") + sample_sheet.to_csv('sample_sheet.csv', index=None) ds.add_param('sample_sheet', 'sample_sheet.csv') ds.logger.info(sample_sheet.to_dict()) diff --git a/.cirro/singlecell_vdjonly/process-form.json b/.cirro/singlecell_vdjonly/process-form.json index 2000dd9..1af374b 100644 --- a/.cirro/singlecell_vdjonly/process-form.json +++ b/.cirro/singlecell_vdjonly/process-form.json @@ -30,6 +30,24 @@ "drop", "hard_stop" ] + }, + "patient_col": { + "type": "string", + "title": "Patient column", + "description": "Sample-sheet column holding the patient/subject ID. GIANA and GLIPH2 pool samples per patient using this; if it is wrong, every sample is treated as its own patient.", + "default": "patient_id" + }, + "giana_dedup_clonotypes": { + "type": "boolean", + "title": "Deduplicate clonotypes before GIANA", + "description": "Collapse to one row per clonotype before GIANA. Patient pooling stacks per-sample rows, so without this GIANA clusters cross-sample duplicates rather than similar sequences.", + "default": true + }, + "tcrdist_radius": { + "type": "integer", + "title": "TCRdist radius", + "description": "Distance threshold for TCRdist3 connected-component clustering.", + "default": 24 } } }, diff --git a/.cirro/singlecell_vdjonly/process-input.json b/.cirro/singlecell_vdjonly/process-input.json index 1bf1131..b61fc4e 100644 --- a/.cirro/singlecell_vdjonly/process-input.json +++ b/.cirro/singlecell_vdjonly/process-input.json @@ -9,5 +9,8 @@ "pseudobulk_qc_min_clones": "$.params.dataset.paramJson.pseudobulk_qc_min_clones", "pseudobulk_qc_min_cells": "$.params.dataset.paramJson.pseudobulk_qc_min_cells", "pseudobulk_qc_mode": "$.params.dataset.paramJson.pseudobulk_qc_mode", + "patient_col": "$.params.dataset.paramJson.patient_col", + "giana_dedup_clonotypes": "$.params.dataset.paramJson.giana_dedup_clonotypes", + "tcrdist_radius": "$.params.dataset.paramJson.tcrdist_radius", "outdir": "$.params.dataset.s3|/data/" } From 4fe3630cd28a3dbadd61dc37ef2afee181ef2f74 Mon Sep 17 00:00:00 2001 From: Syed Shujaat Ali Zaidi Date: Tue, 1 Sep 2026 00:49:27 -0500 Subject: [PATCH 3/7] Add MERGE_VDJ_OBJECT: merged per-cell TCR object, pre- and post-QC The VDJ-only route emits no per-cell object. BULK_TO_EXPORT synthesizes a clonotype-level export whose cells are reconstructed rather than real barcodes (4,372 synthetic rows against 23,764 actual cells on the 8-sample test set). Builds a merged object from VDJ_QC's contig tables at both stages, keeping real barcodes: a Seurat object with 27 TCR metadata columns, the scRepertoire combineTCR object, and flat per-cell and per-sample tables. The pre/post pair makes QC's effect directly measurable. cellranger vdj has no genes, so the Seurat counts assay is a zero placeholder; this is recorded in seu@misc$provenance. Runs on both routes, toggleable via run_merge_vdj_object. Also narrows the run-artifact ignore rule from work*/ to work/ and work_*/ - the wildcard form also matched workflows/, the source directory. --- .gitignore | 4 +- bin/merge_vdj_object.R | 179 ++++++++++++++++++++++++++++ modules/bridges/merge_vdj_object.nf | 61 ++++++++++ nextflow.config | 3 + workflows/tcrtoolkit_sc.nf | 13 ++ 5 files changed, 259 insertions(+), 1 deletion(-) create mode 100755 bin/merge_vdj_object.R create mode 100644 modules/bridges/merge_vdj_object.nf diff --git a/.gitignore b/.gitignore index f1f5157..641377c 100644 --- a/.gitignore +++ b/.gitignore @@ -41,7 +41,9 @@ tmp ## Run artifacts — existing rules only match bare `work` and `results*`, ## so custom -w / --outdir names slip through. -work*/ +# NB: not `work*/` - that also matches workflows/, the source directory. +work/ +work_*/ test_results*/ test_v2port*/ logs/ diff --git a/bin/merge_vdj_object.R b/bin/merge_vdj_object.R new file mode 100755 index 0000000..fbfb043 --- /dev/null +++ b/bin/merge_vdj_object.R @@ -0,0 +1,179 @@ +#!/usr/bin/env Rscript +# +# merge_vdj_object.R +# +# Builds a merged per-cell TCR object from Cell Ranger VDJ contigs, for runs where no +# gene-expression matrix exists (the VDJ-only route, or any `cellranger vdj` output). +# +# IMPORTANT — what this can and cannot be: +# `cellranger vdj` produces receptor sequences per cell barcode and NO genes. A Seurat +# object is fundamentally a counts matrix (genes x cells) with metadata attached, so a +# genuine expression-bearing Seurat cannot be built from VDJ data alone. This script +# therefore emits two things: +# +# 1. _combineTCR.rds scRepertoire::combineTCR() output — the canonical +# VDJ-only merged object. This is the honest artifact and +# is exactly what you hand to combineExpression() later if +# a GEX object turns up. +# +# 2. _seurat.rds A Seurat object carrying one cell per barcode with all +# TCR fields in @meta.data, backed by a PLACEHOLDER assay +# of zeros. Provided only for downstream tools that demand +# the Seurat class. Do NOT normalize, run PCA/UMAP or +# cluster on it — there is no expression to analyse. +# +# 3. _cells.tsv The same per-cell table as a flat TSV. +# +# Usage: +# merge_vdj_object.R --contigs --prefix [--outdir .] + +suppressPackageStartupMessages({ + library(data.table) + library(dplyr) +}) + +# Base-R argument parsing: optparse is not in the single-cell container. +.args <- commandArgs(trailingOnly = TRUE) +.get <- function(flag, default = NA_character_) { + i <- match(flag, .args) + if (!is.na(i) && length(.args) > i) .args[i + 1] else default +} +opt <- list( + contigs = .get("--contigs"), + prefix = .get("--prefix", "vdj"), + outdir = .get("--outdir", "."), + sample_col = .get("--sample-col", "sample") +) +if (is.na(opt$contigs)) stop("--contigs is required") + +dir.create(opt$outdir, recursive = TRUE, showWarnings = FALSE) +msg <- function(...) message("[merge_vdj] ", ...) + +contigs <- as.data.frame(fread(opt$contigs, sep = "\t")) +msg(sprintf("read %s contigs, %d columns", format(nrow(contigs), big.mark = ","), ncol(contigs))) + +scol <- if (opt$sample_col %in% names(contigs)) opt$sample_col else "sample" +stopifnot(scol %in% names(contigs), "barcode" %in% names(contigs)) + +# Cell Ranger barcodes repeat across samples, so the cell key must include the sample. +contigs$cell_id <- paste0(contigs[[scol]], "_", sub("-1$", "", contigs$barcode)) + +# ── 1. scRepertoire combineTCR ──────────────────────────────────────────────── +# combineTCR() parses contig rows positionally and expects Cell Ranger's own column +# names, so feed it the raw frame split by sample rather than a renamed one. +tcr_obj <- NULL +if (requireNamespace("scRepertoire", quietly = TRUE)) { + contig_list <- split(contigs, contigs[[scol]]) + samples <- names(contig_list) + tcr_obj <- tryCatch({ + scRepertoire::combineTCR(contig_list, samples = samples, + removeNA = FALSE, removeMulti = FALSE, filterMulti = FALSE) + }, error = function(e) { + msg("WARN combineTCR() failed: ", conditionMessage(e)); NULL + }) + if (!is.null(tcr_obj)) { + saveRDS(tcr_obj, file.path(opt$outdir, paste0(opt$prefix, "_combineTCR.rds"))) + msg(sprintf("wrote %s_combineTCR.rds (%d samples)", opt$prefix, length(tcr_obj))) + } +} else { + msg("WARN scRepertoire unavailable — skipping combineTCR output") +} + +# ── 2. per-cell table ───────────────────────────────────────────────────────── +# One row per cell. Alpha/beta chains are collapsed separately so a cell with several +# contigs of the same chain keeps the highest-UMI one, matching Cell Ranger convention. +pick <- function(df, want) { + d <- df[df$chain %in% want, , drop = FALSE] + if (!nrow(d)) return(NULL) + if ("umis" %in% names(d)) d <- d[order(-as.numeric(d$umis)), , drop = FALSE] + d[1, , drop = FALSE] +} + +by_cell <- contigs %>% group_split(cell_id) +rows <- lapply(by_cell, function(d) { + a <- pick(d, c("TRA", "TRG")); b <- pick(d, c("TRB", "TRD")) + g <- function(x, col) if (!is.null(x) && col %in% names(x)) as.character(x[[col]]) else NA_character_ + data.frame( + cell_id = d$cell_id[1], + sample = d[[scol]][1], + barcode = sub("-1$", "", d$barcode[1]), + n_contigs = nrow(d), + n_chains = length(unique(d$chain)), + cdr3a = g(a, "cdr3"), cdr3b = g(b, "cdr3"), + cdr3a_nt = g(a, "cdr3_nt"), cdr3b_nt = g(b, "cdr3_nt"), + trav = g(a, "v_gene"), trbv = g(b, "v_gene"), + traj = g(a, "j_gene"), trbj = g(b, "j_gene"), + trac = g(a, "c_gene"), trbc = g(b, "c_gene"), + umis_a = if (!is.null(a) && "umis" %in% names(a)) as.numeric(a$umis) else NA_real_, + umis_b = if (!is.null(b) && "umis" %in% names(b)) as.numeric(b$umis) else NA_real_, + clonotype = if ("raw_clonotype_id" %in% names(d)) as.character(d$raw_clonotype_id[1]) else NA_character_, + multi_alpha = sum(d$chain %in% c("TRA", "TRG")) > 1, + multi_beta = sum(d$chain %in% c("TRB", "TRD")) > 1, + stringsAsFactors = FALSE + ) +}) +cells <- bind_rows(rows) + +# scRepertoire's CTaa convention, so this object joins cleanly to the GEX route's output. +cells$CTaa <- ifelse(is.na(cells$cdr3a) & is.na(cells$cdr3b), NA_character_, + paste0(ifelse(is.na(cells$cdr3a), "NA", paste0("A:", cells$cdr3a)), "|", + ifelse(is.na(cells$cdr3b), "NA", paste0("B:", cells$cdr3b)))) +cells$has_tcr <- !is.na(cells$cdr3a) | !is.na(cells$cdr3b) +cells$paired_tcr <- !is.na(cells$cdr3a) & !is.na(cells$cdr3b) + +# Clone size = cells sharing a CTaa within a sample. +cells <- cells %>% + group_by(sample, CTaa) %>% + mutate(clone_size = ifelse(is.na(CTaa), NA_integer_, dplyr::n())) %>% + ungroup() %>% + as.data.frame() + +fwrite(cells, file.path(opt$outdir, paste0(opt$prefix, "_cells.tsv")), sep = "\t") +msg(sprintf("wrote %s_cells.tsv (%s cells, %d samples, %s paired)", + opt$prefix, format(nrow(cells), big.mark = ","), + length(unique(cells$sample)), format(sum(cells$paired_tcr), big.mark = ","))) + +# ── 3. Seurat object with a placeholder assay ───────────────────────────────── +if (requireNamespace("Seurat", quietly = TRUE)) { + ok <- tryCatch({ + # A Seurat object requires a counts matrix. There is no expression data here, so + # this is a single all-zero feature row purely to satisfy the class contract. The + # information lives entirely in @meta.data. + # Two rows, not one: Seurat 5's Assay5 rejects a single-row layer with + # "Layers must be two-dimensional objects" because it drops to a vector. + m <- Matrix::Matrix(0, nrow = 2, ncol = nrow(cells), sparse = TRUE, + dimnames = list(c("PLACEHOLDER-no-GEX-1", "PLACEHOLDER-no-GEX-2"), + cells$cell_id)) + md <- cells; rownames(md) <- md$cell_id + seu <- Seurat::CreateSeuratObject(counts = m, meta.data = md, + project = opt$prefix, assay = "TCR", + min.cells = 0, min.features = 0) + seu@misc$provenance <- list( + source = normalizePath(opt$contigs), + stage = opt$prefix, + note = paste("VDJ-only: no gene expression. The counts assay is a", + "zero placeholder. Do not normalize, run PCA/UMAP or", + "cluster on this object."), + n_cells = nrow(cells), + n_samples = length(unique(cells$sample)), + created = as.character(Sys.time()) + ) + saveRDS(seu, file.path(opt$outdir, paste0(opt$prefix, "_seurat.rds"))) + msg(sprintf("wrote %s_seurat.rds (%s cells, placeholder assay)", + opt$prefix, format(ncol(seu), big.mark = ","))) + TRUE + }, error = function(e) { msg("WARN Seurat object failed: ", conditionMessage(e)); FALSE }) +} + +# ── 4. summary ──────────────────────────────────────────────────────────────── +summ <- cells %>% + group_by(sample) %>% + summarise(cells = dplyr::n(), + with_tcr = sum(has_tcr), paired = sum(paired_tcr), + pct_paired = round(100 * sum(paired_tcr) / dplyr::n(), 2), + unique_clonotypes = dplyr::n_distinct(CTaa[!is.na(CTaa)]), + multi_chain = sum(multi_alpha | multi_beta), + .groups = "drop") +fwrite(summ, file.path(opt$outdir, paste0(opt$prefix, "_summary.tsv")), sep = "\t") +print(as.data.frame(summ)) +msg("done") diff --git a/modules/bridges/merge_vdj_object.nf b/modules/bridges/merge_vdj_object.nf new file mode 100644 index 0000000..ca15b05 --- /dev/null +++ b/modules/bridges/merge_vdj_object.nf @@ -0,0 +1,61 @@ +/* + * MERGE_VDJ_OBJECT + * + * Builds a merged per-cell TCR object from the VDJ_QC contig tables, at both the pre- and + * post-QC stage. + * + * Why this exists: on the VDJ-only route nothing emits a per-cell object. The pipeline + * pools samples into pseudobulk, and BULK_TO_EXPORT then synthesizes a clonotype-level + * export whose "cells" are reconstructed rather than real barcodes (4,372 synthetic rows + * against 23,764 actual cells on the 8-sample test set). This process keeps real barcode + * identity, so cross-sample clonal questions can be asked at the cell level. + * + * It runs on both routes. With a GEX object TCELL_INTEGRATION already produces the + * authoritative Seurat, but the pre/post-QC pair here is still the cleanest way to see what + * VDJ QC actually removed. + * + * Note the emitted Seurat carries a placeholder counts assay — cellranger vdj has no genes. + * See bin/merge_vdj_object.R's header for what that does and does not permit. + */ +process MERGE_VDJ_OBJECT { + tag "${project_name}" + label 'process_medium' + container "${params.sc_container}" + + publishDir "${params.outdir}/bridge/merged_vdj_object", mode: 'copy', overwrite: true + + input: + path contigs_pre, stageAs: 'contigs_before_qc.tsv' + path contigs_post, stageAs: 'contigs_after_qc.tsv' + // Staged rather than called from ${projectDir}/bin so its contents join the task + // hash; otherwise editing the script leaves -resume reusing the previous output. + path merge_script, stageAs: 'merge_vdj_object.R' + val project_name + + output: + path "pre_qc_*", emit: pre_qc, optional: true + path "post_qc_*", emit: post_qc, optional: true + path "*_cells.tsv", emit: cells, optional: true + + script: + """ + # Each stage is independent: a failure on one must not cost the other, and a 0-byte + # NO_FILE placeholder (either table absent) is skipped rather than treated as data. + # stage:file pairs - the staged names are VDJ_QC's own (before/after), while the + # output prefix is pre/post, so they must be mapped explicitly rather than derived. + for pair in "pre:contigs_before_qc.tsv" "post:contigs_after_qc.tsv"; do + stage="\${pair%%:*}" + f="\${pair##*:}" + if [ ! -s "\$f" ]; then + echo "[merge_vdj] \$f missing or empty — skipping \${stage}-QC object" + continue + fi + Rscript merge_vdj_object.R \\ + --contigs "\$f" \\ + --prefix "\${stage}_qc" \\ + --outdir . \\ + --sample-col "${params.vdj_meta_sample_col ?: 'sample'}" \\ + || echo "[merge_vdj] WARN \${stage}-QC object failed; continuing" + done + """ +} diff --git a/nextflow.config b/nextflow.config index 62e3973..8378344 100644 --- a/nextflow.config +++ b/nextflow.config @@ -124,6 +124,9 @@ params { run_repertoire = true run_master_summary = true run_tcri = true + // Merged per-cell TCR object (pre- and post-QC) built straight from the VDJ contigs. + // Chiefly for the VDJ-only route, where nothing else emits a per-cell object. + run_merge_vdj_object = true // Optional pre-computed TCRi scores; when unset the report computes them itself. tcri_scores_file = null diff --git a/workflows/tcrtoolkit_sc.nf b/workflows/tcrtoolkit_sc.nf index da681e1..f509c7f 100644 --- a/workflows/tcrtoolkit_sc.nf +++ b/workflows/tcrtoolkit_sc.nf @@ -37,6 +37,7 @@ include { CLUSTER_TO_SC_SW } from '../subworkflows/bridges/cluster_to_sc.nf' include { SC_SAMPLE_STATS } from '../modules/bridges/sc_sample_stats.nf' include { BULK_TO_EXPORT } from '../modules/bridges/bulk_to_export.nf' include { CLUSTER_ROLLUP } from '../modules/bridges/cluster_rollup.nf' +include { MERGE_VDJ_OBJECT } from '../modules/bridges/merge_vdj_object.nf' // ── Shared bulk-TCR analysis engine (main/local — unchanged behavior) ───── include { PSEUDOBULK_QC_SW } from '../subworkflows/local/pseudobulk_qc.nf' @@ -96,6 +97,18 @@ workflow TCRTOOLKIT_SC { def vdj_qc_clone_rank_abundance_fig = vdj_qc_out.qc_figures.flatten().filter { f -> f.name == 'clone_rank_abundance.png' }.ifEmpty(nofile) def vdj_qc_multiple_chains_fig = vdj_qc_out.qc_figures.flatten().filter { f -> f.name == 'multiple_chains_by_sample.png' }.ifEmpty(nofile) + // ── Step 1b: merged per-cell TCR object, pre- and post-QC ───────────── + // Real barcodes, all samples pooled. On the VDJ-only route this is the only per-cell + // object produced - BULK_TO_EXPORT's export is clonotype-level with synthesized cells. + if (enabled(params.run_merge_vdj_object)) { + MERGE_VDJ_OBJECT( + pickFile(vdj_qc_out.qc_tables, 'contigs_before_qc.tsv', nofile), + pickFile(vdj_qc_out.qc_tables, 'contigs_after_qc.tsv', nofile), + channel.fromPath("${projectDir}/bin/merge_vdj_object.R", checkIfExists: true), + ch_project_name + ) + } + // ── Step 2: pseudobulk → clonotype table (route-specific source) ────── def sc_samplesheet = nofile if (vdjOnly) { From 2bea78f2eab0df73d219749a007ea63e617e2288 Mon Sep 17 00:00:00 2001 From: Dylan Tamayo <109362047+dltamayo@users.noreply.github.com> Date: Tue, 1 Sep 2026 15:31:52 -0400 Subject: [PATCH 4/7] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- modules/local/compare/giana.nf | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/modules/local/compare/giana.nf b/modules/local/compare/giana.nf index ec9f52d..e0c2636 100644 --- a/modules/local/compare/giana.nf +++ b/modules/local/compare/giana.nf @@ -13,9 +13,7 @@ process GIANA_CALC { path "${patient}_giana.txt", emit: 'giana_output' // path "giana_EncodingMatrix.txt" - script: - def dedup_clonotypes = (params.giana_dedup_clonotypes ?: false) ? 'True' : 'False' - + def dedup_clonotypes = (params.giana_dedup_clonotypes == true) ? 'True' : 'False' """ python3 - < Date: Tue, 1 Sep 2026 17:00:34 -0400 Subject: [PATCH 5/7] Restore script: label in GIANA_CALC, dropped by prior autofix commit The Copilot Autofix commit (2bea78f) that changed the dedup_clonotypes truthiness check from ?: to == true also dropped the script: label above it. Without it, Nextflow's parser stays in output-declaration mode and treats the def statement as an output entry, failing with "Invalid process output" at giana.nf:16 - this broke script compilation for every bulk test touching patient/compare (16 of 17 integration tests failed in PR #99's CI run). Co-Authored-By: Claude Sonnet 5 --- modules/local/compare/giana.nf | 1 + 1 file changed, 1 insertion(+) diff --git a/modules/local/compare/giana.nf b/modules/local/compare/giana.nf index e0c2636..09ede49 100644 --- a/modules/local/compare/giana.nf +++ b/modules/local/compare/giana.nf @@ -13,6 +13,7 @@ process GIANA_CALC { path "${patient}_giana.txt", emit: 'giana_output' // path "giana_EncodingMatrix.txt" + script: def dedup_clonotypes = (params.giana_dedup_clonotypes == true) ? 'True' : 'False' """ python3 - < Date: Tue, 1 Sep 2026 17:24:32 -0500 Subject: [PATCH 6/7] Master Summary: add GIANA / GLIPH2 / TCRdist3 cluster composition The three clustering modules write raw per-patient output and no figures - the bulk notebooks that plot them (template_giana.qmd, template_gliph.qmd) render only in bulk mode, so on the single-cell route GIANA and GLIPH2 computed correctly but nothing displayed them. CLUSTER_ROLLUP now emits five detail tables read straight from that raw output: per-cluster membership for each method, and V-gene usage among clustered clonotypes. Each cluster row carries n_members, n_unique_cdr3 and n_samples, so a cluster whose members are all one CDR3b is visibly duplicate detection rather than similarity. Master Summary gains four figures and a table: cluster size distribution across methods, top GLIPH2 motif clusters, the fraction of clusters spanning more than one sample, the GIANA cluster table, and V-gene usage. 13 -> 17 interactive figures. On the 8-sample test set: 89/89 GLIPH2 motifs and 7/7 GIANA clusters contain more than one distinct CDR3b, and roughly three quarters of clusters from both methods span multiple samples. --- bin/cluster_rollup.py | 72 ++++++++++- modules/bridges/cluster_rollup.nf | 5 + .../MASTER_SUMMARY/Master_Summary_Report.qmd | 117 ++++++++++++++++++ workflows/tcrtoolkit_sc.nf | 7 +- 4 files changed, 194 insertions(+), 7 deletions(-) diff --git a/bin/cluster_rollup.py b/bin/cluster_rollup.py index 76ad1d6..624c027 100755 --- a/bin/cluster_rollup.py +++ b/bin/cluster_rollup.py @@ -138,10 +138,11 @@ def tcrdist_clusters(matrix_paths, radius): except ImportError as exc: print(f"[cluster_rollup] WARN: scipy/numpy unavailable, skipping tcrdist clustering: {exc}", file=sys.stderr) - return 0, {}, 0 + return 0, {}, 0, None per_sample = {} total_clustered = 0 + detail_rows = [] for p in matrix_paths: if not usable(p): @@ -175,11 +176,61 @@ def tcrdist_clusters(matrix_paths, radius): n_comp, labels = connected_components(adj, directed=False) sizes = pd.Series(labels).value_counts() - n_real = int((sizes >= 2).sum()) - per_sample[sample] = n_real - total_clustered += int(sizes[sizes >= 2].sum()) + real = sizes[sizes >= 2] + per_sample[sample] = int(len(real)) + total_clustered += int(real.sum()) + for cid, n in real.items(): + detail_rows.append({"sample": sample, "cluster": int(cid), "n_members": int(n)}) - return sum(per_sample.values()), per_sample, total_clustered + detail = pd.DataFrame(detail_rows).sort_values("n_members", ascending=False) \ + if detail_rows else None + return sum(per_sample.values()), per_sample, total_clustered, detail + + +def write_detail(outdir, name, df, label): + """Detail tables feed the Master Summary's per-method figures.""" + if df is None or not len(df): + return + df.to_csv(os.path.join(outdir, name), sep="\t", index=False) + print(f"[cluster_rollup] wrote {name} ({len(df)} {label})") + + +def cluster_detail(df, cluster_col, cdr3, label): + """ + One row per cluster: size, distinct sequences, and how many samples it spans. + + Cluster identity is scoped by patient (the __group column). GIANA and GLIPH2 both + number/label clusters within a patient, so a bare id collides across patients. + """ + key = cluster_key(df, cluster_col) + if key is None or cdr3 is None: + return None + d = pd.DataFrame({ + "patient": df["__group"].astype(str), + "cluster": df[cluster_col].astype(str), + "_key": key, + "cdr3": cdr3, + "sample": df["sample"].astype(str) if "sample" in df.columns else "NA", + }).dropna(subset=["cdr3"]) + if not len(d): + return None + out = (d.groupby(["_key", "patient", "cluster"], as_index=False) + .agg(n_members=("cdr3", "size"), + n_unique_cdr3=("cdr3", "nunique"), + n_samples=("sample", "nunique"))) + # A cluster whose members are all the same sequence is duplicate detection, not + # similarity - surfacing it makes that visible rather than implied. + out["is_similarity"] = out["n_unique_cdr3"] > 1 + return out.drop(columns="_key").sort_values("n_members", ascending=False) + + +def vgene_usage(df, vcol, label): + if vcol not in df.columns: + return None + d = (df.assign(patient=df["__group"].astype(str), vgene=df[vcol].astype(str)) + .groupby(["patient", "vgene"], as_index=False) + .size().rename(columns={"size": "n_clonotypes"})) + return d.sort_values("n_clonotypes", ascending=False) def first_col(df, candidates): @@ -287,6 +338,10 @@ def main(): ]) clustered["GIANA"] = cluster_series(giana, cdr3, "cluster") counts["GIANA"] = int(n_clusters) + write_detail(args.outdir, "giana_cluster_detail.tsv", + cluster_detail(giana, "cluster", cdr3, "GIANA"), "clusters") + write_detail(args.outdir, "giana_vgene_usage.tsv", + vgene_usage(giana, "TRBV", "GIANA"), "patient/V-gene rows") # ── GLIPH2 ─────────────────────────────────────────────────────────────── if gliph is not None: @@ -305,6 +360,10 @@ def main(): ]) clustered["GLIPH2"] = cluster_series(gliph, cdr3, "tag") counts["GLIPH2"] = int(n_clusters) + write_detail(args.outdir, "gliph2_motif_detail.tsv", + cluster_detail(gliph, "tag", cdr3, "GLIPH2"), "motifs") + write_detail(args.outdir, "gliph2_vgene_usage.tsv", + vgene_usage(gliph, "TRBV", "GLIPH2"), "patient/V-gene rows") # ── TCRdist3 ───────────────────────────────────────────────────────────── # clone_df carries no cluster labels (those come from thresholding the distance @@ -322,9 +381,10 @@ def main(): if export is not None and "tcrdist_cluster" in export.columns: n_tcrdist_clusters = int(export["tcrdist_cluster"].nunique(dropna=True)) elif args.tcrdist_matrix: - n_tcrdist_clusters, per_sample, n_clustered_clones = tcrdist_clusters( + n_tcrdist_clusters, per_sample, n_clustered_clones, td_detail = tcrdist_clusters( args.tcrdist_matrix, args.tcrdist_radius ) + write_detail(args.outdir, "tcrdist_cluster_detail.tsv", td_detail, "clusters") rows = [ ("Samples analysed", tcrd["__group"].nunique()), diff --git a/modules/bridges/cluster_rollup.nf b/modules/bridges/cluster_rollup.nf index 4edd0be..4ed5768 100644 --- a/modules/bridges/cluster_rollup.nf +++ b/modules/bridges/cluster_rollup.nf @@ -38,6 +38,11 @@ process CLUSTER_ROLLUP { path "annotation_giana_summary.tsv", emit: annotation_giana, optional: true path "annotation_gliph2_summary.tsv", emit: annotation_gliph2, optional: true path "annotation_tcrdist3_summary.tsv", emit: annotation_tcrdist3, optional: true + path "giana_cluster_detail.tsv", emit: giana_detail, optional: true + path "gliph2_motif_detail.tsv", emit: gliph2_detail, optional: true + path "tcrdist_cluster_detail.tsv", emit: tcrdist_detail, optional: true + path "giana_vgene_usage.tsv", emit: giana_vgene, optional: true + path "gliph2_vgene_usage.tsv", emit: gliph2_vgene, optional: true script: // nullglob so an empty staging directory expands to nothing rather than a literal diff --git a/modules/scratch/MASTER_SUMMARY/Master_Summary_Report.qmd b/modules/scratch/MASTER_SUMMARY/Master_Summary_Report.qmd index c7ecbc5..e1289d3 100644 --- a/modules/scratch/MASTER_SUMMARY/Master_Summary_Report.qmd +++ b/modules/scratch/MASTER_SUMMARY/Master_Summary_Report.qmd @@ -286,6 +286,11 @@ method_counts <- Tfirst(c("consensus", "rollup"), "method_cluster_counts.tsv") anno_giana <- T("rollup", "annotation_giana_summary.tsv") anno_gliph2 <- T("rollup", "annotation_gliph2_summary.tsv") anno_tcrdist <- T("rollup", "annotation_tcrdist3_summary.tsv") +giana_detail <- T("rollup", "giana_cluster_detail.tsv") +gliph2_detail <- T("rollup", "gliph2_motif_detail.tsv") +tcrdist_detail <- T("rollup", "tcrdist_cluster_detail.tsv") +giana_vgene <- T("rollup", "giana_vgene_usage.tsv") +gliph2_vgene <- T("rollup", "gliph2_vgene_usage.tsv") # ── Repertoire ── rep_rollup <- T("repertoire", "repertoire_summary_rollup.tsv") @@ -781,6 +786,118 @@ duplicates as clusters — inflating cluster counts and producing no output at a patients with a single sample. ::: +## Cluster composition + +Every cluster from each method, sized by membership. `n_unique_cdr3` separates real +sequence similarity from duplicate detection: a cluster whose members are all the same +CDR3β is the same clonotype seen in several samples, not a similarity group. + +```{r} +#| label: cluster-sizes +#| results: asis +det <- list(GIANA = giana_detail, GLIPH2 = gliph2_detail, TCRdist3 = tcrdist_detail) +det <- det[vapply(det, have, logical(1))] +if (!length(det)) { + note("No cluster detail tables available.") +} else { + rows <- lapply(names(det), function(nm) { + d <- det[[nm]] + data.frame(method = nm, + size = as.numeric(d$n_members), + stringsAsFactors = FALSE) + }) + d <- bind_rows(rows) %>% + count(method, size, name = "n_clusters") %>% + transmute(cat = size, val = n_clusters, grp = method, + hover = paste0(method, "
", n_clusters, " cluster(s) of size ", size)) + save_table(d, "master_cluster_size_distribution.tsv") + gbar(d, "Cluster size distribution by method", "Members per cluster", "Clusters") +} +``` + +```{r} +#| label: gliph2-top-motifs +if (!have(gliph2_detail)) { + note("No GLIPH2 motif detail available.") +} else { + save_table(gliph2_detail, "master_gliph2_motif_detail.tsv") + d <- gliph2_detail %>% + slice_max(n_members, n = params$top_n_clusters, with_ties = FALSE) %>% + transmute(cat = paste0(cluster, " (", patient, ")"), + val = n_members, + hover = paste0("motif ", cluster, "
patient ", patient, + "
", n_members, " members, ", n_unique_cdr3, + " distinct CDR3b
", n_samples, " sample(s)")) + hbar(d, paste0("Top ", params$top_n_clusters, " GLIPH2 motif clusters"), + "Members", colour = "#d95f02", height = 560) +} +``` + +```{r} +#| label: cross-sample-reach +#| results: asis +rows <- list() +for (nm in names(det)) { + d <- det[[nm]] + if (!("n_samples" %in% names(d))) next + rows[[nm]] <- data.frame( + method = nm, + frac = mean(as.numeric(d$n_samples) > 1, na.rm = TRUE), + n = nrow(d), stringsAsFactors = FALSE) +} +if (!length(rows)) { + note("No per-cluster sample counts available.") +} else { + d <- bind_rows(rows) + save_table(d, "master_cluster_cross_sample.tsv") + dd <- d %>% transmute(cat = method, val = frac, + hover = paste0(method, "
", + percent(frac, accuracy = 0.1), + " of ", n, " clusters span >1 sample")) + hbar(dd, "Clusters spanning more than one sample", "Fraction of clusters", + height = 380, tickformat = ".0%") +} +``` + +A cluster spanning several samples is the interesting case — the same or similar receptor +recovered independently, rather than confined to one library. + +```{r} +#| label: giana-cluster-table +if (!have(giana_detail)) { + note("No GIANA cluster detail available.") +} else { + save_table(giana_detail, "master_giana_cluster_detail.tsv") + itable(giana_detail, pageSize = 10) +} +``` + +## V-gene usage within clusters + +Germline V-gene composition of the clustered clonotypes, per patient. A motif family +dominated by one V gene is a stronger structural signal than one drawn from many. + +```{r} +#| label: cluster-vgene +#| results: asis +vg <- list(GIANA = giana_vgene, GLIPH2 = gliph2_vgene) +vg <- vg[vapply(vg, have, logical(1))] +if (!length(vg)) { + note("No V-gene usage tables available for the clustering methods.") +} else { + d <- bind_rows(lapply(names(vg), function(nm) + data.frame(method = nm, vg[[nm]], stringsAsFactors = FALSE))) %>% + group_by(method) %>% + slice_max(n_clonotypes, n = params$top_n_genes, with_ties = FALSE) %>% + ungroup() + save_table(d, "master_cluster_vgene_usage.tsv") + dd <- d %>% transmute(cat = vgene, val = n_clonotypes, grp = method, + hover = paste0(method, "
", vgene, "
patient ", patient, + "
", n_clonotypes, " clonotypes")) + gbar(dd, "V-gene usage among clustered clonotypes", "", "Clonotypes", height = 520) +} +``` + ## Per-annotation clustering coverage Fraction of cells clustered by each method, broken down by cell-type annotation. Requires diff --git a/workflows/tcrtoolkit_sc.nf b/workflows/tcrtoolkit_sc.nf index f509c7f..a314766 100644 --- a/workflows/tcrtoolkit_sc.nf +++ b/workflows/tcrtoolkit_sc.nf @@ -251,7 +251,12 @@ workflow TCRTOOLKIT_SC { CLUSTER_ROLLUP.out.method_cluster_counts, CLUSTER_ROLLUP.out.annotation_giana, CLUSTER_ROLLUP.out.annotation_gliph2, - CLUSTER_ROLLUP.out.annotation_tcrdist3) + CLUSTER_ROLLUP.out.annotation_tcrdist3, + CLUSTER_ROLLUP.out.giana_detail, + CLUSTER_ROLLUP.out.gliph2_detail, + CLUSTER_ROLLUP.out.tcrdist_detail, + CLUSTER_ROLLUP.out.giana_vgene, + CLUSTER_ROLLUP.out.gliph2_vgene) .collect().ifEmpty([]) def pseudobulk_tables = channel.empty() From 637af7e8b21aa5c34f86c5b5c2e0c514a5cf9152 Mon Sep 17 00:00:00 2001 From: Syed Shujaat Ali Zaidi Date: Wed, 2 Sep 2026 15:58:44 -0500 Subject: [PATCH 7/7] Flag and report receptors that never matched a GEX cell The full single-cell route pseudobulks from the post-merge export, so a receptor whose barcode did not match a GEX cell is absent from every downstream result - not only the cell-level steps but GIANA, GLIPH2, TCRdist3 and the repertoire statistics. At the 97% match rate of the test set this is immaterial; at a low rate it would silently remove a large share of the data, and nothing reported it. MERGE_VDJ_OBJECT now takes the post-merge export and flags every conserved cell with gex_matched, so the excluded receptors stay identifiable rather than merely present. It runs after TCELL_INTEGRATION for this reason. The Master Summary gains a Receptor retention section reporting the overall match rate with a per-sample breakdown, escalating from a note to a warning below 75%. VDJ-only is unaffected: it performs no GEX join, analyses every receptor, and the flag is NA there. --- bin/merge_vdj_object.R | 30 +++++++- modules/bridges/merge_vdj_object.nf | 4 + .../MASTER_SUMMARY/Master_Summary_Report.qmd | 76 +++++++++++++++++++ modules/scratch/MASTER_SUMMARY/main.nf | 1 + subworkflows/scratch/master_summary.nf | 2 + workflows/tcrtoolkit_sc.nf | 30 +++++--- 6 files changed, 130 insertions(+), 13 deletions(-) diff --git a/bin/merge_vdj_object.R b/bin/merge_vdj_object.R index fbfb043..29b4ace 100755 --- a/bin/merge_vdj_object.R +++ b/bin/merge_vdj_object.R @@ -42,7 +42,11 @@ opt <- list( contigs = .get("--contigs"), prefix = .get("--prefix", "vdj"), outdir = .get("--outdir", "."), - sample_col = .get("--sample-col", "sample") + sample_col = .get("--sample-col", "sample"), + # Optional post-merge GEX export. When supplied, every cell is flagged with whether it + # found a GEX partner - so the receptors dropped from the Full-SC analysis stay + # identifiable in this conserved object rather than just being absent downstream. + gex_export = .get("--gex-export") ) if (is.na(opt$contigs)) stop("--contigs is required") @@ -128,6 +132,28 @@ cells <- cells %>% ungroup() %>% as.data.frame() +# ── GEX match flag ──────────────────────────────────────────────────────────── +# The Full-SC route feeds its pseudobulk from the POST-merge export, so a receptor that +# never matched a GEX barcode is absent from every downstream result. This object keeps it; +# the flag records which ones those are. +cells$gex_matched <- NA +if (!is.na(opt$gex_export) && file.exists(opt$gex_export) && file.size(opt$gex_export) > 0) { + gx <- tryCatch(fread(opt$gex_export, sep = "\t", select = c("cell_id")), + error = function(e) NULL) + if (!is.null(gx) && nrow(gx)) { + # GEX cell_ids carry Cell Ranger's "-1" suffix; this object strips it. Normalise + # both sides before comparing or nothing matches. + gex_ids <- unique(sub("-1$", "", as.character(gx$cell_id))) + cells$gex_matched <- cells$cell_id %in% gex_ids + n_hit <- sum(cells$gex_matched) + msg(sprintf("GEX match: %s of %s cells (%.1f%%) found a GEX partner; %s conserved here only", + format(n_hit, big.mark = ","), format(nrow(cells), big.mark = ","), + 100 * n_hit / nrow(cells), format(nrow(cells) - n_hit, big.mark = ","))) + } else { + msg("WARN could not read --gex-export; gex_matched left NA") + } +} + fwrite(cells, file.path(opt$outdir, paste0(opt$prefix, "_cells.tsv")), sep = "\t") msg(sprintf("wrote %s_cells.tsv (%s cells, %d samples, %s paired)", opt$prefix, format(nrow(cells), big.mark = ","), @@ -173,6 +199,8 @@ summ <- cells %>% pct_paired = round(100 * sum(paired_tcr) / dplyr::n(), 2), unique_clonotypes = dplyr::n_distinct(CTaa[!is.na(CTaa)]), multi_chain = sum(multi_alpha | multi_beta), + gex_matched = if (all(is.na(gex_matched))) NA_integer_ else sum(gex_matched, na.rm = TRUE), + gex_only_here = if (all(is.na(gex_matched))) NA_integer_ else sum(!gex_matched, na.rm = TRUE), .groups = "drop") fwrite(summ, file.path(opt$outdir, paste0(opt$prefix, "_summary.tsv")), sep = "\t") print(as.data.frame(summ)) diff --git a/modules/bridges/merge_vdj_object.nf b/modules/bridges/merge_vdj_object.nf index ca15b05..273cee6 100644 --- a/modules/bridges/merge_vdj_object.nf +++ b/modules/bridges/merge_vdj_object.nf @@ -30,6 +30,9 @@ process MERGE_VDJ_OBJECT { // Staged rather than called from ${projectDir}/bin so its contents join the task // hash; otherwise editing the script leaves -resume reusing the previous output. path merge_script, stageAs: 'merge_vdj_object.R' + // Post-merge GEX export when the full-SC route ran; NO_FILE otherwise. Used only to + // flag which conserved receptors also reached the GEX analysis. + path gex_export, stageAs: 'in_gex_export.tsv' val project_name output: @@ -55,6 +58,7 @@ process MERGE_VDJ_OBJECT { --prefix "\${stage}_qc" \\ --outdir . \\ --sample-col "${params.vdj_meta_sample_col ?: 'sample'}" \\ + --gex-export in_gex_export.tsv \\ || echo "[merge_vdj] WARN \${stage}-QC object failed; continuing" done """ diff --git a/modules/scratch/MASTER_SUMMARY/Master_Summary_Report.qmd b/modules/scratch/MASTER_SUMMARY/Master_Summary_Report.qmd index e1289d3..7f199b1 100644 --- a/modules/scratch/MASTER_SUMMARY/Master_Summary_Report.qmd +++ b/modules/scratch/MASTER_SUMMARY/Master_Summary_Report.qmd @@ -306,6 +306,8 @@ pheno_flux <- T("repertoire", "phenotypic_flux.tsv") tcell_roll <- T("tcell", "tcell_summary_rollup.tsv") tcell_anno <- T("tcell", "tcell_per_annotation_summary.tsv") tcell_topclone <- T("tcell", "tcell_top_clones.tsv") +tcell_overlap <- T("tcell", "barcode_overlap_by_sample.tsv") +merged_vdj <- T("mergedvdj", "post_qc_summary.tsv") conga_roll <- T("conga", "conga_summary_rollup.tsv") conga_clusters <- T("conga", "conga_cluster_summary.tsv") anno_conga <- T("conga", "annotation_conga_summary.tsv") @@ -379,6 +381,80 @@ save_table(ov, "master_overview.tsv") itable(ov, pageSize = 10) ``` +## Receptor retention + +The full single-cell route pseudobulks from the **post-merge** export, so a receptor whose +barcode never matched a GEX cell is absent from every downstream result — not only the +cell-level steps, but GIANA, GLIPH2, TCRdist3 and the repertoire statistics too. Every +receptor is nevertheless conserved in `bridge/merged_vdj_object/`, flagged with whether it +reached the GEX analysis. + +At a high match rate this is immaterial. At a low one it silently removes a large share of +the data, so the rate is reported here rather than left in a module table. + +```{r} +#| label: receptor-retention +#| results: asis +if (!have(tcell_overlap)) { + note(paste("No barcode overlap table — either the VDJ-only route (no GEX join is", + "performed, so every receptor is analysed) or T-cell integration did not run.")) +} else { + d <- tcell_overlap + nb <- intersect(c("tcr_barcodes","n_tcr"), names(d)) + no <- intersect(c("overlap","matched","matched_keys"), names(d)) + if (!length(nb) || !length(no)) { + note("Overlap table present but lacks the expected barcode/overlap columns.") + } else { + tot <- sum(as.numeric(d[[nb[1]]]), na.rm = TRUE) + hit <- sum(as.numeric(d[[no[1]]]), na.rm = TRUE) + rate <- if (tot > 0) hit / tot else NA_real_ + save_table(d, "master_receptor_retention.tsv") + + if (!is.na(rate) && rate < 0.75) { + cat(sprintf(paste0("\n::: {.callout-warning appearance=\"simple\"}\n", + "**%s of %s receptors (%.1f%%) matched a GEX cell.** The remaining %s are ", + "excluded from every downstream result on this route. They are preserved in ", + "`bridge/merged_vdj_object/` with `gex_matched = FALSE`. Below ~75%%, ", + "consider whether barcode conventions are reconciling correctly — the ", + "per-sample rates below will show if it is one sample or all of them.\n:::\n\n"), + format(hit, big.mark=","), format(tot, big.mark=","), 100*rate, + format(tot-hit, big.mark=","))) + } else { + cat(sprintf(paste0("\n::: {.callout-note appearance=\"simple\"}\n", + "**%s of %s receptors (%.1f%%) matched a GEX cell.** The %s unmatched are ", + "excluded from downstream results but preserved in ", + "`bridge/merged_vdj_object/` with `gex_matched = FALSE`.\n:::\n\n"), + format(hit, big.mark=","), format(tot, big.mark=","), 100*rate, + format(tot-hit, big.mark=","))) + } + } +} +``` + +```{r} +#| label: retention-by-sample +if (!have(tcell_overlap)) { + invisible(NULL) +} else { + d <- tcell_overlap + sc <- intersect(c("sample","Sample"), names(d)) + rc <- intersect(c("overlap_rate","rate"), names(d)) + nb <- intersect(c("tcr_barcodes","n_tcr"), names(d)) + no <- intersect(c("overlap","matched","matched_keys"), names(d)) + if (length(sc) && length(nb) && length(no)) { + r <- if (length(rc)) as.numeric(d[[rc[1]]]) else + as.numeric(d[[no[1]]]) / as.numeric(d[[nb[1]]]) + dd <- data.frame(cat = as.character(d[[sc[1]]]), val = r, + hover = paste0(d[[sc[1]]], "
", + format(as.numeric(d[[no[1]]]), big.mark=","), " of ", + format(as.numeric(d[[nb[1]]]), big.mark=","), + " matched
", percent(r, accuracy = 0.1))) + hbar(dd, "Receptors matching a GEX cell, by sample", "Match rate", + height = 400, tickformat = ".0%") + } +} +``` + ## Module coverage Which analysis modules contributed output to this report. `FALSE` for the GEX-dependent diff --git a/modules/scratch/MASTER_SUMMARY/main.nf b/modules/scratch/MASTER_SUMMARY/main.nf index 7309e3b..4b70c2b 100644 --- a/modules/scratch/MASTER_SUMMARY/main.nf +++ b/modules/scratch/MASTER_SUMMARY/main.nf @@ -26,6 +26,7 @@ process MASTER_SUMMARY { path conga_tables, stageAs: 'intables/conga/*' path consensus_tables, stageAs: 'intables/consensus/*' path tcri_tables, stageAs: 'intables/tcri/*' + path mergedvdj_tables, stageAs: 'intables/mergedvdj/*' path qmd val barrier_done diff --git a/subworkflows/scratch/master_summary.nf b/subworkflows/scratch/master_summary.nf index c54bde0..c2b1f0c 100644 --- a/subworkflows/scratch/master_summary.nf +++ b/subworkflows/scratch/master_summary.nf @@ -25,6 +25,7 @@ workflow MASTER_SUMMARY_SW { conga_tables consensus_tables tcri_tables + mergedvdj_tables barrier_done project_name @@ -48,6 +49,7 @@ workflow MASTER_SUMMARY_SW { conga_tables, consensus_tables, tcri_tables, + mergedvdj_tables, ch_notebook, barrier_done, project_name diff --git a/workflows/tcrtoolkit_sc.nf b/workflows/tcrtoolkit_sc.nf index a314766..b6e256a 100644 --- a/workflows/tcrtoolkit_sc.nf +++ b/workflows/tcrtoolkit_sc.nf @@ -97,18 +97,6 @@ workflow TCRTOOLKIT_SC { def vdj_qc_clone_rank_abundance_fig = vdj_qc_out.qc_figures.flatten().filter { f -> f.name == 'clone_rank_abundance.png' }.ifEmpty(nofile) def vdj_qc_multiple_chains_fig = vdj_qc_out.qc_figures.flatten().filter { f -> f.name == 'multiple_chains_by_sample.png' }.ifEmpty(nofile) - // ── Step 1b: merged per-cell TCR object, pre- and post-QC ───────────── - // Real barcodes, all samples pooled. On the VDJ-only route this is the only per-cell - // object produced - BULK_TO_EXPORT's export is clonotype-level with synthesized cells. - if (enabled(params.run_merge_vdj_object)) { - MERGE_VDJ_OBJECT( - pickFile(vdj_qc_out.qc_tables, 'contigs_before_qc.tsv', nofile), - pickFile(vdj_qc_out.qc_tables, 'contigs_after_qc.tsv', nofile), - channel.fromPath("${projectDir}/bin/merge_vdj_object.R", checkIfExists: true), - ch_project_name - ) - } - // ── Step 2: pseudobulk → clonotype table (route-specific source) ────── def sc_samplesheet = nofile if (vdjOnly) { @@ -123,6 +111,23 @@ workflow TCRTOOLKIT_SC { pseudobulk_map = SC_TO_CDR3_SW.out.sample_map } + // ── Step 2b: merged per-cell TCR object, pre- and post-QC ───────────── + // Reads the contig tables directly, so it conserves EVERY receptor - including those + // that never matched a GEX barcode and are therefore absent from the full-SC analysis + // (which pseudobulks from the post-merge export). Runs after TCELL_INTEGRATION so it + // can flag which receptors did match. + mergedvdj_tables = channel.empty() + if (enabled(params.run_merge_vdj_object)) { + MERGE_VDJ_OBJECT( + pickFile(vdj_qc_out.qc_tables, 'contigs_before_qc.tsv', nofile), + pickFile(vdj_qc_out.qc_tables, 'contigs_after_qc.tsv', nofile), + channel.fromPath("${projectDir}/bin/merge_vdj_object.R", checkIfExists: true), + vdjOnly ? channel.fromPath("${projectDir}/assets/NO_FILE") : tcell_out.export_cells, + ch_project_name + ) + mergedvdj_tables = MERGE_VDJ_OBJECT.out.pre_qc.mix(MERGE_VDJ_OBJECT.out.post_qc) + } + // ── Step 3: tcrtoolkit pseudobulk QC gate (both routes) ─────────────── PSEUDOBULK_QC_SW( pseudobulk_map ) @@ -280,6 +285,7 @@ workflow TCRTOOLKIT_SC { conga_tables.flatten().collect().ifEmpty([]), consensus_tables.flatten().collect().ifEmpty([]), tcri_tables.flatten().collect().ifEmpty([]), + mergedvdj_tables.flatten().collect().ifEmpty([]), master_barrier, ch_project_name )