quaqc: QUick ATAC-seq Quality Control

quaqc provides ATAC-seq quality control, read filtering, and genomic track generation, with separate metrics for nuclear, mitochondrial, and plastid reads. It also supports other unspliced DNA-seq experiments, such as ChIP-seq.

Version 1.7 adds optional peak calling and uses +4/-4 as the default Tn5 offset. See NEWS for the full list of changes and Tn5 offsets for compatibility with earlier versions.

Installation

Requires GCC or Clang and GNU Make. Tested on macOS and Linux:

git clone https://github.com/bjmt/quaqc  # Or download latest release
cd quaqc
make release-full
make test  # Make sure quaqc produces the expected outputs (optional)
make install  # Copy quaqc + manual to /usr/local (optional, may require sudo)

See INSTALL for configuration options.

Alternatively, install the version available from Bioconda:

# Example environment setup:
conda create -y --name quaqc -c conda-forge -c bioconda quaqc

Quick start

To see a basic description of all available parameters, run ./quaqc -h. For more detailed help, see the manual page.

Run QC on a coordinate-sorted BAM using the default filters:

quaqc Sample.bam
# -> Sample.quaqc.txt

For an example with peaks, TSSs, and a target region, use the included test data:

./quaqc -v --output-dir . --peaks test/peak.bed --tss test/tss.bed \
    --target-list test/target.bed test/reads.bam
cat ./reads.quaqc.txt

This writes the report to the current directory and uses the peak and TSS files to calculate additional metrics. The example BAM contains reads from a small region of the Arabidopsis genome, so target.bed restricts the analysis and adjusts the effective genome size accordingly. See the example report and metric definitions.

Note that quaqc can also save the final reads passing all filters to a new BAM file via the -S/--keep flag. (This will increase the runtime.)

quaqcr

The companion R package quaqcr can run quaqc and analyse its JSON outputs. See its README for a tutorial. To save JSON alongside the text report:

quaqc --json Sample.quaqc.json Sample.bam

Help

Please open an issue on GitHub.

Citation

Please cite the following journal article if you find this software useful in your research:

Tremblay, B.J.M. and Qüesta, J.I. (2024). quaqc: efficient and quick ATAC-seq quality control and filtering. Bioinformatics 40, btae649. https://doi.org/10.1093/bioinformatics/btae649

Testing

There are two test suites included in test/.

The first (make test) compares the QC report and peak coordinates from test/reads.bam against committed references. It uses standard shell tools and requires no additional bioinformatics software.

The second (bash test/test_suite.sh) is a more comprehensive suite which builds small BAMs from scratch and checks filtering, metrics, output formats, and bug regressions. It requires:

• samtools (for creating and indexing the BAMs)
• python3 (only the standard library, for parsing the JSON output)

Additional use cases

Calling peaks directly (no control needed)

quaqc can call peaks during a QC run using a MACS3-style algorithm without an input/control sample. Add --call-peaks to write a BGZF-compressed ENCODE narrowPeak file per sample:

quaqc -v --call-peaks Sample.bam
# -> Sample.narrowPeak.gz

Briefly, a Tn5-insertion pileup is built by extending each insertion to the size set by --peaks-extsize (150 bp by default), a local Poisson background is estimated from the genome-wide insertion rate and a --peaks-llocal window (10000 bp by default), and the resulting p-values are turned into genome-wide Benjamini-Hochberg q-values using the cutoff set by --peaks-qval (0.05 by default, as in macs3 callpeak). Adding --peaks-qscore will additionally write a -log10(q) bedGraph track. When --peaks is not supplied, the called peaks are also used for FRIP in the report. Peak calling requires a second pass through each BAM.

By default, the effective genome size used for both the background rate and the q-value correction is the effective nuclear genome size reported by quaqc (the summed nuclear chromosome sizes, accounting for --target-list/--blacklist, with organelles excluded). This can be overridden with --peaks-gsize if you would rather supply your own value (e.g. a mappable genome size to match macs3 callpeak -g). In a controlled comparison on a real Arabidopsis ATAC-seq sample (feeding the same quaqc-filtered insertion set to both tools, and giving them the same effective genome size) the calls agree quite closely with macs3 callpeak at the shared default q-value cutoff of 0.05, with a base-pair Jaccard of roughly 0.86-0.88 and the large majority of peaks reciprocally overlapping. This comparison is specific to that sample and matched settings. quaqc uses the genome-wide rate and the large --peaks-llocal window for its background; it does not use a smaller slocal window.

Use --peaks-split FRAC to optionally split peaks at valleys between local pileup maxima. See the manual for all peak-calling options. For a count matrix across samples, provide a shared BED file with --peaks; sample-specific calls cannot be used directly with --quant.

Creating insertion BED files for use by MACS2/3

To call peaks from individual alignment ends, export reads with --bed and pass the result to MACS2/3 as BED input:

quaqc -0 --bed Sample.bam

macs3 callpeak [...] \
    --treatment Sample.bed.gz --format BED \
    --call-summits --keep-dup all \
    --shift -75 --extsize 150 \
    --nomodel --nolambda

Add --bed-ins to export only the 5-prime insertion coordinates in BED3 format, reducing file size. Add --bed-tn5 if Tn5-corrected coordinates are required. Both exports apply quaqc's read filters.

Quantifying reads overlapping peaks from a BED file

Use --quant to count reads in a shared peak set across several BAM files:

quaqc -0 --peaks peaks.bed --quant counts.tsv --quant-ins --quant-tn5 *.bam

This writes a TSV count matrix. --quant-ins counts overlaps at the 5-prime insertion coordinate, and --quant-tn5 applies the Tn5 offset before checking overlaps. Read filters are applied during counting, so a separate filtered BAM is unnecessary.

Creating bedGraph files centered around the 5-prime ends of reads

To visualise insertion density, resize reads around their 5-prime ends and write a bedGraph in one step:

quaqc -0 --bedGraph --bedGraph-qlen 150 Sample.bam

This centres a 150-bp window on each read's 5-prime end. Set --bedGraph-qlen 0 to retain the original alignment coordinates. For single-base insertions with the Tn5 offset applied:

quaqc --bedGraph --bedGraph-tn5 --bedGraph-qlen 1 Sample.bam

See Tn5 offsets for the defaults and how to override them.

Filtering a single-cell ATAC-seq BAM by cell type

Using a set of barcodes, quaqc can subset a BAM by making use of the --rg-list and --keep options. For example, using the CB SAM tag for cell barcodes:

quaqc -0 -A --rg-list barcodes.txt --rg-tag CB --keep --keep-ext .MyCellType.bam Sample.bam

The -0 flag prevents quaqc from creating a QC report, thus making the output BAM the only file created by quaqc. The -A flag can be used if this BAM has already been filtered, and no additional filtering is desired.

Comparing read quality thresholds

Use -n/--target-names to restrict QC to a chromosome when exploring filtering thresholds. This reduces the amount of data scanned per run.

For example, to try out several MAPQ thresholds for an Arabidopsis ATAC-seq sample by only scanning chromosome 1 (which is just '1' for this version of TAIR10):

for i in 5 10 15 20 25 30 ; do
    quaqc -n 1 -q $i -i "MAPQ >= ${i}" -O .quaqc.mapq${i}.txt Sample.bam
done

Each threshold produces a separate report. The -i/--title flag labels each run; -O/--output-ext gives it a distinct output filename.

Quick motif footprints

The TSS pileup can also generate exploratory motif footprints. These are uncorrected for Tn5 sequence bias. Save the output as JSON to recover the pileup for plotting, for example with quaqcr.

In this example code, the --footprint preset will adjust the TSS pileup to produce single base resolution data of Tn5 transposase insertion frequency. The --nfr preset is also used to only consider nucleosome-free reads.

quaqc --nfr --footprint --tss motif.bed --json sample_motif.json Sample.bam

The metrics

See this annotated report for a description of the metrics.

A note on additional metrics

The set of metrics output by quaqc are those that I personally have found useful. If you can think of additional important ones, please create an issue on GitHub and I would be happy to implement them if it seems feasible to do so.

Input preparation and compatibility

Tn5 offsets

In quaqc 1.7, the default Tn5 offset is +4/-4: forward-strand 5-prime coordinates are shifted by +4 bp and reverse-strand coordinates by -4 bp. Earlier versions used +4/-5. To reproduce the earlier offset, add --tn5-rev 5 to your command.

The offsets are configurable with --tn5-fwd and --tn5-rev. Enable shifting for the relevant output using --tss-tn5, --bedGraph-tn5, --bed-tn5, or --quant-tn5; setting an offset alone does not enable it. The --nfr and --footprint presets enable --tss-tn5 automatically.

Information on input BAMs

quaqc is meant to work with unfiltered BAMs produced by sequence alignment tools such as bowtie2. However, some processing steps are recommended and/or required. These include:

• samtools fixmate so that quaqc can compute accurate mate and fragment size metrics for PE experiments (optional, though not having proper read mate data may require using the --use-nomate flag to stop quaqc from filtering all PE reads)
• samtools sort if the reads are not already coordinate sorted (required)
• samtools markdup (or an equivalent command) so that quaqc can compute read duplication metrics (optional)
• samtools index to generate an accompanying index file (optional; quaqc will create one if missing)

My personal recommendation is to run all of these commands simultaneously during sequence alignment. For example:

[bowtie2/etc ... ] \
  | samtools fixmate -m -u - - \
  | samtools sort -u - \
  | samtools markdup --write-index - file.bam

This requires a sequence aligner which can output its results to stdout, as well as being grouped by read name. For further information check out this article.

CRAM

CRAM support depends on the linked HTSlib and reference availability. See INSTALL for optional compression libraries and libcurl support when references must be retrieved remotely.

Long reads

Long-read DNA-seq BAMs have had limited testing. quaqc is designed for short-read data, and its handling of supplementary alignments limits its suitability for long-read QC.

Spliced reads

Spliced alignments, such as those from RNA-seq, are unsupported. Alignment and fragment statistics should not be relied on for these inputs.