Skip to content

Advanced R reference

A practical reference for DNA sequence analysis, reproducibility, and idiomatic R. Use the tabs below to navigate by topic, and the site search to find specific functions or packages.

Read and manipulate sequences with Biostrings

DNAStringSet is the workhorse for handling large collections of sequences without explicit loops. Install via Bioconductor.

# BiocManager::install("Biostrings")
library(Biostrings)

seqs <- readDNAStringSet("sequences.fasta")

# Basic properties
length(seqs)          # number of sequences
width(seqs)           # per-sequence lengths (integer vector)
names(seqs)           # sequence identifiers

# Subsetting
seqs[1:10]
seqs[width(seqs) > 200]   # filter by length

# GC content, complement, reverse complement
letterFrequency(seqs, letters = "GC", as.prob = TRUE)
reverseComplement(seqs)

Bioconductor install

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("Biostrings")

Pattern matching and primer counting

Search for primers, motifs, or adapter sequences across many reads in one call. IUPAC ambiguity codes (M, Y, R, ...) work directly in patterns.

primer <- DNAString("GTGYCAGCMGCCGCGGTAA")

# Count occurrences per read (returns integer vector)
vcountPattern(primer, seqs, max.mismatch = 1)

# Return match positions as IRanges
vmatchPattern(primer, seqs, max.mismatch = 1)

# Trim fixed flanking sequences (e.g. adapters)
trimLRPatterns(
  Lpattern = "AGATCGGAAG",
  subject  = seqs
)

Pairwise and multiple alignment

Align sequences globally or locally; useful for OTU validation or phylogenetic prep.

# Multiple alignment with DECIPHER
# BiocManager::install("DECIPHER")
library(DECIPHER)

aligned <- AlignSeqs(seqs)
BrowseSeqs(aligned)   # interactive viewer in browser

# StaggerAlignment handles variable-length amplicons (ITS, 16S)
aligned <- AlignSeqs(seqs, StaggerAlignment = TRUE)

# Pairwise with Biostrings
pairwiseAlignment(
  seqs[[1]], seqs[[2]],
  type               = "global",
  substitutionMatrix = nucleotideSubstitutionMatrix(match = 1, mismatch = -3)
)

Working with FASTQ and quality scores

Read, filter, and summarise FASTQ files without leaving R. For very large files use FastqStreamer to process in chunks.

library(ShortRead)

fq  <- readFastq("reads.fastq.gz")
sread(fq)    # DNAStringSet of read sequences
quality(fq)  # PhredQuality object

# Per-base quality as numeric matrix (reads x positions)
qmat <- as(quality(fq), "matrix")
rowMeans(qmat)   # mean Q score per read

# Filter low-quality reads
keep     <- rowMeans(qmat) >= 20
fq_clean <- fq[keep]
writeFastq(fq_clean, "reads_clean.fastq.gz")

# Process a large file in chunks
strm <- FastqStreamer("big.fastq.gz", n = 10000)
repeat {
  chunk <- yield(strm)
  if (length(chunk) == 0) break
  # ... process chunk ...
}
close(strm)

Building and filtering phyloseq objects

Combine OTU table, taxonomy, and sample metadata into one object for downstream diversity analysis.

library(phyloseq)

ps <- phyloseq(
  otu_table(otu_mat,  taxa_are_rows = TRUE),
  tax_table(tax_mat),
  sample_data(meta_df)
)

# Always verify orientation before row-wise operations
taxa_are_rows(ps)   # must be TRUE when taxa are rows

# Prevalence filter: keep taxa present in >= 2 samples
prev    <- apply(otu_table(ps), 1, function(x) sum(x > 0))
ps_filt <- prune_taxa(prev >= 2, ps)

# Relative abundance transformation
ps_rel  <- transform_sample_counts(ps, function(x) x / sum(x))

# PCoA ordination
ord <- ordinate(ps_rel, "PCoA", "bray")
plot_ordination(ps_rel, ord, color = "treatment")

Taxa orientation

phyloseq does not enforce a fixed orientation. Always confirm taxa_are_rows(ps) before apply(otu_table(ps), 1, ...). A transposed table produces wrong results silently.

Project startup checklist

Establish structure before writing any analysis code. Consistent layout makes projects navigable by collaborators and by your future self.

usethis::create_project("myproject")  # creates RStudio .Rproj and base dirs
renv::init()                          # package snapshot from the start

# Initialise git in the project root
git init
git add renv.lock .Rprofile .gitignore
git commit -m "initial project setup"

Recommended directory layout:

project/
├── data/          # raw, read-only input files
├── R/             # reusable functions (sourced or packaged)
├── reports/       # .Rmd files
├── results/       # output tables, figures
├── _targets.R     # pipeline definition
├── renv.lock
└── README.md

Tip

Keep data/ read-only. All derived files go in results/. This makes it safe to delete results/ and re-run the pipeline from scratch.

renv: reproducible package libraries

Snapshot exact package versions so analyses run identically months later, on another machine, or in a fresh container.

# One-time setup per project
renv::init()       # initialises renv, creates renv.lock
renv::snapshot()   # record current installed state

# On a collaborator's machine or CI
renv::restore()    # install exact versions from renv.lock

# After updating a package
install.packages("vegan")
renv::snapshot()   # update the lock file

What to commit to git

Commit renv.lock and .Rprofile. Add renv/library/ to .gitignore. The library itself is not portable; only the lock file is.

Seed management and randomness

Any function that uses randomness (clustering, permutation tests, UMAP, subsampling) must have an explicit seed recorded in the script.

set.seed(123)   # set before every stochastic operation

For parallel code, do not rely on set.seed() alone. Use furrr_options(seed = TRUE) instead (see Performance & scaling).

Seed-fishing

Do not choose seeds by running multiple values and keeping whichever gives the most attractive result. Document the seed in the methods section exactly as used.

R Markdown reports with knitr caching

Weave code and prose so every figure is always generated from current data. Use chunk caching carefully: invalidate the cache when the input file changes.

# YAML header (save as analysis.Rmd)
---
title: "16S analysis - project p1208"
date: "`r Sys.Date()`"
output:
  html_document:
    code_folding: hide
    toc: true
    toc_float: true
---

# Global options at the top of the document
knitr::opts_chunk$set(
  echo       = TRUE,
  warning    = FALSE,
  message    = FALSE,
  cache      = TRUE,
  fig.width  = 7,
  fig.height = 5,
  dpi        = 150
)

# Invalidate cache automatically when input file changes
# {r load-data, cache=TRUE, cache.extra=file.info("otu_table.tsv")$mtime}
otu <- read_tsv("otu_table.tsv")

Child documents

Split long analyses into child .Rmd files and include them with `r knitr::knit_child("methods.Rmd")` to keep individual files manageable.

targets pipeline framework

Define the analysis as a directed acyclic graph of targets. Only steps whose inputs have changed are re-run.

# _targets.R (project root)
library(targets)

# Declare packages needed by all targets
tar_option_set(
  packages = c("phyloseq", "vegan", "dplyr", "readr")
)

list(
  tar_target(raw_seqs,  read_fasta("data/seqs.fasta")),
  tar_target(filtered,  filter_length(raw_seqs, min = 200)),
  tar_target(otu_table, run_vsearch(filtered)),
  tar_target(ps,        build_phyloseq(otu_table, "data/meta.tsv")),
  tar_target(fig_pcoa,  plot_pcoa(ps), format = "file")
)

# Running and inspecting the pipeline
tar_make()                    # run outdated targets
tar_visnetwork()              # interactive DAG in browser
tar_read(ps)                  # load a target result into the session
tar_outdated()                # list targets that need updating
tar_meta(fields = error)      # inspect errors from failed targets
tar_destroy()                 # wipe stored results and force a full re-run

Note

tar_option_set(packages = ...) avoids hard-to-debug missing-package errors when targets run in fresh R sessions. Always declare dependencies here rather than relying on the interactive session state.

Session info and provenance

Always append complete session information at the end of reports and scripts. This is the minimum required for a published or archived analysis.

# Concise output (recommended)
sessioninfo::session_info()

# Base R version
sessionInfo()

# Save to file for archiving alongside results
writeLines(
  capture.output(sessioninfo::session_info()),
  "session_info.txt"
)

# Record git commit hash
system("git rev-parse HEAD", intern = TRUE)

In R Markdown

Add a collapsible session info block at the end of every .Rmd using a plain HTML <details> element, which renders fine in knitr HTML output:

<details><summary>Session info</summary>
```{r session-info, echo=FALSE}
sessioninfo::session_info()
```
</details>

Memory-efficient file reading

Many bioinformatics datasets become memory-bound before CPU-bound. Choose the right reader for the file size and access pattern.

# vroom: fast reading of large flat files; uses ALTREP for lazy evaluation
library(vroom)
df <- vroom("large_table.tsv")   # column types inferred; reads are deferred

# arrow: columnar format, no need to load the full file into RAM
library(arrow)

# Single Parquet file
df <- read_parquet("results/otu_table.parquet")

# Entire directory of Parquet files as one virtual dataset
ds <- open_dataset("results/")
ds |>
  filter(depth > 1000) |>
  group_by(sample) |>
  summarise(mean_abund = mean(abundance)) |>
  collect()    # only materialise the filtered result

Convert once, read many times

Convert large TSV outputs to Parquet once and read Parquet in all downstream steps. Parquet is ~3-10x smaller than gzipped TSV and much faster to query.

library(arrow)
otu <- read_tsv("otu_table.tsv")
write_parquet(otu, "results/otu_table.parquet")

Parallel computing

R is single-threaded by default. The future backend plus furrr gives parallel purrr-style iteration with minimal code changes.

library(future)
library(furrr)

# Check how many cores are safely available, then use that number
workers <- parallelly::availableCores()
plan(multisession, workers = workers)

# plan(multicore,  workers = workers)  # fork-based; faster but Linux/macOS only
# plan(cluster, ...)                   # SLURM / HPC: see future.batchtools

# Parallel map: same interface as purrr::map
results <- future_map(
  files,
  process_sample,
  .options = furrr_options(seed = TRUE)   # reproducible parallel RNG
)

# Combine list of data frames into one (modern replacement for future_map_dfr)
dfs <- future_map(samples, load_metadata) |> list_rbind()

# Typed scalar variant still fine
counts <- future_map_int(samples, count_reads)

# Reset to sequential when done
plan(sequential)

When to use each backend:

Backend OS Notes
multisession all separate processes; safe with Bioconductor
multicore Linux / macOS only forks the current process; faster startup
cluster all explicit socket cluster; needed on Windows HPC
future.batchtools HPC submits SLURM / PBS jobs transparently

Bioconductor and forking

Some Bioconductor packages (notably those with C-level global state) are not fork-safe. Prefer multisession when working with Biostrings, ShortRead, or phyloseq in parallel.

Seed handling in parallel

set.seed() before future_map() does not guarantee reproducible results across workers. Always pass .options = furrr_options(seed = TRUE).

Profiling and optimisation

Profile before optimising. Guessing which line is slow is almost always wrong.

# profvis: interactive flame graph in the RStudio viewer
library(profvis)
profvis({
  run_pipeline()
})

# bench::mark: accurate microbenchmark with memory tracking
library(bench)
bench::mark(
  base_apply = apply(mat, 1, mean),
  rowMeans   = rowMeans(mat),
  check      = TRUE     # verify both return the same result
)

# system.time: quick one-off check, no repetition
system.time(run_pipeline())

Typical workflow:

  1. Profile the full pipeline with profvis to find the actual bottleneck.
  2. Benchmark alternative implementations of that function with bench::mark.
  3. Consider data.table, arrow, or parallelisation only after profiling confirms the need.

Tip

bench::mark() runs each expression multiple times and reports median time plus memory allocation. Prefer it over system.time() for anything you plan to optimise.

Efficient wrangling with data.table

For data frames beyond roughly 100 MB, data.table gives 10-100x speedups. dtplyr lets you write dplyr syntax compiled to data.table operations.

library(data.table)

# Read large TSV fast
dt <- fread("large_otu_table.tsv")

# Filter and aggregate without copies
dt[depth > 1000,
   .(mean_abund = mean(abundance)),
   by = .(sample, taxon)]

# dtplyr: dplyr interface, data.table engine
library(dtplyr)
lazy_dt(dt) |>
  filter(depth > 1000) |>
  group_by(sample) |>
  summarise(mean_abund = mean(abundance)) |>
  as_tibble()

Note

Profile with profvis first. Most data frames under ~20 MB need no special handling.

Functional iteration with purrr

Replace fragile for-loops with typed map functions. Errors propagate clearly and the return type is guaranteed.

library(purrr)

# map + list_rbind: modern replacement for the superseded map_dfr()
results <- map(sample_files, function(f) {
  df <- read_tsv(f)
  summarise(df, n = n(), mean_depth = mean(depth))
}) |> list_rbind()

# safely(): catch errors without stopping the loop
safe_read <- safely(read_tsv)
raw     <- map(files, safe_read)
errors  <- keep(raw, ~ !is.null(.x$error))
results <- map(raw, "result") |> compact()

# walk2(): side effects (saving files, writing output)
walk2(plot_list, names(plot_list), \(p, nm)
  ggsave(paste0("results/", nm, ".pdf"), p, width = 7, height = 5))

Tidy eval: writing dplyr-powered functions

Pass column names as arguments to functions using {{ }} (embrace operator) and all_of() for character vectors of names.

library(dplyr)

# {{ }} embraces a single column name passed as a symbol
summarise_by <- function(df, group_col, value_col) {
  df |>
    group_by({{ group_col }}) |>
    summarise(
      mean = mean({{ value_col }}, na.rm = TRUE),
      .groups = "drop"
    )
}
summarise_by(meta, treatment, alpha_div)

# all_of(): select from a character vector of names
diversity_cols <- c("Shannon", "Simpson", "Chao1")
meta |> select(sample_id, all_of(diversity_cols))

# .data pronoun: access columns by string variable
col_name <- "Shannon"
meta |> filter(.data[[col_name]] > 2)

Non-standard evaluation quick reference:

Situation Syntax Example
Single column as argument {{ col }} group_by({{ col }})
Character vector of names all_of(vec) select(all_of(cols))
String variable as column .data[[str]] filter(.data[[col]] > 0)
Inject name on left of := "{name}" := mutate("{nm}" := x + 1)

Data validation with checkmate

Catch bad input data immediately rather than letting silent errors propagate through a long analysis. checkmate gives dramatically better error messages than stopifnot().

library(checkmate)

# Validate function arguments at entry
filter_taxa <- function(ps, min_prev) {
  assert_true(inherits(ps, "phyloseq"))
  assert_number(min_prev, lower = 0)

  prune_taxa(
    apply(otu_table(ps), 1, function(x) sum(x > 0)) >= min_prev,
    ps
  )
}

# Validate a data frame after import
assert_data_frame(meta)
assert_names(names(meta), must.include = c("sample_id", "treatment", "depth"))
assert_numeric(meta$depth, lower = 0, any.missing = FALSE)
assert_subset(meta$treatment, c("control", "treatment"))

assertr for pipeline assertions

checkmate excels at function-argument validation. assertr is better for inline pipeline checks where you want assertions to flow naturally with |>.

library(assertr)
meta |>
  assert(in_set("control", "treatment"), group) |>
  assert(within_bounds(0, Inf), depth) |>
  assert(not_na, sample_id)

Testing with testthat

Write tests for any function you reuse across projects or share with others. Even a handful of tests catches most regressions.

library(testthat)

# tests/testthat/test-filter.R
test_that("filter_length removes short sequences", {
  seqs <- DNAStringSet(c("AAAA", "AAAAAAAAA", "AAAAAAAAAAAAAAAA"))
  out  <- filter_length(seqs, min = 9)
  expect_length(out, 2)
  expect_true(all(width(out) >= 9))
})

test_that("filter_length returns DNAStringSet", {
  seqs <- DNAStringSet(c("AAAA", "AAAAAAAAA"))
  expect_s4_class(filter_length(seqs, min = 5), "DNAStringSet")
})

# Run all tests
testthat::test_dir("tests/testthat/")

Snapshot tests for figures and tables

Use expect_snapshot() or the vdiffr package to test that plots and summary tables do not change unexpectedly between code changes.

Package development basics

When a collection of functions is reused across projects, turning it into a package pays off quickly: documentation is co-located with code, and devtools::load_all() replaces scattered source() calls.

library(usethis)
library(devtools)

# Create package skeleton
usethis::create_package("mybioinf")

# Development cycle
devtools::load_all()    # reload all R/ functions without restarting
devtools::test()        # run testthat suite
devtools::check()       # full R CMD check (CRAN-style)
devtools::document()    # regenerate NAMESPACE and .Rd files from roxygen2

# Document functions with roxygen2 comments directly above each function
#' Filter sequences by minimum length
#'
#' @param seqs A \code{DNAStringSet}.
#' @param min  Minimum length in base pairs (integer).
#' @return A filtered \code{DNAStringSet}.
#' @export
filter_length <- function(seqs, min) {
  seqs[width(seqs) >= min]
}

Debugging

A small toolkit covers most situations. The key habit is to reach for traceback() first, then narrow down with debugonce() before committing to a full browser() session.

# After an error: see the full call stack
traceback()

# tidyverse / rlang errors: richer trace with pipe context
rlang::last_trace()

# Enter an interactive browser at the point of error, one call at a time
debugonce(my_function)   # arms the debugger for the next call only
my_function(args)        # drops into browser() on entry

# Unconditional breakpoint: pause inside a function at a specific line
# Add browser() directly in the function body, then reload with load_all()
filter_taxa <- function(ps, min_prev) {
  browser()              # execution pauses here
  prune_taxa(...)
}

# Pause on any error globally (useful for tracing deep pipelines)
options(error = recover)   # interactive frame selector on error
options(error = NULL)      # restore default behaviour afterwards

Browser commands (type at the Browse[1]> prompt):

Command Action
n next line
s step into function call
f finish current loop or function
c continue to next browser() or end
Q quit debugger
where print call stack

Debugging inside targets

targets catches errors per target. Use tar_meta(fields = error) to read the message, then tar_load(upstream_target) to load inputs into the session and reproduce the failure interactively.

Logging for long-running pipelines

Interactive scripts can use cli for readable progress output. For pipelines that write to file or run on HPC, use logger.

# cli: alerts for key events
library(cli)
cli::cli_alert_info("Loading {length(files)} files")
cli::cli_alert_success("Filtering complete: {nrow(ps)} taxa retained")
cli::cli_alert_warning("Low sequencing depth in sample {sample_id}")

# cli: progress bar for loops over many items
cli::cli_progress_bar("Processing samples", total = length(files))
for (f in files) {
  process_sample(f)
  cli::cli_progress_update()
}
cli::cli_progress_done()

# purrr equivalent: progress bar via .progress argument (purrr >= 1.0)
results <- map(files, process_sample, .progress = "Processing samples")

# logger: structured logging to console and file simultaneously
library(logger)

log_appender(appender_tee("pipeline.log"))   # console + file
log_threshold(INFO)

log_info("Pipeline started: {length(files)} samples")
log_warn("Sample {sid} has depth {depth}, below threshold")
log_error("Failed to read file: {path}")

Tip

Add log_info(paste("targets pipeline run:", Sys.time())) at the top of _targets.R so every tar_make() call is timestamped in the log.