Samtools index fails within nextflow workflow but not outside of it

Question

I am a beginner to Nextflow and using it to setup Chipseq workflow and in the workflow below the samtools_index process throws the following error:

[E::hts_idx_push] Chromosome blocks not continuous
  [E::sam_index] Read 'FCC1L20ACXX:4:1101:1645:2230#0/1' with ref_name='chr17', ref_length=83257441, flags=16, pos=58851543 cannot be indexed
  samtools index: failed to create index for "TCDD_plus_E2_ERa_Rep_1_ChIPSeq_first_1000.Aligned.sorted.bam": Inappropriate ioctl for device

The command executed is:

/Users/siddhaduio.no/Desktop/All_omics_tools/samtools-1.14/samtools index TCDD_plus_E2_ERa_Rep_1_ChIPSeq_first_1000.Aligned.sorted.bam

However, when I index the BAM file outside of nextflow, it runs without any problems. I was wondering what might be the issue here?

#!/usr/bin/env nextflow

nextflow.enable.dsl = 2

params.genome = '/Users/siddhaduio.no/Desktop/All_omics_tools/bowtie2-2.4.4-macos-x86_64/hg38_index/human_hg38_index.*.bt2'
params.reads = '/Volumes/matthews/AHR_ER_project_files/FASTQ/TCDD_plus_E2_ERa_Rep_{1,2,4}_ChIPSeq_first_1000.fastq'
params.results = '/Volumes/matthews/AHR_ER_project_files/FASTQ/results'



/*
* Process to align chip seq reads to genome with bowtie2
*/


process mapping {

publishDir '/Volumes/matthews/AHR_ER_project_files/FASTQ/results', mode: 'copy'

input:

path genomeDir

tuple val(replicateId), path(reads)



output:

tuple val(replicateId), path("${replicateId}.Aligned.sam")


script:

def btw2_index = genomeDir[0].simpleName

"""
/Users/siddhaduio.no/Desktop/All_omics_tools/bowtie2-2.4.4-macos-x86_64/bowtie2 -p 64 -q -x ${btw2_index} -U ${reads.join(',')} -S ${replicateId}.Aligned.sam 2> ${replicateId}.bowtie.log


"""


}

process convert_sam_to_bam {



publishDir '/Volumes/matthews/AHR_ER_project_files/FASTQ/results', mode: 'copy'

input:


tuple val(replicateId), path("${replicateId}.Aligned.sam")



output:

tuple val(replicateId), path("${replicateId}.Aligned.bam"), path("${replicateId}.Aligned.sorted.bam")


script:

"""

/Users/siddhaduio.no/Desktop/All_omics_tools/samtools-1.14/samtools view -S -b -o ${replicateId}.Aligned.bam ${replicateId}.Aligned.sam 
/Users/siddhaduio.no/Desktop/All_omics_tools/samtools-1.14/samtools sort --output-fmt BAM -o ${replicateId}.Aligned.sorted.bam ${replicateId}.Aligned.bam
rm ${replicateId}.Aligned.sam


"""


}


process samtools_index {



publishDir '/Volumes/matthews/AHR_ER_project_files/FASTQ/results', mode: 'copy'

input:


tuple val(replicateId), path("${replicateId}.Aligned.sorted.bam")



output:

tuple val(replicateId), path("${replicateId}.Aligned.sorted.bam.bai")


script:

"""

/Users/siddhaduio.no/Desktop/All_omics_tools/samtools-1.14/samtools index ${replicateId}.Aligned.sorted.bam


"""


}


workflow{

reads_ch= Channel.fromPath(params.reads).map { fastq -> tuple(fastq.simpleName, fastq)}
index_ch= files(params.genome)
reads_ch
index_ch
mapping(index_ch,reads_ch)
convert_sam_to_bam(mapping.out)
samtools_index(convert_sam_to_bam.out)


}

```

Answer 1

You should have also received the following Nextflow warning:

WARN: Input tuple does not match tuple declaration in process `samtools_index`

Specifically, the problem is that the output of your convert_sam_to_bam process is:

tuple val(replicateId), path("${replicateId}.Aligned.bam"), path("${replicateId}.Aligned.sorted.bam")

But the samtools_index process is expecting:

tuple val(replicateId), path("${replicateId}.Aligned.sorted.bam")

---

If you would like to declare multiple output channels, you can use the emit keyword (this is documented under additional options) to define names for the output channel(s). For example:

process convert_sam_to_bam {

    input:
    tuple val(replicateId), path("${replicateId}.Aligned.sam")

    output:
    tuple val(replicateId), path("${replicateId}.Aligned.bam"), emit: aligned_bam
    tuple val(replicateId), path("${replicateId}.Aligned.sorted.bam"), emit: sorted_bam

    script:
    """
    samtools view \\
        -b \\
        -o "${replicateId}.Aligned.bam" \\
        "${replicateId}.Aligned.sam"

    samtools sort \\
        -o "${replicateId}.Aligned.sorted.bam" \\
        "${replicateId}.Aligned.bam"
    """
}

workflow {

    ...

    samtools_index( convert_sam_to_bam.out.sorted_bam )
}

---

Otherwise, either change the convert_sam_to_bam output declaration to:

tuple val(replicateId), path("${replicateId}.Aligned.sorted.bam")

Or use the map operator to get what you want, for example:

workflow {

    ...

    samtools_index(
        convert_sam_to_bam.out.map { rep, aligned, sorted -> tuple( rep, sorted ) }
    )
}

---

Note that a separate indexing process can be avoided using newer versions of samtools with the --write-index option:

To convert a SAM to a compressed BAM using BAI indexing: samtools view --write-index in.sam -o out.bam##idx##out.bam.bai