version 1.0 ########################################################################################## # LRAA Single-Cell RNA-seq Workflow ########################################################################################## # # This workflow performs isoform discovery and/or quantification on single-cell long-read # RNA-seq data with two primary execution modes and multiple optimization paths. # # ===== EXECUTION MODES ===== # # 1. BASIC MODE (single_cell_pipe_mode = "basic") # - Runs initial discovery on the full BAM to generate a preliminary GTF # - Builds single-cell sparse matrices (gene, isoform, splice pattern) from read tracking # - Optionally filters good cells from empty droplets # - Clusters cells using Seurat on gene-level expression # - Outputs: initial GTF, single-cell matrices, cluster assignments, UMAP visualizations # - Use when: You want a quick initial analysis without per-cluster refinement # # 2. CLUSTER-GUIDED MODE (single_cell_pipe_mode = "cluster-guided") # - Performs all basic mode steps PLUS: # - Partitions BAM by cell clusters # - Re-runs LRAA discovery independently on each cluster's reads # - Merges cluster-specific GTFs into a refined final GTF # - Re-quantifies all cells against the final GTF # - Builds final single-cell sparse matrices with improved isoform definitions # - Outputs: refined final GTF, improved single-cell matrices, per-cluster pseudobulk data # - Use when: You want the highest quality isoform definitions tailored to cell types # # ===== QUANTIFICATION-ONLY MODE ===== # # Set quant_only = true to skip isoform discovery and quantify against a known annotation: # - In BASIC mode with quant_only = true: # * Requires: initial_annot_gtf (the annotation to quantify against) # * Skips: Initial discovery phase # * Runs: Quantification, matrix building, cell filtering, clustering # * Outputs: Single-cell matrices and clusters based on the provided annotation # # - In CLUSTER-GUIDED mode with quant_only = true: # * Requires: initial_annot_gtf # * Skips: Both initial discovery AND per-cluster discovery # * Runs: Clustering (if not precomputed), per-cluster quantification only # * Outputs: Single-cell matrices from cluster-guided quantification # # ===== PRECOMPUTED INPUTS (OPTIMIZATION PATHS) ===== # # You can skip expensive computation steps by providing precomputed files from prior runs: # # 1. Skip Initial Discovery: # - Provide: precomputed_init_gtf + precomputed_init_quant_tracking # - Effect: Skips LRAA_init, uses your GTF and tracking files directly # - Use case: Re-run clustering with different parameters without re-discovering isoforms # # 2. Skip Initial Discovery + Clustering: # - Provide: precomputed_init_gtf + precomputed_cluster_assignments_tsv # - Effect: Jumps directly to cluster-guided refinement (if mode = "cluster-guided") # or just outputs the precomputed results (if mode = "basic") # - Use case: Test cluster-guided refinement with different discovery parameters # - NOTE: this is the ONE configuration that emits no init_sc_* matrices, because # no init tracking file exists to build them from. Add # precomputed_init_quant_tracking to get them. # # 3. Kickstart Cluster-Guided Mode: # - Provide: precomputed_init_gtf + precomputed_cluster_assignments_tsv # - Set: single_cell_pipe_mode = "cluster-guided" # - Flow: Uses precomputed GTF as the annotation for cluster-guided discovery # (passed as annot_gtf to guide per-cluster LRAA runs) # - Use case: You have initial results and want refined cluster-specific isoforms # without re-running the expensive initial discovery # # 4. Cluster-Guided Quantification-Only: # - Provide: precomputed_cluster_assignments_tsv + initial_annot_gtf # - Set: single_cell_pipe_mode = "cluster-guided", quant_only = true # - Flow: Skips all discovery, quantifies each cluster against initial_annot_gtf # - Use case: Quantify known isoforms with cluster-aware quantification # # ===== GENE SYMBOL INCORPORATION ===== # # Optionally provide ref_annot_gtf_source_gene_symbols to incorporate gene symbols: # - Runs gffcompare between LRAA GTF and reference annotation # - Updates all single-cell matrices with gene symbols from matched loci # - Only runs when: All required sparse matrices are available (from either basic or # cluster-guided mode, and not skipped due to missing precomputed inputs) # # ===== TYPICAL WORKFLOWS ===== # # Quick exploration: # single_cell_pipe_mode = "basic" # → Get initial GTF, matrices, and clusters quickly # # High-quality production run: # single_cell_pipe_mode = "cluster-guided" # → Full pipeline with cluster-specific isoform refinement # # Iterative refinement: # 1st run: mode = "basic" (save init GTF and cluster assignments) # 2nd run: mode = "cluster-guided" + precomputed files # → Skip expensive initial phase, focus on cluster-specific refinement # # Quantify known annotation: # quant_only = true + initial_annot_gtf # → Fast quantification without discovery # ########################################################################################## import "LRAA.wdl" as LRAA import "subwdls/LRAA-build_sparse_matrices_from_tracking.wdl" as BuildMatrices import "subwdls/LRAA-filter_good_cells.wdl" as FilterCells import "subwdls/LRAA-gene_sparseM_to_seurat_clusters.wdl" as Seurat import "subwdls/Incorporate_gene_symbols.wdl" as GeneSymbols import "LRAA-cell_cluster_guided.wdl" as ClusterGuided # Imported for ONE task, emit_shared_chunk_plan: the run-wide cut plan is selected # here, before the initial pass, and threaded into every phase. The task lives beside # LRAA_quant_by_cluster's comparability argument rather than being duplicated here. import "LRAA_quant_by_cluster.wdl" as QuantByCluster workflow LRAA_singlecell_wf { input { # Core inputs String sample_id File referenceGenome File inputBAM # Optional: guide initial discovery with an annotation File? initial_annot_gtf # Single-cell pipeline mode: 'basic' or 'cluster-guided' String single_cell_pipe_mode = "basic" # Quantification-only mode: skip discovery, requires initial_annot_gtf Boolean quant_only = false # Platform/options Boolean HiFi = false String oversimplify = "chrM" # e.g., "chrM" or "chrM,M" # --polyA_known: trusted cleavage-site BED/GTF (v0.43.0) forwarded to BOTH the init # pass and the per-cluster discovery runs so they waive the internal-priming veto at # the same trusted 3' ends. Default unset => the veto applies to every A-rich read # candidate (a future run may point this at PolyASite 2.0 or the init phase's bed). File? polyA_known String main_chromosomes = "" # if empty, runs without partitioning # Contigs the ONE initial full-library partition extracts concurrently. Defaulted # above 1 because in a single-cell run this step is dead time at the head of the # pipeline on an otherwise idle box: REPORTED on a 188 GB library (1.51 B mapped # reads, 25 contigs, 28-core host), partition_by_chromosome_task ran 27+ minutes # with ONE core busy before any shard started. # # MEASURED on 58.3 M mapped reads across 12 contigs at -@ 4: wall 2:02 at 1 worker, # 1:13.8 at 2, 0:55.8 at 4, 0:56.5 at 6. The KNEE is 4, but the default is 2, on # cost: cpu is derived as workers * (samtools -@ + 1) + 2, so 2 asks for 12 and 4 # asks for 22. Twelve fits a 16-core machine, while 22 forces the next instance # size up -- a 32-core VM billed for the whole task to save the 18 s between # 1:13.8 and 0:55.8. Two captures 1.65x of the available 2.19x for half the # reservation, which is the better trade at every scale measured. # # Raise it to 4 only where the cores are already paid for and idle. Lower it on a # host with fewer than ~14 usable cores, where a 12-core reservation would wait. Int initial_partition_workers = 2 # Contigs each PER-CLUSTER partition extracts concurrently, and the answer differs # by backend, which is why it is separate from initial_partition_workers above. # # 1 is right for LOCAL execution: miniwdl runs the 14-to-32 per-cluster jobs on ONE # host, so widening each one multiplies its reservation by the cluster count and # none of them can be placed until the box drains. That is the regression that cut # this task's cpu to 5 (Partition_data_by_chromosome.wdl:16-21). # # On TERRA each per-cluster job is its own VM, so they are not competing for cores # at all. It stays 1 anyway BY DEFAULT, and this is where the two knobs genuinely # part company: the initial partition is one task whose widening is paid once, # while this one is billed 14 to 32 times over. At 2 that is 14-32 VMs at cpu 12 # instead of 7 to shave seconds off jobs that already run in seconds. Raise it only # if a per-cluster partition is measured as a real share of wall time. Int cluster_partition_workers = 1 String? region # e.g., "chr1:100000-200000"; forces direct mode Boolean rescue_unassigned_reads_via_transcriptome_alignment = true # Turn off 3' end weighting in the EM. This is the outermost entry point, so it # has to reach BOTH the init discovery run below and the cluster-guided # subworkflow, which carries it on to its own two phases. Every layer between # here and LRAA's command line already declared it; only the entry points did # not, which is what made it unreachable from Terra on single-cell runs. Boolean no_weight_reads_by_3prime_agreement = false # The TPM denominator for the INIT pass, and only for it. Unset -- the normal # case -- means LRAA counts the input bam itself, once, with its -F 0x904 # mapped-primary policy. # # Deliberately NOT forwarded to the cluster-guided subworkflow below. Once the # bam is partitioned by cluster, each cluster's denominator is that cluster's own # read total, so each per-cluster LRAA counts the bam it was handed. Passing this # whole-library number down would rescale every cluster's TPM by that cluster's # share of the library, and the result would still look like TPM. Int? num_total_reads # Chunk geometry. Approximate MEGABASES between cut points, and the TOTAL width # of the search window centred on each target cut. LRAA's defaults are 10 and 1, # so a contig shorter than 10 Mb is never cut and chunking degenerates to one # chunk holding the whole thing -- which is what a small test fixture does, # leaving each run only its two strand units to parallelise over. # # The init and the per-cluster runs would each want their own geometry, and no # longer CAN have one: ONE cut plan is emitted below, before any phase runs, and # every phase applies it, so a run has exactly one set of cut positions. Two # geometries would mean a phase whose chunk bounds differ from the plan's, which # ChunkedRun refuses by name and value rather than quietly cutting elsewhere. # # So the *_init pair sets the geometry for the WHOLE run when given, and the plain # pair is the fallback -- the same precedence the init call used to apply to itself. # The init wins because it is the phase that cannot be parallelised any other way: # it is ONE task over every read, and on the 2 Mb test fixture 10 chunks took its # slowest unit from 366 s to 36 s. The per-cluster phases already run side by side, # so finer cuts cost them wall time and nothing else -- measured on the same # fixture, cutting the clusters into 10 took their stage-5 work from 27.2 to 43.7 # min (discovery) and 14.1 to 24.8 min (final quant), because that phase was # already saturating every core. An under-cut init is a serial floor under the # whole run; over-cut clusters are a slower run that is still correct. # # Both pairs unset is the production case and is unchanged: LRAA's own 10 Mb / 1 Mb # defaults then apply to every phase, as they did before one plan governed the run. Float? approx_MB_per_cut Float? approx_MB_per_cut_wiggle_window Float? approx_MB_per_cut_init Float? approx_MB_per_cut_wiggle_window_init # An EXISTING plan to run this pipeline on, instead of selecting one below. # # WHY THIS EXISTS. Cut placement READS a bam -- blocked positions and # spanning-read cost are measured from alignments -- so two runs left to # themselves select DIFFERENT cut positions, and their outputs are then # partitioned differently: boundary overhang drops land in different places and # per-chunk quantities are not comparable row for row. Run A, take its # `shared_chunk_plan` output, pass it here as run B's input, and both runs are # forced onto ONE geometry. That is what makes a de novo run and a ref-guided run # over the same library directly comparable rather than approximately so. # # REUSING ONE PLAN ACROSS ANNOTATION MODES IS SUPPORTED, which is the opposite of # what a reader expects. Nothing compares this plan's annotation to the emitter's # -- there is no digest of it and no check that one was present. The driver checks # the property that actually matters instead: no gene locus in THIS run's own gtf # may be SEVERED by one of the plan's cuts # (ChunkedRun.cut_blocking_annotation_models -- the same blocks_cut predicate the # extractor applies to every chunk boundary). A de novo run carries no annotation # and so has nothing to sever; a ref-guided run is checked against its own gtf on # its own terms. The pair is refused only when a cut really does sever a locus, # which is quantified by nobody because both neighbouring chunks omit it, and # never for a mode disagreement. # # A cut merely INSIDE --margin of a locus is a warning, not a refusal, and this # workflow is why: the plan is emitted before phase 1 and every consumer's model # set differs from the emitter's -- per-cluster discovery reads the init GTF, # final quant the consolidated one -- so a margin shortfall is normal here and has # no remedy once the geometry is fixed. The locus is held whole by one chunk; what # it costs is reads of it that overhang the boundary, which are dropped, counted # and named. The margin remains a hard bar at cut SELECTION, where the wiggle # window can still move the position. # # GEOMETRY ONLY. The TPM denominator and the run mode come from THIS run, so # quant_only, initial_annot_gtf and HiFi remain this run's to set. What IS matched # is refused rather than silently misapplied: a plan whose cut-deciding parameters # disagree with this run's -- so the approx_MB_per_cut pair above must be the one # it was selected with -- or that omits a contig this run processes, or that places # that contig's cuts on a different contig length, fails the run before a single # region is read. # # Supplying this SUPPRESSES the emission below rather than running beside it, and # the supplied file is then what every phase applies: the initial pass, per-cluster # discovery and the final quant. File? chunk_plan # How each LRAA phase divides its work; see LRAA.wdl's `scattering`. The three # phases want different shapes. The initial pass is ONE task over every read, # so by_chunk is the only way to parallelise it. # # PER-CLUSTER DISCOVERY stays by_chromosome, and that default is MEASURED: on the # test fixture, giving those runs 10 chunks took their stage-5 work from 27.2 to # 43.7 min for no gain, because many small clusters running side by side already # saturated every core, and cutting them further only multiplies per-unit fixed # cost. # # ALL THREE phases are by_chromosome, so a run has one scattering shape rather # than three and the phases are comparable to each other. off is the only other # accepted value: validate_scattering in LRAA.wdl refuses by_chunk, so these are # not merely defaults here. # # Everything implementing by_chunk is still in the tree and its arguments still # hold on the backends they were made for -- leaf parallelism buys real # concurrency where leaves land on separate machines, and chunking final quant # bounds a peak memory set by the widest single unit rather than by cluster count, # against a cost of two non-preemptible tasks per cluster. Re-enabling it is the # case statement in validate_scattering and nothing else. String scattering_init = "by_chromosome" String scattering_per_cluster = "by_chromosome" String scattering_final_quant = "by_chromosome" # Optional: reuse outputs from a prior initial discovery run and skip LRAA_init File? precomputed_init_quant_tracking File? precomputed_init_gtf # Optional: provide precomputed cluster assignments to skip steps 1-3 File? precomputed_cluster_assignments_tsv # Single-cell barcode and UMI tags String cell_barcode_tag = "CB" String read_umi_tag = "XM" # Resources and docker (propagated to subcalls where applicable) # Cores per LRAA task: the task's cpu request AND the --cpu_budget it divides across # work units. One value each for direct and chromosome-sharded runs; there is no # second knob to multiply it by. Int cpu = 5 # Used only when main_chromosomes is non-empty and LRAA runs are chromosome-sharded. # OPTIONAL, deliberately. LRAA.wdl resolves this as # select_first([cpuScattered, shard_cpu_computed]), so any concrete value here # wins unconditionally and its contig-length-derived per-shard sizing never # runs -- which is what a default of 5 did on every scattered single-cell run. # Unset lets that estimate apply; set it to pin a value. Int? cpuScattered # The initial whole-library LRAA run, separately from the per-cluster runs. It # is ONE task over every read, so its parallelism is bounded by how many chunks # it was cut into, while the per-cluster runs are many small tasks whose # parallelism comes from running side by side. Those want opposite settings: the # init wants a wide budget it can spend across chunks, the cluster tasks want to # be narrow enough that many fit at once. Unset means "same as cpu/memoryGB", # so nothing changes unless asked. Int? cpuInit Int? memoryGBInit # Optional override for direct LRAA runs. When unset, LRAA.wdl uses max(64 GiB, ceil(1.5 x input BAM GiB)). Int? memoryGB # Optional override for chromosome-sharded LRAA workers only. Has no effect when main_chromosomes is empty. Int? memoryGBPerWorkerScattered # Optional by_chunk overrides: make-chunks, chunk leaves, and merge need # independently sized resources. Unset preserves LRAA.wdl's defaults. Int? chunkMakeChunksCpu Int? chunkMakeChunksMemoryGB Int? chunkCpu Int? chunkMemoryGB Int? chunkMergeCpu Int? chunkMergeMemoryGB # The one-off run-wide chunk-plan task. Cut selection runs one concurrent unit per # contig under --cpu_budget, so cpu IS that budget; the memory number is justified # where these two inputs are declared in LRAA_quant_by_cluster.wdl. Forwarded to # the cluster-guided phase as well, so ONE pair of numbers covers the task # wherever a run ends up emitting it. Int cpu_chunk_plan = 16 Int memoryGB_chunk_plan = 8 # PROVISIONAL -- see LRAA-build_sparse_matrices_from_tracking.wdl. Int memoryGBbuildSparseMatrices = 64 Int memoryGBFilterCells = 32 Int memoryGBSeurat = 32 Int memoryGBmergeGTFs = 32 Int memoryGBquantFinal = 32 Int memoryGBscSparseMatrices = 32 # Build the single-cell matrices per contig shard, alongside quantification, # and merge them; instead of one single-threaded pass over the merged # tracking file. Measured on VILLAGE 2493-NPC (1.78 B tracking rows): the # single pass took 12.48 h at ~100 GiB, the shards 555 s at 3.29 GiB each # -- overlapping quant, so off the critical path -- plus a 0.72 h merge at # 7.22 GiB. Output verified identical to the single pass over 1.87 B # non-zeros, all three levels. # # Requires scattering_init = "by_chromosome", which since v0.32.0 is the # default, so on is what a stock run gets. The excluded case is "off": that # path runs one whole-genome LRAA invocation rather than the per-contig # scatter, so no shard artifacts exist to merge. Such a run does not have to # unset this -- it is forced false below and the library-wide build is used. Boolean sc_sparse_from_shards = true String sparseMatrixCsvEngine = "direct" Int sparseMatrixGzipLevel = 1 Int diskSizeGB = 256 String docker = "us-central1-docker.pkg.dev/methods-dev-lab/lraa/lraa-core:latest" String docker_sc = "us-central1-docker.pkg.dev/methods-dev-lab/lraa/lraa-sc:latest" # Seurat clustering parameters (forwarded to Seurat subworkflow) Int min_cells = 10 Int min_features = 1000 Float percent_mt_max = 20.0 String mt_pattern = "^(MT-|mt-|g:(chrM|MT|M):)" Int npcs = 12 Float resolution = 0.6 Int n_variable_features = 2000 Int seed = 1 # Filter good cells parameters. # # `filter_seed` is the emptyDrops Monte Carlo seed, forwarded to the filter # subworkflow's `seed`. Qualified here only because `seed` above is already # Seurat's; the two are separate knobs on purpose -- changing which cells # Seurat is handed and changing Seurat's own draw are different experiments. Boolean enable_filter_good_cells = true Float fdr_threshold = 0.01 Int? lower_threshold Int filter_seed = 1 File? ref_annot_gtf_source_gene_symbols # used as source for gene symbol assignment at the end. } # Only skip the initial discovery call when the downstream-critical artifacts are provided Boolean has_precomputed_init = defined(precomputed_init_gtf) && (defined(precomputed_init_quant_tracking) || defined(precomputed_cluster_assignments_tsv)) Boolean has_precomputed_clusters = defined(precomputed_cluster_assignments_tsv) Boolean run_initial_phase = !has_precomputed_init Boolean run_clustering_phase = !has_precomputed_clusters Boolean run_cluster_guided = single_cell_pipe_mode == "cluster-guided" # The shard path needs shards, so it needs the initial pass to have RUN and to have # scattered. With scattering_init = "off" there is one whole-genome LRAA invocation # and no per-contig artifacts to merge; with the initial pass skipped for # precomputed inputs there are no artifacts at all. Either way the request is forced # off and the library-wide build over the tracking file is used, rather than a merge # of an empty shard list. Boolean use_sc_sparse_from_shards = sc_sparse_from_shards && scattering_init == "by_chromosome" && run_initial_phase # ONE geometry for the whole run; see the input comments above for the precedence and # for why there cannot be two. Float? approx_MB_per_cut_run = if defined(approx_MB_per_cut_init) then approx_MB_per_cut_init else approx_MB_per_cut Float? approx_MB_per_cut_wiggle_window_run = if defined(approx_MB_per_cut_wiggle_window_init) then approx_MB_per_cut_wiggle_window_init else approx_MB_per_cut_wiggle_window # ONE cut plan for the WHOLE cluster-guided run, emitted BEFORE the initial pass. # # Emitted HERE because this is the earliest point in the pipeline: a plan has to # precede phase 1 to be of any use to it, and all three phases -- initial pass, # per-cluster discovery, final quant -- have to apply the same one. Emitted inside # the final quant instead, which is where it used to live, it served that phase alone # and per-cluster discovery still selected its own cuts: 76 distinct geometries over # 270 extractions on the ref-guided smoke run. # # Selected against the REFERENCE annotation, the only one that exists before any # phase has run, and safe BECAUSE it is that early rather than in spite of it: cut # selection will not place a boundary inside an annotated model, and discovery then # runs INSIDE these chunks, so no novel isoform spanning a cut can be constructed -- # no chunk ever held the reads that would support one. Measured on the two smoke # runs: 0 of 994 consolidated ref-guided models and 0 of 865 de novo models straddle # any of the 9 cut positions. De novo has no reference annotation at all and gets an # unannotated plan: nothing was protected from severing at selection time, and each # consumer checks its own models against the plan's cuts where it applies it. # # Emitted only when the CLUSTER-GUIDED phases run. Those are the many-sibling phases # -- one job per cluster -- and shared geometry is what makes their results # comparable. In basic mode the initial pass is the only LRAA run in the pipeline, so # there is nothing for its cuts to agree with, and a plan would cost rather than buy: # selection would move out of that run's own make_chunks, which OVERLAPS it with # extraction (pylib/ChunkedRun.py's prep pool submits a contig's extractions as soon # as that contig's selection finishes), into a separate task the run then waits on. Boolean want_chunk_plan = run_cluster_guided && (scattering_per_cluster != "off" || scattering_final_quant != "off") # Emitted only when the caller did not bring one. A supplied chunk_plan replaces # this call, exactly as LRAA_quant_by_cluster.wdl gates its own fallback emission # on !defined(internal_chunk_plan): one geometry per run, from one source. if (want_chunk_plan && !defined(chunk_plan)) { call QuantByCluster.emit_shared_chunk_plan as emit_run_chunk_plan { input: inputBAM = inputBAM, referenceGenome = referenceGenome, # The most complete annotation that exists BEFORE any phase runs. With the # initial pass skipped, that is the precomputed init GTF -- the models the # cluster phases then quantify -- rather than the reference, so placement # protects those too. annot_gtf = if run_initial_phase then initial_annot_gtf else precomputed_init_gtf, # Restrict selection to what this run processes; see the task input. main_chromosomes = main_chromosomes, HiFi = HiFi, approx_MB_per_cut = approx_MB_per_cut_run, approx_MB_per_cut_wiggle_window = approx_MB_per_cut_wiggle_window_run, cpu = cpu_chunk_plan, memoryGB = memoryGB_chunk_plan, docker = docker } } # The run's geometry, threaded to every phase below: the caller's when it supplied # one, else the emission above. File? run_chunk_plan = if defined(chunk_plan) then chunk_plan else emit_run_chunk_plan.chunk_plan # The initial pass JOINS that geometry when it runs and is scattered, and it has to: # the models it discovers become the annotation the cluster phases quantify, so a # model found across one of the plan's cuts would be severed by every phase that # applies the plan. With scattering_init=off it cannot -- LRAA.wdl's # validate_scattering refuses a plan for a single whole-genome invocation -- so that # configuration holds the property only empirically, and a model that does straddle a # cut is refused by the consumer that would have severed it rather than quantified # quietly. # # Gated on a plan EXISTING rather than on want_chunk_plan, which is what makes a # supplied plan reach the initial pass in BASIC mode too. Emitting one there would # cost and buy nothing -- the argument above -- but a file the caller already holds # costs no task, and basic mode is precisely the comparison shape this input is for: # one LRAA run per mode, which have to cut identically to be compared. With # chunk_plan unset the two conditions are the same one, since the emission runs # exactly when want_chunk_plan holds and always produces the file. if (defined(run_chunk_plan) && run_initial_phase && scattering_init != "off") { File init_chunk_plan = select_first([run_chunk_plan]) } # 1) Initial transcript discovery + quantification on the full BAM (skipped when precomputed inputs are provided) if (run_initial_phase) { call LRAA.LRAA_wf as LRAA_init { input: sample_id = sample_id, referenceGenome = referenceGenome, inputBAM = inputBAM, annot_gtf = initial_annot_gtf, HiFi = HiFi, oversimplify = oversimplify, polyA_known = polyA_known, main_chromosomes = main_chromosomes, region = region, # THE one full-library chromosome partition in a single-cell run, and the only # place this is raised. The per-cluster LRAA_wf calls in # LRAA-cell_cluster_guided.wdl and LRAA_quant_by_cluster.wdl deliberately do # NOT pass it, so they keep LRAA.wdl's default of 1: they run 14 to 32 times # for seconds each, and a wide cpu reservation there cannot be placed until # the box drains. partition_workers = initial_partition_workers, # single-cell: this workflow never surfaces the normalized splice-graph BAM retain_normalized_splice_graph_bam = false, rescue_unassigned_reads_via_transcriptome_alignment = rescue_unassigned_reads_via_transcriptome_alignment, no_weight_reads_by_3prime_agreement = no_weight_reads_by_3prime_agreement, # Init only. The cluster_guided call below deliberately does not get this -- # see the input's declaration for why. num_total_reads = num_total_reads, cell_barcode_tag = cell_barcode_tag, read_umi_tag = read_umi_tag, scattering = scattering_init, build_sc_sparse_shards = use_sc_sparse_from_shards, docker_sc = docker_sc, # The run's ONE geometry, selected above. The initial pass discovers inside # these chunks, which is what makes the same cuts safe for the cluster phases # that quantify what it found. internal_chunk_plan = init_chunk_plan, approx_MB_per_cut = approx_MB_per_cut_run, approx_MB_per_cut_wiggle_window = approx_MB_per_cut_wiggle_window_run, cpu = select_first([cpuInit, cpu]), cpuScattered = cpuScattered, memoryGB = if defined(memoryGBInit) then memoryGBInit else memoryGB, memoryGBPerWorkerScattered = memoryGBPerWorkerScattered, chunkMakeChunksCpu = chunkMakeChunksCpu, chunkMakeChunksMemoryGB = chunkMakeChunksMemoryGB, chunkCpu = chunkCpu, chunkMemoryGB = chunkMemoryGB, chunkMergeCpu = chunkMergeCpu, chunkMergeMemoryGB = chunkMergeMemoryGB, diskSizeGB = diskSizeGB, docker = docker, quant_only = quant_only } } File? init_quant_expr_file = LRAA_init.mergedQuantExpr File? init_quant_tracking_generated = LRAA_init.mergedQuantTracking File? init_quant_tracking_file = if (!run_initial_phase && defined(precomputed_init_quant_tracking)) then select_first([precomputed_init_quant_tracking]) else init_quant_tracking_generated File? init_gtf_generated = LRAA_init.mergedGTF File? init_gtf_file = if (!run_initial_phase && defined(precomputed_init_gtf)) then select_first([precomputed_init_gtf]) else init_gtf_generated # The initial pseudobulk pass's own read-assignment accounting, in its OWN table. # # Run through the same collator as the cluster phases so both tables parse # identically -- level=chunk rows carrying their chunk_id, then this pass's own # total -- but NEVER pooled with them. The init pass normalizes coverage over the # whole BAM while each cluster normalizes over its own reads, and coverage thinning # is depth-dependent, so the same underlying reads are counted differently by the # two phases and their totals legitimately disagree (measured: ~95k pooled against # ~198k summed per-cluster). A single table holding both would invite a # reconciliation that is not supposed to balance. # # --no_all_clusters because there is one invocation here: the aggregate row would # restate the cluster row it aggregated. if (defined(LRAA_init.mergedReadAssignmentSummary)) { call QuantByCluster.collate_read_assignment_summaries as collate_init_read_assignment_summary { input: summaryFiles = [select_first([LRAA_init.mergedReadAssignmentSummary])], clusterIds = ["pseudobulk"], emitAllClusters = false, outputFilePrefix = sample_id + ".init", # Normalized over the WHOLE bam in one pass, not per cluster. That is the # label that keeps this table from being diffed against a cluster_thinned one. population = "pooled_thinned", docker = docker } } # 2) Single-cell sparse matrices from the INITIAL pass's tracking. # # Built whenever an init tracking artifact EXISTS, independent of whether clustering # runs. This block used to sit inside the run_clustering_phase gate below, which made # the matrices an intermediate of clustering rather than a deliverable of the initial # pass: precomputed_cluster_assignments_tsv says only "do not re-cluster", yet it also # dropped the build, and the init_sc_* outputs came back null with the tracking file # they are derived from sitting right there. # # The two sources are exclusive and between them cover every case that HAS a tracking # file: the streaming merge of the per-contig shard artifacts the initial pass built # alongside quantification, else one library-wide pass over the merged tracking file # -- which is also the path a PRECOMPUTED tracking file takes, since a skipped # initial pass leaves no shards behind. if (!use_sc_sparse_from_shards && defined(init_quant_tracking_file)) { call BuildMatrices.BuildSparseMatricesFromTracking as build_sc_from_init_tracking { input: sample_id = sample_id, tracking_file = select_first([init_quant_tracking_file]), docker = docker_sc, memoryGB = memoryGBbuildSparseMatrices, csv_engine = sparseMatrixCsvEngine, gzip_level = sparseMatrixGzipLevel } } # shared_features stays false here: in basic mode features partition across contigs, # so a feature turning up in two shards is an error worth hearing about. if (use_sc_sparse_from_shards) { call BuildMatrices.merge_sc_shard_sparse as merge_sc_from_shards { input: sample_id = sample_id, # NOT defaulted to []: use_sc_sparse_from_shards now implies the initial pass # ran and scattered, so an absent array is a bug to surface rather than an # empty merge to publish. shard_sparse_tars = select_first([LRAA_init.scShardSparse]), docker = docker_sc, gzip_level = sparseMatrixGzipLevel, shared_features = false } } # Optional now, where these used to be select_first over two branches of a block that # always ran exactly one of them. Null only in the single configuration with no init # tracking at all: precomputed_init_gtf plus precomputed_cluster_assignments_tsv and # no precomputed_init_quant_tracking. Whenever clustering runs a tracking file is # guaranteed to exist -- clusters are not precomputed, so has_precomputed_init # required a tracking file -- which is what makes the select_first uses below safe. File? init_sc_gene_tgz = if defined(merge_sc_from_shards.gene_sparse_dir_tgz) then merge_sc_from_shards.gene_sparse_dir_tgz else build_sc_from_init_tracking.gene_sparse_dir_tgz File? init_sc_isoform_tgz = if defined(merge_sc_from_shards.isoform_sparse_dir_tgz) then merge_sc_from_shards.isoform_sparse_dir_tgz else build_sc_from_init_tracking.isoform_sparse_dir_tgz File? init_sc_splice_tgz = if defined(merge_sc_from_shards.splice_pattern_sparse_dir_tgz) then merge_sc_from_shards.splice_pattern_sparse_dir_tgz else build_sc_from_init_tracking.splice_pattern_sparse_dir_tgz File? init_sc_mapping = if defined(merge_sc_from_shards.mapping_file) then merge_sc_from_shards.mapping_file else build_sc_from_init_tracking.mapping_file # 2.5/3) Empty-droplet filtering and clustering, skipped when cluster assignments # were supplied. Only these steps are gated; the matrices above are not. if (run_clustering_phase) { # 2.5) Filter good cells from the gene-level sparse matrix (optional) if (enable_filter_good_cells) { call FilterCells.FilterGoodCells as filter_good_cells { input: sample_id = sample_id, gene_sparse_tar_gz = select_first([init_sc_gene_tgz]), isoform_sparse_tar_gz = select_first([init_sc_isoform_tgz]), splice_pattern_sparse_tar_gz = select_first([init_sc_splice_tgz]), docker = docker_sc, memoryGB = memoryGBFilterCells, fdr_threshold = fdr_threshold, lower_threshold = lower_threshold, seed = filter_seed } } # 3) Cluster cells from the gene-level sparse matrix (filtered or unfiltered) File gene_sparse_for_clustering = if enable_filter_good_cells then select_first([filter_good_cells.filtered_gene_sparse_tar_gz]) else select_first([init_sc_gene_tgz]) call Seurat.GeneSparseM_To_SeuratClusters as cluster_cells { input: sample_id = sample_id, gene_sparse_tar_gz = gene_sparse_for_clustering, docker = docker_sc, memoryGB = memoryGBSeurat, min_cells = min_cells, min_features = min_features, percent_mt_max = percent_mt_max, mt_pattern = mt_pattern, npcs = npcs, resolution = resolution, n_variable_features = n_variable_features, seed = seed } } # Select cluster assignments: use precomputed if provided, otherwise use generated File? cluster_assignments_generated = cluster_cells.cluster_assignments_tsv File cluster_assignments_file = if (has_precomputed_clusters) then select_first([precomputed_cluster_assignments_tsv]) else select_first([cluster_assignments_generated]) # 4) Cluster-guided reconstruction + final single-cell matrices (only if mode is 'cluster-guided') if (run_cluster_guided) { call ClusterGuided.LRAA_cell_cluster_guided as cluster_guided { input: sample_id = sample_id, referenceGenome = referenceGenome, inputBAM = inputBAM, cell_clusters_info = cluster_assignments_file, annot_gtf = if quant_only then initial_annot_gtf else init_gtf_file, HiFi = HiFi, oversimplify = oversimplify, polyA_known = polyA_known, no_weight_reads_by_3prime_agreement = no_weight_reads_by_3prime_agreement, main_chromosomes = main_chromosomes, # Backend-dependent, and separate from the initial phase's value on purpose: # these partitions run once per cluster, which is 14-32 jobs sharing ONE host # locally and 14-32 independent VMs on Terra. cluster_partition_workers = cluster_partition_workers, cell_barcode_tag = cell_barcode_tag, read_umi_tag = read_umi_tag, scattering = scattering_per_cluster, scattering_final_quant = scattering_final_quant, # The request, not this workflow's resolved value: the cluster-guided # phases are gated by scattering_final_quant, so that workflow applies # the same test against its own setting. sc_sparse_from_shards = sc_sparse_from_shards, # The run's ONE geometry, emitted above before the initial pass. Passed # whenever it exists; that workflow gates it per phase, since a phase running # scattering=off takes no plan. internal_chunk_plan = run_chunk_plan, cpu_chunk_plan = cpu_chunk_plan, memoryGB_chunk_plan = memoryGB_chunk_plan, approx_MB_per_cut = approx_MB_per_cut_run, approx_MB_per_cut_wiggle_window = approx_MB_per_cut_wiggle_window_run, cpu = cpu, cpuScattered = cpuScattered, memoryGB = memoryGB, memoryGBPerWorkerScattered = memoryGBPerWorkerScattered, chunkMakeChunksCpu = chunkMakeChunksCpu, chunkMakeChunksMemoryGB = chunkMakeChunksMemoryGB, chunkCpu = chunkCpu, chunkMemoryGB = chunkMemoryGB, chunkMergeCpu = chunkMergeCpu, chunkMergeMemoryGB = chunkMergeMemoryGB, memoryGBmergeGTFs = memoryGBmergeGTFs, memoryGBquantFinal = memoryGBquantFinal, memoryGBscSparseMatrices = memoryGBscSparseMatrices, sparseMatrixCsvEngine = sparseMatrixCsvEngine, sparseMatrixGzipLevel = sparseMatrixGzipLevel, diskSizeGB = diskSizeGB, docker = docker, docker_sc = docker_sc, quant_only_cluster_guided = quant_only } } # 5) Incorporate gene symbols, from the cluster-guided matrices when that phase # ran, else the initial ones -- filtered if the empty-droplet filter actually # RAN, which is not the same as enable_filter_good_cells being set. # # That distinction is load-bearing now that the initial matrices outlive the # clustering gate. filter_good_cells sits INSIDE that gate, so with cluster # assignments supplied it does not run while the input stays at its default of # true; selecting on the input alone therefore picked a null and took every # selector below with it, leaving ref_annot_gtf_source_gene_symbols silently # inert in exactly the configuration that now has matrices to annotate. # # In quant_only mode the GTF values are naturally undefined, since discovery # produces none. Boolean have_filtered_matrices = run_clustering_phase && enable_filter_good_cells File? gtf_for_symbols = if run_cluster_guided then cluster_guided.LRAA_final_gtf else init_gtf_file File? gene_sparse_for_symbols = if run_cluster_guided then cluster_guided.sc_gene_sparse_tar_gz else (if have_filtered_matrices then filter_good_cells.filtered_gene_sparse_tar_gz else init_sc_gene_tgz) File? isoform_sparse_for_symbols = if run_cluster_guided then cluster_guided.sc_isoform_sparse_tar_gz else (if have_filtered_matrices then filter_good_cells.filtered_isoform_sparse_tar_gz else init_sc_isoform_tgz) File? splice_pattern_sparse_for_symbols = if run_cluster_guided then cluster_guided.sc_splice_pattern_sparse_tar_gz else (if have_filtered_matrices then filter_good_cells.filtered_splice_pattern_sparse_tar_gz else init_sc_splice_tgz) File? mapping_for_symbols = if run_cluster_guided then cluster_guided.sc_gene_transcript_splicehash_mapping else init_sc_mapping if (defined(ref_annot_gtf_source_gene_symbols) && defined(gene_sparse_for_symbols) && defined(isoform_sparse_for_symbols) && defined(splice_pattern_sparse_for_symbols) && defined(mapping_for_symbols)) { call GeneSymbols.Incorporate_gene_symbols as add_gene_symbols { input: sample_id = sample_id, reference_gtf = select_first([ref_annot_gtf_source_gene_symbols]), final_gtf = gtf_for_symbols, final_sc_gene_sparse_tar_gz = select_first([gene_sparse_for_symbols]), final_sc_isoform_sparse_tar_gz = select_first([isoform_sparse_for_symbols]), final_sc_splice_pattern_sparse_tar_gz = select_first([splice_pattern_sparse_for_symbols]), final_sc_gene_transcript_splicehash_mapping = select_first([mapping_for_symbols]), docker = docker } } # 6) Integrated TSS/PolyA sites (cluster-guided mode): the cluster-guided site beds, # supplemented by basic (initial-catalog) sites no cluster-guided site reproduces # within the site-aggregation tolerance. The basic beds come from LRAA_init; when the # initial pass was skipped for a precomputed gtf they are derived here from that gtf, # the same file LRAA_init's beds would have been derived from. if (run_cluster_guided && !defined(LRAA_init.tssBed) && defined(init_gtf_file)) { call LRAA.splice_pattern_collapse as init_sites_from_precomputed { input: lraaGtf = select_first([init_gtf_file]), outputFilePrefix = sample_id + ".init_from_precomputed", docker = docker } } File? basic_TSS_bed = if defined(LRAA_init.tssBed) then LRAA_init.tssBed else init_sites_from_precomputed.tssBed File? basic_PolyA_bed = if defined(LRAA_init.polyaBed) then LRAA_init.polyaBed else init_sites_from_precomputed.polyaBed if (run_cluster_guided && defined(basic_TSS_bed) && defined(basic_PolyA_bed) && defined(cluster_guided.LRAA_final_TSS_bed) && defined(cluster_guided.LRAA_final_PolyA_bed)) { call LRAA.integrate_TSS_PolyA_sites as integrate_sites { input: primary_TSS_bed = select_first([cluster_guided.LRAA_final_TSS_bed]), primary_PolyA_bed = select_first([cluster_guided.LRAA_final_PolyA_bed]), supplement_TSS_bed = select_first([basic_TSS_bed]), supplement_PolyA_bed = select_first([basic_PolyA_bed]), outputFilePrefix = sample_id, docker = docker } } output { # Initial discovery outputs File? init_quant_expr = init_quant_expr_file File? init_quant_tracking = init_quant_tracking_generated File? init_gtf = init_gtf_generated # Derived from init_gtf, and the deliverable set for a BASIC run: no # cluster-guided phase runs there, so the initial catalog is the final one. # Absent when the initial pass was skipped for a precomputed gtf, which was # collapsed by whatever run produced it. File? init_splice_pattern_collapsed_gtf = LRAA_init.splicePatternCollapsedGTF File? init_splice_pattern_collapsed_merge_report = LRAA_init.splicePatternCollapsedMergeReport File? init_splice_pattern_collapsed_gene_conflicts = LRAA_init.splicePatternCollapsedGeneConflicts File? init_TSS_bed = LRAA_init.tssBed File? init_PolyA_bed = LRAA_init.polyaBed # This pass's read-assignment accounting. Reported SEPARATELY from the # cluster-phase tables below and never folded into them; see the call site for # why the two phases' totals are not meant to reconcile. Absent when the initial # pass was skipped for precomputed inputs. File? init_read_assignment_summary = collate_init_read_assignment_summary.collatedSummary # The ONE cut geometry this run applied or chose, and the file a second run is # handed to be forced onto the same cuts. Three sources, in precedence order: # # the chunk_plan input, when the caller supplied one -- every phase applied # exactly that, so it is what a sibling must be given to match; # # else the plan emitted above, which only CLUSTER-GUIDED runs emit (see # want_chunk_plan: basic mode has one LRAA run, so there is nothing for its # cuts to agree with and a separate selection task would cost the overlap # between selection and extraction inside that run's own make_chunks); # # else, in basic mode, the geometry the initial pass's own chunking chose, # GATHERED from its shards rather than re-selected -- LRAA.wdl's # gather_shard_cut_plans over each shard's shard_cut_plan.json. This is what # lets a basic run kickstart a cluster-guided one: its init_gtf becomes # precomputed_init_gtf and its clusters become # precomputed_cluster_assignments_tsv, and this makes its GEOMETRY carry # across the same lineage instead of the derived run picking its own bounds. # # Deliberately not gathered in cluster-guided mode: what that mode emits is # unchanged, so a cluster-guided run whose every phase is scattering=off still # reports no plan rather than the init pass's cuts, which governed only the # init pass. Null also when the initial pass was skipped for precomputed # inputs, or when it ran with no_chunk and so placed no cuts. File? shared_chunk_plan = if defined(chunk_plan) then chunk_plan else if run_cluster_guided then emit_run_chunk_plan.chunk_plan else LRAA_init.gatheredChunkPlan # Initial single-cell matrices and clustering inputs/outputs File? init_sc_gene_sparse_tar_gz = init_sc_gene_tgz File? init_sc_isoform_sparse_tar_gz = init_sc_isoform_tgz File? init_sc_splice_pattern_sparse_tar_gz = init_sc_splice_tgz File? init_sc_gene_transcript_splicehash_mapping = init_sc_mapping # Filter empty droplets File? filtered_gene_sparse_tar_gz = filter_good_cells.filtered_gene_sparse_tar_gz File? filtered_isoform_sparse_tar_gz = filter_good_cells.filtered_isoform_sparse_tar_gz File? filtered_splice_pattern_sparse_tar_gz = filter_good_cells.filtered_splice_pattern_sparse_tar_gz File? good_cell_barcodes = filter_good_cells.good_cell_barcodes File? filtering_summary = filter_good_cells.filtering_summary # Seurat for gene-based filtering, cell clustering, and umap File? seurat_umap_pdf = cluster_cells.umap_pdf File? seurat_umap_with_clusters_tsv = cluster_cells.umap_with_clusters_tsv File? seurat_rds_initial = cluster_cells.seurat_rds_initial File? seurat_rds = cluster_cells.seurat_rds File? seurat_cluster_assignments = cluster_assignments_generated # Final cluster-guided LRAA outputs (main deliverables) File? final_gtf = cluster_guided.LRAA_final_gtf File? final_gtf_tracking = cluster_guided.LRAA_final_gtf_tracking # From the cross-cluster merged catalog, not from any single cluster: gene # assignments are unified by that merge, and the collapse groups by gene_id. File? final_splice_pattern_collapsed_gtf = cluster_guided.LRAA_final_splice_pattern_collapsed_gtf File? final_splice_pattern_collapsed_merge_report = cluster_guided.LRAA_final_splice_pattern_collapsed_merge_report File? final_splice_pattern_collapsed_gene_conflicts = cluster_guided.LRAA_final_splice_pattern_collapsed_gene_conflicts File? final_TSS_bed = cluster_guided.LRAA_final_TSS_bed File? final_PolyA_bed = cluster_guided.LRAA_final_PolyA_bed # Cluster-guided sites plus basic sites outside the aggregation window, each row # labelled by `source` (cluster_guided | basic); counts in the summary. Absent in # basic mode, where the initial catalog's beds are the only ones. File? integrated_TSS_bed = integrate_sites.integratedTssBed File? integrated_PolyA_bed = integrate_sites.integratedPolyaBed File? integrated_sites_summary = integrate_sites.summary File? final_tracking = cluster_guided.LRAA_final_tracking File? final_sc_gene_sparse_tar_gz = cluster_guided.sc_gene_sparse_tar_gz File? final_sc_isoform_sparse_tar_gz = cluster_guided.sc_isoform_sparse_tar_gz File? final_sc_splice_pattern_sparse_tar_gz = cluster_guided.sc_splice_pattern_sparse_tar_gz File? final_sc_gene_transcript_splicehash_mapping = cluster_guided.sc_gene_transcript_splicehash_mapping # Per-cluster read-assignment accounting from the cluster-guided phases, one # table per phase. Each holds level=chunk / level=cluster / level=all_clusters # rows; rows of different level must not be summed. Absent in basic mode, and # the prelim table is additionally absent in quant_only mode. File? prelim_collated_read_assignment_summary = cluster_guided.prelim_collated_read_assignment_summary File? final_collated_read_assignment_summary = cluster_guided.collated_read_assignment_summary # Convenience tarballs and matrices from the cluster-guided phase File? partitioned_cluster_bams_tar = cluster_guided.LRAA_partitioned_cluster_bams_tar File? final_cluster_exprs_tar = cluster_guided.LRAA_final_cluster_exprs_tar File? final_cluster_trackings_tar = cluster_guided.LRAA_final_cluster_trackings_tar # Preliminary intermediate outputs from cluster-guided discovery phase File? prelim_cluster_gtfs = cluster_guided.LRAA_prelim_cluster_gtfs File? prelim_cluster_read_trackings = cluster_guided.LRAA_prelim_cluster_read_trackings File? prelim_cluster_pseudobulk_exprs = cluster_guided.LRAA_prelim_cluster_pseudobulk_exprs File? cluster_gene_counts_matrix = cluster_guided.cluster_gene_counts_matrix File? cluster_gene_TPM_matrix = cluster_guided.cluster_gene_TPM_matrix File? cluster_isoform_counts_matrix = cluster_guided.cluster_isoform_counts_matrix File? cluster_isoform_TPM_matrix = cluster_guided.cluster_isoform_TPM_matrix File? cluster_isoform_counts_forDiffIsoUsage = cluster_guided.cluster_isoform_counts_forDiffIsoUsage File? incl_gene_symbols_gffcompare_tracking = add_gene_symbols.gffcompare_tracking File? incl_gene_symbols_gffcompare_stats = add_gene_symbols.gffcompare_stats File? incl_gene_symbols_updated_gtf = add_gene_symbols.updated_gtf_with_gene_symbols File? incl_gene_symbols_updated_id_mappings = add_gene_symbols.updated_id_mappings File? incl_gene_symbols_gene_sparse_tar_gz = add_gene_symbols.updated_gene_sparse_tar_gz File? incl_gene_symbols_isoform_sparse_tar_gz = add_gene_symbols.updated_isoform_sparse_tar_gz File? incl_gene_symbols_splice_pattern_sparse_tar_gz = add_gene_symbols.updated_splice_pattern_sparse_tar_gz } }