You might have added something like the following to your nextflow.config:
process {
shell = [ '/bin/bash', '-euo', 'pipefail' ]
}
You get an error because zcat has had its destination closed early, and wasn't able to (successfully) write a decompressed file. It has no way of knowing that this was intentional, or if something bad happened (e.g. out of disk space, network error, etc).
If you'd rather not make changes to your shell options, you could instead use a tool that reads the entire input stream. This is what your awk solution does, which is the same as:
zcat test.fastq.gz | awk 'NR<=900000 { print $0 }'
For each input record, AWK tests the current record number (NR) to see if it is less than 900000. If it is, the whole line is printed to stdout. If it's not, try to see if the next line is less 900000. The problem with this solution is that this will be slow if your input FASTQ is really big. So you might be tempted to use:
zcat test.fastq.gz | awk 'NR>900000 { exit }'
But this will work much like the original head solution, and will produce an error with your current shell options. A better solution (which lets you keep your current shell options) is to instead use process substitution to avoid the pipe:
head -n 900000 < <(zcat test.fastq.gz)
Your Nextflow code might look like:
nextflow.enable.dsl=2
params.max_lines = 900000
process TRUNCATE_FASTQ {
input:
tuple val(sample), path(reads)
script:
def (fq1, fq2) = reads
"""
head -n ${params.max_lines} < <(zcat "${fq1}") | gzip > "truncated_${fq1}"
head -n ${params.max_lines} < <(zcat "${fq2}") | gzip > "truncated_${fq2}"
"""
}
workflow {
TRUNCATE_FASTQ( Channel.fromFilePairs('*.{1,2}.fastq.gz') )
}
If you have Python, which has native gzip support, you could also use pysam for this. For example:
params.max_records = 225000
process TRUNCATE_FASTQ {
container 'quay.io/biocontainers/pysam:0.16.0.1--py37hc334e0b_1'
input:
tuple val(sample), path(reads)
script:
def (fq1, fq2) = reads
"""
#!/usr/bin/env python
import gzip
import pysam
def truncate(fn, num):
with pysam.FastxFile(fn) as ifh:
with gzip.open(f'truncated_{fn}', 'wt') as ofh:
for idx, entry in enumerate(ifh):
if idx >= num:
break
print(str(entry), file=ofh)
truncate('${fq1}', ${params.max_records})
truncate('${fq2}', ${params.max_records})
"""
}
Or using Biopython:
params.max_records = 225000
process TRUNCATE_FASTQ {
container 'quay.io/biocontainers/biopython:1.78'
input:
tuple val(sample), path(reads)
script:
def (fq1, fq2) = reads
"""
#!/usr/bin/env python
import gzip
from itertools import islice
from Bio import SeqIO
def xopen(fn, mode='r'):
if fn.endswith('.gz'):
return gzip.open(fn, mode)
else:
return open(fn, mode)
def truncate(fn, num):
with xopen(fn, 'rt') as ifh, xopen(f'truncated_{fn}', 'wt') as ofh:
seqs = islice(SeqIO.parse(ifh, 'fastq'), num)
SeqIO.write(seqs, ofh, 'fastq')
truncate('${fq1}', ${params.max_records})
truncate('${fq2}', ${params.max_records})
"""
}
