BIOINFORMATICS RESEARCH GUIDE

RAD-Seq Population Genomics Workflow: Restriction Site Associated DNA Sequencing

By Xenoibiotech · Last updated: October 2026

Restriction Site Associated DNA Sequencing (RAD-Seq) is a fractional genome sequencing strategy that samples reproducible homologous loci adjacent to specific restriction enzyme cut sites. Using computational pipelines like Stacks 2 and ipyrad, RAD-Seq cost-effectively identifies thousands of SNPs across hundreds of non-model individuals to resolve population structure, gene flow, and phylogenetic lineages.

Key Takeaways

  • RAD-Seq reduces genome complexity by sequencing only DNA fragments flanking specific restriction enzyme recognition motifs.
  • The method enables high-density SNP discovery across large ecological populations without requiring prior reference genome sequence.
  • Upstream demultiplexing via process_radtags rescues barcode reads and filters low-quality sequences prior to locus building.
  • Software suites like Stacks 2 construct homologous loci either de novo (ustacks, cstacks, sstacks) or aligned to a reference genome.
  • Downstream population genetic modules calculate fixation indices (FST), nucleotide diversity (π), and input matrices for demographic inference.

What Is RAD-Seq and Why Is It Preferred in Ecological Population Genetics?

In evolutionary biology, conservation genomics, and molecular ecology, researchers frequently investigate non-model species across large geographic populations (tens to hundreds of individuals). Whole-genome sequencing across hundreds of large-genome organisms (such as conifers, amphibians, or marine teleosts) remains cost-prohibitive. RAD-Seq (Restriction site Associated DNA Sequencing) resolves this economic bottleneck by sampling a consistent, homologous fraction of the genome across all individuals.

By digesting high-molecular-weight genomic DNA with rare-cutting or double-digestion restriction enzymes (e.g., SbfI, PstI, EcoRI) and ligating barcoded sequencing adapters to the cut overhangs, RAD-Seq sequences only the genomic intervals directly adjacent to the restriction cut sites. Because these cut sites are conserved across individuals, researchers sequence the exact same homologous loci repeatedly, generating thousands of informative single nucleotide polymorphisms (SNPs) at a fraction of whole-genome sequencing costs.

How Do You Design Single-Digest vs Double-Digest RAD Experiments?

RAD sequencing encompasses several protocol variants tailored to specific taxonomic genome sizes and required marker densities:

  • Single-Digest RAD (sdRAD): Uses a single rare-cutting restriction enzyme (such as SbfI, an 8-bp cutter), followed by mechanical shearing and size selection. Provides robust, long loci and minimal allele dropout.
  • Double-Digest RAD (ddRAD): Employs two different restriction enzymes (a rare cutter and a frequent cutter, such as PstI and MseI) without mechanical shearing. Library fragment sizes are strictly controlled by dual enzymatic cuts and magnetic bead size selection, lowering laboratory preparation overhead.

Targeted marker density must match the study questions: phylogeography and broad population divergence require 5,000 to 20,000 SNPs, whereas linkage disequilibrium mapping and fine-scale genome-wide association studies (GWAS) require 50,000 to 100,000+ loci.

How Are Barcodes Demultiplexed and Checked for Restriction Overhangs?

Raw multiplexed sequencing reads are demultiplexed into individual sample files using process_radtags in the Stacks software suite.

The software performs three sequential quality checks simultaneously:

  • Barcode Deconvolution: Matches sample-specific in-line barcodes, allowing up to 1 mismatch to rescue high-quality reads with minor sequencing errors.
  • Restriction Site Verification: Confirms the presence of an intact restriction enzyme recognition overhang at the beginning of the read. Reads lacking the correct restriction cut signature represent non-specific shearing artifacts and are purged.
  • Sliding-Window Quality Filtering: Discards reads whose average Phred quality score drops below 20 across a sliding window.

How Does Stacks 2 Build Loci De Novo vs Reference-Guided?

Depending on whether a reference genome exists, analysis diverges into two computational tracks in Stacks 2:

De Novo Locus Building:

  • ustacks: Aligns reads within an individual into matching stacks, building alleles and putative diploid loci based on a maximum nucleotide distance parameter (-M).
  • cstacks: Merges loci across a representative catalog of individuals to create a master catalog of population-wide loci using distance parameter (-n).
  • sstacks and tsv2bam: Matches individual loci back against the catalog and reconstructs paired-end contigs.
  • gstacks: Calls SNPs and genotypes across the population using a maximum-likelihood framework.

Reference-Guided Analysis: When an annotated reference genome is available, demultiplexed reads are mapped directly using BWA-MEM2 or Bowtie2. Coordinate-sorted BAM files are input directly into gstacks, bypassing de novo catalog construction and eliminating paralogous collapse.

How Do You Calculate Population Genetics Statistics (FST, Pi, Heterozygosity)?

The final module of the Stacks pipeline, populations, filters loci and computes population-level summary statistics:

Filtering parameters enforce data completeness across biological groups: for example, requiring a locus to be present in at least 75% of individuals (-r 0.75) across all populations (-p 3), with a minimum minor allele frequency (--min-maf 0.05).

The software outputs classic evolutionary and population metrics: observed and expected heterozygosity (Ho, He), nucleotide diversity (π), inbreeding coefficients (FIS), and pairwise population differentiation indices (FST). Furthermore, populations exports filtered genotype matrices in standard formats (VCF, Structure, PLINK, Treemix, Phylip) for downstream demographic modeling and phylogenetic reconstruction.

What Are the Major Outputs and Common Pitfalls in RAD-Seq?

A full deliverable includes demultiplexing audit summaries, locus coverage depth distributions, filtered VCF files, population genetic statistics tables, PCA and ADMIXTURE population structure plots, and neighbor-joining phylogenetic trees.

Critical analytical pitfalls include:

  • Allele Dropout from Restriction Site Mutations: Mutations within the restriction enzyme recognition site prevent digestion, causing the corresponding allele to drop out entirely and creating artificial homozygous calls. Filtering for locus presence across replicates mitigates dropout bias.
  • Over-Merging Paralogous Loci: Setting the de novo assembly distance parameter -M or -n too high clusters duplicated gene paralogs into single chimeric loci, creating spurious hyper-variable SNPs.
  • Ignoring PCR Duplicates in ddRAD: In ddRAD libraries without unique molecular identifiers (UMIs), identical reads from PCR amplification cannot be distinguished from independent biological fragments, which skews allele frequency estimates.

Apply These Skills in Your Research

Master this computational workflow step by step with our hands-on mentoring, or collaborate with our data science team for production analysis:

Learn this in our 30-day workshop: RAD-Seq Population Genomics Workshop Need this analysis done for your data? Population Genetics & RAD-Seq Services

Frequently Asked Questions

What causes allele dropout in RAD sequencing?

Allele dropout occurs when a single nucleotide mutation alters or destroys the restriction enzyme recognition sequence in one allele. The enzyme cannot cut that DNA molecule, so only the alternative allele is sequenced, falsely presenting a heterozygous individual as homozygous.

Can RAD-seq be performed without a reference genome?

Yes. The Stacks and ipyrad de novo pipelines assemble homologous loci directly from sequence overlaps across individuals without requiring any reference genome, making RAD-seq ideal for non-model organisms.

How does Stacks differentiate biological SNPs from sequencing errors?

Stacks uses a multinomial likelihood model that evaluates read depth and error rates at each nucleotide position, distinguishing genuine heterozygous base calls from random low-frequency sequencing errors.

Further Reading and Authoritative Sources

  • Rochette, N. C., Rivera-Colón, A. G., and Catchen, J. M. (2019). Stacks 2: Analytical methods for paired-end sequencing in RAD-seq. Molecular Ecology, 28(21), 4737–4754. DOI: 10.1111/mec.15201
  • Eaton, D. A., and Overcast, I. (2020). ipyrad: Interactive assembly and analysis of RADseq datasets. Bioinformatics, 36(8), 2592–2594. DOI: 10.1093/bioinformatics/btaa796
  • Andrews, K. R., et al. (2016). Harnessing the power of RADseq for ecological and evolutionary genomics. Nature Reviews Genetics, 17(2), 81–92. DOI: 10.1038/nrg.2016.29

Related Guide: WGS vs WES Comparison Guide