From b77455cf48e5857f9c6a9b577abfa2630a4ab05c Mon Sep 17 00:00:00 2001 From: ch99l Date: Thu, 30 Jul 2026 10:23:20 +0800 Subject: [PATCH 01/24] visium hd implementation --- README.md | 2 +- bin/bambu_discovery.R | 31 ++++ bin/save_counts.R | 20 +-- bin/visium_hd_aggregate_resolution.R | 45 +++++ bin/visium_hd_cluster.R | 59 ++++++ bin/visium_hd_convert_barcode_mappings.py | 36 ++++ bin/visium_hd_convert_tissue_positions.py | 25 +++ bin/visium_hd_spot_bin_mappings.py | 30 ++++ conf/containers.config | 5 +- containers/r/Dockerfile | 4 + lib/Validation.groovy | 80 ++++++++- main.nf | 170 +++++++++++++----- .../clustered_quantification.nf} | 29 ++- .../{ => shared}/construct_read_class.nf | 0 .../bambu/{ => shared}/prepare_annotation.nf | 0 .../standard/single_cell_quantification.nf | 41 +++++ .../{ => standard}/transcript_discovery.nf | 27 ++- modules/bambu/visium_hd/aggregate_bins.nf | 42 +++++ .../visium_hd/spot_level_quantification.nf | 63 +++++++ .../bambu/visium_hd/transcript_discovery.nf | 54 ++++++ .../convert_barcode_mappings.nf | 26 +++ .../convert_tissue_positions.nf | 27 +++ .../filter_barcoded_bam.nf | 27 +++ .../spot_bin_mappings.nf | 27 +++ .../multi_sample_clustering.nf} | 2 +- .../single_sample_clustering.nf} | 2 +- modules/seurat/visium_hd/clustering.nf | 58 ++++++ nextflow.config | 40 +++-- .../{clustering.nf => clustering_standard.nf} | 4 +- subworkflows/prepare_input_visium_hd.nf | 53 ++++++ 30 files changed, 919 insertions(+), 110 deletions(-) create mode 100755 bin/bambu_discovery.R create mode 100755 bin/visium_hd_aggregate_resolution.R create mode 100755 bin/visium_hd_cluster.R create mode 100755 bin/visium_hd_convert_barcode_mappings.py create mode 100755 bin/visium_hd_convert_tissue_positions.py create mode 100755 bin/visium_hd_spot_bin_mappings.py rename modules/bambu/{EM_quant.nf => shared/clustered_quantification.nf} (55%) rename modules/bambu/{ => shared}/construct_read_class.nf (100%) rename modules/bambu/{ => shared}/prepare_annotation.nf (100%) create mode 100644 modules/bambu/standard/single_cell_quantification.nf rename modules/bambu/{ => standard}/transcript_discovery.nf (71%) create mode 100644 modules/bambu/visium_hd/aggregate_bins.nf create mode 100644 modules/bambu/visium_hd/spot_level_quantification.nf create mode 100644 modules/bambu/visium_hd/transcript_discovery.nf create mode 100644 modules/prepare_input_visium_hd/convert_barcode_mappings.nf create mode 100644 modules/prepare_input_visium_hd/convert_tissue_positions.nf create mode 100644 modules/prepare_input_visium_hd/filter_barcoded_bam.nf create mode 100644 modules/prepare_input_visium_hd/spot_bin_mappings.nf rename modules/seurat/{multi_sample.nf => standard/multi_sample_clustering.nf} (95%) rename modules/seurat/{single_sample.nf => standard/single_sample_clustering.nf} (94%) create mode 100644 modules/seurat/visium_hd/clustering.nf rename subworkflows/{clustering.nf => clustering_standard.nf} (71%) create mode 100644 subworkflows/prepare_input_visium_hd.nf diff --git a/README.md b/README.md index 25e7d0a..8b0e4b3 100644 --- a/README.md +++ b/README.md @@ -145,7 +145,7 @@ To configure the executor and container, pass profile types via the `-profile` a - "no_quant": Transcript quantification is not performed - "EM": Performs transcript quantification for each cell/spatial coordinate - "EM_clusters": Performs gene expression-based cell clustering using [Seurat](https://satijalab.org/seurat/), followed by transcript quantification at the cluster level -- `--resolution` [float, default: 0.8]: Seurat clustering resolution +- `--cluster_resolution` [float, default: 0.8]: Seurat clustering resolution ### **Output** All outputs from the pipeline are written to the directory specified by the `--output_dir` parameter. The pipeline produces per-sample alignment files and the combined transcript discovery and quantification results. diff --git a/bin/bambu_discovery.R b/bin/bambu_discovery.R new file mode 100755 index 0000000..37cdae8 --- /dev/null +++ b/bin/bambu_discovery.R @@ -0,0 +1,31 @@ +# Shared discovery sequence for the standard and Visium HD transcript-discovery +# processes. Assumes bambu is loaded by the caller. +# bambuDiscovery() arguments: +# reads named character vector of readClass .rds paths (names = sample +# names), or a single path for a one-sample run. +# annotation the bambu transcript annotation object (a GRangesList). +# genome path to the reference genome FASTA. +# ncore integer, number of cores passed to bambu.singlecell. +# ndr numeric, the novel discovery rate (NDR) threshold for discovery. +# sampleData optional path(s) to spatial / bin metadata joined into colData by +# the quantData stage; NULL (default) attaches none. +# +# Returns a named list with extendedAnno, quantData and seDiscovery (the +# transcript-level unique-counts SE). No file I/O. +bambuDiscovery <- function(reads, annotation, genome, ncore, ndr, sampleData = NULL) { + # Transcript discovery + extendedAnno <- bambu.singlecell(reads = reads, output = "extendedAnnotations", + annotations = annotation, genome = genome, ncore = ncore, + verbose = FALSE, NDR = ndr) + + # Read to transcript assignment + quantData <- bambu.singlecell(reads = reads, output = "quantData", + annotations = extendedAnno, genome = genome, ncore = ncore, + verbose = FALSE, sampleData = sampleData) + + # Quantification without EM + seDiscovery <- bambu.singlecell(reads = quantData, output = "uniqueCounts", + annotations = extendedAnno) + + list(extendedAnno = extendedAnno, quantData = quantData, seDiscovery = seDiscovery) +} diff --git a/bin/save_counts.R b/bin/save_counts.R index eb01f8a..72b12a7 100755 --- a/bin/save_counts.R +++ b/bin/save_counts.R @@ -1,17 +1,17 @@ -# save_counts() - write a SummarizedExperiment as a count directory: one +# saveCounts() - write a SummarizedExperiment as a count directory: one # MTX file per assay, plus a single shared barcodes.tsv.gz / features.tsv.gz, -# and the SummarisedExperiment object. +# and the SummarisedExperiment object. # Assumes DropletUtils and bambu have been loaded -save_counts <- function(se, dir, gene.type = "Gene Expression") { - write10xCounts(dir, assays(se)$counts, version = "3", gene.type = gene.type) +saveCounts <- function(se, dir, geneType = "Gene Expression") { + write10xCounts(dir, assays(se)$counts, version = "3", gene.type = geneType) file.rename(file.path(dir, "matrix.mtx.gz"), file.path(dir, "counts.mtx.gz")) - extra_assays <- setdiff(assayNames(se), "counts") - for (assay_name in extra_assays) { - mat <- as(assays(se)[[assay_name]], "CsparseMatrix") - mtx_path <- file.path(dir, paste0(assay_name, ".mtx")) - Matrix::writeMM(mat, mtx_path) - R.utils::gzip(mtx_path, overwrite = TRUE) + extraAssays <- setdiff(assayNames(se), "counts") + for (assayName in extraAssays) { + mat <- as(assays(se)[[assayName]], "CsparseMatrix") + mtxPath <- file.path(dir, paste0(assayName, ".mtx")) + Matrix::writeMM(mat, mtxPath) + R.utils::gzip(mtxPath, overwrite = TRUE) } saveRDS(se, file.path(dir, paste0("se_", dir, ".rds"))) } diff --git a/bin/visium_hd_aggregate_resolution.R b/bin/visium_hd_aggregate_resolution.R new file mode 100755 index 0000000..5dc5c7a --- /dev/null +++ b/bin/visium_hd_aggregate_resolution.R @@ -0,0 +1,45 @@ +# Assumes SummarizedExperiment, Matrix and dplyr are loaded by the caller. +# A spot is a 2um barcode, i.e. one column of se; a bin is the coarser square +# (e.g. 8um) that several spots fall into. +# Build a bin-resolution SummarizedExperiment from a 2um-resolution SE by summing +# counts across the spots that share each bin. spotMappings gives the spot -> bin +# assignment for one resolution; tissuePositions supplies the new colData, +# alongside the id and sampleName carried over from the 2um SE. +# rowRanges/rowData and metadata are carried over. +aggregateResolution <- function(se, tissuePositions, spotMappings) { + # Every column of se is one spot, so look up the bin each spot falls in. The + # prefixed forms are used so the aggregated columns are named like bambu ids. + matchIdx <- match(colData(se)$id, spotMappings$sample_barcode) + spotToBinMap <- factor(spotMappings$sample_bin[matchIdx]) + + # The bins become the columns of the aggregated matrix, in level order + binIds <- levels(spotToBinMap) # bambu ids, e.g. _s_008um__-1 + spotsPerBin <- t(fac2sparse(spotToBinMap)) # spots x bins, 1 where a spot falls in a bin + + # Sum each bin's spots: (features x spots) %*% (spots x bins) = features x bins + spotCounts <- assays(se)$counts + binCounts <- as(spotCounts %*% spotsPerBin, "CsparseMatrix") + + # build new colData for the aggregated SE object + barcodePerBin <- distinct(spotMappings, sample_bin, barcode = bin) + positionPerBin <- data.frame(sample_bin = binIds) %>% + left_join(barcodePerBin, by = "sample_bin") %>% + left_join(tissuePositions, by = "barcode") %>% + select(-sample_bin) + + sampleName <- unique(colData(se)$sampleName) + binColData <- DataFrame(id = binIds, sampleName = sampleName, positionPerBin, row.names = binIds) + + binSe <- SummarizedExperiment(assays = SimpleList(counts = binCounts), + rowRanges = rowRanges(se), + colData = binColData) + + # Carry the bambu metadata matrices over, summed the same way + incompatiblePerSpot <- metadata(se)$incompatibleCounts + nonuniquePerSpot <- metadata(se)$nonuniqueCounts + + metadata(binSe)$incompatibleCounts <- as(incompatiblePerSpot %*% spotsPerBin, "CsparseMatrix") + metadata(binSe)$nonuniqueCounts <- as(nonuniquePerSpot %*% spotsPerBin, "CsparseMatrix") + metadata(binSe)$seType <- metadata(se)$seType + binSe +} diff --git a/bin/visium_hd_cluster.R b/bin/visium_hd_cluster.R new file mode 100755 index 0000000..ac9e6f0 --- /dev/null +++ b/bin/visium_hd_cluster.R @@ -0,0 +1,59 @@ +# Seurat clustering of a single Visium HD bin resolution, adapted from +# https://satijalab.org/seurat/articles/visiumhd_analysis_vignette +# Assumes Seurat is loaded by the caller, plus SeuratWrappers and Banksy +# when the spatially aware path is used. + +# Seurat stores the map of bin-level barcode (e.g., 8um) -> cluster. To allow Bambu, +# to run quantification, we have to generate the 2um barocde -> cluster map +# as quantData is generated at the 2um resolution +mapSpotsToClusters <- function(binToClustersMap, spotMappings) { + spotToClustersMap <- setNames(unname(binToClustersMap[spotMappings$sample_bin]), spotMappings$sample_barcode) + spotToClustersMap[!is.na(spotToClustersMap)] +} + +# Spatially aware clustering using Banksy (Refer to vignette for more information) +clusterBanksy <- function(object, assay, lambda, kGeom, clusterResolution, npcs = 30) { + # dimx/dimy name the meta.data columns carrying each bin's coordinates + object <- RunBanksy(object, lambda = lambda, assay = assay, slot = "data", + features = "variable", k_geom = kGeom, + dimx = "pxl_col_in_fullres", dimy = "pxl_row_in_fullres", verbose = FALSE) + + DefaultAssay(object) <- "BANKSY" + npcs <- min(npcs, ncol(object) - 1) + object <- RunPCA(object, assay = "BANKSY", reduction.name = "pca.banksy", + features = rownames(object), npcs = npcs, verbose = FALSE) + object <- FindNeighbors(object, reduction = "pca.banksy", dims = 1:npcs, verbose = FALSE) + FindClusters(object, cluster.name = "clusters", resolution = clusterResolution, verbose = FALSE) +} + +# Expression-only clustering (Refer to vignette for more information) +clusterExpression <- function(object, assay, clusterResolution, dims = 15, ncells = 50000) { + object <- FindVariableFeatures(object, verbose = FALSE) + object <- ScaleData(object, verbose = FALSE) + + sketched <- ncol(object) > ncells + if (sketched) { + object <- SketchData(object, ncells = ncells, method = "LeverageScore", sketched.assay = "sketch") + DefaultAssay(object) <- "sketch" + object <- FindVariableFeatures(object, verbose = FALSE) + object <- ScaleData(object, verbose = FALSE) + } + + reduction <- if (sketched) "pca.sketch" else "pca" + npcs <- min(dims, ncol(object) - 1) + object <- RunPCA(object, reduction.name = reduction, npcs = npcs, verbose = FALSE) + object <- FindNeighbors(object, reduction = reduction, dims = 1:npcs, verbose = FALSE) + object <- FindClusters(object, cluster.name = "clusters", resolution = clusterResolution, verbose = FALSE) + + if (sketched) { + # project the sketched labels onto the full set of bins, then adopt them as the + # cluster assignment so both paths emit the same 'clusters' column + object <- ProjectData(object, assay = assay, full.reduction = "full.pca.sketch", + sketched.assay = "sketch", sketched.reduction = "pca.sketch", + dims = 1:npcs, refdata = list(clusters.projected = "clusters")) + DefaultAssay(object) <- assay + object$clusters <- object$clusters.projected + } + + object +} diff --git a/bin/visium_hd_convert_barcode_mappings.py b/bin/visium_hd_convert_barcode_mappings.py new file mode 100755 index 0000000..42c58d4 --- /dev/null +++ b/bin/visium_hd_convert_barcode_mappings.py @@ -0,0 +1,36 @@ +#!/usr/bin/env python3 +import argparse +import sys + +import pandas as pd +import pyarrow.parquet as pq + + +def main(): + parser = argparse.ArgumentParser( + description="Subset a Visium HD barcode_mappings.parquet to a single bin resolution, " + "writing a CSV of the 2um-barcode -> bin assignment with columns 'barcode,bin'" + ) + parser.add_argument("barcode_mappings", help="Path to barcode_mappings.parquet") + parser.add_argument("resolution", help="Bin resolution token, e.g. 008um, 016um") + parser.add_argument("output", help="Output path for the converted CSV") + args = parser.parse_args() + + # read the schema first so only the two needed columns are loaded from the parquet; + # the file has one row per 2um barcode and one column per resolution + schema = pq.read_schema(args.barcode_mappings) + barcode_column = schema.names[0] + bin_column = f"square_{args.resolution}" + if bin_column not in schema.names: + sys.exit( + f"error: column '{bin_column}' not found in {args.barcode_mappings} " + f"(available: {', '.join(schema.names)})" + ) + + barcode_mappings = pd.read_parquet(args.barcode_mappings, columns=[barcode_column, bin_column]) + barcode_mappings = barcode_mappings.rename(columns={barcode_column: "barcode", bin_column: "bin"}) + barcode_mappings.to_csv(args.output, index=False) + + +if __name__ == "__main__": + main() diff --git a/bin/visium_hd_convert_tissue_positions.py b/bin/visium_hd_convert_tissue_positions.py new file mode 100755 index 0000000..8bfe917 --- /dev/null +++ b/bin/visium_hd_convert_tissue_positions.py @@ -0,0 +1,25 @@ +#!/usr/bin/env python3 +import argparse + +import pandas as pd + + +def main(): + parser = argparse.ArgumentParser(description="Convert a Visium HD tissue_positions.parquet to CSV; for the 2um resolution also write the in-tissue barcode list used to filter the BAM") + parser.add_argument("tissue_positions", help="Path to tissue_positions.parquet") + parser.add_argument("resolution", help="Bin resolution token, e.g. 002um, 008um, 016um") + parser.add_argument("csv", help="Output path for the converted CSV") + parser.add_argument("barcodes", help="Output path for the in-tissue barcode list (written only for the 002um resolution)") + args = parser.parse_args() + + tissue_positions = pd.read_parquet(args.tissue_positions) + tissue_positions.to_csv(args.csv, index=False) + + # the BAM carries 2um barcodes, so only the 2um resolution feeds the samtools tissue filter + if args.resolution == "002um": + in_tissue = tissue_positions.loc[tissue_positions["in_tissue"] == 1, "barcode"] + in_tissue.to_csv(args.barcodes, index=False, header=False) + + +if __name__ == "__main__": + main() diff --git a/bin/visium_hd_spot_bin_mappings.py b/bin/visium_hd_spot_bin_mappings.py new file mode 100755 index 0000000..486363d --- /dev/null +++ b/bin/visium_hd_spot_bin_mappings.py @@ -0,0 +1,30 @@ +#!/usr/bin/env python3 +import argparse + +import pandas as pd + + +def main(): + parser = argparse.ArgumentParser(description="Map every 2um spot bambu quantified to its bin at one resolution") + parser.add_argument("barcode_mappings", help="Path to the resolution's barcode mappings CSV (columns: barcode,bin)") + parser.add_argument("barcodes", help="Path to barcodes.tsv.gz from the 2um count directory; one bambu id per column") + parser.add_argument("sample", help="Sample name that bambu prefixes onto each barcode to form the id") + parser.add_argument("output", help="Output path for the spot to bin mapping CSV") + args = parser.parse_args() + + # bambu ids are sampleName_barcode, so strip the prefix to get back to the mapping's barcodes + prefix = f"{args.sample}_" + spot_barcodes = pd.read_csv(args.barcodes, header=None)[0].str.removeprefix(prefix) + + # keep the bins of the spots bambu quantified, then add the prefixed forms so + # downstream modules never have to rebuild a bambu id + barcode_mappings = pd.read_csv(args.barcode_mappings) + spot_mappings = barcode_mappings[barcode_mappings["barcode"].isin(spot_barcodes)].copy() + spot_mappings["sample_barcode"] = prefix + spot_mappings["barcode"] + spot_mappings["sample_bin"] = prefix + spot_mappings["bin"] + + spot_mappings[["sample_barcode", "barcode", "sample_bin", "bin"]].to_csv(args.output, index=False) + + +if __name__ == "__main__": + main() diff --git a/conf/containers.config b/conf/containers.config index 3cac435..83db49c 100644 --- a/conf/containers.config +++ b/conf/containers.config @@ -1,6 +1,7 @@ process { - withLabel: 'r' { container = "ghcr.io/goekelab/bambu-pipe-r:1.0.0" } + withLabel: 'r' { container = "ghcr.io/goekelab/bambu-pipe-r:1.1.0" } withLabel: 'spaceranger' { container = "quay.io/nf-core/spaceranger:9c5e7dc93c32448e" } withLabel: 'minimap2_samtools' { container = "community.wave.seqera.io/library/minimap2_samtools:b09096fc890429ce" } withLabel: 'preprocess' { container = "community.wave.seqera.io/library/chopper_cutadapt_flexiplex_pigz:077c3bc67452482c" } -} \ No newline at end of file + withLabel: 'pyarrow_pandas' { container = "community.wave.seqera.io/library/pip_pandas_pyarrow:e60a5c578ef189b2" } +} diff --git a/containers/r/Dockerfile b/containers/r/Dockerfile index fcdd656..56cf9f1 100644 --- a/containers/r/Dockerfile +++ b/containers/r/Dockerfile @@ -16,10 +16,14 @@ RUN micromamba install -y -n base -c conda-forge -c bioconda \ bioconductor-dropletutils=1.30.0 \ r-seurat=5.4.0 \ r-harmony=2.0.2 \ + bioconductor-banksy=1.4.0 \ r-data.table=1.17.8 \ gxx=15.2.0 \ && micromamba clean --all --yes +# SeuratWrappers provides RunBanksy() and is not packaged on conda +RUN R -e 'devtools::install_github("satijalab/seurat-wrappers", upgrade = "never")' + # Clone bambu single cell feature branch RUN git clone --branch devel_pre_v4 \ https://github.com/GoekeLab/bambu.git /opt/bambu \ diff --git a/lib/Validation.groovy b/lib/Validation.groovy index bc7ae47..2ffb5f9 100644 --- a/lib/Validation.groovy +++ b/lib/Validation.groovy @@ -23,11 +23,87 @@ class Validation { throw new Exception("Invalid params.quantification_mode '${params.quantification_mode}' — must be one of: ${params.valid_quantification_modes.join(', ')}") // Numeric range checks - if (params.resolution <= 0) - throw new Exception("Invalid params.resolution '${params.resolution}' — must be a positive number") + if (params.cluster_resolution <= 0) + throw new Exception("Invalid params.cluster_resolution '${params.cluster_resolution}' — must be a positive number") if (params.ndr != null && (params.ndr < 0 || params.ndr > 1)) throw new Exception("Invalid params.ndr '${params.ndr}' — must be a float between 0 and 1") + + // Visium HD checks + if (params.visium_hd) { + if (params.bins == null) + throw new Exception("params.bins is required when params.visium_hd is true") + + if (!params.bins.exists()) + throw new Exception("params.bins '${params.bins}' does not exist") + + if (params.bins.extension != 'csv') + throw new Exception("params.bins '${params.bins}' must be a CSV file") + + validateVisiumHDBins(params.bins) + + if (params.barcode_mappings == null) + throw new Exception("params.barcode_mappings is required when params.visium_hd is true") + + if (!params.barcode_mappings.exists()) + throw new Exception("params.barcode_mappings '${params.barcode_mappings}' does not exist") + + if (params.barcode_mappings.extension != 'parquet') + throw new Exception("params.barcode_mappings '${params.barcode_mappings}' must be a .parquet file") + + // only the clustering mode reads the clustering bin, so only it has to resolve + if (params.quantification_mode == 'EM_clusters') + validateClusteringBin(params.bins, params.clustering_bin) + } + } + + static def validateVisiumHDRows(rows) { + if (rows.size() != 1) + throw new Exception("Visium HD requires exactly 1 sample in the samplesheet, but found ${rows.size()}") + + def row = rows[0] + ["sample", "path"].each { col -> + if (!row.containsKey(col)) + throw new Exception("Samplesheet is missing a required '${col}' column") + if (!row[col]) + throw new Exception("A row in the samplesheet has an empty '${col}' value") + } + + if (!row.path.endsWith('.bam')) + throw new Exception("Visium HD sample '${row.sample}' must point to a pre-aligned, barcode-tagged BAM file") + } + + static def validateVisiumHDBins(bins) { + def lines = bins.text.readLines().findAll { line -> line.trim() } + if (lines.size() < 2) + throw new Exception("params.bins '${bins}' must have a header row and at least one resolution row") + + def header = lines[0].split(',').collect { col -> col.trim() } + ["resolution", "tissue_positions"].each { col -> + if (!header.contains(col)) + throw new Exception("params.bins '${bins}' is missing a required '${col}' column") + } + + if (!readBinResolutions(bins).contains("2")) + throw new Exception("params.bins '${bins}' must include a row for the native 2 um resolution (resolution = 2)") + } + + static def readBinResolutions(bins) { + def lines = bins.text.readLines().findAll { line -> line.trim() } + def header = lines[0].split(',').collect { col -> col.trim() } + def resIdx = header.indexOf("resolution") + lines[1..-1].collect { line -> line.split(',')[resIdx].trim() } + } + + static def validateClusteringBin(bins, clusteringBin) { + if (clusteringBin == 2) + throw new Exception("Invalid params.clustering_bin '2' — 2 um bins are too sparse to cluster, choose a coarser bin") + + def resolutions = readBinResolutions(bins) + if (!resolutions.contains(clusteringBin.toString())) { + def bin_options = resolutions.findAll { res -> res != "2" } + throw new Exception("params.clustering_bin '${clusteringBin}' is not listed in params.bins '${bins}' — available bins: ${bin_options.join(', ')}") + } } static def validateVisiumSampleCount(samples) { diff --git a/main.nf b/main.nf index 2b45d93..afa53c5 100644 --- a/main.nf +++ b/main.nf @@ -2,16 +2,26 @@ nextflow.enable.types = true -include { DECOMPRESS as DECOMPRESS_GENOME } from './modules/decompress.nf' -include { DECOMPRESS as DECOMPRESS_ANNOTATION } from './modules/decompress.nf' -include { PREPARE_INPUT_STANDARD } from './subworkflows/prepare_input_standard.nf' -include { PREPROCESS_FASTQ } from './modules/preprocess_fastq.nf' -include { ALIGNMENT } from './subworkflows/alignment.nf' -include { BAMBU_CONSTRUCT_READ_CLASS } from './modules/bambu/construct_read_class.nf' -include { BAMBU_PREPARE_ANNOTATION } from './modules/bambu/prepare_annotation.nf' -include { BAMBU_TRANSCRIPT_DISCOVERY } from './modules/bambu/transcript_discovery.nf' -include { CLUSTERING } from './subworkflows/clustering.nf' -include { BAMBU_EM } from './modules/bambu/EM_quant.nf' +// subworkflows +include { PREPARE_INPUT_STANDARD } from './subworkflows/prepare_input_standard.nf' +include { PREPARE_INPUT_VISIUM_HD } from './subworkflows/prepare_input_visium_hd.nf' +include { ALIGNMENT } from './subworkflows/alignment.nf' +include { CLUSTERING } from './subworkflows/clustering_standard.nf' + +// modules +include { DECOMPRESS as DECOMPRESS_GENOME } from './modules/decompress.nf' +include { DECOMPRESS as DECOMPRESS_ANNOTATION } from './modules/decompress.nf' +include { PREPROCESS_FASTQ } from './modules/preprocess_fastq.nf' +include { BAMBU_PREPARE_ANNOTATION } from './modules/bambu/shared/prepare_annotation.nf' +include { BAMBU_CONSTRUCT_READ_CLASS } from './modules/bambu/shared/construct_read_class.nf' +include { BAMBU_TRANSCRIPT_DISCOVERY } from './modules/bambu/standard/transcript_discovery.nf' +include { BAMBU_CLUSTERED_EM } from './modules/bambu/shared/clustered_quantification.nf' +include { BAMBU_EM } from './modules/bambu/standard/single_cell_quantification.nf' +include { BAMBU_TRANSCRIPT_DISCOVERY_VISIUM_HD } from './modules/bambu/visium_hd/transcript_discovery.nf' +include { BAMBU_EM_VISIUM_HD } from './modules/bambu/visium_hd/spot_level_quantification.nf' +include { AGGREGATE_BINS_VISIUM_HD } from './modules/bambu/visium_hd/aggregate_bins.nf' +include { SPOT_BIN_MAPPINGS } from './modules/prepare_input_visium_hd/spot_bin_mappings.nf' +include { SEURAT_VISIUM_HD } from './modules/seurat/visium_hd/clustering.nf' params { input: Path @@ -25,40 +35,34 @@ params { ndr: Float? deduplicate_umis: Boolean quantification_mode: String - resolution: Float + cluster_resolution: Float + visium_hd: Boolean + barcode_mappings: Path? + bins: Path? + clustering_bin: Integer + banksy: Boolean + banksy_lambda: Float + banksy_k_geom: Integer } -workflow { - Validation.validateParams(params, workflow) - - def ndr = params.ndr ?: 'NULL' - - // load reference files - ch_genome = channel.value(params.genome) - ch_annotation = channel.value(params.annotation) - - if (params.genome.extension == 'gz') { - DECOMPRESS_GENOME(ch_genome) - ch_genome = DECOMPRESS_GENOME.out - } +workflow STANDARD { + take: + ch_rows: Channel + ch_genome: Path + ch_annotation: Path + ndr: Float? - if (params.annotation.extension == 'gz') { - DECOMPRESS_ANNOTATION(ch_annotation) - ch_annotation = DECOMPRESS_ANNOTATION.out - } + main: + def ndrArg = ndr != null ? ndr : 'NULL' // load config files ch_barcode_coordinate_config = file("${projectDir}/assets/10x_config/barcode_coordinate_config.csv", checkIfExists: true) ch_adapter_seq_config = file("${projectDir}/assets/10x_config/adapter_seq_config.csv", checkIfExists: true) ch_flank_seq_config = file("${projectDir}/assets/10x_config/flank_seq_config.csv", checkIfExists: true) - // parsing samplesheet csv file - ch_input = channel.of(params.input) - - ch_standard = ch_input.splitCsv(header:true, sep:',') - ch_n_samples = ch_standard.count() + ch_n_samples = ch_rows.count() - PREPARE_INPUT_STANDARD(ch_standard, ch_barcode_coordinate_config) + PREPARE_INPUT_STANDARD(ch_rows, ch_barcode_coordinate_config) // input files are split by type (fastq, bam) ch_input_fastq = PREPARE_INPUT_STANDARD.out.fastq @@ -84,18 +88,94 @@ workflow { def has_spatial = metas.any { meta -> meta.chemistry.startsWith('visium') } // for non-visium samples set the spatial metadata to an empty list (for staging) [samples, paths, metas, has_spatial ? spatial_metadatas : []] } - BAMBU_TRANSCRIPT_DISCOVERY(ch_rds_files_collect, ch_genome, BAMBU_PREPARE_ANNOTATION.out.annotation, ndr) - - if (params.quantification_mode != 'no_quant') { - if (params.quantification_mode == 'EM_clusters') { - CLUSTERING(BAMBU_TRANSCRIPT_DISCOVERY.out.se_gene_counts, ch_n_samples) - ch_clusters = CLUSTERING.out.clusters.map { clusters -> [true, clusters] } // flag to indicate that clustering was performed - } else { - ch_clusters = channel.value([false, []]) // flag to indicate that clustering was not performed - } - BAMBU_EM(BAMBU_TRANSCRIPT_DISCOVERY.out.quant_data, BAMBU_TRANSCRIPT_DISCOVERY.out.extended_annotations, ch_clusters, ch_genome) + BAMBU_TRANSCRIPT_DISCOVERY(ch_rds_files_collect, ch_genome, BAMBU_PREPARE_ANNOTATION.out.annotation, ndrArg) + + // cluster the cells first, then pool each cluster's cells for the EM + if (params.quantification_mode == 'EM_clusters') { + CLUSTERING(BAMBU_TRANSCRIPT_DISCOVERY.out.se_gene_counts, ch_n_samples) + BAMBU_CLUSTERED_EM(CLUSTERING.out.clusters, BAMBU_TRANSCRIPT_DISCOVERY.out.quant_data, BAMBU_TRANSCRIPT_DISCOVERY.out.extended_annotations, ch_genome) + } else if (params.quantification_mode == 'EM') { + BAMBU_EM(BAMBU_TRANSCRIPT_DISCOVERY.out.quant_data, BAMBU_TRANSCRIPT_DISCOVERY.out.extended_annotations, ch_genome) } } +} + +workflow VISIUM_HD { + take: + ch_rows: Channel + ch_genome: Path + ch_annotation: Path + ndr: Float? + + main: + def ndrArg = ndr != null ? ndr : 'NULL' + + PREPARE_INPUT_VISIUM_HD(ch_rows) + BAMBU_PREPARE_ANNOTATION(ch_annotation) + BAMBU_CONSTRUCT_READ_CLASS(PREPARE_INPUT_VISIUM_HD.out.bam, ch_genome, BAMBU_PREPARE_ANNOTATION.out.annotation) + BAMBU_TRANSCRIPT_DISCOVERY_VISIUM_HD(BAMBU_CONSTRUCT_READ_CLASS.out.rds, ch_genome, BAMBU_PREPARE_ANNOTATION.out.annotation, ndrArg, PREPARE_INPUT_VISIUM_HD.out.tissue_positions_002um) + + // perform transcript discovery at 2um first + ch_quant_data = BAMBU_TRANSCRIPT_DISCOVERY_VISIUM_HD.out.quant_data.first() + ch_extended_anno = BAMBU_TRANSCRIPT_DISCOVERY_VISIUM_HD.out.extended_annotations.first() + ch_unique_002um = BAMBU_TRANSCRIPT_DISCOVERY_VISIUM_HD.out.se_unique_002um.first() // unique counts at 2um resolution + ch_barcodes_002um = BAMBU_TRANSCRIPT_DISCOVERY_VISIUM_HD.out.barcodes_002um.first() // list of 2um barcodes (same as the column names) + + // map every 2um spot in the SE to its bin, once per resolution; every module below reads this file + SPOT_BIN_MAPPINGS(PREPARE_INPUT_VISIUM_HD.out.barcode_mappings, ch_barcodes_002um, PREPARE_INPUT_VISIUM_HD.out.sample_name) + + // pair each bin's tissue positions with its spot mappings, keyed on resolution + ch_bins = PREPARE_INPUT_VISIUM_HD.out.tissue_positions_bins.join(SPOT_BIN_MAPPINGS.out.csv) // [resolution, tissue_positions, spot_mappings] + + // aggregate the 2um SEs into each requested bin resolution (e.g., 8um/16um) + AGGREGATE_BINS_VISIUM_HD(ch_bins, ch_unique_002um) + + if (params.quantification_mode == 'EM_clusters') { + def requested_bin = String.format('%03dum', params.clustering_bin) // convert clustering_bin specified as an integer into Spaceranger format + + // perform clustering at the requested resolution only + ch_clustering = AGGREGATE_BINS_VISIUM_HD.out.se_gene_counts + .filter { resolution, _se_gene_counts -> resolution == requested_bin } + .join(SPOT_BIN_MAPPINGS.out.csv) // [resolution, se_gene_counts, spot_mappings] + SEURAT_VISIUM_HD(ch_clustering) + BAMBU_CLUSTERED_EM(SEURAT_VISIUM_HD.out.clusters, ch_quant_data, ch_extended_anno, ch_genome) + } + + if (params.quantification_mode != 'no_quant') { + // run spot level quantification on all resolution + // at 2um resolution, tissue_positions and spot_mappings are not required + ch_resolution_002um = channel.of(['002um', [], []]) + ch_resolutions = ch_resolution_002um.mix(ch_bins) // [resolution, tissue_positions, spot_mappings] + BAMBU_EM_VISIUM_HD(ch_resolutions, ch_quant_data, ch_extended_anno, ch_genome) + } +} + +workflow { + Validation.validateParams(params, workflow) + + // load reference files + ch_genome = channel.value(params.genome) + ch_annotation = channel.value(params.annotation) + + if (params.genome.extension == 'gz') { + DECOMPRESS_GENOME(ch_genome) + ch_genome = DECOMPRESS_GENOME.out + } + + if (params.annotation.extension == 'gz') { + DECOMPRESS_ANNOTATION(ch_annotation) + ch_annotation = DECOMPRESS_ANNOTATION.out + } + + // parsing samplesheet csv file + ch_input = channel.of(params.input) + ch_rows = ch_input.splitCsv(header:true, sep:',') + + if (params.visium_hd) { + VISIUM_HD(ch_rows, ch_genome, ch_annotation, params.ndr) + } else { + STANDARD(ch_rows, ch_genome, ch_annotation, params.ndr) + } channel.topic('versions').collectFile(name: 'software_versions.yml', storeDir: "${params.output_dir}") -} \ No newline at end of file +} diff --git a/modules/bambu/EM_quant.nf b/modules/bambu/shared/clustered_quantification.nf similarity index 55% rename from modules/bambu/EM_quant.nf rename to modules/bambu/shared/clustered_quantification.nf index b10b28e..f907036 100644 --- a/modules/bambu/EM_quant.nf +++ b/modules/bambu/shared/clustered_quantification.nf @@ -1,5 +1,4 @@ -process BAMBU_EM{ - publishDir "$params.output_dir", mode: 'copy', pattern: 'transcript_counts_singlecell' +process BAMBU_CLUSTERED_EM { publishDir "$params.output_dir", mode: 'copy', pattern: 'transcript_counts_clusters' publishDir "$params.output_dir", mode: 'copy', pattern: 'gene_counts_clusters' label "r" @@ -8,15 +7,14 @@ process BAMBU_EM{ label "long" input: + path(clusters) path(quant_data) path(extended_annotation) - tuple val(has_clusters), path(clusters) path(genome) output: - path ('transcript_counts_singlecell'), optional: true - path ('transcript_counts_clusters'), optional: true - path ('gene_counts_clusters'), optional: true + path ('transcript_counts_clusters') + path ('gene_counts_clusters') path "versions.yml", topic: 'versions' script: @@ -28,29 +26,24 @@ process BAMBU_EM{ extendedAnno <- readRDS("$extended_annotation") quantData <- readRDS("$quant_data") - clusters <- if ("$has_clusters" == "true") readRDS("$clusters") else NULL - degBias <- !is.null(clusters) + clusters <- readRDS("$clusters") se <- bambu.singlecell( reads = quantData, - output = if (is.null(clusters)) "EM" else "clusteredEM", + output = "clusteredEM", annotations = extendedAnno, genome = "$genome", ncore = $task.cpus, verbose = FALSE, - opt.em = list(degradationBias = degBias), + opt.em = list(degradationBias = TRUE), clusters = clusters ) - if (is.null(clusters)) { - save_counts(se, "transcript_counts_singlecell", "Transcript Expression") - } else { - save_counts(se, "transcript_counts_clusters", "Transcript Expression") + saveCounts(se, "transcript_counts_clusters", "Transcript Expression") - se_gene <- transcriptToGeneExpression(se) - save_counts(se_gene, "gene_counts_clusters") - } + seGene <- transcriptToGeneExpression(se) + saveCounts(seGene, "gene_counts_clusters") writeLines(c('"${task.process}":', paste0(' R: ', R.Version()\$version.string), paste0(' bambu: ', as.character(packageVersion("bambu")))), "versions.yml") """ -} \ No newline at end of file +} diff --git a/modules/bambu/construct_read_class.nf b/modules/bambu/shared/construct_read_class.nf similarity index 100% rename from modules/bambu/construct_read_class.nf rename to modules/bambu/shared/construct_read_class.nf diff --git a/modules/bambu/prepare_annotation.nf b/modules/bambu/shared/prepare_annotation.nf similarity index 100% rename from modules/bambu/prepare_annotation.nf rename to modules/bambu/shared/prepare_annotation.nf diff --git a/modules/bambu/standard/single_cell_quantification.nf b/modules/bambu/standard/single_cell_quantification.nf new file mode 100644 index 0000000..97df050 --- /dev/null +++ b/modules/bambu/standard/single_cell_quantification.nf @@ -0,0 +1,41 @@ +process BAMBU_EM { + publishDir "$params.output_dir", mode: 'copy', pattern: 'transcript_counts_singlecell' + label "r" + label "low_cpu" + label "high_mem" + label "long" + + input: + path(quant_data) + path(extended_annotation) + path(genome) + + output: + path ('transcript_counts_singlecell') + path "versions.yml", topic: 'versions' + + script: + """ + #!/usr/bin/env Rscript + if ("$params.bambu_path" == "null") { library("bambu") } else { library("devtools"); load_all("$params.bambu_path") } + library(DropletUtils) + source(Sys.which("save_counts.R")) + + extendedAnno <- readRDS("$extended_annotation") + quantData <- readRDS("$quant_data") + + se <- bambu.singlecell( + reads = quantData, + output = "EM", + annotations = extendedAnno, + genome = "$genome", + ncore = $task.cpus, + verbose = FALSE, + opt.em = list(degradationBias = FALSE) + ) + + saveCounts(se, "transcript_counts_singlecell", "Transcript Expression") + + writeLines(c('"${task.process}":', paste0(' R: ', R.Version()\$version.string), paste0(' bambu: ', as.character(packageVersion("bambu")))), "versions.yml") + """ +} diff --git a/modules/bambu/transcript_discovery.nf b/modules/bambu/standard/transcript_discovery.nf similarity index 71% rename from modules/bambu/transcript_discovery.nf rename to modules/bambu/standard/transcript_discovery.nf index 8214740..e89d0d4 100644 --- a/modules/bambu/transcript_discovery.nf +++ b/modules/bambu/standard/transcript_discovery.nf @@ -29,6 +29,7 @@ process BAMBU_TRANSCRIPT_DISCOVERY{ #!/usr/bin/env Rscript if ("$params.bambu_path" == "null") { library("bambu") } else { library("devtools"); load_all("$params.bambu_path") } library(DropletUtils) + source(Sys.which("bambu_discovery.R")) source(Sys.which("save_counts.R")) annotation <- readRDS("$bambu_annotation") @@ -37,28 +38,22 @@ process BAMBU_TRANSCRIPT_DISCOVERY{ sampleData <- strsplit("${spatial_metadata_files.join(',')}", ",")[[1]] chemistry <- setNames(strsplit("${meta.collect { m -> m.chemistry }.join(',')}", ",")[[1]], sampleNames) technology <- setNames(strsplit("${meta.collect { m -> m.technology }.join(',')}", ",")[[1]], sampleNames) - - # Transcript discovery - extendedAnno <- bambu.singlecell(reads = readClassFile, output = "extendedAnnotations", - annotations = annotation, genome = "$genome", ncore = $task.cpus, verbose = FALSE, NDR = $ndr) - saveRDS(extendedAnno, "extended_annotations.rds") - writeToGTF(extendedAnno, "extended_annotations.gtf") - - # Read to transcript assignment sampleData <- if (any(startsWith(chemistry, "visium-v"))) sampleData else NULL # Add spatial metadata for visium samples - quantData <- bambu.singlecell(reads = readClassFile, output = "quantData", - annotations = extendedAnno, genome = "$genome", ncore = $task.cpus, verbose = FALSE, sampleData = sampleData) - saveRDS(quantData, "quant_data.rds") - # Quantification without EM - seDiscovery <- bambu.singlecell(reads = quantData, output = "uniqueCounts", annotations = extendedAnno) + result <- bambuDiscovery(reads = readClassFile, annotation = annotation, genome = "$genome", + ncore = $task.cpus, ndr = $ndr, sampleData = sampleData) + saveRDS(result\$extendedAnno, "extended_annotations.rds") + writeToGTF(result\$extendedAnno, "extended_annotations.gtf") + saveRDS(result\$quantData, "quant_data.rds") + + seDiscovery <- result\$seDiscovery colData(seDiscovery)\$chemistry <- unname(chemistry[colData(seDiscovery)\$sampleName]) # Add chemistry into colData (for subsequent batch correction) colData(seDiscovery)\$technology <- unname(technology[colData(seDiscovery)\$sampleName]) # Add technology into colData (for subsequent batch correction) - save_counts(seDiscovery, "unique_counts", "Transcript Expression") + saveCounts(seDiscovery, "unique_counts", "Transcript Expression") # Generate gene counts SE from unique counts SE - seDiscovery.gene <- transcriptToGeneExpression(seDiscovery) - save_counts(seDiscovery.gene, "gene_counts") + seDiscoveryGene <- transcriptToGeneExpression(seDiscovery) + saveCounts(seDiscoveryGene, "gene_counts") writeLines(c('"${task.process}":', paste0(' R: ', R.Version()\$version.string), paste0(' bambu: ', as.character(packageVersion("bambu")))), "versions.yml") """ diff --git a/modules/bambu/visium_hd/aggregate_bins.nf b/modules/bambu/visium_hd/aggregate_bins.nf new file mode 100644 index 0000000..1996fe9 --- /dev/null +++ b/modules/bambu/visium_hd/aggregate_bins.nf @@ -0,0 +1,42 @@ +process AGGREGATE_BINS_VISIUM_HD { + publishDir "$params.output_dir", mode: 'copy', pattern: 'unique_counts_*' + publishDir "$params.output_dir", mode: 'copy', pattern: 'gene_counts_*' + label "r" + label "low_cpu" + label "medium_mem" + label "short" + + input: + tuple val(resolution), path(tissue_positions), path(spot_mappings) + path(se_unique_002um) + + output: + tuple val(resolution), path("unique_counts_${resolution}"), emit: unique_counts + tuple val(resolution), path("gene_counts_${resolution}"), emit: gene_counts + tuple val(resolution), path("gene_counts_${resolution}/se_gene_counts_${resolution}.rds"), emit: se_gene_counts + path "versions.yml", topic: 'versions' + + script: + """ + #!/usr/bin/env Rscript + if ("$params.bambu_path" == "null") { library("bambu") } else { library("devtools"); load_all("$params.bambu_path") } + library(DropletUtils) + library(Matrix) + library(dplyr) + source(Sys.which("save_counts.R")) + source(Sys.which("visium_hd_aggregate_resolution.R")) + + unique002um <- readRDS("$se_unique_002um") + tissuePositions <- read.csv("$tissue_positions", stringsAsFactors = FALSE) + spotMappings <- read.csv("$spot_mappings", stringsAsFactors = FALSE) + + # aggregate transcript and gene counts + uniqueAgg <- aggregateResolution(unique002um, tissuePositions, spotMappings) + geneAgg <- transcriptToGeneExpression(uniqueAgg) + + saveCounts(uniqueAgg, "unique_counts_$resolution", "Transcript Expression") + saveCounts(geneAgg, "gene_counts_$resolution") + + writeLines(c('"${task.process}":', paste0(' R: ', R.Version()\$version.string), paste0(' bambu: ', as.character(packageVersion("bambu")))), "versions.yml") + """ +} diff --git a/modules/bambu/visium_hd/spot_level_quantification.nf b/modules/bambu/visium_hd/spot_level_quantification.nf new file mode 100644 index 0000000..52b5023 --- /dev/null +++ b/modules/bambu/visium_hd/spot_level_quantification.nf @@ -0,0 +1,63 @@ +process BAMBU_EM_VISIUM_HD { + publishDir "$params.output_dir", mode: 'copy', pattern: 'transcript_counts_*' + label "r" + label "low_cpu" + label "high_mem" + label "long" + + input: + tuple val(resolution), path(tissue_positions), path(spot_mappings) + path(quant_data) + path(extended_annotation) + path(genome) + + output: + path ("transcript_counts_${resolution}") + path "versions.yml", topic: 'versions' + + script: + """ + #!/usr/bin/env Rscript + if ("$params.bambu_path" == "null") { library("bambu") } else { library("devtools"); load_all("$params.bambu_path") } + library(DropletUtils) + library(dplyr) + source(Sys.which("save_counts.R")) + + extendedAnno <- readRDS("$extended_annotation") + quantData <- readRDS("$quant_data") + + # For bin-level quantification, we need to use the clusters argument to collapse the quantData object + # to the desired resolution, since quantData is generated at 2um resolution. + if ("$resolution" == "002um") { + clusters <- NULL + } else { + # bambu prefixes the sample name onto each cluster label, so pass the bare bin + spotMappings <- read.csv("$spot_mappings", stringsAsFactors = FALSE) + clusters <- setNames(spotMappings\$bin, spotMappings\$sample_barcode) + } + + se <- bambu.singlecell( + reads = quantData, + output = if (is.null(clusters)) "EM" else "clusteredEM", + annotations = extendedAnno, + genome = "$genome", + ncore = $task.cpus, + verbose = FALSE, + opt.em = list(degradationBias = FALSE), + clusters = clusters + ) + + # At the bin-level resolutions bambu builds a fresh colData, so we need to reattach the bin metadata. + if (!is.null(clusters)) { + tissuePositions <- read.csv("$tissue_positions", stringsAsFactors = FALSE) + binColData <- as.data.frame(colData(se)) %>% + left_join(tissuePositions, by = c("cluster" = "barcode")) %>% + rename(barcode = cluster) + colData(se) <- DataFrame(binColData, row.names = colnames(se)) + } + + saveCounts(se, "transcript_counts_$resolution", "Transcript Expression") + + writeLines(c('"${task.process}":', paste0(' R: ', R.Version()\$version.string), paste0(' bambu: ', as.character(packageVersion("bambu")))), "versions.yml") + """ +} diff --git a/modules/bambu/visium_hd/transcript_discovery.nf b/modules/bambu/visium_hd/transcript_discovery.nf new file mode 100644 index 0000000..b096507 --- /dev/null +++ b/modules/bambu/visium_hd/transcript_discovery.nf @@ -0,0 +1,54 @@ +process BAMBU_TRANSCRIPT_DISCOVERY_VISIUM_HD { + publishDir "$params.output_dir", mode: 'copy', pattern: 'extended_annotations.gtf' + publishDir "$params.output_dir", mode: 'copy', pattern: 'unique_counts_002um' + publishDir "$params.output_dir", mode: 'copy', pattern: 'gene_counts_002um' + publishDir "$params.output_dir/intermediate_R", mode: 'copy', pattern: '*.rds', enabled: params.save_intermediates + label "r" + label "medium_cpu" + label "high_mem" + label "medium" + + input: + tuple val(sample), path(rds_files), val(meta) + path(genome) + path(bambu_annotation) + val(ndr) + path(tissue_positions) + + output: + path ('quant_data.rds'), emit: quant_data + path ('extended_annotations.rds'), emit: extended_annotations + path ('extended_annotations.gtf') + path ('unique_counts_002um') + path ('gene_counts_002um') + path ('unique_counts_002um/se_unique_counts_002um.rds'), emit: se_unique_002um + path ('unique_counts_002um/barcodes.tsv.gz'), emit: barcodes_002um + path "versions.yml", topic: 'versions' + + script: + """ + #!/usr/bin/env Rscript + if ("$params.bambu_path" == "null") { library("bambu") } else { library("devtools"); load_all("$params.bambu_path") } + library(DropletUtils) + source(Sys.which("bambu_discovery.R")) + source(Sys.which("save_counts.R")) + + annotation <- readRDS("$bambu_annotation") + readClassFile <- setNames("$rds_files", "$sample") + + # the 2um tissue_positions is passed as sampleData so bambu joins the per-spot spatial metadata into colData + result <- bambuDiscovery(reads = readClassFile, annotation = annotation, genome = "$genome", + ncore = $task.cpus, ndr = $ndr, sampleData = "$tissue_positions") + saveRDS(result\$extendedAnno, "extended_annotations.rds") + writeToGTF(result\$extendedAnno, "extended_annotations.gtf") + saveRDS(result\$quantData, "quant_data.rds") + + unique002um <- result\$seDiscovery + gene002um <- transcriptToGeneExpression(unique002um) + + saveCounts(unique002um, "unique_counts_002um", "Transcript Expression") + saveCounts(gene002um, "gene_counts_002um") + + writeLines(c('"${task.process}":', paste0(' R: ', R.Version()\$version.string), paste0(' bambu: ', as.character(packageVersion("bambu")))), "versions.yml") + """ +} diff --git a/modules/prepare_input_visium_hd/convert_barcode_mappings.nf b/modules/prepare_input_visium_hd/convert_barcode_mappings.nf new file mode 100644 index 0000000..422ba04 --- /dev/null +++ b/modules/prepare_input_visium_hd/convert_barcode_mappings.nf @@ -0,0 +1,26 @@ +process CONVERT_BARCODE_MAPPINGS { + publishDir "$params.output_dir/intermediate_visium_hd", mode: 'copy', pattern: '*.csv', enabled: params.save_intermediates + label "pyarrow_pandas" + label "low_cpu" + label "low_mem" + label "short" + + input: + val(resolution) + path(barcode_mappings) + + output: + tuple val(resolution), path("barcode_mappings_${resolution}.csv"), emit: csv + path "versions.yml", topic: 'versions' + + script: + """ + visium_hd_convert_barcode_mappings.py $barcode_mappings ${resolution} barcode_mappings_${resolution}.csv + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + python: \$(python3 --version 2>&1) + pandas: \$(python3 -c 'import pandas; print(pandas.__version__)') + END_VERSIONS + """ +} diff --git a/modules/prepare_input_visium_hd/convert_tissue_positions.nf b/modules/prepare_input_visium_hd/convert_tissue_positions.nf new file mode 100644 index 0000000..580df75 --- /dev/null +++ b/modules/prepare_input_visium_hd/convert_tissue_positions.nf @@ -0,0 +1,27 @@ +process CONVERT_TISSUE_POSITIONS { + publishDir "$params.output_dir/intermediate_visium_hd", mode: 'copy', pattern: '*.csv', enabled: params.save_intermediates + publishDir "$params.output_dir/intermediate_visium_hd", mode: 'copy', pattern: 'barcodes_in_tissue.txt', enabled: params.save_intermediates + label "pyarrow_pandas" + label "low_cpu" + label "low_mem" + label "short" + + input: + tuple val(resolution), path(tissue_positions) + + output: + tuple val(resolution), path("tissue_positions_${resolution}.csv"), emit: csv + path("barcodes_in_tissue.txt"), optional: true, emit: barcodes + path "versions.yml", topic: 'versions' + + script: + """ + visium_hd_convert_tissue_positions.py $tissue_positions ${resolution} tissue_positions_${resolution}.csv barcodes_in_tissue.txt + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + python: \$(python3 --version 2>&1) + pandas: \$(python3 -c 'import pandas; print(pandas.__version__)') + END_VERSIONS + """ +} diff --git a/modules/prepare_input_visium_hd/filter_barcoded_bam.nf b/modules/prepare_input_visium_hd/filter_barcoded_bam.nf new file mode 100644 index 0000000..9b9720b --- /dev/null +++ b/modules/prepare_input_visium_hd/filter_barcoded_bam.nf @@ -0,0 +1,27 @@ +process FILTER_BARCODED_BAM { + publishDir "$params.output_dir/intermediate_bam", mode: 'copy', pattern: '*_filtered.bam*', enabled: params.save_intermediates + label "minimap2_samtools" + label "medium_cpu" + label "medium_mem" + label "medium" + + input: + tuple val(sample), path(bam), val(meta) + path(barcodes) + + output: + tuple val(sample), path("${sample}_filtered.bam"), val(meta), emit: bam + path("${sample}_filtered.bam.bai") + path "versions.yml", topic: 'versions' + + script: + """ + samtools view -@ $task.cpus -D CB:$barcodes -o ${sample}_filtered.bam $bam + samtools index -@ $task.cpus ${sample}_filtered.bam + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + samtools: \$(samtools --version 2>&1 | head -1) + END_VERSIONS + """ +} diff --git a/modules/prepare_input_visium_hd/spot_bin_mappings.nf b/modules/prepare_input_visium_hd/spot_bin_mappings.nf new file mode 100644 index 0000000..8b0a9ac --- /dev/null +++ b/modules/prepare_input_visium_hd/spot_bin_mappings.nf @@ -0,0 +1,27 @@ +process SPOT_BIN_MAPPINGS { + publishDir "$params.output_dir/intermediate_visium_hd", mode: 'copy', pattern: '*.csv.gz', enabled: params.save_intermediates + label "pyarrow_pandas" + label "low_cpu" + label "low_mem" + label "short" + + input: + tuple val(resolution), path(barcode_mappings) + path(barcodes_002um) + val(sample) + + output: + tuple val(resolution), path("spot_bin_mappings_${resolution}.csv.gz"), emit: csv + path "versions.yml", topic: 'versions' + + script: + """ + visium_hd_spot_bin_mappings.py $barcode_mappings $barcodes_002um $sample spot_bin_mappings_${resolution}.csv.gz + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + python: \$(python3 --version 2>&1) + pandas: \$(python3 -c 'import pandas; print(pandas.__version__)') + END_VERSIONS + """ +} diff --git a/modules/seurat/multi_sample.nf b/modules/seurat/standard/multi_sample_clustering.nf similarity index 95% rename from modules/seurat/multi_sample.nf rename to modules/seurat/standard/multi_sample_clustering.nf index 2316cfb..b377b95 100644 --- a/modules/seurat/multi_sample.nf +++ b/modules/seurat/standard/multi_sample_clustering.nf @@ -50,7 +50,7 @@ process SEURAT_MULTI_SAMPLE { dim <- min(dim, ncol(cellMix[["harmony"]])) cellMix <- FindNeighbors(cellMix, reduction = "harmony", dims = 1:dim) - cellMix <- FindClusters(cellMix, resolution = $params.resolution, cluster.name = "harmony_clusters") + cellMix <- FindClusters(cellMix, resolution = $params.cluster_resolution, cluster.name = "harmony_clusters") saveRDS(cellMix, "seurat_obj.rds") clusters <- setNames(paste0("cluster_", cellMix\$harmony_clusters), names(cellMix\$harmony_clusters)) diff --git a/modules/seurat/single_sample.nf b/modules/seurat/standard/single_sample_clustering.nf similarity index 94% rename from modules/seurat/single_sample.nf rename to modules/seurat/standard/single_sample_clustering.nf index 283374c..4730ba8 100644 --- a/modules/seurat/single_sample.nf +++ b/modules/seurat/standard/single_sample_clustering.nf @@ -33,7 +33,7 @@ process SEURAT_SINGLE_SAMPLE { cellMix <- RunPCA(cellMix, features = VariableFeatures(object = cellMix), npcs = npcs) dim <- ifelse(dim >= dim(cellMix@reductions\$pca)[2], dim(cellMix@reductions\$pca)[2], dim) cellMix <- FindNeighbors(cellMix, dims = 1:dim) - cellMix <- FindClusters(cellMix, resolution = $params.resolution, cluster.name = "clusters") + cellMix <- FindClusters(cellMix, resolution = $params.cluster_resolution, cluster.name = "clusters") saveRDS(cellMix, "seurat_obj.rds") diff --git a/modules/seurat/visium_hd/clustering.nf b/modules/seurat/visium_hd/clustering.nf new file mode 100644 index 0000000..1c6f4a2 --- /dev/null +++ b/modules/seurat/visium_hd/clustering.nf @@ -0,0 +1,58 @@ +process SEURAT_VISIUM_HD { + publishDir "$params.output_dir", mode: 'copy', pattern: 'seurat_obj.rds' + label "r" + label "medium_cpu" + label "high_mem" + label "long" + + input: + tuple val(resolution), path(se_gene_counts), path(spot_mappings) + + output: + path ("clusters.rds"), emit: clusters + path ("seurat_obj.rds") + path "versions.yml", topic: 'versions' + + script: + """ + #!/usr/bin/env Rscript + library(Seurat) + library(SummarizedExperiment) + source(Sys.which("visium_hd_cluster.R")) + + banksy <- as.logical("$params.banksy") + assay <- "Spatial.$resolution" + + # the aggregated SE already carries each bin's tissue position in colData + se <- readRDS("$se_gene_counts") + spotMappings <- read.csv("$spot_mappings", stringsAsFactors = FALSE) + object <- CreateSeuratObject(assays(se)\$counts, assay = assay, meta.data = as.data.frame(colData(se))) + + DefaultAssay(object) <- assay + object <- NormalizeData(object, verbose = FALSE) + + # Cluster spatially with Banksy, or on gene expression alone + if (banksy) { + library(SeuratWrappers) + library(Banksy) + object <- clusterBanksy(object, assay, $params.banksy_lambda, $params.banksy_k_geom, $params.cluster_resolution) + } else { + object <- clusterExpression(object, assay, $params.cluster_resolution) + } + + # bambu quantifies at 2um, so map the bin level labels down to the 2um spots + binToClustersMap <- setNames(paste0("cluster_", object\$clusters), names(object\$clusters)) + spotToClustersMap <- mapSpotsToClusters(binToClustersMap, spotMappings) + + # Save the 2um cluster assignment, plus the Seurat object carrying the bin level labels + saveRDS(spotToClustersMap, "clusters.rds") + saveRDS(object, "seurat_obj.rds") + + writeLines(c( + '"${task.process}":', + paste0(' R: ', R.Version()\$version.string), + paste0(' seurat: ', as.character(packageVersion("Seurat"))), + paste0(' banksy: ', as.character(packageVersion("Banksy"))) + ), "versions.yml") + """ +} diff --git a/nextflow.config b/nextflow.config index 5c8bf0a..8099667 100644 --- a/nextflow.config +++ b/nextflow.config @@ -12,33 +12,49 @@ params { input = null // Path to samplesheet .csv file genome = null // Path to .fa or .fasta file annotation = null // Path to .gtf or .gff file - + // Optional: Output directory output_dir = "output" // Path to output directory /* Optional: Samplesheet settings (Non Visium HD samples only) Note: Use this if all samples share the same chemistry/technology - */ - chemistry = null // Examples: "10x3v2", "10x3v3", "10x5v2", "visium-v1" - technology = null // Options: "ONT", "PacBio" - + */ + chemistry = null // Examples: "10x3v2", "10x3v3", "10x5v2", "visium-v1" + technology = null // Options: "ONT", "PacBio" + // Optional: Stop after alignment and save BAM files only bam_only = false // boolean // Optional: Q-score filtering - qscore_filtering = true // boolean + qscore_filtering = true // boolean // Optional: Bambu parameters - ndr = null // null or float + ndr = null // null or float deduplicate_umis = true // boolean // Optional: Quantification mode quantification_mode = "EM_clusters" // Options: "no_quant", "EM", "EM_clusters" // Optional: Seurat clustering - resolution = 0.8 // float - + cluster_resolution = 0.8 // float + + /* + Optional: Visium HD + Note: Set --visium_hd true for a single-sample run starting from a pre-aligned, barcode-tagged BAM file. + The samplesheet only needs 'sample' and 'path' columns — 'chemistry' and 'technology' are not used. + Counts are produced at 2um and at every resolution listed in the bins samplesheet. + */ + visium_hd = false // boolean + bins = null // CSV samplesheet of resolutions (columns: resolution,tissue_positions); must include a 2 row (native base) plus any bins; required when --visium_hd true + barcode_mappings = null // Path to barcode_mappings.parquet (required when --visium_hd true) + + // Optional: Visium HD clustering + clustering_bin = 8 // Integer bin size to cluster at; must be listed in the bins samplesheet + banksy = true // boolean; spatially aware clustering + banksy_lambda = 0.8 // float; 0.8 segments tissue domains, 0.2 gives spatially informed cell types + banksy_k_geom = 50 // integer; spatial neighbours per bin + } includeConfig 'conf/base.config' @@ -50,7 +66,7 @@ profiles { singularity { singularity.enabled = true singularity.autoMounts = true - docker.enabled = false + docker.enabled = false singularity.runOptions = "--bind ${env('PWD')}" singularity.envWhitelist = "JAVA_HOME" } @@ -88,7 +104,7 @@ profiles { includeConfig 'conf/smoke_test.config' -// Output an html timeline report +// Output an html timeline report timeline { enabled = true; file = "${params.output_dir}/pipeline_info/execution_timeline.html" } // Output resource and runtime reports for a workflow run @@ -98,4 +114,4 @@ report { enabled = true; file = "${params.output_dir}/pipeline_info/execution_ trace { enabled = true; file = "${params.output_dir}/pipeline_info/execution_trace.txt" } // Produce a workflow diagram -dag { enabled = true; file = "${params.output_dir}/pipeline_info/pipeline_dag.svg" } \ No newline at end of file +dag { enabled = true; file = "${params.output_dir}/pipeline_info/pipeline_dag.svg" } diff --git a/subworkflows/clustering.nf b/subworkflows/clustering_standard.nf similarity index 71% rename from subworkflows/clustering.nf rename to subworkflows/clustering_standard.nf index 08a3979..8436ed8 100644 --- a/subworkflows/clustering.nf +++ b/subworkflows/clustering_standard.nf @@ -1,5 +1,5 @@ -include { SEURAT_SINGLE_SAMPLE } from '../modules/seurat/single_sample.nf' -include { SEURAT_MULTI_SAMPLE } from '../modules/seurat/multi_sample.nf' +include { SEURAT_SINGLE_SAMPLE } from '../modules/seurat/standard/single_sample_clustering.nf' +include { SEURAT_MULTI_SAMPLE } from '../modules/seurat/standard/multi_sample_clustering.nf' workflow CLUSTERING { take: diff --git a/subworkflows/prepare_input_visium_hd.nf b/subworkflows/prepare_input_visium_hd.nf new file mode 100644 index 0000000..8c47ae3 --- /dev/null +++ b/subworkflows/prepare_input_visium_hd.nf @@ -0,0 +1,53 @@ +include { FILTER_BARCODED_BAM } from '../modules/prepare_input_visium_hd/filter_barcoded_bam.nf' +include { CONVERT_BARCODE_MAPPINGS } from '../modules/prepare_input_visium_hd/convert_barcode_mappings.nf' +include { CONVERT_TISSUE_POSITIONS } from '../modules/prepare_input_visium_hd/convert_tissue_positions.nf' + +workflow PREPARE_INPUT_VISIUM_HD { + take: + ch_rows // raw samplesheet rows + + main: + // Visium HD: single sample, starting from a pre-aligned, barcode-tagged BAM file + ch_rows.collect(flat: false).map { rows -> Validation.validateVisiumHDRows(rows) } + + ch_sample = ch_rows.map { row -> + def sample_path = file(row.path, checkIfExists: true) + def meta = [chemistry: 'visium-hd', technology: 'NA'] + [row.sample, sample_path, meta] + } + + // parse the bins samplesheet into one [resolution, tissue_positions] tuple per row (includes mandatory 2um base) + ch_resolutions = channel.of(params.bins) + .splitCsv(header: true, sep: ',') + .map { row -> + // Convert integer values in the resolution column of the CSV into Spaceranger resolution, e.g. 8 -> "008um" + def resolution = String.format('%03dum', row.resolution as Integer) + def tissue_positions = file(row.tissue_positions, checkIfExists: true) + [resolution, tissue_positions] + } + ch_barcode_mappings = channel.value(params.barcode_mappings) + + // convert every resolution's tissue positions parquet to CSV; the 2um task also extracts the in-tissue barcode list + CONVERT_TISSUE_POSITIONS(ch_resolutions) + + // the 2um in-tissue barcodes are used to filter out-of-tissue reads from the BAM file + FILTER_BARCODED_BAM(ch_sample, CONVERT_TISSUE_POSITIONS.out.barcodes) + + // convert the barcode mappings parquet to one 'barcode,bin' CSV per bin resolution level (e.g., 8um/16um) + ch_bin_resolutions = ch_resolutions + .map { resolution, _tissue_positions -> resolution } + .filter { resolution -> resolution != '002um' } + CONVERT_BARCODE_MAPPINGS(ch_bin_resolutions, ch_barcode_mappings) + + // extract the 2um tissue position separately since transcript discovery is performed on 2um resolution + // first, before aggregating the SE object to the other lower resolutions + ch_tissue_positions_002um = CONVERT_TISSUE_POSITIONS.out.csv.filter { resolution, _csv -> resolution == '002um' }.map { _resolution, csv -> csv } + ch_tissue_positions_bins = CONVERT_TISSUE_POSITIONS.out.csv.filter { resolution, _csv -> resolution != '002um' } + + emit: + bam = FILTER_BARCODED_BAM.out.bam + sample_name = ch_sample.map { sample, _path, _meta -> sample }.first() + tissue_positions_002um = ch_tissue_positions_002um + tissue_positions_bins = ch_tissue_positions_bins + barcode_mappings = CONVERT_BARCODE_MAPPINGS.out.csv +} From 27196faa2315a6c6d1be6613ecbb5ab0bb5871d6 Mon Sep 17 00:00:00 2001 From: ch99l Date: Thu, 30 Jul 2026 10:30:20 +0800 Subject: [PATCH 02/24] update Dockerfile --- containers/r/Dockerfile | 1 - 1 file changed, 1 deletion(-) diff --git a/containers/r/Dockerfile b/containers/r/Dockerfile index 56cf9f1..9759a33 100644 --- a/containers/r/Dockerfile +++ b/containers/r/Dockerfile @@ -16,7 +16,6 @@ RUN micromamba install -y -n base -c conda-forge -c bioconda \ bioconductor-dropletutils=1.30.0 \ r-seurat=5.4.0 \ r-harmony=2.0.2 \ - bioconductor-banksy=1.4.0 \ r-data.table=1.17.8 \ gxx=15.2.0 \ && micromamba clean --all --yes From 937a68380da5a93d10ad57b8bb2c46ed15711bdc Mon Sep 17 00:00:00 2001 From: ch99l Date: Thu, 30 Jul 2026 10:39:35 +0800 Subject: [PATCH 03/24] update Dockerfile --- containers/r/Dockerfile | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/containers/r/Dockerfile b/containers/r/Dockerfile index 9759a33..7a248e3 100644 --- a/containers/r/Dockerfile +++ b/containers/r/Dockerfile @@ -11,6 +11,7 @@ ARG MAMBA_DOCKERFILE_ACTIVATE=1 RUN micromamba install -y -n base -c conda-forge -c bioconda \ r-base=4.5.3 \ r-devtools=2.5.2 \ + r-remotes=2.5.0 \ r-biocmanager=1.30.27 \ bioconductor-bambu=3.12.1 \ bioconductor-dropletutils=1.30.0 \ @@ -20,8 +21,8 @@ RUN micromamba install -y -n base -c conda-forge -c bioconda \ gxx=15.2.0 \ && micromamba clean --all --yes -# SeuratWrappers provides RunBanksy() and is not packaged on conda -RUN R -e 'devtools::install_github("satijalab/seurat-wrappers", upgrade = "never")' +# Install Seurat Wrappers +RUN R -e 'remotes::install_github("satijalab/seurat-wrappers", upgrade = "never")' # Clone bambu single cell feature branch RUN git clone --branch devel_pre_v4 \ From e95955345783941b4d3bd94f0b908df4729cdde5 Mon Sep 17 00:00:00 2001 From: ch99l Date: Thu, 30 Jul 2026 11:05:57 +0800 Subject: [PATCH 04/24] update file structure --- main.nf | 2 +- .../standard}/extract_barcodes.nf | 0 .../standard}/extract_spatial_coordinates.nf | 0 .../visium_hd}/convert_barcode_mappings.nf | 0 .../visium_hd}/convert_tissue_positions.nf | 0 .../visium_hd}/filter_barcoded_bam.nf | 0 .../visium_hd}/spot_bin_mappings.nf | 0 subworkflows/prepare_input_standard.nf | 4 ++-- subworkflows/prepare_input_visium_hd.nf | 6 +++--- 9 files changed, 6 insertions(+), 6 deletions(-) rename modules/{prepare_input_standard => prepare_input/standard}/extract_barcodes.nf (100%) rename modules/{prepare_input_standard => prepare_input/standard}/extract_spatial_coordinates.nf (100%) rename modules/{prepare_input_visium_hd => prepare_input/visium_hd}/convert_barcode_mappings.nf (100%) rename modules/{prepare_input_visium_hd => prepare_input/visium_hd}/convert_tissue_positions.nf (100%) rename modules/{prepare_input_visium_hd => prepare_input/visium_hd}/filter_barcoded_bam.nf (100%) rename modules/{prepare_input_visium_hd => prepare_input/visium_hd}/spot_bin_mappings.nf (100%) diff --git a/main.nf b/main.nf index afa53c5..f895ebf 100644 --- a/main.nf +++ b/main.nf @@ -20,7 +20,7 @@ include { BAMBU_EM } from './modules/bambu/standard include { BAMBU_TRANSCRIPT_DISCOVERY_VISIUM_HD } from './modules/bambu/visium_hd/transcript_discovery.nf' include { BAMBU_EM_VISIUM_HD } from './modules/bambu/visium_hd/spot_level_quantification.nf' include { AGGREGATE_BINS_VISIUM_HD } from './modules/bambu/visium_hd/aggregate_bins.nf' -include { SPOT_BIN_MAPPINGS } from './modules/prepare_input_visium_hd/spot_bin_mappings.nf' +include { SPOT_BIN_MAPPINGS } from './modules/prepare_input/visium_hd/spot_bin_mappings.nf' include { SEURAT_VISIUM_HD } from './modules/seurat/visium_hd/clustering.nf' params { diff --git a/modules/prepare_input_standard/extract_barcodes.nf b/modules/prepare_input/standard/extract_barcodes.nf similarity index 100% rename from modules/prepare_input_standard/extract_barcodes.nf rename to modules/prepare_input/standard/extract_barcodes.nf diff --git a/modules/prepare_input_standard/extract_spatial_coordinates.nf b/modules/prepare_input/standard/extract_spatial_coordinates.nf similarity index 100% rename from modules/prepare_input_standard/extract_spatial_coordinates.nf rename to modules/prepare_input/standard/extract_spatial_coordinates.nf diff --git a/modules/prepare_input_visium_hd/convert_barcode_mappings.nf b/modules/prepare_input/visium_hd/convert_barcode_mappings.nf similarity index 100% rename from modules/prepare_input_visium_hd/convert_barcode_mappings.nf rename to modules/prepare_input/visium_hd/convert_barcode_mappings.nf diff --git a/modules/prepare_input_visium_hd/convert_tissue_positions.nf b/modules/prepare_input/visium_hd/convert_tissue_positions.nf similarity index 100% rename from modules/prepare_input_visium_hd/convert_tissue_positions.nf rename to modules/prepare_input/visium_hd/convert_tissue_positions.nf diff --git a/modules/prepare_input_visium_hd/filter_barcoded_bam.nf b/modules/prepare_input/visium_hd/filter_barcoded_bam.nf similarity index 100% rename from modules/prepare_input_visium_hd/filter_barcoded_bam.nf rename to modules/prepare_input/visium_hd/filter_barcoded_bam.nf diff --git a/modules/prepare_input_visium_hd/spot_bin_mappings.nf b/modules/prepare_input/visium_hd/spot_bin_mappings.nf similarity index 100% rename from modules/prepare_input_visium_hd/spot_bin_mappings.nf rename to modules/prepare_input/visium_hd/spot_bin_mappings.nf diff --git a/subworkflows/prepare_input_standard.nf b/subworkflows/prepare_input_standard.nf index 4275965..8e01e07 100644 --- a/subworkflows/prepare_input_standard.nf +++ b/subworkflows/prepare_input_standard.nf @@ -1,5 +1,5 @@ -include { EXTRACT_10X_BARCODES } from '../modules/prepare_input_standard/extract_barcodes.nf' -include { EXTRACT_10X_SPATIAL_COORDINATES } from '../modules/prepare_input_standard/extract_spatial_coordinates.nf' +include { EXTRACT_10X_BARCODES } from '../modules/prepare_input/standard/extract_barcodes.nf' +include { EXTRACT_10X_SPATIAL_COORDINATES } from '../modules/prepare_input/standard/extract_spatial_coordinates.nf' workflow PREPARE_INPUT_STANDARD { take: diff --git a/subworkflows/prepare_input_visium_hd.nf b/subworkflows/prepare_input_visium_hd.nf index 8c47ae3..4866c7c 100644 --- a/subworkflows/prepare_input_visium_hd.nf +++ b/subworkflows/prepare_input_visium_hd.nf @@ -1,6 +1,6 @@ -include { FILTER_BARCODED_BAM } from '../modules/prepare_input_visium_hd/filter_barcoded_bam.nf' -include { CONVERT_BARCODE_MAPPINGS } from '../modules/prepare_input_visium_hd/convert_barcode_mappings.nf' -include { CONVERT_TISSUE_POSITIONS } from '../modules/prepare_input_visium_hd/convert_tissue_positions.nf' +include { FILTER_BARCODED_BAM } from '../modules/prepare_input/visium_hd/filter_barcoded_bam.nf' +include { CONVERT_BARCODE_MAPPINGS } from '../modules/prepare_input/visium_hd/convert_barcode_mappings.nf' +include { CONVERT_TISSUE_POSITIONS } from '../modules/prepare_input/visium_hd/convert_tissue_positions.nf' workflow PREPARE_INPUT_VISIUM_HD { take: From b13afa5b420a357be2512c5c287b56209ed513d5 Mon Sep 17 00:00:00 2001 From: ch99l Date: Thu, 30 Jul 2026 14:29:34 +0800 Subject: [PATCH 05/24] bug fix and refactor --- main.nf | 50 +++++++++---------- ...ion.nf => cluster_level_quantification.nf} | 2 +- .../standard/single_cell_quantification.nf | 2 +- .../visium_hd/spot_level_quantification.nf | 2 +- 4 files changed, 27 insertions(+), 29 deletions(-) rename modules/bambu/shared/{clustered_quantification.nf => cluster_level_quantification.nf} (97%) diff --git a/main.nf b/main.nf index f895ebf..48ce15d 100644 --- a/main.nf +++ b/main.nf @@ -3,25 +3,25 @@ nextflow.enable.types = true // subworkflows -include { PREPARE_INPUT_STANDARD } from './subworkflows/prepare_input_standard.nf' -include { PREPARE_INPUT_VISIUM_HD } from './subworkflows/prepare_input_visium_hd.nf' -include { ALIGNMENT } from './subworkflows/alignment.nf' -include { CLUSTERING } from './subworkflows/clustering_standard.nf' +include { PREPARE_INPUT_STANDARD } from './subworkflows/prepare_input_standard.nf' +include { PREPARE_INPUT_VISIUM_HD } from './subworkflows/prepare_input_visium_hd.nf' +include { ALIGNMENT } from './subworkflows/alignment.nf' +include { CLUSTERING } from './subworkflows/clustering_standard.nf' // modules -include { DECOMPRESS as DECOMPRESS_GENOME } from './modules/decompress.nf' -include { DECOMPRESS as DECOMPRESS_ANNOTATION } from './modules/decompress.nf' -include { PREPROCESS_FASTQ } from './modules/preprocess_fastq.nf' -include { BAMBU_PREPARE_ANNOTATION } from './modules/bambu/shared/prepare_annotation.nf' -include { BAMBU_CONSTRUCT_READ_CLASS } from './modules/bambu/shared/construct_read_class.nf' -include { BAMBU_TRANSCRIPT_DISCOVERY } from './modules/bambu/standard/transcript_discovery.nf' -include { BAMBU_CLUSTERED_EM } from './modules/bambu/shared/clustered_quantification.nf' -include { BAMBU_EM } from './modules/bambu/standard/single_cell_quantification.nf' -include { BAMBU_TRANSCRIPT_DISCOVERY_VISIUM_HD } from './modules/bambu/visium_hd/transcript_discovery.nf' -include { BAMBU_EM_VISIUM_HD } from './modules/bambu/visium_hd/spot_level_quantification.nf' -include { AGGREGATE_BINS_VISIUM_HD } from './modules/bambu/visium_hd/aggregate_bins.nf' -include { SPOT_BIN_MAPPINGS } from './modules/prepare_input/visium_hd/spot_bin_mappings.nf' -include { SEURAT_VISIUM_HD } from './modules/seurat/visium_hd/clustering.nf' +include { DECOMPRESS as DECOMPRESS_GENOME } from './modules/decompress.nf' +include { DECOMPRESS as DECOMPRESS_ANNOTATION } from './modules/decompress.nf' +include { PREPROCESS_FASTQ } from './modules/preprocess_fastq.nf' +include { BAMBU_PREPARE_ANNOTATION } from './modules/bambu/shared/prepare_annotation.nf' +include { BAMBU_CONSTRUCT_READ_CLASS } from './modules/bambu/shared/construct_read_class.nf' +include { BAMBU_CLUSTER_LEVEL_QUANTIFICATION } from './modules/bambu/shared/cluster_level_quantification.nf' +include { BAMBU_TRANSCRIPT_DISCOVERY } from './modules/bambu/standard/transcript_discovery.nf' +include { BAMBU_SINGLE_CELL_QUANTIFICATION } from './modules/bambu/standard/single_cell_quantification.nf' +include { BAMBU_TRANSCRIPT_DISCOVERY_VISIUM_HD } from './modules/bambu/visium_hd/transcript_discovery.nf' +include { BAMBU_SPOT_LEVEL_QUANTIFICATION } from './modules/bambu/visium_hd/spot_level_quantification.nf' +include { AGGREGATE_BINS_VISIUM_HD } from './modules/bambu/visium_hd/aggregate_bins.nf' +include { SPOT_BIN_MAPPINGS } from './modules/prepare_input/visium_hd/spot_bin_mappings.nf' +include { SEURAT_VISIUM_HD } from './modules/seurat/visium_hd/clustering.nf' params { input: Path @@ -93,9 +93,9 @@ workflow STANDARD { // cluster the cells first, then pool each cluster's cells for the EM if (params.quantification_mode == 'EM_clusters') { CLUSTERING(BAMBU_TRANSCRIPT_DISCOVERY.out.se_gene_counts, ch_n_samples) - BAMBU_CLUSTERED_EM(CLUSTERING.out.clusters, BAMBU_TRANSCRIPT_DISCOVERY.out.quant_data, BAMBU_TRANSCRIPT_DISCOVERY.out.extended_annotations, ch_genome) + BAMBU_CLUSTER_LEVEL_QUANTIFICATION(CLUSTERING.out.clusters, BAMBU_TRANSCRIPT_DISCOVERY.out.quant_data, BAMBU_TRANSCRIPT_DISCOVERY.out.extended_annotations, ch_genome) } else if (params.quantification_mode == 'EM') { - BAMBU_EM(BAMBU_TRANSCRIPT_DISCOVERY.out.quant_data, BAMBU_TRANSCRIPT_DISCOVERY.out.extended_annotations, ch_genome) + BAMBU_SINGLE_CELL_QUANTIFICATION(BAMBU_TRANSCRIPT_DISCOVERY.out.quant_data, BAMBU_TRANSCRIPT_DISCOVERY.out.extended_annotations, ch_genome) } } } @@ -132,21 +132,19 @@ workflow VISIUM_HD { if (params.quantification_mode == 'EM_clusters') { def requested_bin = String.format('%03dum', params.clustering_bin) // convert clustering_bin specified as an integer into Spaceranger format - // perform clustering at the requested resolution only ch_clustering = AGGREGATE_BINS_VISIUM_HD.out.se_gene_counts .filter { resolution, _se_gene_counts -> resolution == requested_bin } .join(SPOT_BIN_MAPPINGS.out.csv) // [resolution, se_gene_counts, spot_mappings] SEURAT_VISIUM_HD(ch_clustering) - BAMBU_CLUSTERED_EM(SEURAT_VISIUM_HD.out.clusters, ch_quant_data, ch_extended_anno, ch_genome) - } + BAMBU_CLUSTER_LEVEL_QUANTIFICATION(SEURAT_VISIUM_HD.out.clusters, ch_quant_data, ch_extended_anno, ch_genome) - if (params.quantification_mode != 'no_quant') { + } else if (params.quantification_mode == 'EM') { // run spot level quantification on all resolution // at 2um resolution, tissue_positions and spot_mappings are not required - ch_resolution_002um = channel.of(['002um', [], []]) - ch_resolutions = ch_resolution_002um.mix(ch_bins) // [resolution, tissue_positions, spot_mappings] - BAMBU_EM_VISIUM_HD(ch_resolutions, ch_quant_data, ch_extended_anno, ch_genome) + ch_resolution_002um = channel.of(['002um', [], []]) + ch_resolutions = ch_resolution_002um.concat(ch_bins) // [resolution, tissue_positions, spot_mappings] + BAMBU_SPOT_LEVEL_QUANTIFICATION(ch_resolutions, ch_quant_data, ch_extended_anno, ch_genome) } } diff --git a/modules/bambu/shared/clustered_quantification.nf b/modules/bambu/shared/cluster_level_quantification.nf similarity index 97% rename from modules/bambu/shared/clustered_quantification.nf rename to modules/bambu/shared/cluster_level_quantification.nf index f907036..a77c65f 100644 --- a/modules/bambu/shared/clustered_quantification.nf +++ b/modules/bambu/shared/cluster_level_quantification.nf @@ -1,4 +1,4 @@ -process BAMBU_CLUSTERED_EM { +process BAMBU_CLUSTER_LEVEL_QUANTIFICATION { publishDir "$params.output_dir", mode: 'copy', pattern: 'transcript_counts_clusters' publishDir "$params.output_dir", mode: 'copy', pattern: 'gene_counts_clusters' label "r" diff --git a/modules/bambu/standard/single_cell_quantification.nf b/modules/bambu/standard/single_cell_quantification.nf index 97df050..7a713f0 100644 --- a/modules/bambu/standard/single_cell_quantification.nf +++ b/modules/bambu/standard/single_cell_quantification.nf @@ -1,4 +1,4 @@ -process BAMBU_EM { +process BAMBU_SINGLE_CELL_QUANTIFICATION { publishDir "$params.output_dir", mode: 'copy', pattern: 'transcript_counts_singlecell' label "r" label "low_cpu" diff --git a/modules/bambu/visium_hd/spot_level_quantification.nf b/modules/bambu/visium_hd/spot_level_quantification.nf index 52b5023..9d2cb8a 100644 --- a/modules/bambu/visium_hd/spot_level_quantification.nf +++ b/modules/bambu/visium_hd/spot_level_quantification.nf @@ -1,4 +1,4 @@ -process BAMBU_EM_VISIUM_HD { +process BAMBU_SPOT_LEVEL_QUANTIFICATION { publishDir "$params.output_dir", mode: 'copy', pattern: 'transcript_counts_*' label "r" label "low_cpu" From 027aaf092d95b4aa403c5397c7ed028f6df368f9 Mon Sep 17 00:00:00 2001 From: ch99l Date: Thu, 30 Jul 2026 15:40:38 +0800 Subject: [PATCH 06/24] fix Dockerfile --- containers/r/Dockerfile | 2 ++ 1 file changed, 2 insertions(+) diff --git a/containers/r/Dockerfile b/containers/r/Dockerfile index 7a248e3..ddbd039 100644 --- a/containers/r/Dockerfile +++ b/containers/r/Dockerfile @@ -18,6 +18,8 @@ RUN micromamba install -y -n base -c conda-forge -c bioconda \ r-seurat=5.4.0 \ r-harmony=2.0.2 \ r-data.table=1.17.8 \ + r-magick \ + r-leidenalg \ gxx=15.2.0 \ && micromamba clean --all --yes From 781caa890d3969d6f3c8ef7b58adb15f718d22ab Mon Sep 17 00:00:00 2001 From: ch99l Date: Thu, 30 Jul 2026 17:14:35 +0800 Subject: [PATCH 07/24] update --- bin/visium_hd_aggregate_resolution.R | 22 ++++++++++++++----- modules/bambu/visium_hd/aggregate_bins.nf | 4 ++-- .../visium_hd/spot_level_quantification.nf | 2 +- .../bambu/visium_hd/transcript_discovery.nf | 4 ++-- 4 files changed, 21 insertions(+), 11 deletions(-) diff --git a/bin/visium_hd_aggregate_resolution.R b/bin/visium_hd_aggregate_resolution.R index 5dc5c7a..753d846 100755 --- a/bin/visium_hd_aggregate_resolution.R +++ b/bin/visium_hd_aggregate_resolution.R @@ -20,15 +20,25 @@ aggregateResolution <- function(se, tissuePositions, spotMappings) { spotCounts <- assays(se)$counts binCounts <- as(spotCounts %*% spotsPerBin, "CsparseMatrix") - # build new colData for the aggregated SE object - barcodePerBin <- distinct(spotMappings, sample_bin, barcode = bin) - positionPerBin <- data.frame(sample_bin = binIds) %>% - left_join(barcodePerBin, by = "sample_bin") %>% - left_join(tissuePositions, by = "barcode") %>% + ## build new colData for the aggregated SE object + + # bin id -> barcode (without sampleName prefix), e.g. sample_s_008um_00247_00090-1 -> s_008um_00247_00090-1 + binIdToBarcodeMap <- distinct(spotMappings, sample_bin, barcode = bin) + + # bin id -> the 2um barcodes the bin contains + binToSpotsMap <- data.frame(sample_bin = as.character(spotToBinMap), barcode = colData(se)$barcode) %>% + group_by(sample_bin) %>% + summarise(barcodes = list(barcode), .groups = "drop") + + # creates dataframe containing one row per bin, each row contains metadata information (e.g., barcode, in_tissue, array and pixel coordinates) + metadataPerBin <- data.frame(sample_bin = binIds) %>% + left_join(binIdToBarcodeMap, by = "sample_bin") %>% + left_join(binToSpotsMap, by = "sample_bin") %>% + left_join(tissuePositions, by = "barcode") %>% select(-sample_bin) sampleName <- unique(colData(se)$sampleName) - binColData <- DataFrame(id = binIds, sampleName = sampleName, positionPerBin, row.names = binIds) + binColData <- DataFrame(id = binIds, sampleName = sampleName, metadataPerBin, row.names = binIds) binSe <- SummarizedExperiment(assays = SimpleList(counts = binCounts), rowRanges = rowRanges(se), diff --git a/modules/bambu/visium_hd/aggregate_bins.nf b/modules/bambu/visium_hd/aggregate_bins.nf index 1996fe9..6042d62 100644 --- a/modules/bambu/visium_hd/aggregate_bins.nf +++ b/modules/bambu/visium_hd/aggregate_bins.nf @@ -1,6 +1,6 @@ process AGGREGATE_BINS_VISIUM_HD { - publishDir "$params.output_dir", mode: 'copy', pattern: 'unique_counts_*' - publishDir "$params.output_dir", mode: 'copy', pattern: 'gene_counts_*' + publishDir "$params.output_dir/unique_counts", mode: 'copy', pattern: 'unique_counts_*' + publishDir "$params.output_dir/gene_counts", mode: 'copy', pattern: 'gene_counts_*' label "r" label "low_cpu" label "medium_mem" diff --git a/modules/bambu/visium_hd/spot_level_quantification.nf b/modules/bambu/visium_hd/spot_level_quantification.nf index 9d2cb8a..6050c1b 100644 --- a/modules/bambu/visium_hd/spot_level_quantification.nf +++ b/modules/bambu/visium_hd/spot_level_quantification.nf @@ -1,5 +1,5 @@ process BAMBU_SPOT_LEVEL_QUANTIFICATION { - publishDir "$params.output_dir", mode: 'copy', pattern: 'transcript_counts_*' + publishDir "$params.output_dir/transcript_counts", mode: 'copy', pattern: 'transcript_counts_*' label "r" label "low_cpu" label "high_mem" diff --git a/modules/bambu/visium_hd/transcript_discovery.nf b/modules/bambu/visium_hd/transcript_discovery.nf index b096507..b3b596e 100644 --- a/modules/bambu/visium_hd/transcript_discovery.nf +++ b/modules/bambu/visium_hd/transcript_discovery.nf @@ -1,7 +1,7 @@ process BAMBU_TRANSCRIPT_DISCOVERY_VISIUM_HD { publishDir "$params.output_dir", mode: 'copy', pattern: 'extended_annotations.gtf' - publishDir "$params.output_dir", mode: 'copy', pattern: 'unique_counts_002um' - publishDir "$params.output_dir", mode: 'copy', pattern: 'gene_counts_002um' + publishDir "$params.output_dir/unique_counts", mode: 'copy', pattern: 'unique_counts_002um' + publishDir "$params.output_dir/gene_counts", mode: 'copy', pattern: 'gene_counts_002um' publishDir "$params.output_dir/intermediate_R", mode: 'copy', pattern: '*.rds', enabled: params.save_intermediates label "r" label "medium_cpu" From eae76cc8dd0786b91b54e91f6b4b3a1321528379 Mon Sep 17 00:00:00 2001 From: Chin Hao <60098604+ch99l@users.noreply.github.com> Date: Thu, 30 Jul 2026 21:01:42 +0800 Subject: [PATCH 08/24] docker refactor (#39) --- .github/workflows/build_container.yml | 158 +++++++++++------- bin/create_seurat_object.R | 14 ++ conf/base.config | 2 +- conf/containers.config | 3 +- containers/{r => bambu}/Dockerfile | 14 +- containers/seurat/Dockerfile | 13 ++ main.nf | 9 +- .../shared/cluster_level_quantification.nf | 2 +- modules/bambu/shared/construct_read_class.nf | 2 +- modules/bambu/shared/prepare_annotation.nf | 2 +- .../standard/single_cell_quantification.nf | 2 +- .../bambu/standard/transcript_discovery.nf | 10 +- modules/bambu/visium_hd/aggregate_bins.nf | 8 +- .../visium_hd/spot_level_quantification.nf | 2 +- .../bambu/visium_hd/transcript_discovery.nf | 2 +- .../standard/multi_sample_clustering.nf | 29 ++-- .../standard/single_sample_clustering.nf | 22 +-- modules/seurat/visium_hd/clustering.nf | 15 +- subworkflows/clustering_standard.nf | 12 +- 19 files changed, 194 insertions(+), 127 deletions(-) create mode 100755 bin/create_seurat_object.R rename containers/{r => bambu}/Dockerfile (71%) create mode 100644 containers/seurat/Dockerfile diff --git a/.github/workflows/build_container.yml b/.github/workflows/build_container.yml index 1b97a6a..6907417 100644 --- a/.github/workflows/build_container.yml +++ b/.github/workflows/build_container.yml @@ -1,13 +1,16 @@ -name: Build R container +name: Build R containers -# The image tag is read from conf/containers.config, so that single line is both what the -# pipeline pulls and what this workflow publishes -- the two cannot drift. -# Runs on pushes to a feature branch, whenever the container changes. +# The image tag of the container is read from conf/containers.config, so that the tag is +# what the pipeline pulls and what this workflow publishes. Each container is built independently, +# so updating one tag in conf/containers.config rebuilds only that image: +# feature branch - rebuilds its tag on every push, unless main or devel is pinned to that tag +# main or devel - builds only a tag that has never been published, so most merges are a no-op +# manual run - 'force' rebuilds and overwrites the tag on any branch on: push: paths: - - 'containers/r/**' + - 'containers/**' - 'conf/containers.config' - '.github/workflows/build_container.yml' workflow_dispatch: @@ -22,81 +25,116 @@ concurrency: cancel-in-progress: true jobs: - build: - name: Build and publish + discover: + name: Select containers to build runs-on: ubuntu-latest permissions: contents: read - packages: write + packages: read + outputs: + containers: ${{ steps.select.outputs.containers }} steps: - uses: actions/checkout@v4 - - name: Resolve image from conf/containers.config - id: image - run: | - IMAGE=$(sed -n "s|.*withLabel: *'r'.*container *= *\"\([^\"]*\)\".*|\1|p" conf/containers.config) - if [ -z "$IMAGE" ]; then - echo "::error::could not parse the 'r' container from conf/containers.config" - exit 1 - fi - case "$IMAGE" in - ghcr.io/goekelab/*) ;; - *) echo "::error::refusing to push outside ghcr.io/goekelab: $IMAGE"; exit 1 ;; - esac - echo "image=$IMAGE" >> "$GITHUB_OUTPUT" - echo "Resolved container: $IMAGE" - - - name: Refuse to overwrite a tag main or devel is pinned to - if: github.ref_name != 'main' && github.ref_name != 'devel' + - uses: docker/login-action@v3 + with: + registry: ghcr.io + username: ${{ github.actor }} + password: ${{ secrets.GITHUB_TOKEN }} + + - name: Select containers to build + id: select run: | - IMAGE="${{ steps.image.outputs.image }}" - for BASE in main devel; do - git fetch --no-tags --depth=1 origin "$BASE" 2>/dev/null || continue - PINNED=$(git show "FETCH_HEAD:conf/containers.config" 2>/dev/null \ - | sed -n "s|.*withLabel: *'r'.*container *= *\"\([^\"]*\)\".*|\1|p") - if [ "$IMAGE" = "$PINNED" ]; then - echo "::error::$BASE is pinned to $IMAGE -- bump the tag in conf/containers.config before rebuilding it" + RELEASED=false + case "${{ github.ref_name }}" in main|devel) RELEASED=true ;; esac + + # Each base branch is fetched into its own ref so both stay readable in the loop + BASES="" + for BASE in devel main; do + if git fetch --no-tags --depth=1 origin "$BASE:refs/base/$BASE" 2>/dev/null; then + BASES="$BASES refs/base/$BASE" + fi + done + + # devel is at or ahead of main, so it is the baseline for "did this container change" + DIFF_BASE=refs/base/devel + + # Every containers/