Enable the use of CSI bam index

[JoseEspinosa]

Description of feature

See nf-core/atacseq#327 and nf-core/atacseq#328

[a1ultima]

@JoseEspinosa

We have successfully patched in a new parameter:

--bam_csi_index = true

Solved by

making three changes:

1. conf/modules.config: add the following to the bottom of modules.config

process {
    // wildcard ext.args: matches all SAMTOOLS_INDEX subtasks w/i the scope of other tasks in the modules.config
    withName: /.*SAMTOOLS_INDEX$/ {
        ext.args = params.bam_csi_index ? '-c' : ''
    }
}

2. modules/nf-core/samtools/index/main.nf: updated SAMTOOLS_INDEX main.nf, to handle CSI vs. standard (.bai) index options

/*
 *  honours a new Boolean pipeline parameter `params.bam_csi_index`
 *    – when true, passes `--csi` to samtools index
 *    – when false / unset, behaves exactly as before (creates .bai)
*/
process SAMTOOLS_INDEX {
    tag "$meta.id"
    label 'process_low'

    conda "${moduleDir}/environment.yml"
    container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
        'https://depot.galaxyproject.org/singularity/samtools:1.20--h50ea8bc_0' :
        'biocontainers/samtools:1.20--h50ea8bc_0' }"

    input:
        tuple val(meta), path(bam)

    output:
        tuple val(meta), path("*.bai"), optional:true, emit:bai
        tuple val(meta), path("*.csi"), optional:true, emit:csi
        tuple val(meta), path("*.crai"), optional:true, emit:crai
        path  "versions.yml",                  emit:versions

    when:
        task.ext.when == null || task.ext.when

    script:
        /*
         * build the flag list:
         *   • task.ext.args still works for power-users
         *   • --csi added iff params.bam_csi_index is true
         */
        def extraArgs = []
        if (params.bam_csi_index)  extraArgs << '--csi'
        if (task.ext.args)         extraArgs << task.ext.args
        def joinedArgs = extraArgs.join(' ')

        """
        samtools index -@ ${task.cpus - 1} ${joinedArgs} ${bam}

        cat <<-END_VERSIONS > versions.yml
        "${task.process}":
            samtools: \$(samtools --version | head -n1 | awk '{print \$2}')
        END_VERSIONS
        """
    stub:
        /*
         * replicate real-run behaviour so `.nextflow -stub-run` works
         */
        def makeIndex = { Path inFile ->
            String ext
            if (inFile.extension == 'cram') {
                ext = 'crai'
            } else if (params.bam_csi_index || (task.ext.args ?: '').contains('-c')) {
                ext = 'csi'
            } else {
                ext = 'bai'
            }
            file("${inFile}.${ext}").text = ''
        }

        makeIndex(bam)

        """
        cat <<-END_VERSIONS > versions.yml
        "${task.process}":
            samtools: stub
        END_VERSIONS
        """
}

3. nextflow_schema.json:

...
                "bam_csi_index": {
                "type"       : "boolean",
                "description": "Create a CSI index instead of BAI for BAM files (required for genomes with >2 Gb contigs).",
                "fa_icon"    : "fas fa-ruler-combined",
                "default"    : false,
                "help_text"  : "When true the SAMTOOLS_INDEX process will call `samtools index -c` to produce `.csi` files and downstream tasks will expect that extension."
                }
...

Usage

Similar to nf-core/rnaseq, we can now add the following boolean option in nextflow.config:

params{
...
    bam_csi_index              = true  // true: for large contig-genomes makes .csi index // or false: to use default .bai index
...
}

Or as a CLI flag:

nextflow run main.nf \
...
--bam_csi_index true \

[JoseEspinosa]

Hi @a1ultima, thanks a lot for the detailed proposal!
Would you be interested in opening a PR to incorporate these changes into the pipeline? If not, no worries, just let me know and I’ll take care of it. 😄