Genetic maps and recombination
The crossing functions — cross(),
selfcross() and double_haploid() — need to
know how often two markers recombine. This vignette explains what they
need, how synthetic_map() supplies it when you have nothing
better, and how to replace it with a real map.
Why base pairs are not enough
Marker files almost always carry physical positions in base pairs. Meiosis does not care about those. What matters is genetic distance in centiMorgans: 1 cM means a 1% chance of a crossover between two loci per generation.
The two are not proportional. Recombination is heavily suppressed near the centromere and elevated toward the chromosome ends, so two markers 10 Mb apart can be tightly linked in one region and freely recombining in another.
simplePHENOTYPES keeps the genetic map in the cm column,
alongside pos:
head(geno[, 1:5])
# snp allele chr pos cm
# 1 ss196422159 A/G 1 379844 0.0000000
# 2 ss196422171 G/A 1 613257 0.2482394
# 3 ss196422173 G/A 1 659354 0.2972609
# 4 ss196422186 G/C 1 992572 0.6515820
# 5 ss196500940 G/A 1 2044264 1.7694445
# 6 ss196422234 A/G 1 2044555 1.7697537If cm is missing, as_population() stops
rather than guessing:
no_map <- geno
no_map$cm <- NA
as_population(no_map, individuals = c("33-16", "38-11"))
# Error:
# ! The genetic map (`cm`) is all NA, so recombination distances are unknown and meiosis cannot be simulated. Build one from the physical positions with synthetic_map(), e.g.
# geno$cm <- synthetic_map(geno$chr, geno$pos)synthetic_map()
synthetic_map() turns physical positions into a
plausible centiMorgan map.
cm <- synthetic_map(geno$chr, geno$pos)
aggregate(cm ~ chr, data = data.frame(chr = geno$chr, cm = cm), FUN = max)
# chr cm
# 1 1 219.3573
# 2 2 172.2582
# 3 3 168.4068
# 4 4 175.4795
# 5 5 158.9096
# 6 6 123.2703
# 7 7 128.5558
# 8 8 128.1356
# 9 9 114.2056
# 10 10 108.5956This is a model, not a measurement. It is not derived from any published linkage map, and its values should not be reported as estimated recombination distances. Its purpose is to give a simulation a realistic recombination landscape when no real map is at hand. If you have a real map, use it (see below).
What the model does
Each marker gets a local recombination rate that dips around the centromere:
\[r(p) = 1 - s \exp\left(-\tfrac{1}{2}\left(\frac{p - c}{wL}\right)^2\right)\]
with centromere c, chromosome span L,
suppression strength s and relative width w.
That rate is integrated along the chromosome and rescaled to the target
length, so cM always increases with bp.
The consequence is the flat middle section of the curve — a large stretch of physical genome that contributes very little genetic distance:
on1 <- geno$chr == 1
plot(geno$pos[on1] / 1e6, geno$cm[on1], type = "l", lwd = 2,
xlab = "Physical position, chromosome 1 (Mb)", ylab = "Genetic position (cM)")
lines(geno$pos[on1] / 1e6, synthetic_map(geno$chr, geno$pos, suppression = 0)[on1],
lwd = 2, lty = 2, col = "grey50")
legend("topleft", c("centromere-suppressed", "uniform"),
lwd = 2, lty = c(1, 2), col = c("black", "grey50"), bty = "n")
Local rate, in cM per Mb, makes the same point more directly:
pos1 <- geno$pos[on1] / 1e6
cm1 <- geno$cm[on1]
bin <- cut(pos1, breaks = seq(0, 310, by = 20), labels = FALSE)
along <- data.frame(pos = pos1, cm = cm1, bin = bin)
span <- aggregate(cbind(cm, pos) ~ bin, data = along,
FUN = function(x) diff(range(x)))
span$cM_per_Mb <- round(span$cm / span$pos, 2)
span[, c("bin", "cM_per_Mb")]
# bin cM_per_Mb
# 1 1 1.06
# 2 2 1.02
# 3 3 0.94
# 4 4 0.79
# 5 5 0.56
# 6 6 0.33
# 7 7 0.19
# 8 8 0.20
# 9 9 0.36
# 10 10 0.60
# 11 11 0.81
# 12 12 0.95
# 13 13 1.03
# 14 14 1.06
# 15 15 1.07Rates drop to roughly 0.2 cM/Mb in the middle against about 1 cM/Mb on the arms.
Tuning it
| Argument | Effect |
|---|---|
cm_per_mb |
overall scale; default 0.73 gives ~1,500 cM for this genome |
total_cm |
set chromosome lengths directly instead of by rate |
suppression |
0 = uniform map, closer to 1 = stronger centromeric suppression |
width |
how much of the chromosome the suppressed region covers |
centromere |
centromere positions in bp; defaults to each chromosome’s midpoint |
# A uniform map: cM strictly proportional to bp
flat <- synthetic_map(geno$chr, geno$pos, suppression = 0)
# Fixed lengths, one per chromosome
fixed <- synthetic_map(geno$chr, geno$pos, total_cm = 150)
aggregate(cm ~ chr, data = data.frame(chr = geno$chr, cm = fixed), FUN = max)
# chr cm
# 1 1 150
# 2 2 150
# 3 3 150
# 4 4 150
# 5 5 150
# 6 6 150
# 7 7 150
# 8 8 150
# 9 9 150
# 10 10 150
# Real centromere positions, if you know them
limits <- aggregate(pos ~ chr, data = geno,
FUN = function(p) min(p) + 0.45 * diff(range(p)))
tuned <- synthetic_map(geno$chr, geno$pos, centromere = limits$pos)Using a real map
If you have measured genetic positions, use them — just put them in
cm, in centiMorgans, sorted within each chromosome:
The map must satisfy three conditions, which
as_population() checks:
- no missing values,
- non-decreasing within each chromosome,
- markers ordered by position.
Co-located markers may share a cM value; that correctly means no recombination between them.
How the map changes results
The map is not decoration — it determines how much recombination
happens. The expected number of crossovers on a chromosome equals its
length in Morgans (under the default gamma interference model and under
interference = "poisson" alike), so a 150 cM chromosome
averages 1.5 crossovers per meiosis.
Compare linkage under the default map with linkage under a map ten times shorter, where crossovers are rare and markers stay linked much further:
long_map <- geno
short_map <- geno
short_map$cm <- geno$cm / 10
linkage <- function(g, seed) {
pop <- as_population(g, individuals = c("33-16", "38-11"))
f1 <- cross(pop[1], pop[2], n = 1, seed = seed)
dh <- double_haploid(f1, n = 100, seed = seed + 1)
on1 <- g$chr == 1
p1 <- dosages(pop[1])[on1, 1]; p2 <- dosages(pop[2])[on1, 1]
inf <- which(p1 != p2 & p1 != 0 & p2 != 0)
from_p1 <- dosages(dh)[on1, ][inf, ] == p1[inf]
# recombination between the first informative locus and one ~100 loci away
mean(from_p1[100, ] != from_p1[1, ])
}
c(default = linkage(long_map, 1), ten_times_shorter = linkage(short_map, 1))
# default ten_times_shorter
# 0.18 0.04Shortening the map keeps the same markers linked far more tightly. If your simulated population shows implausibly long haplotype blocks, or breaks up far faster than your species should, the map is the first thing to check.
Session information
sessionInfo()
# R version 4.6.1 (2026-06-24)
# Platform: x86_64-pc-linux-gnu
# Running under: Ubuntu 26.04.1 LTS
#
# Matrix products: default
# BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
# LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
#
# locale:
# [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
# [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
# [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
# [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
# [9] LC_ADDRESS=C LC_TELEPHONE=C
# [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#
# time zone: Etc/UTC
# tzcode source: system (glibc)
#
# attached base packages:
# [1] stats graphics grDevices utils datasets methods base
#
# other attached packages:
# [1] simplePHENOTYPES_2.0.0.9003 rmarkdown_2.32
#
# loaded via a namespace (and not attached):
# [1] digest_0.6.39 R6_2.6.1 fastmap_1.2.0
# [4] xfun_0.61 maketools_1.3.2 cachem_1.1.0
# [7] knitr_1.52 RhpcBLASctl_0.23-42 htmltools_0.5.9
# [10] buildtools_1.0.0 lifecycle_1.0.5 cli_3.6.6
# [13] gdsfmt_1.49.8 SNPRelate_1.47.0 sass_0.4.10
# [16] data.table_1.18.6.1 jquerylib_0.1.4 compiler_4.6.1
# [19] sys_3.4.3 tools_4.6.1 evaluate_1.0.5
# [22] bslib_0.12.0 yaml_2.3.12 otel_0.2.0
# [25] jsonlite_2.0.0 rlang_1.3.0