Also works with any unspliced DNA-seq experiment! Compatible with plant genomes!!
Requires gcc/clang and GNU Make, tested with macOS and Linux. Basic install:
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.
To see a basic description of all available parameters, run ./quaqc -h
. For
more detailed help, see the manual page.
quaqc comes built-in with sensible defaults, meaning few parameters must be specified for regular usage. (The hardcoded defaults can be changed before compilation by editing quaqc.h.) To run quaqc using the 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 command changes the location of the output QC report to the current directory,
and uses the peak and TSS files to generate additional stats. Since this example
BAM is just a subset of reads contained within a small region of the Arabidopsis
genome, the target.bed
file restricts quaqc to only consider that region of the genome (and thus adjust
how it calculates the stats). See a copy of the output
report here.
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.)
For those who prefer to analyze their data with R, a companion package
quaqcr
is available. This package contains
functions for both executing quaqc from within R as well as analyzing all of
its outputs made available when outputting the results in JSON format. See
the README for a detailed tutorial.
Please open an issue on GitHub.
No article about quaqc is yet published. Until then, please cite this repository if you find this software useful in your research.
Tremblay, B.J.M. (2024). quaqc: Quick ATAC-seq QC (Version 0.1) [Computer software]. https://github.com/bjmt/quaqc
quaqc is fairly performant, but it can still require some patience for decently
sized BAMs (e.g., about one minute for 80 million reads from Arabidopsis ATAC-seq).
To substantially speed up quaqc to iterate through different filtering
thresholds, make use of the -n/--target-names
to restrict quaqc to only
look at single chromosomes. This is often enough to get an idea of the quality of
the data while getting a several fold speed up (e.g., restricting quaqc to the
first chromosome of Arabidopsis takes just over 10 seconds for the previously
mentioned example).
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
Now, six different MAPQ thresholds have been tested in the time it would take
to process the entire BAM file once across all chromosomes. (The -i/--title
flag is used to add a unique title to each run.)
The TSS pileup functionality of quaqc can be used to instead generate motif footprints. Note that these are only useful for testing, as no correction of base composition is peformed to reduce Tn5 transposase insertion biases. The output must be saved in JSON format to recover the footprint data, which can then be plotted with an external program.
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
See this annotated report for a description of the 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.
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.
quaqc only works with indexed BAMs, though actually indexing them before using quaqc is optional. If no index can be found then quaqc will attempt to create the index.
As far as I can tell CRAM files work fine as well, as long as the reference is available locally. See the notes in the installation instructions about compiling quaqc with libcurl if you need to work with CRAMs depending on remote references.
Some quick testing shows that quaqc seems to work with DNA-seq long read BAMs too. However if you try it out and see obvious errors in the output, please open an issue and I will do my best to fix the problem. Overall this is probably the wrong tool to use to run QC on long reads. Keep in mind quaqc will not properly handle supplementary alignments.
I have not tested quaqc with BAMs containing spliced reads (e.g. RNA-seq), but I expect most of the statistics (such as alignment and fragment lengths) would be incorrectly calculated. Even the number of reads may be wrongly counted, since quaqc is not set up to properly handle supplementary alignments.