Pipelines: Nextflow and SLURM#

Two ways to run migec over more than one sample. They answer different questions: Nextflow when the pipeline continues past the consensus into alignment, calling or AIRR; SLURM when you want three commands over a cohort and nothing else.

Nextflow#

integrations/nextflow/ is an nf-core-style local module set with a runnable entry point.

main.nf                                  --mode consensus | ctdna | airr

modules/migec/checkout/main.nf           reads     -> tagged FASTQ + QC tables
modules/migec/refine/main.nf             tagged    -> corrected FASTQ + barcode table + cells
modules/migec/assemble/main.nf           corrected -> one consensus per molecule
modules/migec/plot/main.nf               the tables -> SVG figures (needs gnuplot)

modules/downstream/align/main.nf         consensus -> BAM, tags carried and checked
modules/downstream/callvariants/main.nf  BAM       -> VCF (LoFreq or Mutect2)
modules/downstream/arda/main.nf          consensus -> AIRR clonotypes

subworkflows/migec/main.nf               the three migec stages chained
subworkflows/migec_ctdna/main.nf         + align + call: rare somatic variants
subworkflows/migec_airr/main.nf          + arda: immune repertoires
nextflow.config                          defaults, every one read with getOrDefault

nextflow run integrations/nextflow -profile docker \
    --mode ctdna --input 'data/*_R{1,2}.fq.gz' --preset tso500 \
    --fasta ref.fa --outdir results/

nextflow run integrations/nextflow --mode airr --input 'data/*_R{1,2}.fq.gz' --preset migec

Or include a subworkflow in a pipeline of your own:

include { MIGEC } from './integrations/nextflow/subworkflows/migec/main'

workflow {
    ch_reads = Channel.fromFilePairs(params.input)
        .map { id, files -> [ [ id: id, preset: '10x-v2', payload_mate: 2 ], files ] }

    MIGEC(ch_reads, file(params.cell_whitelist ?: 'NO_FILE'))

    MIGEC.out.consensus.view()
}

Per-sample keys in meta win over the params.* defaults, so one run can mix chemistries: bc_pattern, preset, read_structure, read_structure2, barcodes, max_offset, payload_mate, expect_cells, rt_error, contig, fast, aligner, caller, species.

What the downstream modules encode#

Two rules, both from Variant calling: which caller, and what it can possibly see:

Collapse first, then align once. Aligning raw reads and grouping on (position, UMI) is the other order in use, and it costs one alignment per read rather than per molecule, with the aligner seeing uncorrected sequence.

A standard variant caller, never a UMI-aware one. After assemble a caller’s depth already is a molecule count. UMI-VarCal and UMIErrorCorrect group and consense themselves, so they replace assemble rather than following it.

Never: do not set a family-size filter downstream of assemble. Every family has size 1 by construction, so --min-family-size 3 discards the entire library and reports zero variants without an error.

The align module checks its own output: if no MI:Z: tag reaches the BAM it exits non-zero and names the flag each aligner needs (-y for minimap2 and minibwa map, -C for bwa mem). Without that check the failure is silent and surfaces much later as an untagged BAM.

SLURM#

integrations/slurm/ is two sbatch templates and a sample sheet, for the case where the deliverable is the consensus and a workflow engine would be ceremony.

# one sample
sbatch --export=ALL,R1=s1_R1.fq.gz,R2=s1_R2.fq.gz,SAMPLE=s1,PRESET=10x-v2,PAYLOAD_MATE=2 \
       integrations/slurm/migec_sample.sbatch

# a cohort, one array task per row
sbatch --array=1-$(($(wc -l < samples.tsv) - 1)) \
       integrations/slurm/migec_array.sbatch samples.tsv

Both run without SLURM as ordinary bash scripts – every SLURM variable has a fallback – which is how they are tested and how a layout should be checked before a cohort is queued:

R1=s1.fq.gz SAMPLE=s1 BC_PATTERN='0:12' bash integrations/slurm/migec_sample.sbatch

Note: array task 1 is the first data row. The header is skipped rather than counted, so the range is 1-(rows - 1). An --array=0-N runs a task that reads the header as a sample and fails somewhere confusing instead of at the sheet.

Sizing the request: checkout scales with reads (give it cores), refine with distinct barcodes (give it memory; table_bytes in its JSON sizes the next run), assemble with one bucket at a time (cores; peak memory is set by the bucket count, not the library). 16 cores and 32 GB covers a typical targeted or single-cell library.

Never: -t changes the wall clock and nothing else. Every stage is byte-identical at any thread count, so a retry on a different node or an escalating --requeue cannot produce a result that disagrees with the first attempt. That is what makes automatic retries safe here, and tests/benchmark/ asserts it.

Why three processes#

The stages have different shapes, and one process forces one answer for all three:

stage	scales with	threads	memory
`checkout`	reads	1,548,835 reads/s at 16	chunk-bounded, plus the UMI counters
`refine`	distinct barcodes	1,554,156 reads/s at 16	the barcode table, ~96 B each
`assemble`	reads, then buckets	2,470,928 reads/s at 16	one bucket per worker

Splitting them means a failed assemble resumes without re-running the demultiplex, and each stage gets the label and the retry that fits it. refine carries process_high_memory because its memory is set by the number of distinct barcodes and by nothing else: a 200 GB shallow run and a 200 GB deep one need wildly different amounts of it, and FASTQ size predicts neither.

Retries are safe. Every stage’s output is byte-identical at any thread count (Speed and memory), so an attempt that gets 4 cpus and an attempt that gets 32 produce the same bytes – which is what makes errorStrategy 'retry' with escalating resources sound here rather than merely convenient.

process {
    withLabel: process_high        { cpus = 16; memory = 32.GB }
    withLabel: process_medium      { cpus = 8;  memory = 16.GB }
    withLabel: process_high_memory { memory = { 64.GB * task.attempt } }
}

The three things that go wrong#

The barcode read is not always the payload read. On 10x, R1 is 26 nt of cell barcode and UMI and nothing else. checkout is given both mates; the later stages must then run on the mate that carries cDNA. payload_mate: 2 says so. Assuming R1 produces empty consensuses on every droplet chemistry, and nothing in the run reports it as an error.

A per-sample `false` is not an absent value. Groovy’s ?: treats false as absent, so [ id: 's1', contig: false ] against params.migec_contig = true would silently mean its opposite. The modules use meta.containsKey('contig') ? meta.contig : params... for every boolean, which is the one direction a per-sample override exists to make possible.

``–rt-error`` names a chemistry, not a number. rt (1e-4, caps at Q40) for anything with a reverse transcription step, medium (1e-5) for an ordinary polymerase and no RT, high (1e-6) for a proofreading one. It is the one-molecule floor – 10x’s Q60 requires two UMIs to agree, and combining molecules is arda’s job, downstream of this. See The RT/PCR floor — the quality no consensus can beat.

Smoke tests#

Every process has a stub: block, so nextflow run -stub-run walks the whole graph without migec installed – the test to run before the real one.

params.migec_limit_read and params.migec_limit_umi pass through to the stages that take them (--limit-read on all three, --limit-umi on refine and assemble), which gets an answer out of a 400 GB run in a minute.

Warning

A limit is not a sample, and every limited run says so in its own report. Use subsample – a smaller library that is still a library when the output has to behave like the library.

Note

The -stub-run and the real run have not been executed on the machine these docs were written on – nextflow is not installed there. The modules are reviewed against the nf-core module spec, not verified by a pipeline run; treat the first run in your own pipeline as the verification.

The SLURM templates are the opposite case: they have been run end to end here, as ordinary bash, producing a consensus and its figures, and their error paths were exercised too. Between the two, the migec commands themselves are covered either way – both call the same three.