fraq: A high-throughput extensible toolkit for processing fastq data

fraq is a high-throughput extensible toolkit for processing fastq data. The goal of this package is to empower users to quickly build out programmatic ‘kernels’ to define any FASTQ processing task they may need. fraq then takes those kernels and handles I/O, compression and multithreading. It builds on Intel TBB’s flow graph to orchestrate concurrency and data processing; throughput can be as fast as compression and disk speed allow.

The package ships with a suite of predefined ‘kernels’ for common FASTQ tasks, detailed in this vignette.

Why use fraq?

Quick start

Install

if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}
BiocManager::install("fraq")

Pre-built functions

fraq ships with a collection of ready-to-go kernels that cover common preprocessing steps:

Walkthroughs with synthetic data

For illustration purposes, we create a small synthetic dataset.

set.seed(314156)
example_dir <- file.path(tempdir(), "fraq_function_examples")
dir.create(example_dir, showWarnings = FALSE, recursive = TRUE)
R1 <- file.path(example_dir, "example_R1.fastq")
R2 <- file.path(example_dir, "example_R2.fastq")
generate_random_fastq(c(R1, R2), n_reads = 2000, read_length = 150)
example_reads <- c(R1, R2)

Summarization

fraq_summary() rolls up QC tables for the R1/R2 fastq pairs.

library(fraq)
# summarize quality metrics
qc <- fraq_summary(c(R1, R2))
Quick QC overview from fraq_summary().

Filtering

Various filtering operations are illustrated below. Here we downsample to reduce dataset size and then discard mates whose mean PHRED drops below 22 or accumulate more than 6 bases with PHRED < 18 (roughly 10% of reads will be filtered).

filter_dir <- file.path(example_dir, "filtering")
dir.create(filter_dir, showWarnings = FALSE)

downsampled_reads <- file.path(
    filter_dir,
    c("example_R1_ds.fastq.gz", "example_R2_ds.fastq.gz")
)
fraq_downsample(example_reads, downsampled_reads, amount = 0.50, nthreads = 2L)

quality_reads <- file.path(
    filter_dir,
    c("example_R1_q20.fastq.gz", "example_R2_q20.fastq.gz")
)
fraq_quality_filter(
    input = downsampled_reads,
    output = quality_reads,
    min_mean_quality = 22,
    max_low_q_bases = 6L,
    low_q_threshold = 18L,
    nthreads = 2L
)

filtered_stats <- fraq_summary(quality_reads, nthreads = 2L)$basic_stats_R1
filtered_stats

##   total_sequences total_bases seq_len_min seq_len_mean seq_len_max gc_percent
## 1             970      145500         150          150         150   50.00619

Splitting

Splitting workflows either direct reads into barcode-specific files, chunk long runs into bite-sized batches, or simply inspect barcode usage.

split_dir <- file.path(example_dir, "splitting")
if (dir.exists(split_dir)) unlink(split_dir, recursive = TRUE)
dir.create(split_dir, showWarnings = FALSE)
barcode_set <- c("ACGTAA", "TTGGCC")
demux_patterns <- file.path(
    split_dir,
    c("R1_{barcode}.fastq.gz", "R2_{barcode}.fastq.gz")
)

Barcode-guided demultiplexing

Demux looks for barcode/adapter/primer sequence at the start of the first given fastq file.

fraq_demux(
    input = example_reads,
    output_format = demux_patterns,
    barcodes = barcode_set,
    max_distance = 1L,
    nthreads = 2L
)
basename(sort(
    list.files(split_dir, pattern = "R1_.*\\.fastq.gz$", full.names = TRUE)
))

## [1] "R1_ACGTAA.fastq.gz"   "R1_NO_MATCH.fastq.gz" "R1_TTGGCC.fastq.gz"

Fixed-size chunking

fraq_chunk splits reads into fixed-size batches with incremental file names.

chunk_prefix <- file.path(split_dir, "chunk_demo")
fraq_chunk(
    input = example_reads[1],
    output_prefix = chunk_prefix,
    output_suffix = "gz",
    chunk_size = 200,
    nthreads = 2L
)
basename(sort(
    list.files(
        split_dir,
        pattern = "chunk_demo_chunk.+\\.fastq.gz$",
        full.names = TRUE
    )
))

##  [1] "chunk_demo_chunk0.fastq.gz" "chunk_demo_chunk1.fastq.gz"
##  [3] "chunk_demo_chunk2.fastq.gz" "chunk_demo_chunk3.fastq.gz"
##  [5] "chunk_demo_chunk4.fastq.gz" "chunk_demo_chunk5.fastq.gz"
##  [7] "chunk_demo_chunk6.fastq.gz" "chunk_demo_chunk7.fastq.gz"
##  [9] "chunk_demo_chunk8.fastq.gz" "chunk_demo_chunk9.fastq.gz"

Barcode counting

Barcode counting looks for sequence substrings anywhere in the fastq reads and outputs a data frame of counts.

fraq_count_barcodes(
    input = example_reads,
    barcodes = barcode_set,
    max_distance = 0L,
    allow_revcomp = FALSE,
    nthreads = 2L
)

##       barcode count
## 1    NO_MATCH  1734
## 2      TTGGCC   132
## 3      ACGTAA   126
## 4 MULTI_MATCH     8

Modification

Format conversion, adapter trimming, consensus merging, etc.

mod_dir <- file.path(example_dir, "modification")
if (dir.exists(mod_dir)) unlink(mod_dir, recursive = TRUE)
dir.create(mod_dir, showWarnings = FALSE)

Convert between formats

Re-wrap the paired-end files in Zstandard-compressed FASTQ.

converted_fastq <- file.path(
    mod_dir,
    c("example_R1.fastq.zst", "example_R2.fastq.zst")
)
fraq_convert(example_reads, converted_fastq, nthreads = 2L)

Concatenate files

Combine the converted shards into a single gzip-compressed stream.

concatenated <- file.path(mod_dir, "example_all.fastq.gz")
fraq_concat(converted_fastq, concatenated, nthreads = 2L)

Merge overlapping pairs

Generate consensus single-end reads while keeping optional unmerged outputs.

merged_reads <- file.path(mod_dir, "example_merged.fastq")
unmerged_reads <- file.path(
    mod_dir,
    c("example_unmerged_R1.fastq", "example_unmerged_R2.fastq")
)
fraq_merge_pairs(
    input = example_reads,
    output_merged = merged_reads,
    output_unmerged = unmerged_reads,
    min_overlap = 20L,
    max_mismatch_rate = 0.05,
    nthreads = 2L
)

## $merged_reads
## [1] 648
## 
## $unmerged_reads
## [1] 1352
## 
## $mean_insert_size
## [1] 267.7454
## 
## $sd_insert_size
## [1] 10.25771
## 
## $mean_overlap
## [1] 32.25463
## 
## $mean_mismatch_rate
## [1] 0

Trim adapters

Remove adapter prefixes from R1 (and optionally drop untrimmed reads).

trimmed_reads <- file.path(
    mod_dir,
    c("example_trimmed_R1.fastq", "example_trimmed_R2.fastq")
)
fraq_trim_adapters(
    input = example_reads,
    output = trimmed_reads,
    adapters = "ACTAC",
    max_distance = 1L,
    filter_untrimmed = FALSE,
    nthreads = 2L
)

##      adapter count
## 1 NO_ADAPTER  1973
## 2      ACTAC    27

Supported formats

fraq chooses formats from file names, so the extension you supply controls how data is decoded and encoded.

Extension system: flow graph overview

The fraq_run() pipeline wires a TBB flow graph so that IO and data flow happen concurrently with data processing. A high-level view looks like this:

Each block of fastq is read from Primary Reader (for the first mate R1) and Secondary Readers (for R2 or additional fastq files). The key node in this graph is the Process Kernel. It takes fastq records (i.e. R1 and R2 reads), processes or alters them, decides whether to keep them or not (filtering) and outputs any number of fastq records (demux and splitting). This simple pattern naturally supports lots of different fastq processing operations and can be customized.

For R kernels, the Process Kernel step uses a different execution path, described below.

Extension system: writing a custom kernel in R

You can build a custom kernel with fraq_run_r(). The R function runs on the main R thread rather than inside a TBB graph node, so it is safe for it to work with R objects and call back into R. When nthreads > 1, IO happens on background threads while the R kernel stays on the main R thread. When nthreads = 1, or when fraq detects that it is running inside a forked R process, fraq uses a serial path that does not construct a TBB graph.

Blocks are delivered to the R kernel in increasing block-index order, and the index vector within each call is increasing.

An R kernel is called as kernel(reads, index):

Vectorized R kernels can still perform well on large FASTQ datasets because they operate on full blocks rather than one read at a time.

input_paths <- c("input_R1.fastq.gz", "input_R2.fastq.gz")
output_paths <- c("even_R1.fastq.gz", "even_R2.fastq.gz")

even_read_kernel <- function(reads, index) {
    keep <- index %% 2 == 0
    filtered_read1 <- reads$read1[keep, , drop = FALSE]
    filtered_read2 <- reads$read2[keep, , drop = FALSE]

    output <- list()
    output[[output_paths[1]]] <- filtered_read1
    output[[output_paths[2]]] <- filtered_read2
    output
}

fraq_run_r(
    input_paths,
    even_read_kernel,
    nthreads = 2L
)

Do not use parallel::mclapply() inside a fraq_run_r() kernel. On Unix-like systems it forks the R process, and forking while fraq has active background IO and compression threads can deadlock or crash.

Extension system: writing a custom kernel with Rcpp

If the prebuilt kernels are insufficient, you can write your own via an Rcpp script. You supply a lambda to fraq::run, and the runtime handles all batching, IO, and parallelism. More information can be found in ?fraq_rcpp_template.

Below is an example that keeps only reads whose GC fraction falls in a window (default 35-65%).

// [[Rcpp::depends(fraq)]]
#include <Rcpp.h>
#include <fraq.h>

double calc_gc_content(const std::string &s) {
    double gc = 0.0;
    for(char c : s) { if(c == 'G' || c == 'C') gc += 1.0 ; }
    return gc / (double) s.size();
}

// [[Rcpp::export(rng=false)]]
void fraq_gc_filter(std::vector<std::string> input,
                    std::vector<std::string> output,
                    double gc_min = 0.35, double gc_max = 0.65) {
  
    auto gc_filter_kernel = [&](fraq::input_t reads, size_t read_index)
        -> fraq::output_t {
    for(auto & read : reads) {
        double gc = calc_gc_content(read.seq);
        if(gc < gc_min || gc > gc_max) return {};
    }
    return fraq::zip(output, std::move(reads));
    };
  
    fraq::FraqRunConfig cfg;
    cfg.zstd_compress_level = 5;
    int nthreads = 4;
    fraq::run(input, gc_filter_kernel, nthreads, cfg);
}

All fraq classes are transparent structures with no private members, built on standard library types.

Compile via Rcpp::sourceCpp() then in R you can call your custom kernel as a normal function. The extensions on the output paths decide output format automatically.

input  <- c("sample_R1.fastq.gz", "sample_R2.fastq.gz")
output <- c("filtered_R1.fastq.zst", "filtered_R2.fastq.zst")
fraq_gc_filter(input, output, gc_min = 0.30, gc_max = 0.70)

Some important tips when building custom kernels

Streaming with named pipes

You can use named pipes (Linux/Mac only - Windows is not supported) to stream input and output directly into other command line programs.

Below is an example using fraq_downsample on input fastqs (random fastqs in this example) and streaming the output directly to bwa-mem2.

downsample_fifo.R

library(fraq)
generate_random_fastq("R1.fastq")
generate_random_fastq("R2.fastq")
fraq_downsample(input=c("R1.fastq", "R2.fastq"),
                output=c("ds_R1.fastq.fifo", "ds_R2.fastq.fifo"),
                amount = 0.25, nthreads = 5L)

In bash (Linux/macOS), create the named pipes first before any operations:

HG38_REF=/path/to/hg38.fa.gz
mkfifo ds_R1.fastq.fifo ds_R2.fastq.fifo
Rscript downsample_fifo.R &
bwa-mem2 mem -t 8 $HG38_REF ds_R1.fastq.fifo ds_R2.fastq.fifo > output.sam

Windows users should stick with regular files or platform-specific piping; .fifo paths are not available there.

Tuning and threading

Global knobs live behind fraq_options():

# Inspect the current block size.
fraq_options("blocksize")

# Shrink batches when running small tests.
fraq_options("blocksize", 16384L)

# Tune compression levels for new outputs.
fraq_options("zstd_compress_level", 6L)
fraq_options("gzip_compress_level", 4L)

Most kernels accept nthreads. With nthreads = 1, fraq uses a serial path. With nthreads > 1, fraq caps the TBB scheduler to the requested parallelism. If fraq detects that it is running inside a forked R process, it forces nthreads = 1 to avoid using TBB after fork().

FRAQ file format

The .fraq container stores FASTQ reads in independent blocks so that IO and block-level compression can run concurrently. Each block holds up to 65,535 reads (fraq_options("blocksize")) and stores block-level info such as zstd compression, name-prefix factoring, and optional nucleotide bit-packing.

FRAQ takes inspiration from the Nucleotide Archival Format (NAF) by concatenating nucleotide and quality payloads before compressing them with zstd (improving compression efficiency), following the approach described by Kryukov et al. (2019). FRAQ differs by block compressing the stream, which enables multithreaded compression, streaming and tailoring the layout specifically to the FASTQ format instead of being more general.

Specification

Reference: Kryukov, Kirill, et al. “Nucleotide Archival Format (NAF) enables efficient lossless reference-free compression of DNA sequences.” Bioinformatics 35.19 (2019): 3826-3828.