3
$\begingroup$

I am trying to run my nextflow pipeline, and have gotten this error:
samtools sort: failed to read header from "-"

I'm not sure why this error is occurring.

#!/usr/bin/env nextflow
nextflow.enable.dsl=2
import groovy.io.*

params.t = 8

process joinFastq {     
input:
path(ds_directory)

output:
tuple path("joined_${ds_directory.baseName}.fastq"), val("$ds_directory.baseName")

script:
"""
gunzip -c ${ds_directory}/*/fastq_pass/*.fastq.gz > 
 joined_${ds_directory.baseName}.fastq
"""

}
 process minimap2 {
cpus 8
memory "20 GB"
input:
tuple path(fastq), val(baseName)

output:
tuple path("*.bam"), val(baseName)

  script:
   """
  minimap2 -ax splice -uf -k14 -t 8 ${params.genome} ${fastq} > mapped_${baseName}.sam
    samtools view -h -b mapped_${baseName}.sam > ${baseName}.bam
   samtools sort -o ${baseName}.bam > ${baseName}_sorted.bam
   """
   }

I have tried applying the -h flag to samtools sort in the script

   samtools sort -h -o ${baseName}.bam > ${baseName}_sorted.bam

and the same error is produced. I'm not sure why the header flag didn't solve this issue.

$\endgroup$

1 Answer 1

3
$\begingroup$

This is not a Nextflow error. The problem is that without an input file, samtools sort tries to read from stdin. Using a recent samtools, you can however coordinate sort the SAM and write a sorted BAM using:

samtools sort -o "${baseName}.bam" "mapped_${baseName}.sam"

You may have been intending to pipe the output to samtools sort, which would avoid writing large SAM files and is usually preferable. A minimal example might look like:

params.reads = './*.fastq.gz'
params.reference = './gencode.v43.transcripts.fa'
process minimap2 {

    tag { sample }

    cpus 8
    memory 20.GB

    input:
    tuple val(sample), path(fastq)
    path ref_fasta

    output:
    tuple val(sample), path("${sample}.bam")

    script:
    def cpus = task.cpus > 1 ? task.cpus - 1 : task.cpus

    """
    minimap2 \\
        -ax splice \\
        -uf \\
        -k14 \\
        -t ${cpus} \\
        "${ref_fasta}" \\
        "${fastq}" |
    samtools sort \\
        -o "${sample}.bam" 
    """
}
workflow {

    reads = Channel.fromFilePairs( params.reads, size: 1 )
    reference = file( params.reference )

    minimap2( reads, reference )
}

Also, be sure to avoid using absolute paths in your workflow scripts. For your workflow to be portable, processes should only access files in the working directory. The example above avoids this, by using a value channel, to stage the reference FASTA file into the working directory.

In case anyone is wondering, the ternary is to ensure that we do not oversubscribe the node(s) that these jobs land on. Otherwise, each job could use (up to) 900% (800% + 100% = 900%) CPU.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.