MACS2 “0 reads have been read” error after alignmentSieve in ATAC-seq Nextflow pipeline (single-end)
Hi all,
I am developing an ATAC-seq analysis pipeline for single-end reads using Nextflow and am running into an error at the peak calling step with MACS2.
The error is:
ERROR ~ Error executing process > 'CALL_PEAKS (2)'
Caused by:
Process `CALL_PEAKS (2)` terminated with an error exit status (1)
...
Command output:
0
Command error:
INFO ...: #1 total tags in treatment: 0
...
ZeroDivisionError: float division by zero
My pipeline is:
- Download with SRA tools. Do QC with FastQC, trimmomatic, MultiQC, and samtools flagstat.
- Alignment using Bowtie2 for single end reads
- Remove mitochondrial reads with grep
- Blacklist removal with bedtools intersect -v
- Tn5 shift correction using alignmentSieve
- MACS2 peak calling, using callpeak with -f BAM --nomodel --extsize 147 for single-end ATAC-seq
The reference genome is mm10/GRCm38.
The relevant part of the workflow is:
REMOVE_BLACKLIST(BOWTIE2_ALIGN.out.filtered_bam, params.blacklist)
SHIFT_ALIGN(REMOVE_BLACKLIST.out)
CALL_PEAKS(SHIFT_ALIGN.out)
The relevant processes are:
process REMOVE_BLACKLIST {
conda 'envs/bedtools_env.yml'
input:
tuple val(run), val(biosample), val(samplename), val(library_layout), val(library_source), val(experiment), path(filtered_bam), path(filtered_index)
path blacklist
output:
tuple val(run), val(biosample), val(samplename), val(library_layout), val(library_source), val(experiment), path("${run}.blacklist_removed.sorted.bam"), path("${run}.blacklist_removed.sorted.bam.bai")
script:
"""
bedtools intersect -v -abam ${filtered_bam} -b ${blacklist} > ${run}.blacklist_removed.bam
# Sort & index the new BAM
samtools sort -o ${run}.blacklist_removed.sorted.bam ${run}.blacklist_removed.bam
samtools index ${run}.blacklist_removed.sorted.bam
rm ${run}.blacklist_removed.bam
samtools view -c ${run}.blacklist_removed.sorted.bam
"""
}
process SHIFT_ALIGN {
conda 'envs/deeptools_env.yml'
input:
tuple val(run), val(biosample), val(samplename), val(library_layout), val(library_source), val(experiment), path(filtered_bam), path(filtered_index)
output:
tuple val(run), val(biosample), val(samplename), val(library_layout), val(library_source), val(experiment), path("${run}_shifted.bam"), path("${run}_shifted.bam.bai")
script:
"""
alignmentSieve -b ${filtered_bam} -o ${run}_shifted.bam --ATACshift --minMappingQuality 10 --filterMetrics log.txt
samtools index ${run}_shifted.bam
"""
}
process CALL_PEAKS {
conda 'envs/macs2_env.yml'
publishDir params.outdir, mode: "copy", pattern: ['*.narrowPeak','*.bed','*.xls','*.bdg']
input:
tuple val(run), val(biosample), val(samplename), val(library_layout), val(library_source), val(experiment), path(shifted_blacklist_removed_bam), path(shifted_blacklist_removed_bam_index)
output:
tuple val(run), val(biosample), val(samplename), val(library_layout), val(library_source), val(experiment),
path("${run}_peaks.narrowPeak"),
path("${run}_summits.bed"),
path("${run}_peaks.xls"),
path("${run}_treat_pileup.bdg")
script:
"""
samtools view -c ${shifted_blacklist_removed_bam}
macs2 callpeak \
-t ${shifted_blacklist_removed_bam} \
-f BAM \
-g mm \
--nomodel \
--keep-dup auto \
--extsize 147 \
-n ${run}
"""
}
My package versions are:
deeptools=3.5.6
samtools=1.20
macs2=2.2.9.1
bedtools=2.31.1
I have narrowed down the problem to alignmentSieve when trying to apply the ATAC-seq shift. Output BAM files from the remove blacklist process have millions of reads, but the input BAM files to the call peaks process appear to have 0 reads. The result is when I try to call MACS2, MACS2 sees the input BAM files have 0 reads.
I am not sure why alignmentSieve outputs BAM files with 0 reads. I checked the documentation, and inputs/outputs between processes seem to line up.
• 3,220 views
•
link
written
by
DdogBoss •
1 answer
If the --shift or --ATACshift options are used, then only properly-paired reads will be used.
--ATACshift does not work on single end reads. The solution may be to forego shifting the alignments.
• 0 views
•
link
written
by
DdogBoss •
Log in to answer this question.
More posts like this
-
Nextflow process fails with “Invalid method invocation call” despite correct tuple
written by DdogBoss •
-
Nextflow: How to format input tuple for STAR_ALIGN process with STAR index directory
written by DdogBoss •
-
Nextflow process not executing
written by myoui3122010 •
-
Nextflow bwa-mem2 pipeline/channel question
written by megaphonecat •
-
Nextflow, transform multiple outputs of one process to paired outputs and use them as input to the…
written by MolGeek •
-
Nextflow output from different processes mixed up as input in later processes
written by Eliveri •
-
BWA alignment failing after adding nextflow.config file?
written by Eliveri •
-
how to pass trimmed fastP output to bwa in Nextflow
written by rohitsatyam102
-
Help in writing NGS pipeline (or: the road from fasta to vcf)
written by Hervé •
-
How to retrieve BAM file from NCBI SRA
written by mufernando •
Can you share the
log.txtfile from the alignmentSieve output?Hi, Arup.
Thank you for commenting. Here are the log.txt files for the samples that made it through alignmentSieve:
So, it does look like alignmentSieve still does produce BAM files with reads.
Manually inspect the bam files whether they have content. If so, run the peak calling manually outside of nextflow to debug. Do that for the troublesome sample and post the full logs. Debugging should be always from most lowlevel (manual execution) to most highlevel (inside wrappers/cotnainers/workflow managers).
Hi, ATpoint.
I checked the header. This is what the header in the work directory of the call peak process looks like for one BAM file:
I also tried running call peaks outside of the nextflow pipeline, and I still run into the same error:
I have also tried removing the blacklist after shifting alignments and calling peaks, but I still run into the same Zero Division error with calling peaks.
In terms of removing the blacklist after shifting alignments, this was the workflow:
The bam file is empty. Now you have to dig back why. Is it the alignment that failed, or any filtering or the sieving.
I figured out the problem and it was --ATACshift and alignmentSieve. According to the documentation,
--ATACshift does not work on single end reads. alignmentSieve documents read statistics before filtering, which is confusing.
The
alignmentSievestep is definitely generating a BAM file. I would suggest looking into the corresponding work directory of the run.