Detecting candidate chromosomal inversions with radr

Scope

The analytical starting points are Li and Ralph (2019) for local PCA, Faria et al. (2019) and Wellenreuther and Bernatchez (2018) for inversion biology and interpretation, and empirical work combining genotype groups, LD, recombination, diversity, and structural evidence. These references should be read before treating a local-structure signal as an inversion.

detect_inversions() is a screening tool for finding genomic regions whose local population structure is unusual relative to the rest of a GDS dataset. It is designed for diploid, biallelic SNP data, including carefully filtered RADseq and similar reduced-representation datasets.

The function detects candidate inversion-associated haploblocks. SNP data can show the expected population-genomic footprint of an inversion, but they do not directly demonstrate reversed chromosome orientation or locate physical breakpoints. Use long reads, split reads, discordant read pairs, linkage maps, comparative assemblies, or cytogenetics for structural confirmation.

Caution: do not scan unfiltered or batch-confounded data. Structured missingness, uneven depth, allele dropout, paralogs, close relatives, and plate, lane, library, sequencing-run, or genotype-calling effects can create local PCA clusters, heterozygosity differences, and extended LD that closely resemble an inversion signal. Perform dataset-level quality control first, avoid confounding biological groups with technical batches, and cross-tabulate every candidate genotype group against the available technical metadata. detect_inversions() is not a substitute for data filtering or batch assessment.

Screen before LD pruning. Extended LD is part of the signal being investigated. Run the first inversion screen before LD pruning, then use the inferred regions and arrangements in explicit sensitivity analyses.

How `radr` builds on previous inversion work

radr::detect_inversions() builds on analytical principles explored by Li and Ralph (2019), who formalised local PCA for detecting changes in population structure along a genome; Huang et al. (2020), who demonstrated that candidate inversions can be recovered from reduced-representation SNP data; and Pearse et al. (2019) and Akopyan et al. (2025), who showed the value of combining regional genotype groups, linkage disequilibrium, recombination, diversity, and structural evidence when interpreting inversion-associated haploblocks.

The function is an independent implementation rather than a wrapper around another inversion package. It was designed to make these transferable ideas practical for the quality-controlled RADseq and similar datasets already used in radr. In particular, detect_inversions():

reads genotypes directly from GDS;
constructs every window within a chromosome, linkage group, or scaffold;
offers fixed-SNP, fixed-base-pair, and experimental LD-scaled windows;
combines local covariance PCA with contiguous-window candidate detection;
evaluates regional PCA groups quantitatively instead of treating a requested three-group k-means solution as evidence by itself;
summarises heterozygosity, overall and within-group LD, candidate boundaries, flanking windows, internal transitions, and window-size sensitivity;
returns ordered putative arrangement genotypes, relative call confidence, arrangement-specific whitelists, and homokaryotype-only sensitivity sets;
accepts centromere, low-recombination, repeat, and assembly-gap annotations without allowing them to determine the candidate calls;
explicitly records the consequences of missing-data filtering and temporary mean imputation for RADseq interpretation; and
writes reproducible tables and publication-ready ggplot2 diagnostics to a standard results folder while returning all results for further exploration.

This integrated workflow is the main advantage of using radr: it moves from a broad local-structure scan to an auditable candidate evidence table without claiming more resolution than the marker data provide. The output is intended to guide the next experiment, not replace it. Linkage mapping can test recombination suppression; haplotagging or other linked-read approaches can phase the alternative haplotypes and refine structural hypotheses; and ONT or PacBio long reads, breakpoint PCR, comparative assemblies, or cytogenetics can provide physical confirmation. Haplotagging has recovered long population haplotypes and inversion-associated barcode-sharing patterns at scale (Meier et al. 2021), but exact or repetitive breakpoints may still require continuous long reads or a junction-specific assay.

What is a chromosomal inversion?

An inversion occurs when a chromosome segment is reversed. Individuals may carry two copies of one arrangement, one copy of each arrangement, or two copies of the alternative arrangement. These are often described as two homokaryotypes and one heterokaryotype.

Recombination is commonly reduced between alternative arrangements in heterokaryotypes. Consequently, alleles across a large interval can remain associated as a haplotype. This can protect combinations of locally adaptive alleles from being broken apart by recombination, especially when populations exchange migrants. Inversions can therefore maintain genomic differentiation even when much of the genome shows weak population structure.

Three-panel diagram showing a chromosome segment reversing orientation, formation of an inversion loop in a heterokaryotype, and the linkage disequilibrium, local PCA groups, and possible population FST peak detectable from SNP data.

Not every inversion is adaptive, and not every long differentiated haploblock is an inversion. Selection, centromeres, low-recombination regions, recent admixture, family structure, paralogs, assembly errors, marker-density changes, and technical batches can generate partially similar patterns.

Types of chromosomal inversion and how to distinguish them

The classical distinction depends on the relationship between the inverted segment and the centromere, not on the size of the candidate or the number of genotype clusters.

Three-panel diagram showing paracentric-compatible, pericentric-compatible, and undetermined candidate inversion classifications relative to the centromere.

The labels in this schematic deliberately describe compatibility, not a structurally confirmed inversion type. Classification requires credible breakpoints and an independently supported centromere position. A threshold-defined candidate core should remain undetermined when either is uncertain. For a more detailed illustration of inversion-loop pairing and meiotic products, see the figure by Bérénice Bougas in Wellenreuther and Bernatchez (2018), https://doi.org/10.1016/j.tree.2018.04.002.

Inversion type	Structural definition	Classical meiotic consequence in a heterokaryotype	Evidence needed for classification
Paracentric	Both breakpoints occur on the same chromosome arm, so the centromere is outside the inverted segment.	A crossover within the inversion loop can produce dicentric and acentric recombinant chromatids. Their reduced recovery contributes to apparent recombination suppression.	Two validated breakpoints on the same annotated chromosome arm.
Pericentric	The breakpoints occur on opposite chromosome arms, so the inverted segment includes the centromere.	A crossover within the inversion loop can produce recombinant chromatids carrying duplications and deletions. Reduced recovery can suppress observed recombination.	Validated breakpoints on opposite sides of a reliably annotated centromere.

These are descriptions of simple inversions. Real chromosomes can instead contain nested, overlapping, adjacent, or serial inversions, and an inversion may be accompanied by duplications, deletions, translocations, repetitive sequence, or assembly differences. A broad population-genomic signal can therefore represent one simple inversion, several linked rearrangements, or a larger structurally complex haploblock.

What evidence separates the types?

Classification normally requires all of the following:

A chromosome-level reference framework. The candidate must be placed on an assembled chromosome with credible chromosome-arm and centromere annotation. The largest observed SNP position is not sufficient.
Breakpoint evidence. Split reads, discordant read pairs, long reads, breakpoint-spanning assemblies, optical maps, or breakpoint PCR should identify the two orientation changes. Repetitive breakpoint regions may require long continuous reads or assembly comparison.
Centromere context. The validated breakpoints must then be compared with centromere coordinates. Breakpoints on one arm support a paracentric classification; breakpoints on opposite arms support a pericentric one.
Independent consistency checks. Linkage maps, pedigrees, cytogenetics, synteny with another assembly, and patterns of dosage or read depth can test the proposed structure and expose more complex alternatives.

Centromere annotations are uncertain in many non-model fish assemblies. Merely overlapping a putative centromeric region, showing suppressed recombination, or changing chromosome-arm proportions is supportive context rather than decisive classification. The reference assembly and coordinate version should always be reported.

Local PCA, three arrangement-like clusters, heterozygosity, LD, and windowed population differentiation can identify an inversion-associated haploblock, but they do not reveal chromosome orientation and cannot by themselves distinguish a paracentric inversion from a pericentric inversion. Until both breakpoints and the centromere relationship are established, report the result as a candidate region or putative inversion rather than assigning an inversion type.

Evidence expected from SNP data

A convincing candidate normally combines several lines of evidence:

Local structure: contiguous windows show a population structure that differs from the genomic background.
Three genotype groups: regional PCA may show two homokaryotype groups and an intermediate heterokaryotype group.
Heterozygosity: the intermediate group is often more heterozygous within the region than the two outer groups.
Extended LD: SNPs across the region show elevated linkage disequilibrium.
Genomic continuity: the signal spans adjacent markers and does not merely follow an isolated outlier SNP.
Technical independence: the signal does not track plates, lanes, libraries, missingness, coverage, or other processing variables.
Population differentiation: when biologically defined populations differ in arrangement frequency, a broad windowed F_{ST} peak overlaps the local structure signal and helps prioritise the region for follow-up.

None of these criteria alone proves an inversion. For example, LD may be uneven inside an old inversion because gene conversion and double crossovers permit some gene flux between arrangements.

Lessons from empirical fish examples

Published fish studies illustrate both the value and the limits of an inversion screen. detect_inversions() does not reproduce their complete analyses: several relied on whole-genome sequence, linkage maps, pedigrees, long reads, or phenotype data that are not present in a typical RADseq GDS. Instead, the function uses the parts of their reasoning that transfer to marker data and reports candidates for further investigation.

Rainbow trout: a structurally validated double inversion

Pearse et al. (2019) characterized a roughly 55-Mb double-inversion supergene on rainbow trout chromosome Omy05. Linkage mapping showed almost complete recombination suppression in heterokaryotypic parents, while comparative genomics and long-read assembly helped resolve the structure. The two arrangements were associated with migratory tendency, with effects that depended on sex and dominance.

This example motivates several choices in detect_inversions():

evidence should extend across a contiguous haploblock rather than a single strongly associated marker;
regional genotype groups, LD, heterozygosity, and boundary changes should be interpreted together;
one broad signal may contain adjacent or nested rearrangements and should not automatically be described as one simple inversion; and
association with sex, migration, or another phenotype can support biological relevance, but it does not establish chromosome orientation.

The outer candidate coordinates returned by detect_inversions() therefore describe a marker-supported interval. Internal changes in window scores, LD, or clustering may justify examining more than one rearrangement inside that interval.

Chinook salmon: a major-effect locus is not necessarily an inversion

Thompson et al. (2020) showed that variation near GREB1L is strongly associated with adult migration timing in Chinook salmon. This is an important counterexample for interpreting a genomic peak: a locus with a large phenotypic effect, strong differentiation, or extended haplotypes is not by itself evidence of a chromosomal inversion.

Consequently, detect_inversions() does not use phenotype association or an isolated differentiation peak as sufficient evidence. A narrow selected region may be biologically important without producing the extended local-PCA, three-group, heterozygosity, and LD pattern expected from an inversion haploblock.

Atlantic silverside: distinguish inversions from centromeres

Akopyan et al. (2025) compared 168 Atlantic silverside genomes from four populations using whole-genome variation and recombination maps. Large, abruptly bounded differentiation haploblocks coincided with known inversions and showed the characteristic three tight regional-PCA clusters. Narrower differentiation peaks frequently coincided with putative centromeres. Regional PCA in centromeric regions could also reflect reduced recombination, but the individuals were more dispersed and retained more haplotype variation than in the inversion regions.

The study also demonstrates why relative differentiation alone can mislead. Elevated F_{ST} near centromeres resulted partly from low within-population diversity, whereas absolute sequence divergence (d_{XY}) was elevated in the large inversions but reduced near centromeres. For radr, the transferable diagnostic lessons are to:

compare broad, abrupt and contiguous haploblocks with narrow or gradual peaks;
examine the compactness and separation of the three inferred genotype groups, not merely the existence of a regional PCA pattern;
compare candidates with known or putative centromeres, recombination maps, marker density, repeats, and assembly gaps; and
avoid calling a low-recombination region an inversion without independent structural or linkage evidence.

Sparse, ascertained RADseq markers generally do not support the same robust windowed d_{XY} analysis as whole-genome sequence. detect_inversions() therefore does not manufacture an absolute-divergence statistic from insufficient data. Where dense sequence data and defensible population groups are available, F_{ST}, diversity, and d_{XY} provide valuable downstream validation alongside the radr candidate scan.

Cross-reference candidates with chromosome-resolved F_{ST}

Population differentiation is most informative when it is plotted along each linkage group rather than reduced to one genome-wide value. The complementary assigner::fst_WC84() function calculates Weir and Cockerham’s (1984) F_{ST} from independently defined populations and, when CHROM and numeric POS metadata are available, returns:

fst.linkage.groups, with one summary per chromosome or linkage group;
fst.windows, with WC84 variance components combined in physical windows; and
fst.genome.plot, a chromosome-concatenated Manhattan-style figure.

The physical window size should reflect marker density and the expected scale of the candidate. For sparse RADseq data, a 25-kb window may contain too few markers; larger windows and a sensitivity analysis across several window sizes are usually more defensible. Use the same genome assembly, chromosome labels, position metadata, sample filters, and marker filters in both analyses.

library(ggplot2)

strata <- readr::read_tsv("strata.tsv", show_col_types = FALSE)

fst <- assigner::fst_WC84(
  data = genome,
  strata = strata,
  linkage.group = TRUE,
  window.size = 1e6,
  window.method = "fixed",
  filename = "inversion_fst_crosscheck"
)

# Whole-genome overview supplied by assigner
fst$fst.genome.plot

# Linkage-group view with radr candidate intervals shaded
candidate_intervals <- inv$candidates |>
  dplyr::transmute(
    CHROM = as.character(chromosome),
    xmin = start,
    xmax = end
  )

ggplot(fst$fst.windows, aes(WINDOW_MID, FST_WC84)) +
  geom_rect(
    data = candidate_intervals,
    aes(xmin = xmin, xmax = xmax, ymin = -Inf, ymax = Inf),
    inherit.aes = FALSE,
    fill = "#E69F00",
    alpha = 0.18
  ) +
  geom_line(colour = "grey45") +
  geom_point(aes(colour = OUTLIER), size = 1.4) +
  facet_wrap(~ CHROM, scales = "free_x") +
  scale_x_continuous(
    labels = scales::label_number(scale = 1e-6, suffix = " Mb")
  ) +
  labs(
    x = "Genomic position",
    y = expression("Windowed WC84 " * F[ST]),
    colour = "Exploratory outlier"
  ) +
  theme_bw()

An overlapping, spatially coherent F_{ST} peak strengthens the case that a candidate haploblock contributes to differentiation among the sampled populations. Its absence does not reject an inversion: both arrangements may occur at similar frequencies among populations. Conversely, an F_{ST} peak can arise from local selection, a centromere, low recombination, reduced diversity, marker ascertainment, or technical structure without an inversion.

Population strata must be specified independently of the inversion scan. Calculating F_{ST} between the regional PCA clusters inferred from the same candidate genotypes is circular and should not be presented as independent support. Where whole-genome sequence is available, compare relative differentiation with d_{XY}, within-population diversity, recombination, and structural evidence, as illustrated by Akopyan et al. (2025).

Rainbow trout across populations: inversion signals are context-dependent

Campbell, Anderson, Garza, and Pearse (2021) compared independently landlocked rainbow trout populations with an anadromous source population. All landlocked populations showed an increased frequency of the large Omy05 inversion, while three of four showed an increased frequency of the Omy20 inversion. Outside these major rearrangements, however, the genomic response to similar selection was much less parallel.

This work reinforces two principles used in detect_inversions(). First, a long inversion-associated haplotype can represent standing adaptive variation that changes frequency repeatedly across populations. Second, neither the strength of a genomic signal nor its association with migration should be assumed to be constant across the species range. Candidate clusters should therefore be cross-tabulated against population and environment, and biological associations should be tested rather than assigned from the candidate region alone.

How the local analysis is organised

Windows are created within chromosome, linkage-group, or scaffold boundaries. If LG1 ends with 40 unused SNPs and LG2 follows in the GDS, those 40 SNPs are never combined with the first 60 SNPs of LG2. This prevents an artificial window spanning unrelated genomic regions.

The analysis is not one PCA per linkage group. Instead:

each linkage group is divided into SNP windows;
a separate local covariance PCA is calculated for every window;
each window is represented by its leading covariance axes;
distances between window representations are calculated;
multidimensional scaling (MDS) places similar windows together and unusual windows farther from the genomic background;
contiguous unusual windows are joined into candidate regions;
a regional PCA, clustering, heterozygosity summary, and LD summary are calculated for each candidate.

The window calculations are independent. parallel.core can distribute them across separate R workers. Each worker opens its own read-only GDS connection; workers do not share one SeqArray connection. A progress bar reports completed windows in sequential and parallel mode. Because every worker constructs a sample covariance matrix, more workers also require more memory. Start with two or four workers rather than automatically using every available core.

When several linkage groups are scanned, their window summaries are compared in the same MDS analysis. The returned table retains the linkage-group label, so chromosome-position plots should normally be faceted by linkage group. A regional PCA belongs to one candidate on one linkage group and therefore does not itself require an LG facet.

As a complementary diagnostic, chromosome.pca = TRUE also calculates one independent PCA for every chromosome or linkage group. These chromosome-wide PCAs do not select candidates. After a candidate has been inferred, its AA, AB, and BB labels are projected as colours onto every panel. This makes it easy to ask whether the same individual grouping is localized to the candidate linkage group or recurs throughout the genome.

Run a scan from GDS

The examples are not evaluated because they refer to user files.

library(radr)

genome <- genometranslator::read_genome("filtered_dataset.gds")

sample_metadata <- readr::read_tsv("maintained_sample_metadata.tsv")

# Stage 1: ordinary discovery scan with the defaults
screen <- detect_inversions(
  data = genome,
  strata = sample_metadata,
  parallel.core = 4
)

screen$candidates
screen$chromosome.lengths

# Stage 2: a more detailed genome-wide scan
inv <- detect_inversions(
  data = genome,
  strata = sample_metadata,
  window.snps = 100,
  step.snps = 100,
  sensitivity.window.snps = c(100, 250, 500, 1000),
  min.call.rate = 0.90,
  min.candidate.windows = 2,
  stability.replicates = 100,
  parallel.core = 4,
  chromosome.pca = TRUE,
  chromosome.pca.max.snps = 2000,
  return.ld = TRUE
)

inv
inv$candidates
inv$arrangement.genotypes
inv$homokaryotype.whitelist
inv$homokaryotype.all.candidates
inv$path.folder

The first call is the normal starting point. It scans the filtered genome before LD pruning and reports candidate regions without requiring the user to choose a large set of tuning values. The second call adds SNP-window sensitivity, assignment stability, retained LD matrices, and more parallel workers. These extra diagnostics are valuable, but they need not be paid for before the basic signal has been located.

strata requires INDIVIDUALS. Its rows act as a sample whitelist and its other columns remain descriptive metadata. Include useful variables such as STRATA, sequencing batch, library, lane, plate, extraction method, caller, sampling year, or other project-specific factors. The function compares these variables with regional PC1 and inferred arrangements, but does not copy them into or use them to modify the GDS.

Window-size sensitivity and SNP-resampling stability are optional because both add computation. sensitivity.window.snps = NULL and stability.replicates = 0 skip them. A rapid first screen can use those defaults; a focused rerun of promising chromosomes can then request both. Chromosome-wide PCA also adds computation. Set chromosome.pca = FALSE for a minimal screen, or reduce chromosome.pca.max.snps when chromosomes contain very dense whole-genome data. SNPs are selected at evenly distributed marker indices within each chromosome so the cap preserves chromosome-wide coverage.

Known centromeres, assembly gaps, or other low-recombination annotations can be supplied without allowing those annotations to determine which windows are selected:

known_regions <- data.frame(
  chromosome = c("LG5", "LG14", "LG14"),
  start = c(22000000, 18000000, 41000000),
  end = c(26000000, 23000000, 42500000),
  type = c("putative_centromere", "low_recombination", "assembly_gap")
)

inv <- detect_inversions(
  data = genome,
  window.snps = 100,
  step.snps = 50,
  known.regions = known_regions
)

The known_region_overlap column lists overlapping annotation types. An overlap is a warning for interpretation, not evidence for or against an inversion, and it does not change candidate selection or the evidence score.

Chromosome size and candidate extent

When the GDS retains the original VCF contig dictionary, chromosome lengths are read automatically. Candidate tables then report candidate_span_bp, chromosome_length_bp, chromosome_fraction, chromosome_percent, and the left and right flanking lengths. The complete length provenance is returned in inv$chromosome.lengths and written to chromosome_length_context.tsv.

If declared lengths are absent from the GDS, provide either a named vector or a table with CHROM and LENGTH columns:

lengths <- readr::read_tsv("chromosome_lengths.tsv")

inv <- detect_inversions(
  data = genome,
  chromosome.lengths = lengths
)

A reference FASTA path or its existing .fai index can instead be supplied with reference.genome. The function reads the index rather than scanning the complete FASTA:

inv <- detect_inversions(
  data = genome,
  reference.genome = "genome.fasta"
)

Without declared sequence lengths, the largest observed marker position is used only as a labelled fallback. It underestimates chromosome length and can therefore overestimate the percentage occupied by a candidate.

To investigate a known chromosome-14 signal without letting other linkage groups define the background:

chr14 <- detect_inversions(
  data = genome,
  strata = sample_metadata,
  chromosome = "14",
  window.snps = 100,
  step.snps = 50,
  min.call.rate = 0.90,
  min.candidate.windows = 2
)

Once this focused call confirms the signal, compare independent primary window definitions. These calls genuinely recall the candidate boundaries:

chr14_50_25 <- detect_inversions(
  data = genome, strata = sample_metadata, chromosome = "14",
  window.snps = 50, step.snps = 25, parallel.core = 4
)

chr14_100_50 <- detect_inversions(
  data = genome, strata = sample_metadata, chromosome = "14",
  window.snps = 100, step.snps = 50, parallel.core = 4
)

chr14_250_125 <- detect_inversions(
  data = genome, strata = sample_metadata, chromosome = "14",
  window.snps = 250, step.snps = 125, parallel.core = 4
)

Compare candidate start, end, span, chromosome percentage, and reciprocal overlap across these runs. The shared interval is a defensible stable candidate core; the union describes boundary uncertainty. Do not select only the run that produces the cleanest or widest interval.

sensitivity.window.snps is complementary. It summarizes PC1 variance and LD at additional scales within one run, but it does not independently rerun candidate selection at every size:

chr14_sensitivity <- detect_inversions(
  data = genome,
  strata = sample_metadata,
  chromosome = "14",
  window.snps = 100,
  step.snps = 50,
  sensitivity.window.snps = c(50, 100, 250, 500),
  stability.replicates = 100,
  parallel.core = 4
)

After boundary sensitivity, repeat the candidate analysis with stricter call-rate thresholds and alternative sample sets, particularly after removing close relatives or suspect sequencing batches. Only then should stable marker boundaries be compared with recombination maps, structural-variant callers, long reads, read-pair evidence, assemblies, or breakpoint assays.

Overlapping windows (step.snps < window.snps) can help localise transitions, but adjacent results are then strongly dependent. Coordinates should not be reported as precise breakpoints merely because a score changes between two overlapping windows.

For physical windows instead of equal-SNP windows, supply window.bp. This is particularly useful for checking whether a candidate is robust to variation in RAD marker density:

chr14_mb <- detect_inversions(
  data = genome,
  chromosome = "14",
  window.bp = 1e6,
  step.bp = 5e5,
  min.window.snps = 20
)

Experimental LD-scaled windows adapt their number of markers to local adjacent LD. They help assess whether a signal depends on an arbitrary fixed marker count, but they do not estimate an inversion breakpoint or a formal LD block:

chr14_ld <- detect_inversions(
  data = genome,
  chromosome = "14",
  window.method = "ld",
  ld.window.threshold = 0.10,
  ld.window.min.snps = 50,
  ld.window.max.snps = 500
)

Plot window evidence by linkage group

detect_inversions() automatically saves the standard window-score, MDS, LD, call-rate, regional-PCA, heterozygosity, and candidate-LD figures. It also returns the ggplot2 objects in inv$output.files$plots, so they can be changed without rerunning the genomic calculations.

library(ggplot2)

ggplot(
  inv$windows,
  aes(x = (start + end) / 2, y = robust_score,
      colour = candidate_window)
) +
  geom_point() +
  geom_hline(
    yintercept = inv$settings$score.threshold,
    linetype = 2
  ) +
  facet_wrap(~ chromosome, scales = "free_x") +
  scale_x_continuous(labels = scales::label_number(scale = 1e-6,
                                                    suffix = " Mb")) +
  labs(
    x = "Window midpoint",
    y = "Robust local-structure score",
    colour = "Candidate"
  ) +
  theme_bw()

This is the appropriate LG-faceted view: window score versus genomic position. The MDS coordinates can also be inspected, but faceting them by linkage group answers a different question because MDS axes represent similarity among windows rather than physical position.

ggplot(inv$windows, aes(MDS1, MDS2, colour = chromosome,
                        shape = candidate_window)) +
  geom_point(size = 2) +
  theme_bw()

Inspect candidate genotype groups

Each candidate has a matching entry in inv$diagnostics. A regional PCA plot shows individuals, not windows:

candidate1 <- inv$diagnostics[[1]]

ggplot(candidate1$scores, aes(PC1, PC2, colour = arrangement)) +
  geom_point(size = 2, alpha = 0.8) +
  labs(colour = "Putative arrangement genotype") +
  theme_bw()

candidate1$cluster_summary
candidate1$cluster_models
candidate1$assignment_stability
candidate1$metadata_audit
candidate1$metadata_contingency
candidate1$scores |>
  dplyr::select(individual, arrangement, arrangement_dosage,
                arrangement_confidence, assignment_stability)

When sample metadata are supplied, each candidate also writes a *_metadata_by_arrangement.tsv table and one regional-PCA figure per usable metadata variable. These outputs show directly whether inferred arrangements are concentrated in particular populations, projects, plates, libraries, or lanes. A strong association produces a console warning. This is a confounder audit, not a correction: strata does not residualize or otherwise adjust the local PCA. Where biological groups and wet-laboratory batches are confounded, repeat the scan in balanced subsets or within independently replicated groups before interpreting the candidate as an inversion.

Compare candidate groups across linkage groups

The chromosome-PCA overview uses the arrangement calls from one candidate only as colours. It does not rerun clustering or force three groups independently on every linkage group. Every panel is calculated from genotypes on that linkage group, using at most chromosome.pca.max.snps evenly distributed SNPs.

inv$chromosome.pca$summary
inv$output.files$plots$`INV-CAND-1_chromosome_pca`

Separation that is strongest on the candidate linkage group and weak elsewhere supports a localized haploblock interpretation. Similar separation across many linkage groups instead suggests genome-wide population structure, family structure, admixture, or a technical batch effect that needs to be resolved before describing the region as a putative inversion.

The reverse is not a formal rejection test. A short candidate interval can be diluted by thousands of unrelated SNPs in a whole-linkage-group PCA, so weak or absent separation in the candidate linkage-group panel does not invalidate a clear regional PCA. Interpret this overview together with the window scan, regional PCA, heterozygosity, LD, missingness, depth, sample metadata, and window-size sensitivity.

For candidates with quantitative three-cluster support, AA, AB, and BB mean the low, middle, and high regional-PC1 clusters. They are convenient putative arrangement genotypes, not validated breakpoint genotypes. AA and BB do not identify which arrangement is reference, ancestral, derived, or physically inverted. Reversing the arbitrary sign of PC1 can exchange the two outer labels without changing the biological result.

Candidates without quantitative three-cluster support use neutral Group 1, Group 2, and Group 3 labels. Their numeric cluster IDs remain in exported tables, but arrangement dosage and homokaryotype classification are withheld. Relative confidence describes separation from the next algorithmic cluster; it is not a posterior probability. For supported candidates, the result folder contains one whitelist per arrangement, one homokaryotype-only whitelist per candidate, and a whitelist of individuals classified as a homokaryotype across every candidate.

The homokaryotype sets are useful for asking whether a downstream signal persists after excluding putative heterokaryotypes. Haplotype scans need further care: selscan, for example, does not accept missing genotype or haplotype data. Do not convert missing calls silently into a homozygous arrangement.

Three clusters are biologically suggestive only when the central PCA cluster also has the expected heterozygosity and the pattern is not explained by population, family, or batch. Two clusters may occur when one arrangement is rare or absent, while continuous PCA scores may indicate ordinary population structure rather than a polymorphic inversion.

The requested cluster.k controls the reported arrangement calls, but cluster_models compares one-, two-, and three-cluster descriptions of PC1 using an approximate Gaussian BIC. recommended_cluster_k in the candidate table records the lowest-BIC model. This check prevents a three-group story from being accepted only because k-means was asked for three centres.

When stability.replicates > 0, regional SNPs are resampled and the PCA and clustering are repeated. assignment_stability is the proportion of repeated calls that agree after accounting for the arbitrary direction of PC1. Low overall stability or unstable individuals near cluster boundaries should make arrangement labels provisional.

The LD diagnostic uses the upper triangle for all individuals and the lower triangle for the more common outer PCA cluster, treated provisionally as a homokaryotype. This comparison is useful because mixing arrangements can create strong LD even when LD is lower within one arrangement. It is not independent validation: the subset was inferred from the same candidate-region genotypes. The group-specific values are available in candidate1$ld_summary.

Coverage, allele balance, and missingness artifacts

Technical summaries are generated from whatever compatible count information was retained in the GDS, not from the original filename or input format. The function recognizes standard DP and biallelic AD nodes and genometranslator genotype-metadata fields named READ_DEPTH, ALLELE_REF_DEPTH, and ALLELE_ALT_DEPTH. These may originate from a VCF, DArT two-row count file, or another count-based format.

candidate1$coverage_source
candidate1$technical_summary

candidate1$scores |>
  dplyr::select(individual, arrangement, call_rate, mean_depth,
                mean_heterozygote_allele_balance)

inv$output.files$plots$`INV-CAND-1_coverage_call_rate`

A depth shift among arrangement classes may reflect a real structural variant, but it can also indicate paralogy, a copy-number variant, collapsed repeats, or mapping bias. Allele balance far from 0.5 among called heterozygotes is another warning. These measurements help distinguish a local genotype signal from a coverage artifact; they do not prove which biological mechanism produced it.

Diagnostic SNP loadings and arrangement differentiation

Regional PC1 loadings rank the SNPs contributing most strongly to the local structure. They can guide diagnostic-panel exploration, but marker selection and evaluation must use separate samples or resampling to avoid upward bias.

candidate1$pc1_loadings |>
  dplyr::slice_min(loading_rank, n = 25)

candidate1$arrangement_differentiation
inv$output.files$plots$`INV-CAND-1_marker_loadings`

The arrangement table reports mean absolute allele-frequency differences and pairwise Hudson FST. Hudson FST is used only as a compact description of allele-frequency separation between inferred classes. AA, AB, and BB are not natural populations, and the same regional SNPs helped define the classes. Consequently, this FST is partly circular and is not independent evidence for an inversion. Population differentiation should instead be estimated from independently defined populations with an estimator such as assigner::fst_WC84().

Dxy is not calculated among AA, AB, and BB. In particular, AB is an inferred heterokaryotype rather than an independently sampled evolutionary lineage. Population-level Dxy belongs in a later comparison among genuine populations or lineages.

Mean imputation: what it does

Within each window, SNPs below min.call.rate are excluded. At each remaining SNP, a missing genotype is temporarily replaced by the observed mean allele dosage for that SNP. The value is used only to make covariance PCA matrices complete; the GDS is not changed. After marker centring, an imputed value contributes zero deviation at that SNP. It therefore supplies no evidence that the individual carries either arrangement.

Only covariance PCA uses these imputed values. LD is calculated from observed genotypes using pairwise-complete correlations; missing LD genotypes are not mean-imputed.

This is intentionally simple and transparent. It prevents PCA from discarding every individual with one missing call, but it is not a correction for biased missingness.

Why RADseq missingness needs special attention

RADseq missingness is often structured rather than random. Relevant causes include:

polymorphism at restriction sites and consequent allele dropout;
uneven sequencing depth among samples or libraries;
lane, plate, library-preparation, or genotyping batches;
population differences in reference divergence or mapping quality;
separate marker discovery or genotype calling among sample groups;
paralogous loci and collapsed repetitive regions;
marker density that varies with restriction sites, assembly quality, or filtering.

If one population or batch is missing preferentially at a set of linked SNPs, mean imputation pulls those samples toward the local PCA centre. Depending on the pattern, this can weaken a biological cluster, manufacture an intermediate group, or make a batch-specific genomic region look unusual.

For every candidate, perform sensitivity analyses:

plot marker and individual call rate, depth, and heterozygosity;
cross-tabulate inferred clusters against population, plate, lane, library, family, and sampling date;
repeat the scan with stricter min.call.rate values;
repeat after excluding the most incomplete samples;
compare fixed-SNP and fixed-base-pair windows when marker density is uneven;
inspect whether the signal remains when close relatives are removed;
verify that marker density and assembly gaps do not define the candidate boundaries.

A robust biological signal should not disappear under every reasonable QC choice, and its clusters should not be synonymous with one technical batch.

Results folder and customised plots

Each call creates a dated detect_inversions folder following the standard radr workflow. It contains the recorded call, TSV tables, and PNG and PDF figures. Timestamped folders prevent a normal run from overwriting an earlier analysis. The returned object also contains the window table, candidate table, PCA scores, cluster summaries, LD summaries, and plot objects. A standard plot can therefore be adjusted and saved under an additional name:

p <- inv$output.files$plots$window_scores +
  labs(title = "Candidate inversion scan after batch QC")

ggsave(file.path(inv$path.folder, "window_scores_annotated.png"), p,
       width = 10, height = 7, dpi = 300)

Use save.plots = FALSE when only tables and returned objects are wanted. The results folder and reproducibility tables are still created.

Add context to genome scans

Use genome_scan_context() to place scan peaks beside marker density, missingness, heterozygosity, allele frequency, LD, and candidate-region annotations:

context <- radr::genome_scan_context(
  data = genome,
  inversion.regions = inv,
  scan.statistics = bayescan_or_fst_table,
  window.snps = 250
)

context$context

For each regional peak, ask:

Does it overlap low recombination or a centromere?
Is marker density unusual?
Does missingness increase in the same region?
Is the signal restricted to one caller, mapping strategy, sequencing batch, or reference assembly?
Does it remain after excluding putative inversion heterozygotes?
Is it supported by a method based on a different genomic signature?

These questions should also guide analyses in assigner and grouper. Population structure should be explored with the full marker set, with candidate inversion regions excluded, and with arrangement genotypes included explicitly. Otherwise a large haploblock can dominate clustering and make one regional feature look like genome-wide structure.

Interpreting coordinates

The candidate start and end values are the outer marker coordinates of the selected windows. They depend on window size, step size, marker density, missingness filters, and the reference assembly. Report them as an approximate inversion-associated interval. With RADseq data, there may be a substantial gap between the last marker outside the signal and the first marker inside it.

A defensible report might state:

A candidate inversion-associated haploblock was detected from approximately X to Y Mb on LG14, supported by concordant local-PCA, regional genotype-group, heterozygosity, and LD patterns. These coordinates describe the marker-defined interval and not validated structural breakpoints.

Interpreting the candidate evidence table

The candidate table combines continuous diagnostics with a deliberately simple screening grade. Important columns include:

cluster_separation: the smallest gap between adjacent regional-PC1 cluster centres, divided by the pooled within-cluster standard deviation;
cluster_compactness: the proportion of regional PC1 variation explained by the inferred clusters;
recommended_cluster_k: the one-, two-, or three-cluster PC1 model with the lowest approximate BIC;
assignment_stability: agreement of arrangement calls across SNP-resampling replicates, when requested;
smallest_cluster_n and smallest_cluster_frequency: protection against treating a few outlying individuals as an inversion arrangement;
middle_heterozygosity_excess: middle-cluster heterozygosity minus the mean of the two outer clusters;
regional_mean_ld_r2, homokaryotype_mean_ld_r2, and heterokaryotype_mean_ld_r2: LD across all individuals and within inferred genotype groups;
flanking_mean_ld_r2 and boundary_contrast: comparison with the immediate noncandidate windows;
internal_transition_max: the largest score change between candidate windows, useful for finding complex, adjacent, or nested signals; and
known_region_overlap: any user-supplied centromere, low-recombination, assembly-gap, or other annotation intersecting the candidate;
candidate_class: a cautious interpretation such as putative inversion-associated haploblock, unresolved haploblock, low-recombination or centromeric candidate, technical or assembly-associated region, or another structural-variation-associated haploblock; and
alternative_explanations: reminders that local PCA can reflect a centromere, low recombination, assembly or mapping problems, introgression, population-specific missingness, or another structural variant.

three_cluster_evidence is not merely a record that k-means was run with cluster.k = 3. It additionally requires all three inferred groups to contain at least three samples, the smallest group to represent at least 5% of samples, and adequate separation relative to within-group spread.

evidence_score assigns one point for each of five observations: quantitative three-cluster support, positive middle-cluster heterozygosity excess, positive boundary contrast, regional LD above flanking LD, and continuity across at least two windows. Scores of 0–2, 3–4, and 5 are labelled weak, moderate, and strong, respectively. This grade prioritises candidates for inspection; it is not a posterior probability and does not prove an inversion.

Recommended workflow

Use the scan as one part of a staged analysis:

complete sample, marker, batch, depth, and missingness QC;
scan chromosome-specific windows for unusual local structure;
inspect contiguous candidates rather than isolated outlier windows;
evaluate regional PCA groups, heterozygosity, and LD jointly;
test sensitivity to filtering, window size, step size, and relatedness;
compare candidates with centromeres, assembly gaps, repeats, and gene annotations;
use linkage maps or haplotagging to test recombination and phase the alternative arrangements; and
seek breakpoint-spanning long reads, junction PCR, comparative assembly, or cytogenetic evidence before calling a structurally confirmed inversion.

For genome scans and population-structure analyses, report three deliberate views whenever the data permit:

the complete genome;
a collinear sensitivity dataset excluding candidate inversion-associated or major low-recombination regions; and
candidate-region analyses using arrangement frequencies, heterozygosity, LD, differentiation, environmental associations, and homokaryotype-only sensitivity datasets.

Inversions should be detected, annotated, and analysed explicitly. They should neither be ignored nor automatically removed. Use “candidate”, “putative”, or “inversion-associated haploblock” until structural evidence justifies a more definitive term.

References

Akopyan M, Tigano A, Jacobs A, Wilder AP, Therkildsen NO (2025) Genetic differentiation is constrained to chromosomal inversions and putative centromeres in locally adapted populations with higher gene flow. Molecular Biology and Evolution, 42, msaf092.

Campbell MA, Anderson EC, Garza JC, Pearse DE (2021) Polygenic basis and the role of genome duplication in adaptation to similar selective environments. Journal of Heredity, 112, 614-625.

Bhatia G, Patterson N, Sankararaman S, Price AL (2013) Estimating and interpreting FST: the impact of rare variants. Genome Research, 23, 1514-1521.

Faria R, Johannesson K, Butlin RK, Westram AM (2019) Evolving inversions. Trends in Ecology & Evolution, 34, 239-248.

Gosselin T, Anderson EC, Bradbury I (2020) assigner: assignment analysis with GBS/RAD data using R. R package. https://doi.org/10.5281/zenodo.592677.

Hoffmann AA, Rieseberg LH (2008) Revisiting the impact of inversions in evolution: from population genetic markers to drivers of adaptive shifts and speciation? Annual Review of Ecology, Evolution, and Systematics, 39, 21-42.

Hudson RR, Slatkin M, Maddison WP (1992) Estimation of levels of gene flow from DNA sequence data. Genetics, 132, 583-589.

Huang K, Andrew RL, Owens GL, Ostevik KL, Rieseberg LH (2020) Multiple chromosomal inversions contribute to adaptive divergence of a dune sunflower ecotype. Molecular Ecology, 29, 2535-2549.

Kess T, Bentzen P, Lehnert SJ, Sylvester EVA, Lien S, Kent MP, Sinclair-Waters M, Morris CJ, Regular P, Fairweather R, Bradbury IR (2019) A migration-associated supergene reveals loss of biocomplexity in Atlantic cod. Science Advances, 5, eaav2461.

Li H, Ralph P (2019) Local PCA shows how the effect of population structure differs along the genome. Genetics, 211, 289-304.

Meier JI, Salazar PA, Kučka M, Davies RW, Dréau A, Aldás I, Power OB, Nadeau NJ, Bridle JR, Rolian C, Barton NH, McMillan WO, Jiggins CD, Chan YF (2021) Haplotype tagging reveals parallel formation of hybrid races in two butterfly species. Proceedings of the National Academy of Sciences of the United States of America, 118, e2015005118.

Mérot C (2020) Making the most of population genomic data to understand the importance of chromosomal inversions for adaptation and speciation. Molecular Ecology, 29, 2513-2516.

Pearse DE, Barson NJ, Nome T, Gao G, Campbell MA, Abadía-Cardoso A, Anderson EC, Rundio DE, Williams TH, Naish KA, Moen T, Liu S, Kent M, Moser M, Minkley DR, Rondeau EB, Brieuc MSO, Sandve SR, Miller MR, Cedillo L, Baruch K, Hernandez AG, Ben-Zvi G, Shem-Tov D, Barad O, Kuzishchin K, Garza JC, Lindley ST, Koop BF, Thorgaard GH, Palti Y, Lien S (2019) Sex-dependent dominance maintains migration supergene in rainbow trout. Nature Ecology & Evolution, 3, 1731-1742.

Thompson NF, Anderson EC, Clemento AJ, Campbell MA, Pearse DE, Hearsey JW, Kinziger AP, Garza JC (2020) A complex phenotype in salmon controlled by a simple change in migratory timing. Science, 370, 609-613.

Weir BS, Cockerham CC (1984) Estimating F-statistics for the analysis of population structure. Evolution, 38, 1358-1370.

Wellenreuther M, Bernatchez L (2018) Eco-evolutionary genomics of chromosomal inversions. Trends in Ecology & Evolution, 33, 427-440. https://doi.org/10.1016/j.tree.2018.04.002