#!/usr/bin/env nextflow
/*
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    OwnGenome — whole genome, from reads to report
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    Each samplesheet row starts from FASTQ (trimmed, aligned and called
    here), from a BAM or CRAM (called here), or from a VCF with an optional
    BAM or CRAM from any other caller (e.g. nf-core/sarek). From there it runs
    pharmacogenomics, variant annotation, clinical screening, BAM analysis,
    structural variant calling, and consolidated reporting.

    https://github.com/GeiserX/OwnGenome
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/

nextflow.enable.dsl = 2
// Module-scoped tools in modules/local/<module>/resources/usr/bin (the MoChA
// image needs a ps for the task metrics).
nextflow.enable.moduleBinaries = true

/*
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    IMPORT WORKFLOWS
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/

include { UPSTREAM     } from './workflows/upstream'
include { PGX          } from './workflows/pgx'
include { ANNOTATION   } from './workflows/annotation'
include { CLINICAL     } from './workflows/clinical'
include { BAM_ANALYSIS } from './workflows/bam_analysis'
include { SV           } from './workflows/sv'
include { REPORTING    } from './workflows/reporting'
include { FAMILY       } from './workflows/family'
include { VCF_PRECHECK } from './modules/local/vcf_precheck/main'
include { CRAM_TO_BAM  } from './modules/local/cram_archive/main'
include { CRAM_ARCHIVE } from './modules/local/cram_archive/main'
include { TRAITS_SITES } from './modules/local/traits/main'
include { TRAITS_PHASE } from './modules/local/traits/main'
include { TRAITS       } from './modules/local/traits/main'
include { BLOODTYPE_REGIONS; BLOODTYPE_DEPTH; BLOODTYPE } from './modules/local/bloodtype/main'
include { ARCHAIC_GENOTYPES; ARCHAIC_IBDMIX; ARCHAIC_REPORT } from './modules/local/archaic/main'
include { ANCIENT_SITES; ANCIENT_GENO; ANCIENT_PCA; ANCIENT_REPORT } from './modules/local/ancient/main'
include { LPA_REGIONS; LPA_DEPTH; LPA } from './modules/local/lpa/main'
include { VNTYPER } from './modules/local/vntyper/main'
include { LOCITYPER } from './modules/local/locityper/main'
include { VIRAL_READS; VIRAL_REPORT } from './modules/local/viral/main'
include { EXOMISER_PREPARE; EXOMISER; EXOMISER_REPORT } from './modules/local/exomiser/main'
include { CHIP_PREPARE; CHIP_MUTECT2; CHIP_ANNOTATE; CHIP_REPORT } from './modules/local/chip/main'
include { MEI_DEPTH; MEI; MEI_REPORT } from './modules/local/mei/main'
include { MOCHA_PREP; MOCHA_PHASE; MOCHA_CALL; MOCHA_REPORT } from './modules/local/mocha/main'

/*
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    MAIN WORKFLOW
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/

workflow {

    // ─── Validate inputs ────────────────────────────────────────────────
    // Every parameter must be one nextflow_schema.json declares: Nextflow
    // takes any --name as a new parameter, so a misspelt option (--outdirr,
    // --max-cpus) would be ignored with no word. --pedigree is read with
    // containsKey and has no default; config_profile_* are the test
    // profiles' labels.
    def schema_params = new groovy.json.JsonSlurper()
        .parseText(file("${projectDir}/nextflow_schema.json").text)['definitions']
        .values().collectMany { group -> (group['properties'] ?: [:]).keySet() as List }
    def known_params = schema_params + ['pedigree', 'config_profile_name', 'config_profile_description']
    def unknown_params = params.keySet().findAll { k -> !known_params.contains(k) }.sort()
    if (unknown_params) {
        error "unknown parameter${unknown_params.size() > 1 ? 's' : ''} " +
              "${unknown_params.collect { k -> '--' + k }.join(', ')}: not an option of this pipeline " +
              "(nextflow_schema.json lists them; options take underscores, as --max_cpus). Check the spelling."
    }

    if (!params.input) {
        error "Please provide a samplesheet with --input <samplesheet.csv>"
    }

    if (!params.reference) {
        error "Please provide a reference FASTA with --reference <path/to/GRCh38.fasta>"
    }
    if (!(params.sex_check in ['fail', 'warn'])) {
        error "--sex_check must be 'fail' or 'warn', got '${params.sex_check}'."
    }
    if (!(params.freemix_warn instanceof Number) || params.freemix_warn < 0 || params.freemix_warn >= 1) {
        error "--freemix_warn must be a fraction from 0 up to (not including) 1 (FREEMIX above it is reported " +
              "as possible contamination), got '${params.freemix_warn}'."
    }
    // A bare --prs_keep_ambiguous can reach here as the string 'true'.
    if (!(params.prs_keep_ambiguous.toString() in ['true', 'false'])) {
        error "--prs_keep_ambiguous takes no value (or true / false), got '${params.prs_keep_ambiguous}'."
    }
    if (params.reference.endsWith('.gz')) {
        error "--reference ${params.reference} is compressed. The tools here need a plain FASTA with a .fai index: " +
              "decompress it (gunzip, or bgzip -d) and run 'samtools faidx' on the result."
    }

    def tools_list = params.tools ? params.tools.split(',').collect { it.trim() }.findAll { it } : []

    // Every name the workflows gate on. A name outside this list would be
    // ignored silently (e.g. 'clinvar_screen' instead of 'clinvar').
    def known_tools = [
        'pharmcat', 'cpic', 'clinvar', 'pypgx',
        'vep', 'vcfanno', 'slivar', 'clinical_filter',
        'cpsr', 'roh', 'prs', 'ancestry', 'ancestry_fine', 'mito_haplogroup',
        'hla_typing', 'expansion_hunter', 'stranger', 'telomere_hunter', 'mosdepth', 'mito_variants', 'cyrius',
        'manta', 'delly', 'cnvpytor', 'duphold', 'annotsv', 'survivor_merge',
        'sample_qc', 'cram_archive', 'parascopy', 'paraphase', 'y_haplogroup', 'html_report', 'multiqc',
        'consumer_export',
        'bloodtype',
        'lpa',
        'vntyper',
        'traits',
        'archaic',
        'ancient',
        'y_str',
        'viral',
        'exomiser',
        'chip',
        'mei',
        'mocha',
        'locityper',
    ]
    def unknown_tools = tools_list.findAll { !known_tools.contains(it) }
    if (unknown_tools) {
        error "unknown tool${unknown_tools.size() > 1 ? 's' : ''} ${unknown_tools.join(', ')} in --tools. " +
              "Known tools: ${known_tools.join(', ')}."
    }

    // exomiser (step 48) runs when --hpo names the symptoms; naming it in
    // --tools without them stops the run.
    if (params.hpo && !tools_list.contains('exomiser')) {
        tools_list << 'exomiser'
    }
    if (tools_list.contains('exomiser')) {
        def hpo_terms = params.hpo ? params.hpo.toString().split(/[,\s]+/).findAll { it } : []
        if (!hpo_terms) {
            error "Tool 'exomiser' ranks variants against a person's symptoms: give them with --hpo as HPO ids " +
                  "(--hpo HP:0001250,HP:0001263; https://hpo.jax.org)."
        }
        def bad_terms = hpo_terms.findAll { !(it ==~ /HP:[0-9]{7}/) }
        if (bad_terms) {
            error "--hpo: ${bad_terms.join(', ')} ${bad_terms.size() > 1 ? 'are not HPO ids' : 'is not an HPO id'}: " +
                  "an HPO id is HP: and seven digits, e.g. HP:0001250 (Seizure)."
        }
        if (!(params.exomiser_assembly in ['hg38', 'hg19'])) {
            error "--exomiser_assembly must be hg38 or hg19, got '${params.exomiser_assembly}'."
        }
    }

    // Fail-fast: stop when enabled tools lack required databases

    def db_requirements = [
        ['vep',              'vep_cache',         '--vep_cache'],
        ['cpsr',             'pcgr_data',         '--pcgr_data'],
        ['cpsr',             'vep_cache_cpsr',    '--vep_cache_cpsr'],
        ['expansion_hunter', 'expansion_catalog', '--expansion_catalog'],
        ['hla_typing',       'hla_dat',           '--hla_dat'],
        ['hla_typing',       'hla_genes',         '--hla_genes'],
        ['clinvar',          'clinvar',           '--clinvar'],
        ['clinvar',          'clinvar_index',     '--clinvar_index'],
        ['pypgx',            'pypgx_bundle',      '--pypgx_bundle'],
        ['annotsv',          'annotsv_annotations', '--annotsv_annotations'],
        ['cnvpytor',         'cnvpytor_resources', '--cnvpytor_resources'],
        ['sample_qc',        'somalier_sites',     '--somalier_sites'],
        ['sample_qc',        'verifybamid2_panel', '--verifybamid2_panel'],
        ['cyrius',           'cyrius_install',     '--cyrius_install'],
        ['parascopy',        'parascopy_data',     '--parascopy_data'],
        ['y_haplogroup',     'yleaf_data',         '--yleaf_data'],
        ['ancestry_fine',    'ancestry_fine_panel', '--ancestry_fine_panel'],
        ['consumer_export',  'consumer_sites',     '--consumer_sites'],
        ['archaic',          'archaic_data',       '--archaic_data'],
        ['ancient',          'ancient_panel',      '--ancient_panel'],
        ['viral',            'viral_ref',          '--viral_ref'],
        ['exomiser',         'exomiser_data',      '--exomiser_data'],
        ['chip',             'chip_gff',           '--chip_gff'],
        ['mei',              'megane_kmers',       '--megane_kmers'],
        ['mocha',            'mocha_data',         '--mocha_data'],
        ['locityper',        'locityper_data',     '--locityper_data'],
    ]

    db_requirements.each { tool, param_name, flag ->
        if (tools_list.contains(tool) && !params[param_name]) {
            error "Tool '${tool}' is enabled in --tools but ${flag} is not set. " +
                  "Either provide ${flag} or remove '${tool}' from --tools."
        }
    }

    // --slivar_bin is accepted for one release so an old command line still
    // starts; slivar now runs from the pinned image (SLIVAR_IMAGE).
    if (params.slivar_bin) {
        log.warn "--slivar_bin is deprecated and ignored: slivar now runs from the pinned image " +
                 "(SLIVAR_IMAGE in versions.env). Remove the option; the next release drops it."
    }

    // KIR is a second T1K pass of hla_typing, against IPD-KIR
    if (params.kir && !tools_list.contains('hla_typing')) {
        error "--kir types the KIR genes in a second T1K pass of hla_typing: add 'hla_typing' to --tools."
    }
    if (params.kir && !params.kir_dat) {
        error "--kir needs --kir_dat, IPD-KIR's kir.dat (scripts/setup.sh --kir-data installs it)."
    }
    // --hla_qc: OptiType, a second HLA caller, checks T1K's HLA-A and HLA-B.
    // A warning, not a stop: the test_all profile turns it on, and runs of
    // that profile with a short --tools list must still start.
    if (params.hla_qc && !tools_list.contains('hla_typing')) {
        log.warn "--hla_qc does nothing without hla_typing in --tools: OptiType runs as a second caller beside T1K."
    }
    if (tools_list.contains('parascopy') && !(params.parascopy_population in ['AFR', 'AMR', 'EAS', 'EUR', 'SAS'])) {
        error "--parascopy_population must be one of AFR, AMR, EAS, EUR, SAS (the 1000 Genomes models), " +
              "got '${params.parascopy_population}'."
    }

    if (params.paraphase_regions && !(params.paraphase_regions.toString() ==~ /[A-Za-z0-9_]+(,[A-Za-z0-9_]+)*/)) {
        error "--paraphase_regions must be Paraphase region names separated by commas (e.g. smn1,hba), " +
              "got '${params.paraphase_regions}'."
    }

    // cpic requires pharmcat (it parses PharmCAT JSON output)
    if (tools_list.contains('cpic') && !tools_list.contains('pharmcat')) {
        error "Tool 'cpic' is enabled in --tools but 'pharmcat' is not. " +
              "CPIC lookup requires PharmCAT JSON output — add 'pharmcat' to --tools or remove 'cpic'."
    }

    // stranger requires expansion_hunter (it annotates EH STR VCF output)
    if (tools_list.contains('stranger') && !tools_list.contains('expansion_hunter')) {
        error "Tool 'stranger' is enabled in --tools but 'expansion_hunter' is not. " +
              "Stranger annotates ExpansionHunter VCF output — add 'expansion_hunter' to --tools or remove 'stranger'."
    }

    // survivor_merge keeps calls seen by two or more callers; with one caller
    // it would write a header-only consensus and report success.
    def sv_callers = ['manta', 'delly', 'cnvpytor'].findAll { tools_list.contains(it) }
    if (tools_list.contains('survivor_merge') && sv_callers.size() < 2) {
        error "Tool 'survivor_merge' needs at least two SV callers in --tools (manta, delly, cnvpytor); " +
              "got ${sv_callers ? sv_callers.join(', ') : 'none'}. Add another caller or remove 'survivor_merge'."
    }

    // ClinVar: paired inputs required together
    if (params.clinvar && !params.clinvar_index) {
        error "When --clinvar is provided, --clinvar_index must also be provided."
    }
    if (!params.clinvar && params.clinvar_index) {
        error "When --clinvar_index is provided, --clinvar must also be provided."
    }

    // ─── Parse samplesheet ──────────────────────────────────────────────
    // Columns: sample, then one starting point per row:
    //   fastq_1,fastq_2  reads: trimmed, aligned and called here
    //   bam,bam_index    a BAM without a VCF: called here
    //   cram,crai        a CRAM instead of the BAM (read with --reference,
    //                    which must be the FASTA it was written with)
    //   vcf,vcf_index    a VCF from any caller, with bam,bam_index or
    //                    cram,crai optional, and gvcf,gvcf_index optional:
    //                    the gVCF of the same calls (step 03 writes one), read
    //                    by PharmCAT and PRS as DEEPVARIANT's would be
    // and sex (male or female), required on every row that is called.
    // Rows are read and checked here, before any task starts, so a bad row
    // stops the run at once.
    def samplesheet_rows = file(params.input, checkIfExists: true).splitCsv(header: true, strip: true)
    def seen_samples = [] as Set
    samplesheet_rows.each { row ->
        if (!row.sample) {
            error "Samplesheet row without a 'sample' value. Columns found: ${row.keySet()}"
        }
        // Sanitize sample ID — used in shell commands, file paths, and HTML output
        if (!(row.sample ==~ /^[a-zA-Z0-9._-]+$/)) {
            error "Sample name '${row.sample}' contains invalid characters. Use only a-z, A-Z, 0-9, '.', '_', '-'"
        }
        // CPSR, bin/pgx_parse.py and bin/collect_summary.py take the id as an
        // argument value, and argparse reads '-x' as an option. (CPSR's 3 to
        // 40 characters are handled by bin/cpsr_sample_id.)
        if (row.sample.startsWith('-')) {
            error "Sample name '${row.sample}' starts with '-', which tools read as an option. Start it with a letter or digit."
        }
        // Sample ids name the output directory and key every per-sample join
        if (!seen_samples.add(row.sample)) {
            error "Sample '${row.sample}' appears more than once in ${params.input}. Each sample needs exactly one row."
        }
        // Pairs go together
        [['fastq_1', 'fastq_2'], ['bam', 'bam_index'], ['cram', 'crai'], ['vcf', 'vcf_index'], ['gvcf', 'gvcf_index']].each { a, b ->
            if (row[a] && !row[b]) {
                error "Sample '${row.sample}': '${a}' provided without '${b}'. Both are required together."
            }
            if (!row[a] && row[b]) {
                error "Sample '${row.sample}': '${b}' provided without '${a}'. Both are required together."
            }
        }
        if (row.fastq_1 && (row.bam || row.cram || row.vcf)) {
            error "Sample '${row.sample}': a row starts from FASTQ or from a BAM, CRAM or VCF, not both. " +
                  "Remove the fastq columns to use the BAM, CRAM or VCF, or the other columns to align the reads."
        }
        if (row.gvcf && !row.vcf) {
            error "Sample '${row.sample}': a gvcf goes with the vcf it belongs to. A row without a VCF is " +
                  "called here, and DEEPVARIANT writes its own gVCF: remove the gvcf columns, or add the vcf."
        }
        if (row.bam && row.cram) {
            error "Sample '${row.sample}': the row has both a BAM and a CRAM. Keep one of them."
        }
        if (row.cram && !(row.cram ==~ /.*\.cram$/)) {
            error "Sample '${row.sample}': the cram column names ${row.cram}, which does not end in .cram. " +
                  "A BAM goes in the bam column."
        }
        if (!row.fastq_1 && !row.bam && !row.cram && !row.vcf) {
            error "Sample '${row.sample}': the row has no input. Give fastq_1 and fastq_2, or bam and bam_index, " +
                  "or cram and crai, or vcf and vcf_index (with an optional bam and bam_index, or cram and crai)."
        }
        // Sex sets the chrX/chrY ploidy of DeepVariant and of ExpansionHunter
        // (whose default is female), and INDEXCOV checks it against the BAM
        if (row.sex && !(row.sex.toLowerCase() in ['male', 'female'])) {
            error "Sample '${row.sample}': sex '${row.sex}' is not recognised. Use 'male' or 'female'."
        }
        if ((row.fastq_1 || ((row.bam || row.cram) && !row.vcf)) && !row.sex) {
            error "Sample '${row.sample}': DeepVariant calls this row and needs the sample's sex: a male sample " +
                  "is called haploid on chrX and chrY outside the pseudoautosomal regions. Add a 'sex' column " +
                  "(male or female)."
        }
        if (tools_list.contains('expansion_hunter') && (row.bam || row.cram) && !row.sex) {
            error "Sample '${row.sample}': expansion_hunter needs the sample's sex to genotype chrX loci. " +
                  "Add a 'sex' column (male or female) to the samplesheet, or remove 'expansion_hunter' from --tools."
        }
    }

    // ─── Pedigree (opt-in: --pedigree) ──────────────────────────────────
    // The family workflow runs only when --pedigree is given, and only on
    // the samples of the samplesheet the pedigree names: two or more, each
    // with a gVCF (DEEPVARIANT's, or the row's gvcf column). --pedigree may
    // be a families file (assets/families.example.tsv), a pedigree with a
    // leading family column: one family runs, the one --family_name names
    // (needed when the file holds more than one). Here only the family and
    // sample columns are read; FAMILY_PREPARE (bin/family_report.py prepare)
    // checks the rest before the joint call starts. Read with containsKey,
    // so a run without the option is not told about an undefined parameter.
    def pedigree   = params.containsKey('pedigree') ? params.pedigree : null
    // The family's folder under --outdir: --family_name, else 'family'
    def family_dir = params.family_name ?: 'family'
    def family_ids = []
    def couple_ids = []
    if (pedigree) {
        def ped_lines = file(pedigree, checkIfExists: true).readLines().findAll { l -> l.trim() }
        def ped_cols  = ped_lines ? ped_lines[0].replaceFirst(/^#/, '').split('\t').collect { c -> c.trim().toLowerCase() } : []
        def si        = ped_cols.indexOf('sample')
        if (si < 0) {
            error "--pedigree ${pedigree}: the first line must be a header with a 'sample' column " +
                  "(sample, sex, father, mother, partner, tab separated, after an optional family column; see " +
                  "assets/pedigree.example.tsv and assets/families.example.tsv)."
        }
        if (!(family_dir ==~ /^[a-zA-Z0-9._-]+$/) || family_dir.startsWith('-')) {
            error "--family_name '${family_dir}' names the family's folder: use only a-z, A-Z, 0-9, '.', '_', '-', " +
                  "and do not start it with '-'."
        }
        def ped_rows = ped_lines.drop(1).findAll { l -> !l.startsWith('#') }
            .collect { l -> l.split('\t', -1).collect { c -> c.trim() } }
        def fi = ped_cols.indexOf('family')
        if (fi >= 0) {
            def families = ped_rows.collect { v -> fi < v.size() ? v[fi] : '' }.findAll { f -> f && f != '.' }.unique().sort()
            if (params.family_name && !families.contains(params.family_name)) {
                error "--family_name ${params.family_name}: ${pedigree} has no such family; its families are: " +
                      "${families ? families.join(', ') : 'none'}."
            }
            if (!params.family_name && families.size() > 1) {
                error "--pedigree ${pedigree} holds ${families.size()} families (${families.join(', ')}): name the " +
                      "one to run with --family_name."
            }
            def chosen = params.family_name ?: (families ? families[0] : null)
            if (chosen) {
                ped_rows = ped_rows.findAll { v -> fi < v.size() && v[fi] == chosen }
            }
        }
        def ped_samples = ped_rows.collect { v -> si < v.size() ? v[si] : '' }
        family_ids = ped_samples.findAll { s -> seen_samples.contains(s) }.unique()
        // Partners (the partner column) among those samples: the couple
        // carrier screen (step 39) runs inside FAMILY when there are two.
        def pi = ped_cols.indexOf('partner')
        couple_ids = ped_rows
            .findAll { v -> pi >= 0 && si < v.size() && pi < v.size() && family_ids.contains(v[si]) && family_ids.contains(v[pi]) }
            .collectMany { v -> [v[si], v[pi]] }
            .unique()
            .sort()
        if (family_ids.size() < 2) {
            error "--pedigree ${pedigree} names ${family_ids.size()} sample${family_ids.size() == 1 ? '' : 's'} of " +
                  "${params.input} (${family_ids ? family_ids.join(', ') : 'none'}): the family analysis needs two " +
                  "or more. Add the other members' rows to the samplesheet, or leave --pedigree out."
        }
        if (seen_samples.contains(family_dir)) {
            error "A sample is named '${family_dir}', the folder of <outdir> where --pedigree writes the family " +
                  "results. Rename the sample, or give the family another name with --family_name."
        }
        if (!params.somalier_sites) {
            error "--pedigree checks the declared relationships with somalier and needs --somalier_sites " +
                  "(scripts/setup.sh --sample-qc-data installs reference/somalier/sites.hg38.vcf.gz)."
        }
        samplesheet_rows.findAll { row -> family_ids.contains(row.sample) && row.vcf && !row.gvcf }.each { row ->
            error "Sample '${row.sample}' is in --pedigree, and the family joint call reads each member's gVCF, but " +
                  "its row gives a vcf without a gvcf. Add the gvcf and gvcf_index columns (step 03 and DEEPVARIANT " +
                  "write <sample>.g.vcf.gz), or give the row's BAM or FASTQ instead so DeepVariant calls it here."
        }
    }

    ch_rows = Channel
        .fromList(samplesheet_rows)
        .map { row -> [[id: row.sample, sex: row.sex ? row.sex.toLowerCase() : null], row] }

    ch_fastq = ch_rows
        .filter { meta, row -> row.fastq_1 }
        .map { meta, row -> [meta, file(row.fastq_1, checkIfExists: true), file(row.fastq_2, checkIfExists: true)] }
    ch_bam_call = ch_rows
        .filter { meta, row -> row.bam && !row.vcf }
        .map { meta, row -> [meta, file(row.bam, checkIfExists: true), file(row.bam_index, checkIfExists: true)] }
    ch_bam_given = ch_rows
        .filter { meta, row -> row.bam && row.vcf }
        .map { meta, row -> [meta, file(row.bam, checkIfExists: true), file(row.bam_index, checkIfExists: true)] }
    ch_vcf_given = ch_rows
        .filter { meta, row -> row.vcf }
        .map { meta, row -> [meta, file(row.vcf, checkIfExists: true), file(row.vcf_index, checkIfExists: true)] }
    ch_gvcf_given = ch_rows
        .filter { meta, row -> row.gvcf }
        .map { meta, row -> [meta, file(row.gvcf, checkIfExists: true), file(row.gvcf_index, checkIfExists: true)] }

    // ─── Reference genome ───────────────────────────────────────────────
    ch_reference      = Channel.value(file(params.reference, checkIfExists: true))
    ch_reference_fai  = Channel.value(file("${params.reference}.fai", checkIfExists: true))
    ch_par_bed        = Channel.value(file("${projectDir}/assets/par_grch38.bed", checkIfExists: true))

    // ─── CRAM rows ──────────────────────────────────────────────────────
    // CRAM_TO_BAM writes each CRAM out as a BAM in the work directory, checked
    // against the CRAM (samtools flagstat), and the row goes on as a BAM row:
    // every BAM step then reads it as it reads any BAM. The decoding needs the
    // reference the CRAM was written with; samtools stops on another one.
    def cram_ids = samplesheet_rows.findAll { row -> row.cram }.collect { row -> row.sample } as Set
    def cram_vcf_ids = samplesheet_rows.findAll { row -> row.cram && row.vcf }.collect { row -> row.sample } as Set
    ch_cram = ch_rows
        .filter { meta, row -> row.cram }
        .map { meta, row -> [meta, file(row.cram, checkIfExists: true), file(row.crai, checkIfExists: true)] }
    CRAM_TO_BAM(ch_cram, ch_reference, ch_reference_fai)
    ch_bam_call  = ch_bam_call.mix(CRAM_TO_BAM.out.bam.filter { meta, bam, bai -> !cram_vcf_ids.contains(meta.id) })
    ch_bam_given = ch_bam_given.mix(CRAM_TO_BAM.out.bam.filter { meta, bam, bai -> cram_vcf_ids.contains(meta.id) })

    // ═══════════════════════════════════════════════════════════════════
    // WORKFLOW 0: UPSTREAM — FASTQ to BAM, sex check, BAM to VCF and gVCF
    // ═══════════════════════════════════════════════════════════════════
    UPSTREAM(
        ch_fastq,
        ch_bam_call,
        ch_bam_given,
        ch_reference,
        ch_reference_fai,
        ch_par_bed
    )

    // Every sample's VCF: the one given, or the one DeepVariant called; and
    // its gVCF, when the samplesheet gives one or DeepVariant wrote one
    ch_vcf_input = ch_vcf_given.mix(UPSTREAM.out.vcf)
    ch_gvcf      = ch_gvcf_given.mix(UPSTREAM.out.gvcf)

    // ─── Input check ────────────────────────────────────────────────────
    // VCF_PRECHECK reads each VCF once before any analysis. Three problems
    // stop the run, with the fix in the message:
    //   - more than one sample column (a joint-called family VCF): the steps
    //     would mix people. VCF_PRECHECK itself stops on it, with the sample
    //     count, the first names and the `bcftools view -s` command that keeps
    //     one sample;
    //   - no contig is chr-named (1, MT): the mito haplogroup comes out empty
    //     and chrX leaks into the ROH summary, both with exit 0;
    //   - a gVCF (or a file named like one) with pharmcat selected: PharmCAT
    //     refuses it.
    // A file not on GRCh38 stops in VCF_PRECHECK too: by its ##contig lengths,
    // or, when the header gives none, by the REF bases of its first records,
    // compared with --reference.
    // Downstream filters keep FILTER=PASS only. A VCF with no PASS record
    // (FILTER '.' everywhere) would give zero hits in every step, so stop,
    // or with --allow_unfiltered use a copy where '.' reads as PASS.
    VCF_PRECHECK(ch_vcf_input, ch_reference, ch_reference_fai)

    ch_vcf_checked = ch_vcf_input
        .map { meta, vcf, idx -> [meta.id, meta, vcf, idx] }
        .join(VCF_PRECHECK.out.status.map { meta, status, counts, contig_style, contigs_seen, gvcf ->
            [meta.id, status, counts, contig_style, contigs_seen, gvcf]
        })
        .map { id, meta, vcf, idx, status, counts, contig_style, contigs_seen, gvcf ->
            if (contig_style == 'other') {
                error "Sample '${id}': no contig in ${vcf.name} is named the chr way (it has ${contigs_seen}...). " +
                      "The pipeline needs GRCh38 contigs named chr1 to chr22, chrX, chrY and chrM; without them the " +
                      "mito haplogroup comes out empty and the ROH summary is wrong. Rename them, then point the " +
                      "samplesheet at the new file:\n" +
                      '    for c in $(seq 1 22) X Y; do echo "$c chr$c"; done > chr_map.txt\n' +
                      '    echo "MT chrM" >> chr_map.txt\n' +
                      "    bcftools annotate --rename-chrs chr_map.txt -Oz -o ${id}.chr.vcf.gz ${vcf.name}\n" +
                      "    bcftools index -t ${id}.chr.vcf.gz\n" +
                      "Contigs outside the map (unplaced scaffolds) keep their names, and the ClinVar screen leaves " +
                      "them out. A GRCh37 file needs more than a rename: see docs/vcf-first.md."
            }
            if (gvcf == 'blocks' && tools_list.contains('pharmcat')) {
                error "Sample '${id}': ${vcf.name} is a gVCF (it has reference blocks: ALT <*>, <NON_REF> or '.' " +
                      "with INFO/END, or a ##GVCFBlock header line), and PharmCAT refuses a gVCF. Remove the " +
                      "reference blocks and give the file a name without .g.vcf or .genomic.vcf (docs/vcf-first.md " +
                      "has the commands), or run without pharmcat and cpic."
            }
            if (gvcf == 'name' && tools_list.contains('pharmcat')) {
                error "Sample '${id}': ${vcf.name} has no reference blocks, but its name contains .g.vcf or " +
                      ".genomic.vcf, and PharmCAT refuses any file named so. Rename the file and its index, for " +
                      "example to ${id}.vcf.gz and ${id}.vcf.gz.tbi, and point the samplesheet at them."
            }
            [id, meta, vcf, idx, status, counts]
        }
        .branch { id, meta, vcf, idx, status, counts ->
            pass:       status == 'pass'
            unfiltered: status == 'unfiltered'
            no_pass:    true
        }

    ch_vcf = ch_vcf_checked.pass
        .map { id, meta, vcf, idx, status, counts -> [meta, vcf, idx] }
        .mix(
            ch_vcf_checked.unfiltered
                .map { id, meta, vcf, idx, status, counts -> [id, meta, vcf.name, counts] }
                .join(VCF_PRECHECK.out.relaxed.map { meta, vcf, idx -> [meta.id, vcf, idx] })
                .map { id, meta, name, counts, vcf, idx ->
                    log.warn "Sample '${id}': no PASS record in ${name} (${counts}); --allow_unfiltered is set, " +
                             "so records with FILTER '.' are used as PASS."
                    [meta, vcf, idx]
                },
            ch_vcf_checked.no_pass
                .map { id, meta, vcf, idx, status, counts ->
                    def remedy = params.allow_unfiltered
                        ? "Filter the VCF with your caller's recommended filters; --allow_unfiltered cannot " +
                          "help here, because no record has FILTER '.'."
                        : "Filter the VCF with your caller's recommended filters, or rerun with " +
                          "--allow_unfiltered to treat FILTER '.' as PASS."
                    error "Sample '${id}': no record in ${vcf.name} has FILTER=PASS (${counts}). " +
                          "Every PASS-only step would report zero hits. ${remedy}"
                }
        )

    // Every BAM, given or aligned here, once INDEXCOV has checked its sex
    ch_bam = UPSTREAM.out.bam

    // ─── Sequence dictionary ────────────────────────────────────────────
    // Only MITO_VARIANTS and CHIP_MUTECT2 (GATK) read the sequence dictionary
    ch_reference_dict = Channel.value([])
    if (tools_list.contains('mito_variants') || tools_list.contains('chip')) {
        if (!(params.reference ==~ /.*\.(fasta|fa|fna)$/)) {
            error "mito_variants and chip need the reference's .dict next to it, found by replacing the FASTA extension; " +
                  "--reference ${params.reference} does not end in .fasta, .fa or .fna."
        }
        ch_reference_dict = Channel.value(
            file(params.reference.replaceAll(/\.(fasta|fa|fna)$/, '.dict'), checkIfExists: true)
        )
    }

    // ─── Optional reference databases ───────────────────────────────────
    // Empty list [] = "no file" — standard Nextflow pattern for optional path inputs.
    // Processes check truthiness (e.g., `if (myfile)`) to skip absent databases.

    // ClinVar
    ch_clinvar       = params.clinvar       ? Channel.value(file(params.clinvar, checkIfExists: true))       : Channel.value([])
    ch_clinvar_index = params.clinvar_index  ? Channel.value(file(params.clinvar_index, checkIfExists: true)) : Channel.value([])

    // VEP caches
    ch_vep_cache      = Channel.value(params.vep_cache      ? file(params.vep_cache, checkIfExists: true)      : [])
    ch_vep_cache_cpsr = Channel.value(params.vep_cache_cpsr ? file(params.vep_cache_cpsr, checkIfExists: true) : [])

    // PCGR/CPSR data bundle
    ch_pcgr_data = Channel.value(params.pcgr_data ? file(params.pcgr_data, checkIfExists: true) : [])

    // PyPGx bundle
    ch_pypgx_bundle = Channel.value(params.pypgx_bundle ? file(params.pypgx_bundle, checkIfExists: true) : [])

    // Annotation score databases (CADD, SpliceAI, REVEL, AlphaMissense)
    ch_cadd_snv             = Channel.value(params.cadd_snv             ? file(params.cadd_snv, checkIfExists: true)             : [])
    ch_cadd_snv_index       = Channel.value(params.cadd_snv_index       ? file(params.cadd_snv_index, checkIfExists: true)       : [])
    ch_cadd_indel           = Channel.value(params.cadd_indel           ? file(params.cadd_indel, checkIfExists: true)           : [])
    ch_cadd_indel_index     = Channel.value(params.cadd_indel_index     ? file(params.cadd_indel_index, checkIfExists: true)     : [])
    ch_spliceai_snv         = Channel.value(params.spliceai_snv         ? file(params.spliceai_snv, checkIfExists: true)         : [])
    ch_spliceai_snv_index   = Channel.value(params.spliceai_snv_index   ? file(params.spliceai_snv_index, checkIfExists: true)   : [])
    ch_spliceai_indel       = Channel.value(params.spliceai_indel       ? file(params.spliceai_indel, checkIfExists: true)       : [])
    ch_spliceai_indel_index = Channel.value(params.spliceai_indel_index ? file(params.spliceai_indel_index, checkIfExists: true) : [])
    ch_revel                = Channel.value(params.revel                ? file(params.revel, checkIfExists: true)                : [])
    ch_revel_index          = Channel.value(params.revel_index          ? file(params.revel_index, checkIfExists: true)          : [])
    ch_alphamissense        = Channel.value(params.alphamissense        ? file(params.alphamissense, checkIfExists: true)        : [])
    ch_alphamissense_index  = Channel.value(params.alphamissense_index  ? file(params.alphamissense_index, checkIfExists: true)  : [])
    ch_gnomad_constraint    = Channel.value(params.gnomad_constraint    ? file(params.gnomad_constraint, checkIfExists: true)    : [])

    // PGS scoring & ancestry reference
    ch_pgs_scoring  = Channel.value(params.pgs_scoring  ? file(params.pgs_scoring, checkIfExists: true)  : [])
    ch_ancestry_ref = Channel.value(params.ancestry_ref ? file(params.ancestry_ref, checkIfExists: true) : [])
    // The panel's GRCh38 SNVs, which setup.sh --ancestry-panel writes beside
    // it: PRS_SCORE_SITES genotypes them from the gVCF, so the projection
    // sees the sites where the sample matches the reference too.
    def ancestry_sites = params.ancestry_ref ? file(params.ancestry_ref.toString().replaceFirst(/\.tar\.zst$/, '') + '_GRCh38_sites.tsv') : null
    if (ancestry_sites && !ancestry_sites.exists() && !workflow.stubRun) {
        error "ERROR: --ancestry_ref ${params.ancestry_ref} has no site list beside it (${ancestry_sites}); run scripts/setup.sh --ancestry-panel <genome_dir> again."
    }
    ch_ancestry_sites = Channel.value(ancestry_sites && ancestry_sites.exists() ? ancestry_sites : [])
    // The catalog's trait label of each PGS id (scripts/ci/check-pgs-labels.sh checks them).
    ch_pgs_labels = Channel.value(file("${projectDir}/assets/pgs_scores.tsv", checkIfExists: true))
    if (params.pgsc_calc && !file("${params.pgsc_calc}/main.nf").exists()) {
        error "ERROR: --pgsc_calc ${params.pgsc_calc} holds no main.nf; it is the pgsc_calc checkout setup.sh makes (tools/pgsc_calc-<release>)."
    }

    // Step 26b's panel folder (setup.sh --ancestry-fine-panel): its model/DONE
    // is written last, so a folder without it is an install that stopped.
    ch_ancestry_fine_panel = Channel.value(params.ancestry_fine_panel ? file(params.ancestry_fine_panel, checkIfExists: true) : [])
    if (params.ancestry_fine_panel && !workflow.stubRun && !file("${params.ancestry_fine_panel}/model/DONE").exists()) {
        error "ERROR: --ancestry_fine_panel ${params.ancestry_fine_panel} is not a finished panel (no model/DONE); run scripts/setup.sh --ancestry-fine-panel <genome_dir> again."
    }

    // Step 44's panel folder (setup.sh --ancient-panel): its DONE is written
    // last, once its axes were checked, so a folder without it is an install
    // that stopped.
    ch_ancient_panel = Channel.value(params.ancient_panel ? file(params.ancient_panel, checkIfExists: true) : [])
    if (params.ancient_panel && !workflow.stubRun && !file("${params.ancient_panel}/DONE").exists()) {
        error "ERROR: --ancient_panel ${params.ancient_panel} is not a finished panel (no DONE); run scripts/setup.sh --ancient-panel <genome_dir> again."
    }

    // ExpansionHunter variant catalog
    ch_expansion_catalog = Channel.value(params.expansion_catalog ? file(params.expansion_catalog, checkIfExists: true) : [])

    // HLA reference database (IPD-IMGT/HLA hla.dat) and the gene annotation
    // T1K takes the genes' GRCh38 coordinates from
    ch_hla_dat   = Channel.value(params.hla_dat   ? file(params.hla_dat, checkIfExists: true)   : [])
    ch_hla_genes = Channel.value(params.hla_genes ? file(params.hla_genes, checkIfExists: true) : [])

    // AnnotSV annotation directory (the biocontainer ships no annotation data)
    ch_annotsv_annotations = Channel.value(params.annotsv_annotations ? file(params.annotsv_annotations, checkIfExists: true) : [])

    // Delly exclude map (regions skipped by delly sr -x)
    ch_delly_exclude = Channel.value(params.delly_exclude ? file(params.delly_exclude, checkIfExists: true) : [])

    // Manta call regions (configManta.py --callRegions): a bgzipped BED with its .tbi beside it
    ch_manta_call_regions       = Channel.value(params.manta_call_regions ? file(params.manta_call_regions, checkIfExists: true) : [])
    ch_manta_call_regions_index = Channel.value(params.manta_call_regions ? file("${params.manta_call_regions}.tbi", checkIfExists: true) : [])

    // UCSC GRCh38 chromosome bands for TelomereHunter (-b)
    ch_cytoband = Channel.value(params.cytoband ? file(params.cytoband, checkIfExists: true) : [])

    // Sample identity and contamination (sample_qc): somalier's sites VCF and
    // the folder of VerifyBamID2's marker panel (.UD, .mu, .bed)
    ch_somalier_sites     = Channel.value(params.somalier_sites     ? file(params.somalier_sites, checkIfExists: true)     : [])
    ch_verifybamid2_panel = Channel.value(params.verifybamid2_panel ? file(params.verifybamid2_panel, checkIfExists: true) : [])

    // Opt-in steps: Cyrius's install (setup.sh --cyrius), IPD-KIR's kir.dat
    // (--kir), Parascopy's homology table and models, and its background
    // windows for a BAM that covers part of the genome
    ch_cyrius_install = Channel.value(params.cyrius_install ? file(params.cyrius_install, checkIfExists: true) : [])
    ch_kir_dat        = Channel.value(params.kir_dat ? file(params.kir_dat, checkIfExists: true) : [])
    ch_parascopy_data = Channel.value(params.parascopy_data ? file(params.parascopy_data, checkIfExists: true) : [])
    ch_parascopy_bed  = Channel.value(params.parascopy_depth_bed ? file(params.parascopy_depth_bed, checkIfExists: true) : [])
    ch_yleaf_data     = Channel.value(params.yleaf_data ? file(params.yleaf_data, checkIfExists: true) : [])
    // GENCODE's exons (setup.sh: reference/gencode.v<release>.basic.exons.bed):
    // the cDNA read check of CLINVAR_SCREEN's hits and of CPSR's P/LP rows,
    // and the couple screen's deletions over an exon
    ch_gencode_exons  = Channel.value(params.gencode_exons ? file(params.gencode_exons, checkIfExists: true) : [])
    // Step 20's annotation of the chrM calls (MITO_VARIANTS_ANNOTATE): the
    // whole ClinVar VCF with chr names, its .tbi beside it (setup.sh installs
    // clinvar/clinvar_chr.vcf.gz), and gnomAD v3.1's chrM sites (setup.sh
    // installs gnomad_mito/, GNOMAD_MITO_URL in versions.env). Without them
    // the table names the missing file and leaves its columns '.'.
    ch_mito_clinvar       = Channel.value(params.mito_clinvar ? file(params.mito_clinvar, checkIfExists: true) : [])
    ch_mito_clinvar_index = Channel.value(params.mito_clinvar ? file("${params.mito_clinvar}.tbi", checkIfExists: true) : [])
    ch_gnomad_mito        = Channel.value(params.gnomad_mito ? file(params.gnomad_mito, checkIfExists: true) : [])

    // ═══════════════════════════════════════════════════════════════════
    // WORKFLOW 1: BAM_ANALYSIS — HLA (and KIR), STR, telomere, coverage,
    // mito, SMN1/SMN2, paralog genes (HiFi), sample QC, Y haplogroup. First, because PGX reads the HLA types.
    // ═══════════════════════════════════════════════════════════════════
    BAM_ANALYSIS(
        ch_bam,
        ch_reference,
        ch_reference_fai,
        ch_reference_dict,
        ch_expansion_catalog,
        ch_hla_dat,
        ch_hla_genes,
        ch_cytoband,
        ch_somalier_sites,
        ch_verifybamid2_panel,
        ch_kir_dat,
        ch_parascopy_data,
        ch_parascopy_bed,
        ch_yleaf_data,
        ch_vcf,
        ch_mito_clinvar,
        ch_mito_clinvar_index,
        ch_gnomad_mito
    )

    // ═══════════════════════════════════════════════════════════════════
    // WORKFLOW 2: PGX — Pharmacogenomics & ClinVar screening; PharmCAT gets
    // the HLA types, an agreed CYP2D6 call and MT-RNR1 from the chrM calls
    // as outside calls
    // ═══════════════════════════════════════════════════════════════════
    // Opt-in (--medications): the user's own list of the drugs taken, never a
    // default. CPIC_LOOKUP checks it for inhibitors of CYP2D6, CYP2C19, CYP2C9
    // and CYP3A4/5.
    ch_medications = params.medications
        ? Channel.value(file(params.medications, checkIfExists: true))
        : Channel.value([])
    ch_inhibitors  = Channel.value(file("${projectDir}/assets/pgx_inhibitors.tsv", checkIfExists: true))
    // The drugs a pathogenic ClinVar hit in a disease gene makes risky (POLG with valproate).
    ch_contraindications = Channel.value(file("${projectDir}/assets/pgx_contraindications.tsv", checkIfExists: true))
    PGX(
        ch_vcf,
        ch_reference,
        ch_reference_fai,
        ch_clinvar,
        ch_clinvar_index,
        ch_bam,
        ch_pypgx_bundle,
        ch_gvcf,
        BAM_ANALYSIS.out.hla_alleles,
        ch_cyrius_install,
        ch_gencode_exons,
        BAM_ANALYSIS.out.mito_vcf,
        BAM_ANALYSIS.out.hla_optitype,
        ch_medications,
        ch_inhibitors,
        ch_contraindications
    )

    // ═══════════════════════════════════════════════════════════════════
    // WORKFLOW 3: ANNOTATION — VEP → vcfanno → slivar / clinical_filter
    // ═══════════════════════════════════════════════════════════════════
    ANNOTATION(
        ch_vcf,
        ch_reference,
        ch_reference_fai,
        ch_vep_cache,
        ch_cadd_snv,
        ch_cadd_snv_index,
        ch_cadd_indel,
        ch_cadd_indel_index,
        ch_spliceai_snv,
        ch_spliceai_snv_index,
        ch_spliceai_indel,
        ch_spliceai_indel_index,
        ch_revel,
        ch_revel_index,
        ch_alphamissense,
        ch_alphamissense_index,
        ch_gnomad_constraint,
        ch_clinvar,
        ch_clinvar_index,
        ch_bam
    )

    // ═══════════════════════════════════════════════════════════════════
    // WORKFLOW 4: CLINICAL — CPSR, ROH, PRS and ancestry (pgsc_calc), fine-grained ancestry, mito
    // haplogroup (from the Mutect2 chrM calls of BAM_ANALYSIS when they exist)
    // ═══════════════════════════════════════════════════════════════════
    CLINICAL(
        ch_vcf,
        ch_pcgr_data,
        ch_vep_cache_cpsr,
        ch_pgs_scoring,
        ch_ancestry_ref,
        ch_ancestry_sites,
        ch_pgs_labels,
        ch_gvcf,
        ch_reference,
        ch_reference_fai,
        BAM_ANALYSIS.out.mito_vcf,
        ch_ancestry_fine_panel,
        ch_bam,
        ch_gencode_exons
    )

    // ─── Traits panel (traits, step 42) ─────────────────────────────────
    // The genotype at each position of assets/traits.tsv, read from the gVCF
    // only: a sample without one is skipped with a warning, because its
    // variant-only VCF cannot tell a reference match from a position
    // without reads.
    ch_traits          = Channel.empty()
    ch_genotype_rules  = Channel.empty()   // [meta, <id>_genotype_rules.tsv], for the report
    ch_traits_versions = Channel.empty()
    if (tools_list.contains('traits')) {
        ch_traits_table = Channel.value(file("${projectDir}/assets/traits.tsv", checkIfExists: true))
        ch_vcf
            .map { meta, _vcf, _idx -> [meta.id, meta] }
            .join(ch_gvcf.map { meta, _gvcf, _gidx -> [meta.id, true] }, remainder: true)
            .filter { _id, meta, has_gvcf -> meta != null && !has_gvcf }
            .subscribe { id, _meta, _has_gvcf ->
                log.warn "traits skipped for ${id}: it reads the gVCF, and the samplesheet row has none " +
                         "(add the gvcf and gvcf_index columns, or start from the BAM or FASTQ)."
            }
        // assets/genotype_rules.tsv reads several sites together (APOE,
        // TAS2R38); TRAITS_PHASE phases the rule sites that need it from the
        // BAM, and a sample without a BAM gets [] for the phase.
        ch_traits_rules = Channel.value(file("${projectDir}/assets/genotype_rules.tsv", checkIfExists: true))
        TRAITS_SITES(ch_gvcf, ch_traits_table, ch_traits_rules, ch_reference, ch_reference_fai)
        ch_traits_genotypes = TRAITS_SITES.out.genotypes.map { meta, g -> [meta.id, meta, g] }
        TRAITS_PHASE(
            ch_traits_genotypes
                .join(ch_bam.map { meta, bam, bai -> [meta.id, bam, bai] })
                .map { _id, meta, g, bam, bai -> [meta, g, bam, bai] },
            ch_traits_rules, ch_reference, ch_reference_fai
        )
        TRAITS(
            ch_traits_genotypes
                .join(TRAITS_PHASE.out.phased.map { meta, f -> [meta.id, f] }, remainder: true)
                .filter { row -> row[1] != null }
                .map { row -> [row[1], row[2], row[3] ?: []] },
            ch_traits_table, ch_traits_rules
        )
        ch_traits          = TRAITS.out.traits
        ch_genotype_rules  = TRAITS.out.rules
        ch_traits_versions = TRAITS_SITES.out.versions.mix(TRAITS_PHASE.out.versions, TRAITS.out.versions)
    }

    // ─── Ancient populations (opt-in: ancient, step 44) ─────────────────
    // The sample projected onto the panel's present-day West Eurasian PCs,
    // beside the AADR's ancient people, from the gVCF only: a sample without
    // one is skipped with a warning. Equivalent to scripts/44-ancient.sh.
    ch_ancient          = Channel.empty()
    ch_ancient_info     = Channel.empty()
    ch_ancient_svg      = Channel.empty()
    ch_ancient_versions = Channel.empty()
    if (tools_list.contains('ancient')) {
        ch_vcf
            .map { meta, _vcf, _idx -> [meta.id, meta] }
            .join(ch_gvcf.map { meta, _gvcf, _gidx -> [meta.id, true] }, remainder: true)
            .filter { _id, meta, has_gvcf -> meta != null && !has_gvcf }
            .subscribe { id, _meta, _has_gvcf ->
                log.warn "ancient skipped for ${id}: it reads the gVCF, and the samplesheet row has none " +
                         "(add the gvcf and gvcf_index columns, or start from the BAM or FASTQ)."
            }
        ANCIENT_SITES(ch_gvcf, ch_ancient_panel, ch_reference, ch_reference_fai)
        ANCIENT_GENO(ANCIENT_SITES.out.genotypes, ch_ancient_panel)
        ANCIENT_PCA(ANCIENT_GENO.out.merged, ch_ancient_panel)
        ANCIENT_REPORT(ANCIENT_PCA.out.evec, ch_ancient_panel)
        ch_ancient          = ANCIENT_REPORT.out.table
        ch_ancient_info     = ANCIENT_REPORT.out.info
        ch_ancient_svg      = ANCIENT_REPORT.out.svg
        ch_ancient_versions = ANCIENT_SITES.out.versions.mix(ANCIENT_GENO.out.versions, ANCIENT_PCA.out.versions,
                                                             ANCIENT_REPORT.out.versions)
    }

    // ─── CRAM archive (opt-in: cram_archive) ────────────────────────────
    // The CRAM of every BAM, checked against it, published beside it; a row
    // given as CRAM already has one. The BAM is never deleted here.
    ch_cram_versions = Channel.empty()
    if (tools_list.contains('cram_archive')) {
        CRAM_ARCHIVE(ch_bam.filter { meta, bam, bai -> !cram_ids.contains(meta.id) }, ch_reference, ch_reference_fai)
        ch_cram_versions = CRAM_ARCHIVE.out.versions
    }

    // ═══════════════════════════════════════════════════════════════════
    // WORKFLOW 5: SV — Structural variant calling & annotation (opt-in)
    // ═══════════════════════════════════════════════════════════════════
    SV(
        ch_bam,
        ch_reference,
        ch_reference_fai,
        ch_delly_exclude,
        ch_annotsv_annotations,
        ch_manta_call_regions,
        ch_manta_call_regions_index
    )

    // ═══════════════════════════════════════════════════════════════════
    // WORKFLOW 5b: FAMILY — joint call, relationship check, phasing by
    // inheritance (opt-in: --pedigree naming two or more samples), and the
    // couple carrier screen when it marks two of them as partners. After SV,
    // because the couple screen reads the partners' Delly calls.
    // ═══════════════════════════════════════════════════════════════════
    ch_family_versions = Channel.empty()
    ch_couple_tables   = Channel.empty()
    if (family_ids) {
        if (couple_ids && !tools_list.contains('clinvar')) {
            log.warn "--pedigree marks ${couple_ids.join(' and ')} as partners, but clinvar is not in --tools: " +
                     "the couple carrier screen reads each partner's ClinVar hits, so it does not run."
        }
        def ch_couple_outputs = couple_ids && tools_list.contains('clinvar') ? Channel.empty().mix(
            PGX.out.clinvar_dir, SV.out.delly_vcf, BAM_ANALYSIS.out.expansion_vcf, BAM_ANALYSIS.out.smn_copy_number,
            BAM_ANALYSIS.out.cn_intervals, SV.out.cnvpytor_vcf
        ).filter { meta, f -> couple_ids.contains(meta.id) }.map { meta, f -> f } : Channel.empty()
        FAMILY(
            ch_gvcf.filter { meta, g, idx -> family_ids.contains(meta.id) },
            ANNOTATION.out.clinical_vcf.filter { meta, f -> family_ids.contains(meta.id) },
            Channel.value(file(pedigree, checkIfExists: true)),
            // The partners with ClinVar hits this run (family_plan.tsv's couple row)
            Channel.value(couple_ids && tools_list.contains('clinvar') ? couple_ids : []),
            ch_somalier_sites,
            ch_reference,
            ch_reference_fai,
            ch_par_bed,
            ch_couple_outputs,
            Channel.value(file("${projectDir}/assets/carrier_genes.tsv", checkIfExists: true)),
            ch_gencode_exons,
            // The shared segments phase against step 26b's panel: no panel, no task
            params.ancestry_fine_panel ? ch_ancestry_fine_panel : Channel.empty()
        )
        ch_family_versions = FAMILY.out.versions
        // <id>_couple_carrier.tsv: the id is the name without its 19-character suffix
        ch_couple_tables   = FAMILY.out.couple_tables.flatten()
            .map { f -> [[id: f.name.substring(0, f.name.length() - 19)], f] }
    }

    // ─── Blood groups (opt-in: bloodtype) ───────────────────────────────
    // The antigens of assets/bloodtype_markers.tsv from the VCF, and RhD from
    // the BAM's read depth. A sample without a BAM gets [] for the depth: its
    // RhD is not determined. Equivalent to scripts/41-bloodtype.sh.
    ch_bloodtype          = Channel.empty()
    ch_bloodtype_versions = Channel.empty()
    if (tools_list.contains('bloodtype')) {
        ch_bloodtype_markers = Channel.value(file("${projectDir}/assets/bloodtype_markers.tsv", checkIfExists: true))
        BLOODTYPE_REGIONS()
        BLOODTYPE_DEPTH(ch_bam, ch_bloodtype_markers, BLOODTYPE_REGIONS.out.bed)
        BLOODTYPE(
            ch_vcf
                .map { meta, vcf, idx -> [meta.id, meta, vcf] }
                .join(BLOODTYPE_DEPTH.out.depth.map { meta, f -> [meta.id, f] }, remainder: true)
                .filter { row -> row[1] != null }
                .map { row -> [row[1], row[2], row[3] ?: []] },
            ch_bloodtype_markers
        )
        ch_bloodtype          = BLOODTYPE.out.table
        ch_bloodtype_versions = BLOODTYPE_REGIONS.out.versions.mix(BLOODTYPE_DEPTH.out.versions, BLOODTYPE.out.versions)
    }

    // ─── Archaic ancestry (opt-in: archaic, step 43) ────────────────────
    // IBDmix against the Altai and Vindija Neanderthals and the Altai
    // Denisovan, from the VCF and the table setup.sh --archaic builds; the
    // genes of the segments come from --hla_genes (GENCODE's gene lines).
    // Equivalent to scripts/43-archaic.sh.
    ch_archaic          = Channel.empty()
    ch_archaic_versions = Channel.empty()
    if (tools_list.contains('archaic')) {
        if (!workflow.stubRun && !file("${params.archaic_data}/DONE").exists()) {
            error "ERROR: --archaic_data ${params.archaic_data} is not a finished archaic table (no DONE); run scripts/setup.sh --archaic <genome_dir> again."
        }
        if (!(params.archaic_pop in ['ALL', 'AFR', 'AMR', 'EAS', 'EUR', 'SAS'])) {
            error "--archaic_pop must be ALL, AFR, AMR, EAS, EUR or SAS, got '${params.archaic_pop}'."
        }
        ARCHAIC_GENOTYPES(
            ch_vcf.map { meta, vcf, idx -> [meta, vcf] },
            Channel.value(file(params.archaic_data, checkIfExists: true))
        )
        ARCHAIC_IBDMIX(ARCHAIC_GENOTYPES.out.genotypes)
        ARCHAIC_REPORT(ARCHAIC_IBDMIX.out.segments, ch_hla_genes)
        ch_archaic          = ARCHAIC_REPORT.out.summary
        ch_archaic_versions = ARCHAIC_GENOTYPES.out.versions.mix(ARCHAIC_IBDMIX.out.versions, ARCHAIC_REPORT.out.versions)
    }

    // ─── Lp(a): LPA KIV-2 repeat count (opt-in: lpa) ────────────────────
    // Read depth over the KIV-2 exons against single-copy LPA exons, for each
    // sample with a BAM. Equivalent to scripts/46-lpa.sh.
    ch_lpa          = Channel.empty()
    ch_lpa_versions = Channel.empty()
    if (tools_list.contains('lpa')) {
        LPA_REGIONS()
        LPA_DEPTH(ch_bam, LPA_REGIONS.out.bed)
        LPA(LPA_DEPTH.out.depth)
        ch_lpa          = LPA.out.table
        ch_lpa_versions = LPA_REGIONS.out.versions.mix(LPA_DEPTH.out.versions, LPA.out.versions)
    }

    // ─── MUC1 VNTR frameshifts (opt-in: vntyper, step 52) ───────────────
    // VNtyper on each sample with a BAM: the ADTKD-MUC1 frameshifts no
    // standard caller sees. Equivalent to scripts/52-vntyper.sh.
    ch_vntyper          = Channel.empty()
    ch_vntyper_versions = Channel.empty()
    if (tools_list.contains('vntyper')) {
        VNTYPER(ch_bam)
        ch_vntyper          = VNTYPER.out.table
        ch_vntyper_versions = VNTYPER.out.versions
    }

    // ─── Viral DNA in the unmapped reads (opt-in: viral, step 47) ──────
    // The BAM's unmapped and chrEBV reads realigned to the viral genomes of
    // --viral_ref (setup.sh --viral-data). A sample without a BAM has no
    // unmapped reads to read and is left out. Equivalent to scripts/47-viral.sh.
    ch_viral          = Channel.empty()
    ch_viral_versions = Channel.empty()
    if (tools_list.contains('viral')) {
        ch_viral_ref = Channel.value(file(params.viral_ref, checkIfExists: true))
        VIRAL_READS(ch_bam, ch_viral_ref)
        VIRAL_REPORT(VIRAL_READS.out.coverage, ch_viral_ref)
        ch_viral          = VIRAL_REPORT.out.table
        ch_viral_versions = VIRAL_READS.out.versions.mix(VIRAL_REPORT.out.versions)
    }

    // ─── Symptom-driven ranking (exomiser, step 48: --hpo) ──────────────
    // Exomiser on the VCF of the one sample the --hpo terms describe
    // (--exomiser_proband, or the only sample). Equivalent to
    // scripts/48-exomiser.sh without its family mode, which is bash-only.
    ch_exomiser          = Channel.empty()
    ch_exomiser_terms    = Channel.empty()
    ch_exomiser_versions = Channel.empty()
    if (tools_list.contains('exomiser')) {
        def proband = params.exomiser_proband ?: (seen_samples.size() == 1 ? seen_samples.first() : null)
        if (!proband) {
            error "--hpo describes one person and the samplesheet has ${seen_samples.size()} samples: name the one " +
                  "with --exomiser_proband (${seen_samples.join(', ')})."
        }
        if (!seen_samples.contains(proband)) {
            error "--exomiser_proband ${proband} is not a sample of ${params.input} (${seen_samples.join(', ')})."
        }
        ch_exomiser_data = Channel.value(file(params.exomiser_data, checkIfExists: true))
        EXOMISER_PREPARE(
            ch_vcf.filter { meta, _vcf, _idx -> meta.id == proband }.map { meta, vcf, _idx -> [meta, vcf] },
            Channel.value(params.hpo.toString().split(/[,\s]+/).findAll { it }.join(',')),
            ch_exomiser_data
        )
        EXOMISER(EXOMISER_PREPARE.out.analysis, ch_exomiser_data)
        EXOMISER_REPORT(
            EXOMISER.out.variants.join(EXOMISER_PREPARE.out.terms)
        )
        ch_exomiser          = EXOMISER_REPORT.out.top
        ch_exomiser_terms    = EXOMISER_REPORT.out.terms
        ch_exomiser_versions = EXOMISER_PREPARE.out.versions.mix(EXOMISER.out.versions, EXOMISER_REPORT.out.versions)
    }

    // ─── Clonal haematopoiesis (opt-in: chip, step 49) ───────────────────
    // Mutect2 tumour-only on the coding exons of the CHIP genes'
    // MANE Select transcripts, from the BAM; a sample without one is
    // skipped with a warning. Equivalent to scripts/49-chip.sh.
    ch_chip          = Channel.empty()
    ch_chip_depth    = Channel.empty()
    ch_chip_versions = Channel.empty()
    if (tools_list.contains('chip')) {
        ch_chip_genes = Channel.value(file("${projectDir}/assets/chip_genes.tsv", checkIfExists: true))
        // gnomAD's AF-only germline resource and a panel of normals, each with its .tbi beside it
        def germ = params.chip_germline_resource
        def pon  = params.chip_pon
        def chip_germline       = germ ? file(germ, checkIfExists: true) : []
        def chip_germline_index = germ ? file("${germ}.tbi", checkIfExists: true) : []
        def chip_pon            = pon ? file(pon, checkIfExists: true) : []
        def chip_pon_index      = pon ? file("${pon}.tbi", checkIfExists: true) : []
        ch_vcf
            .map { meta, _vcf, _idx -> [meta.id, meta] }
            .join(ch_bam.map { meta, _bam, _bai -> [meta.id, true] }, remainder: true)
            .filter { _id, meta, has_bam -> meta != null && !has_bam }
            .subscribe { id, _meta, _has_bam ->
                log.warn "chip skipped for ${id}: Mutect2 reads the BAM, and the samplesheet row has none."
            }
        CHIP_PREPARE(Channel.value(file(params.chip_gff, checkIfExists: true)), ch_chip_genes, ch_reference_fai)
        CHIP_MUTECT2(ch_bam, CHIP_PREPARE.out.bed, ch_reference, ch_reference_fai, ch_reference_dict,
                     chip_germline, chip_germline_index, chip_pon, chip_pon_index)
        CHIP_ANNOTATE(
            CHIP_MUTECT2.out.vcf.map { meta, vcf -> [meta.id, meta, vcf] }
                .join(ch_bam.map { meta, bam, bai -> [meta.id, bam, bai] })
                .map { _id, meta, vcf, bam, bai -> [meta, vcf, bam, bai] },
            CHIP_PREPARE.out.gff, CHIP_PREPARE.out.bed, ch_reference, ch_reference_fai
        )
        CHIP_REPORT(CHIP_ANNOTATE.out.tables, ch_chip_genes, CHIP_PREPARE.out.gff, CHIP_PREPARE.out.bed)
        ch_chip          = CHIP_REPORT.out.chip
        ch_chip_depth    = CHIP_REPORT.out.depth
        ch_chip_versions = CHIP_PREPARE.out.versions.mix(CHIP_MUTECT2.out.versions, CHIP_ANNOTATE.out.versions,
                                                         CHIP_REPORT.out.versions)
    }

    // ─── Mobile element insertions (opt-in: mei, step 50) ──────────────
    // MEGAnE on each sample's BAM with the k-mer set of --megane_kmers
    // (setup.sh --megane-kmers; its .mi count file beside it), then the calls
    // and the known-site catalog of assets/mei_sites.tsv. With
    // --annotsv_annotations, MEI_REPORT runs AnnotSV on the calls first.
    // A sample without a BAM is left out. Equivalent to scripts/50-mei.sh.
    ch_mei          = Channel.empty()
    ch_mei_versions = Channel.empty()
    if (tools_list.contains('mei')) {
        def mk = file(params.megane_kmers, checkIfExists: true)
        ch_megane_kmers = Channel.value([mk, file(mk.toString().replaceFirst(/\.mk$/, '.mi'), checkIfExists: true)])
        ch_mei_sites = Channel.value(file("${projectDir}/assets/mei_sites.tsv", checkIfExists: true))
        MEI_DEPTH(ch_bam, ch_mei_sites)
        MEI(ch_bam.join(MEI_DEPTH.out.depth.map { meta, mainchr, depth -> [meta, mainchr] }),
            ch_reference, ch_reference_fai, ch_megane_kmers)
        MEI_REPORT(MEI.out.vcf.join(MEI_DEPTH.out.depth.map { meta, mainchr, depth -> [meta, depth] }), ch_mei_sites,
                   ch_annotsv_annotations)
        ch_mei          = MEI_REPORT.out.tables
        ch_mei_versions = MEI_DEPTH.out.versions.mix(MEI.out.versions, MEI_REPORT.out.versions)
    }

    // ─── Mosaic chromosomal alterations (opt-in: mocha, step 51) ───────
    // The VCF's heterozygous SNVs phased with Beagle against the panel of
    // --mocha_data (setup.sh --mocha-data), then MoChA on chr1 to chr22 and
    // the calls above the floor. Equivalent to scripts/51-mocha.sh.
    ch_mocha          = Channel.empty()
    ch_mocha_summary  = Channel.empty()
    ch_mocha_versions = Channel.empty()
    if (tools_list.contains('mocha')) {
        ch_mocha_data = Channel.value(file(params.mocha_data, checkIfExists: true))
        MOCHA_PREP(ch_vcf, ch_reference, ch_reference_fai)
        MOCHA_PHASE(MOCHA_PREP.out.phase_in, ch_mocha_data)
        MOCHA_CALL(
            MOCHA_PREP.out.prep.map { meta, bcf, csi -> [meta.id, meta, bcf, csi] }
                .join(MOCHA_PHASE.out.phased.map { meta, d -> [meta.id, d] })
                .map { id, meta, bcf, csi, d -> [meta, bcf, csi, d] }
        )
        MOCHA_REPORT(MOCHA_CALL.out.tables)
        ch_mocha          = MOCHA_REPORT.out.tables.map { meta, calls, summary -> [meta, calls] }
        ch_mocha_summary  = MOCHA_REPORT.out.tables.map { meta, calls, summary -> [meta, summary] }
        ch_mocha_versions = MOCHA_PREP.out.versions.mix(MOCHA_PHASE.out.versions, MOCHA_CALL.out.versions,
                                                        MOCHA_REPORT.out.versions)
    }

    // ─── Locityper second opinion at hard loci (opt-in: locityper, step 53)
    // Locityper on each sample's BAM with the database of --locityper_data
    // (setup.sh --locityper-data, built from the same reference). A sample
    // without a BAM is left out. Equivalent to scripts/53-locityper.sh.
    ch_locityper          = Channel.empty()
    ch_locityper_versions = Channel.empty()
    if (tools_list.contains('locityper')) {
        LOCITYPER(ch_bam, ch_reference, ch_reference_fai,
                  Channel.value(file(params.locityper_data, checkIfExists: true)))
        ch_locityper          = LOCITYPER.out.table
        ch_locityper_versions = LOCITYPER.out.versions
    }

    // ═══════════════════════════════════════════════════════════════════
    // WORKFLOW 6: REPORTING — HTML report & MultiQC
    // ═══════════════════════════════════════════════════════════════════

    // Per-sample report inputs: the sample's VCF, and the output files of the
    // steps that ran for it, as one list (HTML_REPORT links each where
    // bin/collect_summary.py reads it). remainder: true keeps a sample a step
    // did not run for, or whose task failed and was ignored; that step's slot
    // is null and drops out of the list. The join also makes the report wait
    // for every selected step it shows, the SV consensus merge included.
    // The ACMG SF rules, the carrier list, the gene callability table, the
    // confirmed mtDNA variants and the two paralog catalogs (step 35's loci,
    // step 35b's regions) the report applies (bin/ is the only repository
    // folder a task sees).
    def report_assets = ["acmg_sf.tsv", "carrier_genes.tsv", "callability_genes.tsv", "mito_confirmed.tsv",
                         "paralog_loci.tsv", "paraphase_genes.tsv"].collect { n -> file("${projectDir}/assets/${n}", checkIfExists: true) }
    ch_report_inputs = ch_vcf
        .map { meta, vcf, idx -> [meta.id, meta, vcf] }
        .join(PGX.out.clinvar_dir.map             { meta, f -> [meta.id, f] }, remainder: true)
        .join(PGX.out.clinvar_cdna.map            { meta, f -> [meta.id, f] }, remainder: true)
        .join(PGX.out.clinvar_cdna_scan.map       { meta, f -> [meta.id, f] }, remainder: true)
        .join(PGX.out.clinvar_deep_intronic.map   { meta, f -> [meta.id, f] }, remainder: true)
        .join(PGX.out.pharmcat_json.map           { meta, f -> [meta.id, f] }, remainder: true)
        .join(PGX.out.pharmcat_html.map           { meta, f -> [meta.id, f] }, remainder: true)
        .join(PGX.out.cpic_phenotypes.map         { meta, f -> [meta.id, f] }, remainder: true)
        .join(PGX.out.cpic_recommendations.map    { meta, f -> [meta.id, f] }, remainder: true)
        .join(ANNOTATION.out.vep_vcf.map          { meta, f -> [meta.id, f] }, remainder: true)
        .join(ANNOTATION.out.vep_vcf_index.map    { meta, f -> [meta.id, f] }, remainder: true)
        .join(ANNOTATION.out.clinical_vcf.map     { meta, f -> [meta.id, f] }, remainder: true)
        .join(ANNOTATION.out.clinical_summary.map { meta, f -> [meta.id, f] }, remainder: true)
        .join(ANNOTATION.out.slivar_vcf.map       { meta, f -> [meta.id, f] }, remainder: true)
        .join(ANNOTATION.out.slivar_summary.map   { meta, f -> [meta.id, f] }, remainder: true)
        .join(ANNOTATION.out.comphet_pairs.map    { meta, f -> [meta.id, f] }, remainder: true)
        .join(CLINICAL.out.cpsr_html.map          { meta, f -> [meta.id, f] }, remainder: true)
        .join(CLINICAL.out.cpsr_tsv.map           { meta, f -> [meta.id, f] }, remainder: true)
        .join(CLINICAL.out.cpsr_cdna.map          { meta, f -> [meta.id, f] }, remainder: true)
        .join(UPSTREAM.out.read_count.map         { meta, f -> [meta.id, f] }, remainder: true)
        .join(CLINICAL.out.roh_regions.map        { meta, f -> [meta.id, f] }, remainder: true)
        .join(CLINICAL.out.haplogroup.map         { meta, f -> [meta.id, f] }, remainder: true)
        .join(CLINICAL.out.haplocheck.map         { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.y_haplogroup.map   { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.y_tree.map         { meta, tree, _path -> [meta.id, tree] }, remainder: true)
        .join(BAM_ANALYSIS.out.y_tree.map         { meta, _tree, path -> [meta.id, path] }, remainder: true)
        .join(BAM_ANALYSIS.out.y_str.map          { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.paraphase_cn.map   { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.paraphase_vars.map { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.smn_copy_number.map { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.hba_copy_number.map { meta, f -> [meta.id, f] }, remainder: true)
        // step 35's other catalog tables, a list per sample (flattened below)
        .join(BAM_ANALYSIS.out.parascopy_tables.map { meta, fs -> [meta.id, fs] }, remainder: true)
        .join(CLINICAL.out.prs_scores.map         { meta, f -> [meta.id, f] }, remainder: true)
        .join(CLINICAL.out.ancestry_results.map   { meta, f -> [meta.id, f] }, remainder: true)
        .join(CLINICAL.out.ancestry_fine.map      { meta, f -> [meta.id, f] }, remainder: true)
        .join(CLINICAL.out.ancestry_fine_png.map  { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.coverage.map       { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.callability.map    { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.cn_intervals.map   { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.sample_qc.map      { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.telomere_results.map { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.mito_vcf.map       { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.mito_annotated.map { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.mito_deletions.map { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.mito_copies.map    { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.hla_alleles.map    { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.hla_optitype.map   { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.hla_concordance.map { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.kir_genotype.map   { meta, f, _release -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.expansion_vcf.map  { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.stranger_vcf.map   { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.rfc1_motif.map     { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.str_motifs.map     { meta, f -> [meta.id, f] }, remainder: true)
        .join(BAM_ANALYSIS.out.eh_catalog.map     { meta, f -> [meta.id, f] }, remainder: true)
        .join(SV.out.manta_vcf.map                { meta, f -> [meta.id, f] }, remainder: true)
        .join(SV.out.delly_vcf.map                { meta, f -> [meta.id, f] }, remainder: true)
        .join(SV.out.cnvpytor_calls.map           { meta, f -> [meta.id, f] }, remainder: true)
        .join(SV.out.merged_sv.map                { meta, f -> [meta.id, f] }, remainder: true)
        .join(SV.out.annotsv_class4_5.map         { meta, f -> [meta.id, f] }, remainder: true)
        .join(PGX.out.pypgx_summary.map           { meta, f -> [meta.id, f] }, remainder: true)
        .join(PGX.out.depth_check.map             { meta, f -> [meta.id, f] }, remainder: true)
        .join(PGX.out.cyrius_results.map          { meta, f -> [meta.id, f] }, remainder: true)
        .join(PGX.out.pgx_consensus.map           { meta, f -> [meta.id, f] }, remainder: true)
        .join(PGX.out.cpic_comparison.map         { meta, f -> [meta.id, f] }, remainder: true)
        .join(UPSTREAM.out.indexcov.map           { meta, f -> [meta.id, f] }, remainder: true)
        .join(CLINICAL.out.consumer_export.map    { meta, f -> [meta.id, f] }, remainder: true)
        .join(ch_traits.map                       { meta, f -> [meta.id, f] }, remainder: true)
        .join(ch_genotype_rules.map               { meta, f -> [meta.id, f] }, remainder: true)
        .join(ch_bloodtype.map                    { meta, f -> [meta.id, f] }, remainder: true)
        .join(ch_archaic.map                      { meta, f -> [meta.id, f] }, remainder: true)
        .join(ch_ancient.map                      { meta, f -> [meta.id, f] }, remainder: true)
        .join(ch_ancient_info.map                 { meta, f -> [meta.id, f] }, remainder: true)
        .join(ch_ancient_svg.map                  { meta, f -> [meta.id, f] }, remainder: true)
        .join(ch_lpa.map                          { meta, f -> [meta.id, f] }, remainder: true)
        .join(ch_vntyper.map                      { meta, f -> [meta.id, f] }, remainder: true)
        .join(ch_viral.map                        { meta, f -> [meta.id, f] }, remainder: true)
        .join(ch_exomiser.map                     { meta, f -> [meta.id, f] }, remainder: true)
        .join(ch_exomiser_terms.map               { meta, f -> [meta.id, f] }, remainder: true)
        .join(ch_chip.map                         { meta, f -> [meta.id, f] }, remainder: true)
        .join(ch_chip_depth.map                   { meta, f -> [meta.id, f] }, remainder: true)
        .join(ch_mei.map                          { meta, calls, sites -> [meta.id, calls] }, remainder: true)
        .join(ch_mei.map                          { meta, calls, sites -> [meta.id, sites] }, remainder: true)
        .join(ch_mocha.map                        { meta, f -> [meta.id, f] }, remainder: true)
        .join(ch_mocha_summary.map                { meta, f -> [meta.id, f] }, remainder: true)
        .join(ch_locityper.map                    { meta, f -> [meta.id, f] }, remainder: true)
        .join(ch_couple_tables.map                { meta, f -> [meta.id, f] }, remainder: true)
        .filter { items -> items[1] != null }
        .map { items -> [items[1], items[2], items[3..-1].collectMany { f -> (f instanceof List ? f : [f]) }.findAll { f -> f != null } + report_assets] }

    // QC files for MultiQC (mosdepth summaries)
    ch_multiqc_files = BAM_ANALYSIS.out.coverage.map { meta, f -> f }

    REPORTING(
        ch_report_inputs,
        ch_multiqc_files
    )

    // ─── Software versions ──────────────────────────────────────────────
    // Every task writes a versions.yml whose first line is the tag of the
    // image it ran in; one file collects them, one block per process.
    Channel.empty()
        .mix(
            UPSTREAM.out.versions,
            VCF_PRECHECK.out.versions,
            PGX.out.versions,
            ANNOTATION.out.versions,
            CLINICAL.out.versions,
            BAM_ANALYSIS.out.versions,
            SV.out.versions,
            REPORTING.out.versions,
            CRAM_TO_BAM.out.versions,
            ch_cram_versions,
            ch_family_versions,
            ch_traits_versions,
            ch_bloodtype_versions,
            ch_archaic_versions,
            ch_ancient_versions,
            ch_lpa_versions,
            ch_vntyper_versions,
            ch_viral_versions,
            ch_chip_versions,
            ch_mei_versions,
            ch_exomiser_versions,
            ch_mocha_versions,
            ch_locityper_versions
        )
        .map { f -> f.text }
        .unique()
        .collectFile(name: 'software_versions.yml', storeDir: "${params.outdir}/pipeline_info", sort: true)

    // ─── Completion handler ─────────────────────────────────────────────
    // Nextflow runs the handler with the script binding's variable map as its
    // delegate, and a map answers null for a name it lacks: `workflow` inside
    // the handler was null ("Cannot get property 'success' on null object" on
    // every run under 25.10), so the completion message never printed. Local
    // variables are resolved where the closure is written, so the handler
    // reads only these.
    //
    // A process with a failed task is named, at the end and in
    // failed_tasks.tsv in the launch directory (process, then 'ignored' for a
    // report-only tool nextflow.config lets fail, or 'failed'): run-all.sh
    // marks that step failed. A run with an ignored task still exits 0, so it
    // does not say "completed successfully".
    //
    // process_status.tsv, beside it, has one row per process of this run,
    // whatever the run's end: how many of its tasks succeeded, came from the
    // cache (-resume, or storeDir), failed for good (not retried, not ignored),
    // failed and were ignored, and were left unfinished (pending, submitted,
    // running or aborted when the run stopped). A process that got no task has
    // a row of zeros. run-all.sh reads it after a failed run, so a step whose
    // tasks completed or were cached is recorded 'ok' and only a step the stop
    // cut short is 'not finished'. It is written last: the lines above do not
    // depend on it.
    def run_info    = workflow
    def run_log     = log
    def outdir      = params.outdir
    def tasks_file  = file("${workflow.launchDir}/failed_tasks.tsv")
    def status_file = file("${workflow.launchDir}/process_status.tsv")
    //
    // A published index (.bai, .crai, .csi, .tbi) older than its data file
    // makes samtools and htslib warn "The index file is older than the data
    // file" on every later query. publishDir copies a task's files in
    // parallel, so the small index is often in place before its BAM or VCF
    // is. Nextflow runs this handler after the last copy has finished. It
    // sets each such index to its data file's time (or now, when later) when
    // the index is a plain file with one link that this run wrote. Any other
    // older index is left as it is and named in a warning: written before
    // this run, beside a data file this run replaced, it may index other
    // data; a hard link shares its time with the task's file in work/, and
    // touching that would make the next -resume run its downstream tasks
    // again. The search covers each sample's folder and the family's under
    // --outdir, six levels down, and skips the folder that holds the work
    // directory (run-all.sh puts it in <sample>/nextflow/). An index named
    // sample.bai or sample.crai belongs to sample.bam or sample.cram when
    // there is no sample.bai's data file.
    def index_ext   = ['.bai', '.crai', '.csi', '.tbi']
    def index_roots = (samplesheet_rows.collect { row -> row.sample } + [params.family_name ?: 'family']).unique()
        .collect { d -> file("${params.outdir}/${d}") }
    def index_skip  = workflow.workDir.toAbsolutePath().normalize().toString() + '/'
    run_info.onComplete {
        try {
            def start_ms = run_info.start.toInstant().toEpochMilli()
            def kept = []
            index_roots
                .findAll { d -> d.getFileSystem() == java.nio.file.FileSystems.getDefault() && d.isDirectory() }
                .collectMany { d -> java.nio.file.Files.list(d).withCloseable { s -> s.iterator().toList() } }
                .findAll { c -> !index_skip.startsWith(c.toAbsolutePath().normalize().toString() + '/') }
                .each { c ->
                    def found = []
                    try {
                        found = java.nio.file.Files.walk(c, 5).withCloseable { s -> s.iterator().toList() }
                    }
                    catch( e ) {
                        run_log.warn "Could not read ${c} to check its index times: ${e}"
                    }
                    found.findAll { p -> index_ext.any { x -> p.getFileName().toString().endsWith(x) } }.each { idx ->
                        def n = idx.getFileName().toString()
                        def data = idx.resolveSibling(n.substring(0, n.lastIndexOf('.')))
                        if (!java.nio.file.Files.exists(data) && (n.endsWith('.bai') || n.endsWith('.crai'))) {
                            data = idx.resolveSibling(n.endsWith('.bai') ? n.substring(0, n.length() - 4) + '.bam' : n.substring(0, n.length() - 5) + '.cram')
                        }
                        if (java.nio.file.Files.exists(data)) {
                            def t_idx  = java.nio.file.Files.getLastModifiedTime(idx)
                            def t_data = java.nio.file.Files.getLastModifiedTime(data)
                            def links = 0
                            try {
                                links = java.nio.file.Files.isRegularFile(idx, java.nio.file.LinkOption.NOFOLLOW_LINKS)
                                    ? java.nio.file.Files.getAttribute(idx, 'unix:nlink', java.nio.file.LinkOption.NOFOLLOW_LINKS) : 0
                            }
                            catch( _e ) {
                                links = 0
                            }
                            if (t_idx < t_data && links == 1 && t_idx.toMillis() >= start_ms) {
                                def now = java.nio.file.attribute.FileTime.fromMillis(System.currentTimeMillis())
                                java.nio.file.Files.setLastModifiedTime(idx, t_data > now ? t_data : now)
                            } else if (t_idx < t_data) {
                                kept << idx.toString()
                            }
                        }
                    }
                }
            if (kept) {
                run_log.warn "Index files older than their data file, left as they are: " +
                    "${kept.take(10).join(', ')}${kept.size() > 10 ? " and ${kept.size() - 10} more" : ''}. " +
                    "Each may index an earlier version of its file, and samtools and htslib warn on every query of it. " +
                    "Index the file again (samtools index for a BAM or CRAM, bcftools index for a VCF), " +
                    "or touch the index when the file was only copied."
            }
        }
        catch( e ) {
            run_log.warn "Could not check the times of the published index files: ${e}"
        }
        def ignored = run_info.stats.processes.findAll { r -> r.ignored > 0 }.collect { r -> r.name }
        def failed  = run_info.stats.processes.findAll { r -> r.failed - r.ignored - r.retries > 0 }.collect { r -> r.name }
        if (ignored || failed) {
            tasks_file.text = 'process\tstatus\n' +
                ignored.collect { n -> "${n}\tignored\n" }.join('') + failed.collect { n -> "${n}\tfailed\n" }.join('')
        } else if (tasks_file.exists()) {
            tasks_file.delete()   // an earlier run's
        }
        if (run_info.success && !ignored) {
            run_log.info ""
            run_log.info "Pipeline completed successfully!"
            run_log.info "Results: ${outdir}"
            run_log.info ""
        } else if (run_info.success) {
            run_log.warn "Pipeline completed, but report-only tools failed and were skipped: ${ignored.join(', ')}. " +
                         "Their results are missing; see .nextflow.log. Results: ${outdir}"
        } else {
            run_log.error "Pipeline failed${failed ? ' in ' + failed.join(', ') : ''}. Check .nextflow.log for details."
        }
        status_file.text = 'process\tsucceeded\tcached\tfailed\tignored\tunfinished\n' +
            run_info.stats.processes.collect { r ->
                "${r.name}\t${r.succeeded}\t${r.cached + r.stored}\t${r.failed - r.ignored - r.retries}\t${r.ignored}\t" +
                    "${r.pending + r.submitted + r.running + r.aborted}\n"
            }.join('')
    }
}