<p>I haven't seen that particular error before, but I'd first make sure your pipeline errors aren't being masked by your shell. Consider adding the following to your <code>nextflow.config</code>:</p>
<pre><code>process {

  shell = [ '/bin/bash', '-euo', 'pipefail' ]
}
</code></pre>
<p>IIRC, the default shell (as provided by Nextflow) does not include the pipefail option for <a href="https://www.gnu.org/savannah-checkouts/gnu/bash/manual/bash.html#Pipelines" rel="nofollow noreferrer">pipelines</a>.</p>
<hr />
<p>The only other thing I can think of is to make sure your reference FASTA (and BWA index files) are localized in the workDir. You might find the intermittent (filesystem?) errors maybe go away even if you are staging using symlinks.</p>
<pre><code>params.ref_fasta = './bwa_index/ref.fasta'

align_to_ref = Channel
    .fromPath( [ params.ref_fasta, &quot;${params.ref_fasta}.*&quot; ] )
    .collect()


process bwa_mem_proc {

    tag { &quot;<span class="math-container">${sample_name}:$</span>{fastq.name}&quot; }

    label 'bwa_mem_proc'

    input:
    tuple val(sample_name), path(fastq), path(header)
    path align_to_ref

    output:
    tuple val(sample_name), path(&quot;${fastq.baseName}.unsorted.mini.bam&quot;)

    script:
    def task_cpus = task.cpus &gt; 1 ? task.cpus - 1 : task.cpus

    if( params.output_format == &quot;bam&quot; ) {
        &quot;&quot;&quot;
        bwa mem \\
            -p \\
            -t <span class="math-container">${task_cpus} \\
            -M \\
            -C \\
            -H &lt;(grep "^@RG" "$</span>{header}&quot;) \\
            &quot;<span class="math-container">${align_to_ref.first()}" \\
            "$</span>{fastq}&quot; |
        samtools view \\
            -1 \\
            -o &quot;<span class="math-container">${fastq.baseName}.unsorted.mini.bam" \\
            -
        """
    } else {
        """
        bwa mem \\
            -p \\
            -t ${task_cpus} \\
            -M \\
            -C \\
            -H &lt;(grep "^@RG" "${header}") \\
            "${align_to_ref.first()}" \\
            "${fastq}" |
        samtools view \\
            -C \\
            -T "${align_to_ref.first()}" \\
            -o "$</span>{fastq.baseName}.unsorted.mini.bam&quot; \\
            -
        &quot;&quot;&quot;
    }
}
</code></pre>
<p>Note that in the first case (where <code>params.output_format == &quot;bam&quot;</code>), the reference FASTA file itself is not required by BWA MEM. Only the BWA index files are needed here.</p>

CC BY-SA 4.0
2022-02-23T16:26:39.037Z
Indeed! You were right again. Turns out `bwa` was running out of memory but that error wasn't floating to the top.
Mark Ebbert
CC BY-SA 4.0
2022-02-23T01:39:25.750Z
2022-02-23T01:39:25.750Z
Steve
<h1>Update</h1>
<p>The issue was that <code>bwa</code> was running out of memory and failing, but that error wasn't floating to the top (see @Steve's answer, below). I was getting an error from <code>samtools</code>. I should have been able to figure that out, but hopefully this post saves someone else the trouble.</p>
<h1>Original post</h1>
<p>I'm experiencing a strange intermittent parsing error when piping <code>bwa mem</code> to <code>samtools view</code>. It only seems to happen when the command runs within a <code>Nextflow</code> pipeline (though it seems unlikely that <code>Nextflow</code> has anything to do with it). If I <code>cd</code> to the directory and run <code>bash .command.sh</code>, it runs to completion without issue.</p>
<h1>Nextflow process</h1>
<p>Here's the <code>Nextflow</code> process I'm running (code credit to @Steve, though modified to handle <code>.bam</code> or <code>.cram</code>). The portion that's running in this case</p>
<pre><code>process bwa_mem_proc {

    tag { &quot;<span class="math-container">${sample_name}:$</span>{fastq.name}&quot; }

    label 'bwa_mem_proc'

    input:
    tuple val(sample_name), path(fastq), path(header)

    output:
    tuple val(sample_name), path(&quot;${fastq.baseName}.unsorted.mini.bam&quot;)

    script:
    def task_cpus = task.cpus &gt; 1 ? task.cpus - 1 : task.cpus

    &quot;&quot;&quot;
    if [ &quot;<span class="math-container">$params.output_format" = "bam" ]; then
    bwa mem \\
        -p \\
        -t $</span>{task_cpus} \\
            -M \\
            -C \\
            -H &lt;(grep &quot;^@RG&quot; &quot;<span class="math-container">${header}") \\
        "$</span>{params.align_to_ref}&quot; \\
            &quot;<span class="math-container">${fastq}" |
    samtools view \\
        -1 \\
        -o "$</span>{fastq.baseName}.unsorted.mini.bam&quot; \\
            -
    else
        bwa mem \\
            -p \\
            -t <span class="math-container">${task_cpus} \\
        -M \\
        -C \\
        -H &lt;(grep "^@RG" "$</span>{header}&quot;) \\
            &quot;<span class="math-container">${params.align_to_ref}" \\
        "$</span>{fastq}&quot; |
        samtools view \\
            -C \\
            -T params.align_to_ref \\
            -o &quot;${fastq.baseName}.unsorted.mini.bam&quot; \\
            -
    fi
    &quot;&quot;&quot;
}
</code></pre>
<h1>Exact error</h1>
<p>Here's an example of the intermittent errors I'm getting. I would have assumed that some read (e.g., read # <code>2_102_429</code> in the <code>.fastq</code>) is malformed, but as I mentioned, running <code>bash .command.sh</code> in the working directory for the failed run finishes without issue.</p>
<pre><code>Command error:
  [M::mem_pestat] skip orientation FF
  [M::mem_process_seqs] Processed 700000 reads in 109.002 CPU sec, 16.411 real sec
  [M::process] 0 single-end sequences; 700000 paired-end sequences
  [M::process] read 700000 sequences (70000000 bp)...
  [M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (24, 301945, 4, 10)
  [M::mem_pestat] analyzing insert size distribution for orientation FF...
  [M::mem_pestat] (25, 50, 75) percentile: (118, 180, 399)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 961)
  [M::mem_pestat] mean and std.dev: (179.38, 124.36)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 1242)
  [M::mem_pestat] analyzing insert size distribution for orientation FR...
  [M::mem_pestat] (25, 50, 75) percentile: (247, 319, 420)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 766)
  [M::mem_pestat] mean and std.dev: (342.90, 127.88)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 939)
  [M::mem_pestat] skip orientation RF as there are not enough pairs
  [M::mem_pestat] analyzing insert size distribution for orientation RR...
  [M::mem_pestat] (25, 50, 75) percentile: (326, 764, 1452)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 3704)
  [M::mem_pestat] mean and std.dev: (893.30, 530.98)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 4830)
  [M::mem_pestat] skip orientation FF
  [M::mem_pestat] skip orientation RR
  [M::mem_process_seqs] Processed 700000 reads in 97.823 CPU sec, 14.292 real sec
  [M::process] 0 single-end sequences; 700000 paired-end sequences
  [M::process] read 700000 sequences (70000000 bp)...
  [M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (14, 196658, 9, 15)
  [M::mem_pestat] analyzing insert size distribution for orientation FF...
  [M::mem_pestat] (25, 50, 75) percentile: (206, 422, 1021)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 2651)
  [M::mem_pestat] mean and std.dev: (506.92, 458.91)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 3466)
  [M::mem_pestat] analyzing insert size distribution for orientation FR...
  [M::mem_pestat] (25, 50, 75) percentile: (248, 322, 429)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 791)
  [M::mem_pestat] mean and std.dev: (345.17, 129.22)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 972)
  [M::mem_pestat] skip orientation RF as there are not enough pairs
  [M::mem_pestat] analyzing insert size distribution for orientation RR...
  [M::mem_pestat] (25, 50, 75) percentile: (485, 777, 1315)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 2975)
  [M::mem_pestat] mean and std.dev: (961.60, 569.27)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 3805)
  [M::mem_pestat] skip orientation FF
  [M::mem_pestat] skip orientation RR
  [M::mem_process_seqs] Processed 700000 reads in 275.897 CPU sec, 39.910 real sec
  [M::process] 0 single-end sequences; 700000 paired-end sequences
  [M::process] read 700000 sequences (70000000 bp)...
  [W::sam_read1_sam] Parse error at line 2102429
  samtools view: error reading file &quot;-&quot;
</code></pre>
<p>A previous sample complained about an <code>@SQ</code> header line missing the <code>LN:</code> tag. It's like the pipe gets interrupted and lines get truncated, or something.</p>
<pre><code>  [E::sam_hrecs_update_hashes] Header includes @SQ line &quot;KI270330&quot; with no LN: tag
  samtools view: failed to add PG line to the header
</code></pre>
<p>Anyone seen this before?</p>

<phylogenetics><samtools><bwa><nextflow>
BWA mem | samtools view: Intermittent parsing error
# BWA mem | samtools view: Intermittent parsing error

(Asked by: username_0 on 2022-02-22)

Update
The issue was that bwa was running out of memory and failing, but that error wasn't floating to the top (see @Steve's answer, below). I was getting an error from samtools. I should have been able to figure that out, but hopefully this post saves someone else the trouble.
Original post
I'm experiencing a strange intermittent parsing error when piping bwa mem to samtools view. It only seems to happen when the command runs within a Nextflow pipeline (though it seems unlikely that Nextflow has anything to do with it). If I cd to the directory and run bash .command.sh, it runs to completion without issue.
Nextflow process
Here's the Nextflow process I'm running (code credit to @Steve, though modified to handle .bam or .cram). The portion that's running in this case
process bwa_mem_proc {

    tag { "${sample_name}:${fastq.name}" }

    label 'bwa_mem_proc'

    input:
    tuple val(sample_name), path(fastq), path(header)

    output:
    tuple val(sample_name), path("${fastq.baseName}.unsorted.mini.bam")

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

    """
    if [ "$params.output_format" = "bam" ]; then
    bwa mem \\
        -p \\
        -t ${task_cpus} \\
            -M \\
            -C \\
            -H <(grep "^@RG" "${header}") \\
        "${params.align_to_ref}" \\
            "${fastq}" |
    samtools view \\
        -1 \\
        -o "${fastq.baseName}.unsorted.mini.bam" \\
            -
    else
        bwa mem \\
            -p \\
            -t ${task_cpus} \\
        -M \\
        -C \\
        -H <(grep "^@RG" "${header}") \\
            "${params.align_to_ref}" \\
        "${fastq}" |
        samtools view \\
            -C \\
            -T params.align_to_ref \\
            -o "${fastq.baseName}.unsorted.mini.bam" \\
            -
    fi
    """
}

Exact error
Here's an example of the intermittent errors I'm getting. I would have assumed that some read (e.g., read # 2_102_429 in the .fastq) is malformed, but as I mentioned, running bash .command.sh in the working directory for the failed run finishes without issue.
Command error:
  [M::mem_pestat] skip orientation FF
  [M::mem_process_seqs] Processed 700000 reads in 109.002 CPU sec, 16.411 real sec
  [M::process] 0 single-end sequences; 700000 paired-end sequences
  [M::process] read 700000 sequences (70000000 bp)...
  [M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (24, 301945, 4, 10)
  [M::mem_pestat] analyzing insert size distribution for orientation FF...
  [M::mem_pestat] (25, 50, 75) percentile: (118, 180, 399)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 961)
  [M::mem_pestat] mean and std.dev: (179.38, 124.36)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 1242)
  [M::mem_pestat] analyzing insert size distribution for orientation FR...
  [M::mem_pestat] (25, 50, 75) percentile: (247, 319, 420)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 766)
  [M::mem_pestat] mean and std.dev: (342.90, 127.88)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 939)
  [M::mem_pestat] skip orientation RF as there are not enough pairs
  [M::mem_pestat] analyzing insert size distribution for orientation RR...
  [M::mem_pestat] (25, 50, 75) percentile: (326, 764, 1452)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 3704)
  [M::mem_pestat] mean and std.dev: (893.30, 530.98)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 4830)
  [M::mem_pestat] skip orientation FF
  [M::mem_pestat] skip orientation RR
  [M::mem_process_seqs] Processed 700000 reads in 97.823 CPU sec, 14.292 real sec
  [M::process] 0 single-end sequences; 700000 paired-end sequences
  [M::process] read 700000 sequences (70000000 bp)...
  [M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (14, 196658, 9, 15)
  [M::mem_pestat] analyzing insert size distribution for orientation FF...
  [M::mem_pestat] (25, 50, 75) percentile: (206, 422, 1021)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 2651)
  [M::mem_pestat] mean and std.dev: (506.92, 458.91)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 3466)
  [M::mem_pestat] analyzing insert size distribution for orientation FR...
  [M::mem_pestat] (25, 50, 75) percentile: (248, 322, 429)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 791)
  [M::mem_pestat] mean and std.dev: (345.17, 129.22)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 972)
  [M::mem_pestat] skip orientation RF as there are not enough pairs
  [M::mem_pestat] analyzing insert size distribution for orientation RR...
  [M::mem_pestat] (25, 50, 75) percentile: (485, 777, 1315)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 2975)
  [M::mem_pestat] mean and std.dev: (961.60, 569.27)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 3805)
  [M::mem_pestat] skip orientation FF
  [M::mem_pestat] skip orientation RR
  [M::mem_process_seqs] Processed 700000 reads in 275.897 CPU sec, 39.910 real sec
  [M::process] 0 single-end sequences; 700000 paired-end sequences
  [M::process] read 700000 sequences (70000000 bp)...
  [W::sam_read1_sam] Parse error at line 2102429
  samtools view: error reading file "-"

A previous sample complained about an @SQ header line missing the LN: tag. It's like the pipe gets interrupted and lines get truncated, or something.
  [E::sam_hrecs_update_hashes] Header includes @SQ line "KI270330" with no LN: tag
  samtools view: failed to add PG line to the header

Anyone seen this before?

**Tags:** phylogenetics, samtools, bwa, nextflow

---

## 1 answers

### Answer by username_1 on 2022-02-23. Score: 2 | Accepted answer

I haven't seen that particular error before, but I'd first make sure your pipeline errors aren't being masked by your shell. Consider adding the following to your nextflow.config:
process {

  shell = [ '/bin/bash', '-euo', 'pipefail' ]
}

IIRC, the default shell (as provided by Nextflow) does not include the pipefail option for pipelines.

The only other thing I can think of is to make sure your reference FASTA (and BWA index files) are localized in the workDir. You might find the intermittent (filesystem?) errors maybe go away even if you are staging using symlinks.
params.ref_fasta = './bwa_index/ref.fasta'

align_to_ref = Channel
    .fromPath( [ params.ref_fasta, "${params.ref_fasta}.*" ] )
    .collect()


process bwa_mem_proc {

    tag { "${sample_name}:${fastq.name}" }

    label 'bwa_mem_proc'

    input:
    tuple val(sample_name), path(fastq), path(header)
    path align_to_ref

    output:
    tuple val(sample_name), path("${fastq.baseName}.unsorted.mini.bam")

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

    if( params.output_format == "bam" ) {
        """
        bwa mem \\
            -p \\
            -t ${task_cpus} \\
            -M \\
            -C \\
            -H <(grep "^@RG" "${header}") \\
            "${align_to_ref.first()}" \\
            "${fastq}" |
        samtools view \\
            -1 \\
            -o "${fastq.baseName}.unsorted.mini.bam" \\
            -
        """
    } else {
        """
        bwa mem \\
            -p \\
            -t ${task_cpus} \\
            -M \\
            -C \\
            -H <(grep "^@RG" "${header}") \\
            "${align_to_ref.first()}" \\
            "${fastq}" |
        samtools view \\
            -C \\
            -T "${align_to_ref.first()}" \\
            -o "${fastq.baseName}.unsorted.mini.bam" \\
            -
        """
    }
}

Note that in the first case (where params.output_format == "bam"), the reference FASTA file itself is not required by BWA MEM. Only the BWA index files are needed here.

> **Comments:**
> **username_0**: Indeed! You were right again. Turns out `bwa` was running out of memory but that error wasn't floating to the top.  

---

