Demographic models from the joint SFS in Goodman’s mouse lemur

0. Technical details

• You need two R packages, ggplot2 and patchwork. If you do not have them yet: install.packages(c("ggplot2", "patchwork")).
• Navigate to an easy-to-find directory in R with setwd("path/to/your/directory"), the same way as last time.
• Download the data from drive and download the data/ folder inside your local directory.
• fastsimcoal2 itself is not run in R. It is a command-line program for Linux, Windows, and older (Intel) Macs. We will run it on the cluster, but you might want to visualize some data in R on your own laptop.
• The data and an example SLURM submission script are on the cluster at /home/fmt019/work/mwcgls/sfs. I suggest you copy that whole directory to where you are working on the cluster. The files are small so it will be reasonable in this case.
• Skip directly to section 7 to submit a job to the SLURM scheduler, then come back and continue from the beginning

1. The dataset and the question

Microcebus lehilahytsara is a small, forest-dwelling primate found in the rainforests of eastern Madagascar and in isolated forest patches of the Central Highlands. Because it needs forest, its population history should serve as a record of forest history: when populations split, whether they stayed connected, and whether they shrank. Tiley, van Elst et al. (2022) used RADseq data to ask whether that history is dominated by ancient climate change (for example the Last Glacial Maximum, LGM, about 19–26.5 thousand years ago) or by human arrival (about 2 thousand years ago).

We will work with one of the paper’s three model sets: north vs. south. Ambavala is a rainforest population north of the Mangoro rivier. Ambatovy is rainforest population south of the Mangoro. Tsinjoarivo lies along the eastern escarpment of the Central Highlands. The phylogenetic analyses gave us an expected relationship of:

\[((\text{Tsinjoarivo}, \text{Ambatovy}), \text{Ambavala})\]

What a phylogeny cannot tell us alone is if populations exchanged migrants after they split, when that gene flow stopped, and whether population sizes changed recently. Those are demographic questions, and we answer them by comparing the observed joint site frequency spectrum with what different demographic models predict. Fastsimcoal2 (fsc) accomplishes this by simulating sequence evolution under the coalescent, for a user-defined model, and compares those expectations to the observation to find the least worst fit. Because we are depending on simulation, this means we have to perform many replicates to feel confident about our model selection and parameter estimates.

Two caveats about this data:

• The spectra were estimated from genotype likelihoods in ANGSD rather than from called genotypes. That is useful for low-coverage RADseq, and it is why the counts you will see are not whole numbers: they are expected numbers of sites. There are other approaches to genotyping, as you have covered, and tools to create sfs from genotype calls, but we start from the published data to save time.
• The samples are small: 5 Tsinjoarivo, 9 Ambatovy and 7 Ambavala lemurs. Recall from lecture that \(\mathbb E[T_{MRCA}] = 2(1 - 1/n)\): ten gene copies already capture 90% of the depth of the tree. Small samples are less of a problem for demographic inference than you might expect, at least for older events. Recent events are a different story: they leave their mark in the rarest variants, and small samples see few of those — one reason the timing of the recent size change in section 9 is so uncertain.

2. Reading a joint SFS

Every model folder in data/northVSsouth/ (models 0 to 10) contains the same three data files plus the two model files we use in meet in section 6. The data files *.obs are included in every directory because of a quirk in how fsc expects file names and looks for the input files. Here is an example of the inputs:

data/northVSsouth/10/
  10.tpl                    the model (template)
  10.est                    the parameters to estimate and their search ranges
  10_jointMAFpop1_0.obs     joint SFS: population 1 (rows) x population 0 (columns)
  10_jointMAFpop2_0.obs     joint SFS: population 2 (rows) x population 0 (columns)
  10_jointMAFpop2_1.obs     joint SFS: population 2 (rows) x population 1 (columns)

fastsimcoal2 is strict about file names. MAF means minor allele frequency, so these are folded spectra, and the population numbers refer to the order the populations are listed in the .tpl file: 0 = Tsinjoarivo, 1 = Ambatovy, 2 = Ambavala.

library(ggplot2)
library(patchwork)

# The file path is set up a little unclear here to try and make the script universal for everyone. Feel free to edit as need be to make the file inputs more clear.
DATA <- file.path("data", "northVSsouth")

# number of gene copies sampled per population (diploid individuals x 2),
# in the same order as the .tpl file
n_genes <- c(Tsinjoarivo = 10, Ambatovy = 18, Ambavala = 14)

Before we write any code, look at the top-left corner of one file in a text editor:

1 observations
    d0_0    d0_1    d0_2    d0_3    ...
d1_0    29026526.216652 43663.425033    23733.469008    ...
d1_1    148701.589575   20632.093979    11685.910229    ...

The first line says how many spectra are in the file. The second is a header of column labels. Every row after that starts with a label like d1_0 (“deme 1, count 0”) followed by numbers. So the cell in row d1_i and column d0_j is the number of sites at which the minor allele is carried by \(i\) gene copies in Ambatovy and \(j\) gene copies in Tsinjoarivo.

We could use read.table(..., header = TRUE), but we can also read the rows ourselves and then check the dimensions against the sample sizes. A header is only a set of labels, and fastsimcoal2 ignores it.

# Read one fastsimcoal2 joint SFS file into a matrix.
read_obs <- function(file) {
  lines <- readLines(file)

  # Keep only the data rows. They are the lines that start with "d" and a number
  # (d1_0, d1_1, ...). The "1 observations" line and the header are dropped.
  lines <- lines[grepl("^d[0-9]", lines)]

  rows <- list()
  for (k in seq_along(lines)) {
    pieces <- strsplit(lines[k], "[[:space:]]+")[[1]]   # split on tabs/spaces
    rows[[k]] <- as.numeric(pieces[-1])                  # drop the label, keep the counts
  }
  m <- do.call(rbind, rows)        # stack the rows into one matrix

  # Name rows and columns by the allele count they stand for: 0, 1, 2, ...
  rownames(m) <- 0:(nrow(m) - 1)
  colnames(m) <- 0:(ncol(m) - 1)
  m
}

ambat_tsin  <- read_obs(file.path(DATA, "10", "10_jointMAFpop1_0.obs"))   # rows Ambatovy,  cols Tsinjoarivo
ambav_tsin  <- read_obs(file.path(DATA, "10", "10_jointMAFpop2_0.obs"))   # rows Ambavala,  cols Tsinjoarivo
ambav_ambat <- read_obs(file.path(DATA, "10", "10_jointMAFpop2_1.obs"))   # rows Ambavala,  cols Ambatovy

Q1: Before you run the next chunk, predict the dimensions of each matrix from the sample sizes. Why is each dimension one more than the number of gene copies?

dim(ambat_tsin)
dim(ambav_tsin)
dim(ambav_ambat)

# If any of these fail, the rows and columns do not mean what we think they mean
stopifnot(nrow(ambat_tsin)  == n_genes["Ambatovy"] + 1, ncol(ambat_tsin)  == n_genes["Tsinjoarivo"] + 1)
stopifnot(nrow(ambav_tsin)  == n_genes["Ambavala"] + 1, ncol(ambav_tsin)  == n_genes["Tsinjoarivo"] + 1)
stopifnot(nrow(ambav_ambat) == n_genes["Ambavala"] + 1, ncol(ambav_ambat) == n_genes["Ambatovy"] + 1)

Now the most important cell in the whole file: the top-left one, [0, 0]. These are sites where neither population carries the minor allele — monomorphic sites.

total_sites <- sum(ambat_tsin)
monomorphic <- ambat_tsin["0", "0"]
snps        <- total_sites - monomorphic

cat("total sites:      ", round(total_sites), "\n")
cat("monomorphic sites:", round(monomorphic), "\n")
cat("SNPs:             ", round(snps), sprintf("(%.1f%% of sites)", 100 * snps / total_sites), "\n")

You should get 29,661,758 sites, of which 635,232 are SNPs (2.1%). Why keep 29 million boring sites? Because of \(\theta = 4N_e\mu\). The shape of the spectrum — how SNPs are spread across the cells — fixes the relative sizes and times of a model. The number of SNPs compared with the number of sites that could have mutated and did not gives \(\theta\) per site. And with a known mutation rate (\(\mu = 1.52\times10^{-8}\) per site per generation, we can rescale \(\theta\) into an estimated number of individuals and estimated divergence times in years; section 10 does that conversion. Without monomorphic sites, we could still select a model (such as post-divergence gene flow), but not interpret the divergence times and population sizes.

3. Looking at the data

A joint SFS is a two-dimensional histogram, so we draw it as a heatmap. Two functions that you can reuse in R for reading and plotting joint SFS are provided here. One turns the matrix into a table that ggplot can use, one draws it.

# Turn a joint SFS matrix into a data frame with one row per cell.
jsfs_to_df <- function(m) {
  df <- expand.grid(row_count = 0:(nrow(m) - 1), col_count = 0:(ncol(m) - 1))
  # as.vector() reads a matrix column by column -- the same order expand.grid() used
  df$sites <- as.vector(m)
  df
}

# Axis ticks at whole allele counts: every 2 copies, or every 5 for larger samples
count_breaks <- function(lims) {
  top <- floor(lims[2])
  if (top > 12) seq(0, top, by = 5) else seq(0, top, by = 2)
}

# Draw a joint SFS as a heatmap of log10(number of sites).
# The log is used, otherwise variation in many cells would be diffcult to see because a few bins can be very large
plot_jsfs <- function(df, limits = NULL) {
  df$sites[df$row_count == 0 & df$col_count == 0] <- NA   # hide the monomorphic cell: it would swamp the colours
  df$sites[df$sites < 1] <- NA                            # cells with (almost) no sites are left white
  ggplot(df, aes(x = col_count, y = row_count, fill = log10(sites))) +
    geom_tile() +
    scale_fill_viridis_c(na.value = "white", limits = limits, name = "log10\nsites") +
    scale_x_continuous(breaks = count_breaks) +
    scale_y_continuous(breaks = count_breaks) +
    theme_bw(base_size = 11) +
    theme(panel.grid = element_blank())
}

plot_jsfs(jsfs_to_df(ambat_tsin)) +
  coord_fixed() +
  labs(x = "minor-allele copies in Tsinjoarivo (of 10)",
       y = "minor-allele copies in Ambatovy (of 18)",
       title = "Ambatovy x Tsinjoarivo")

The observed joint SFS for Ambatovy and Tsinjoarivo. Each cell is the number of sites with \(i\) minor-allele copies in Ambatovy and \(j\) in Tsinjoarivo, on a log scale. The monomorphic cell at [0, 0] is masked to not skew the color scale.

Three things to notice:

1. The bottom row and the left column are the brightest. These are private variants, segregating in one population but absent from the other sample. Most of them are rare (low counts), just as \(\mathbb E[\xi_i] = \theta/i\) predicts for each population on its own.
2. The interior is far from empty. Hundreds of thousands of sites carry the minor allele in both populations. These shared variants were: 1) already segregating in the common ancestor or 2) carried across by gene flow. Among sites with the same total number of copies, the most common cells are those where the allele has a similar frequency in both populations. Note that this is frequency, not count: 9 copies out of 18 matches 5 out of 10.
3. The upper-right triangle is empty. That is because we are using the folded SFS. Even if you have a decent outgroup, determining the ancestral allele is very difficult in practice. So, best to use the folded (minor) SFS from the start. ANGSD determines the minor allele, by using all 28 gene copies together for the two populations. Other programs might determine the minor allele by using all populations, not just the pair being folded. Thus, you might need to adjust the fsc commands from this exercise for your own work depending on how the SFS are estimated and folded.

Now all three pairs side by side, on a shared colour scale:

shared_limits <- c(0, 5.5)     # log10 sites, the same for every panel

p1 <- plot_jsfs(jsfs_to_df(ambat_tsin), shared_limits) +
  labs(x = "Tsinjoarivo", y = "Ambatovy", title = "Ambatovy x Tsinjoarivo")
p2 <- plot_jsfs(jsfs_to_df(ambav_tsin), shared_limits) +
  labs(x = "Tsinjoarivo", y = "Ambavala", title = "Ambavala x Tsinjoarivo")
p3 <- plot_jsfs(jsfs_to_df(ambav_ambat), shared_limits) +
  labs(x = "Ambatovy", y = "Ambavala", title = "Ambavala x Ambatovy")

p1 + p2 + p3 + plot_layout(guides = "collect")

The three observed pairwise joint spectra. fastsimcoal2 fits all three at once.

The 1D SFS is hiding inside the 2D one

Add up a joint SFS across one population and you get a one-dimensional spectrum of the other — its margin. For an unfolded joint SFS the margin is that population’s ordinary SFS. Ours was folded over both populations together, so the margin counts the allele that is minor overall, which can be the majority allele within one population. Folding the margin once more, within the population, fixes that. Then we can compute everything from lecture: \(\pi\), \(\theta_W\) and Tajima’s \(D\).

# Fold a 1D SFS. x[1] is the count of sites with 0 copies, x[2] with 1 copy, ..., x[n+1] with n copies.
# Returns one entry per minor-allele count 1 .. n-1 (entries above n/2 stay 0), which is the
# input TajimaD() expects. Because pi weights count i by i(n - i), and i(n - i) is the same for
# i and n - i, pi and theta_W come out exactly right from a folded spectrum.
fold_1d <- function(x) {
  n <- length(x) - 1
  folded <- numeric(n - 1)
  for (i in 1:(n - 1)) {
    minor <- min(i, n - i)
    folded[minor] <- folded[minor] + x[i + 1]
  }
  folded
}

# Tajima's D from an SFS -- the same function as in lecture
TajimaD <- function(sfs) {              # sfs = c(#singletons, #doubletons, ...)
  n  <- length(sfs) + 1                 # number of sampled chromosomes
  ss <- sum(sfs)                        # segregating sites
  a1 <- sum(1 / seq_len(n - 1)); a2 <- sum(1 / seq_len(n - 1)^2)
  b1 <- (n + 1) / (3 * (n - 1)); b2 <- 2 * (n^2 + n + 3) / (9 * n * (n - 1))
  c1 <- b1 - 1/a1; c2 <- b2 - (n + 2)/(a1 * n) + a2/a1^2
  e1 <- c1 / a1;  e2 <- c2 / (a1^2 + a2)
  theta_pi <- sum(sapply(seq_along(sfs), function(i) i*(n-i)*sfs[i])) / choose(n, 2)
  theta_w  <- ss / a1
  (theta_pi - theta_w) / sqrt(e1*ss + e2*ss*(ss-1))
}

# pi, theta_W (per site) and Tajima's D for one population, from a margin
summarise_population <- function(margin) {
  n   <- length(margin) - 1
  sfs <- fold_1d(margin)
  i   <- 1:(n - 1)
  data.frame(n_genes    = n,
             SNPs       = round(sum(sfs)),
             pi         = sum(i * (n - i) * sfs) / choose(n, 2) / sum(margin),
             theta_W    = sum(sfs) / sum(1 / i) / sum(margin),
             tajimas_D  = TajimaD(sfs))
}

one_d <- rbind(
  cbind(population = "Tsinjoarivo", summarise_population(colSums(ambat_tsin))),   # add up over Ambatovy
  cbind(population = "Ambatovy",    summarise_population(rowSums(ambat_tsin))),   # add up over Tsinjoarivo
  cbind(population = "Ambavala",    summarise_population(rowSums(ambav_ambat)))   # add up over Ambatovy
)
one_d$pi      <- signif(one_d$pi, 3)
one_d$theta_W <- signif(one_d$theta_W, 3)
one_d$tajimas_D <- round(one_d$tajimas_D, 2)
one_d

You should get \(\pi \approx 0.004\)–\(0.005\) per site in all three populations, and Tajima’s \(D\) of +0.11 for Tsinjoarivo, −0.43 for Ambatovy and −0.20 for Ambavala. (Ambavala’s value shifts a little, to −0.17, if you take its margin from the Ambavala × Tsinjoarivo spectrum instead. Each pairwise spectrum was estimated from a slightly different set of sites.)

Do not hold these numbers up against the ±2 rule from lecture. Tajima’s variance assumes that all the SNPs share one genealogy, with no recombination between them. Across many SNPs we do not rely on hard cutoffs, but aim to interpret the relative differeneces between and within populations.

Q2: Rank the three populations by Tajima’s \(D\). Using the “Interpreting the sign” table from lecture, which population looks most like it is inside a contraction right now? The best model (section 9) will say that both Tsinjoarivo and Ambavala crashed recently. Why can Tajima’s \(D\) alone not tell you that Ambavala crashed? What does a joint SFS have that a single \(D\) value does not?

4. What shapes a joint SFS?

Before fitting models to the lemur data, it helps to see what different histories do to a joint SFS. Below are four simulated histories for two populations, A and B, with the coalescent simulator msprime. Each population has 10 diploids (20 gene copies) and \(N_e = 10{,}000\):

scenario	what happened (forward in time)
recent split	A and B split 2,000 generations ago (0.1 coalescent units)
old split	A and B split 40,000 generations ago (2 coalescent units)
old split + gene flow	as above, but A and B kept exchanging migrants (\(m = 5\times10^{-5}\) per generation, so \(4N_em = 2\))
old split + crash in B	as the old split, but B crashed to \(N_e = 500\) 400 generations ago and is still small

The simulation code can be provided for those interested, but to keep the time, we can have a look at the results for now:

Joint spectra simulated under four histories, unfolded (top) and folded (bottom).

Read the top row first, because unfolded spectra are easier to interpret:

• Recent split: shared variants sit close to the diagonal. There has been little time for drift, so an allele at frequency 30% in A is at about 30% in B. (Rare private variants still line the axes, as in every panel.)
• Old split: the diagonal empties and mass moves to the axes (private variants) and into the corners — top-left and bottom-right are sites fixed in one population and absent from the other. These are fixed differences.
• Gene flow: migrants carry alleles across, so the area between the axes and the diagonal fills back in, even though the split is just as old.
• Crash in B: drift in a small population quickly loses rare alleles, so B’s axis (the left column) thins out, and what remains is pushed toward intermediate and high frequencies in B. This is the two-dimensional version of “a contraction you are still inside” from lecture.

Folding (bottom row) keeps most of this information. It mirrors the upper-right triangle onto the lower-left: we can no longer tell “derived allele common in both populations” from “derived allele rare in both”, but the contrast between the axes and the interior survives, and so do the fixed differences (here they land exactly on the folding line).

Q3: The recent split and the old split with gene flow both create segregating sites shared between populations. Are there some differences between the spectra though that might be useful for developing your models for fastsimcoal2?

5. Three populations, three pairwise spectra

With three populations there are three pairs, and fastsimcoal2 fits all three pairwise spectra at once. To see why fitting all three spectra at once matters, we simulated the north-vs-south tree under three of the paper’s hypotheses. All three use the paper’s estimated split times and its population sizes from before the recent crash (section 9), with no crash at all, and they differ only in gene flow:

• model 0 — no gene flow at all
• model 1 — gene flow between Ambavala and the ancestor of the two southern populations, which stops when the southern populations split
• model 7 — gene flow between Ambatovy and Ambavala, from that split to the present

Model 1 uses the paper’s estimated ancestral migration rates. Model 7’s rate (\(5\times10^{-5}\) per generation each way) is simply chosen to make the effect visible.

Again, we will not simulate the spectra here, but we can make some visual comparisons

The observed spectra (top row) and the spectra predicted by three hypotheses of gene flow.

What fastsimcoal2 will do for us is take simulated SFS under models like the ones specified above and use its composite likelihood function to help estimate parameters and help us decide which is the most likely demographic scenario!

Why “composite” likelihood

For any model, fastsimcoal2 computes the probability of each cell of each spectrum and multiplies them together over cells and over the three pairs. That treats things as independent that are not. Every site is counted once in each of the three pairwise spectra, and sites within a RAD locus are linked. Multiplying probabilities anyway gives a composite (or pseudo-) likelihood. It is still a good tool for finding the best parameters, but it overstates how certain we are. Keep that in mind for section 8.

6. Writing the model for fastsimcoal2

fastsimcoal2 needs two files per model. The template (.tpl) describes the history, with parameter names where numbers would go. The estimation file (.est) says which of those parameters to estimate and where to start searching. Here is model 0, the simplest:

//Number of population samples (demes)
3
//Population effective sizes (number of genes)
Theta_Tsin
Theta_Ambat
Theta_Ambav
//Sample sizes
10
18
14
//Growth rates : negative growth implies population expansion
0
0
0
//Number of migration matrices : 0 implies no migration between demes
0
//historical event: time, source, sink, migrants, new size, new growth rate, migr. matrix
2 historical event
Tau_NodeS 1 0 1 Theta_NodeS 0 0 absoluteResize
Tau_NodeR 2 0 1 Theta_NodeR 0 0 absoluteResize
//Number of independent loci [chromosome]
1 0
//Per chromosome: Number of linkage blocks
1
//per Block: data type, num loci, rec. rate and mut rate + optional parameters
FREQ 1 0 1.52e-8 OUTEXP

Line by line:

• 3 demes, listed in a fixed order: 0 = Tsinjoarivo, 1 = Ambatovy, 2 = Ambavala. This order is why the data file for Ambatovy × Tsinjoarivo is called _jointMAFpop1_0.obs.
• Population sizes are names (Theta_Tsin, …) that fsc fills in. Despite the name, they are not \(\theta\); they are sizes in gene copies (haploid), so twice the number of diploid individuals.
• Sample sizes are also in gene copies: 10, 18 and 14.
• Historical events are the heart of the model, and they run backward in time, just like the coalescent. Tau_NodeS 1 0 1 Theta_NodeS 0 0 absoluteResize reads: at time Tau_NodeS generations ago, move a proportion 1 (all) of the lineages from deme 1 (Ambatovy) into deme 0 (Tsinjoarivo); set deme 0’s size to Theta_NodeS; set its growth rate to 0; use migration matrix 0 from now on. Without the absoluteResize keyword, the fifth number would be a ratio by which the current size is multiplied. Backward in time, a split is a merge. Deme 0 now stands for the ancestor of the two southern populations (“node S”). The second event merges Ambavala into that ancestor at the root (“node R”).
• FREQ 1 0 1.52e-8 OUTEXP: the data are site frequencies, with mutation rate \(1.52\times10^{-8}\); OUTEXP asks fsc to write out the expected spectra.

And the matching .est file:

[PARAMETERS]
//#isInt? #name #dist. #min #max
//all Ns are in number of haploid individuals (gene copies)
1 Theta_Tsin unif 10 100000 output
1 Theta_Ambat unif 10 100000 output
1 Theta_Ambav unif 10 100000 output
1 Theta_NodeS unif 10 100000 output
1 Theta_NodeR unif 10 100000 output
1 Tau_NodeS unif 10 50000 output
1 Tau_NodeR logunif Tau_NodeS 200000 output paramInRange  // root is never younger than node S
[COMPLEX PARAMETERS]

• The seven lines under [PARAMETERS] are the free parameters: five sizes and two times. The first column says whether the value is a whole number; then come a distribution (unif or logunif) and a range. The ranges are where each search starts: fsc can wander above an ordinary upper limit unless a parameter is marked bounded. Several of the paper’s size estimates are well over 100,000.
• Tau_NodeR has a range that refers to another parameter: from Tau_NodeS up to 200,000 generations, flagged with paramInRange. That is how the file guarantees the root is never younger than node S. Such ranges are hard limits, so the upper one is set generously. After a fit, it is worth checking that no estimate is sitting right at a limit.
• [COMPLEX PARAMETERS] are calculated from the free ones. Model 0 has none, but the section header must still be there; fsc28 crashes without it.
• Models with gene flow add free parameters such as NM_12: a scaled migration parameter (a population size times a migration rate). It is converted to the per-generation rate that fsc actually uses with MIG_12 = NM_12/Theta_Ambav. This is the same kind of reparameterization as in the fsc manual’s own examples, and it mostly makes sensible search ranges easier to write. In the .tpl, those rates go into migration matrices, where the entry in row \(i\), column \(j\) is the probability that a lineage moves from deme \(i\) to deme \(j\) per generation, again backward in time.

NoteSoftware versions change input files

The paper used fastsimcoal 2.6. Its files had no absoluteResize or paramInRange, which arrived in version 2.7. Sizes were set as ratios computed in [COMPLEX PARAMETERS] (Resize_NodeS = Theta_NodeS/Theta_Tsin), the root time was built as Tau_NodeR = Tau_NodeS + BL_S, and the files carried a [RULES] section that version 2.7 dropped. The files in data/ use the newer syntax. They describe the same eleven models with the same number of free parameters, and fsc28 simulates identical spectra from the old and new versions given the same parameter values (we checked). But fsc has improved its optimization algorithm since, so your best likelihoods will not match Table S12. Whenever you reuse someone’s input files, check which version of the program they were written for.

The 11 hypotheses

The paper built its models step by step: first test which pairs exchanged genes, then refine when, then add a recent size change. Every folder holds one model. The number of free parameters, \(k\), is the number of lines under [PARAMETERS]. Have a look at some of the input *.est files and see if you can count the parameters.

model	ancestor S ↔︎ Ambavala	Ambatovy ↔︎ Ambavala	Tsinjoarivo ↔︎ Ambatovy	recent size change	\(k\)
0					7
1	✓				9
2	✓	✓			11
3	✓	✓	✓		13
4	✓		✓		11
5			✓		9
6		✓	✓		11
7		✓			9
8	✓	✓ stops before the present			12
9	✓	✓ starts again recently			12
10	✓	✓ stops before the present		✓ all three	16

You can visualize the 11 possible models below as well:

The eleven north-vs-south models. Blue arrows are gene flow; a dashed line marks when gene flow stops (models 8 and 10) or resumes (model 9); the red line in model 10 is the recent size change in all three populations.

Note that models 8 and 9 have the same parameters and differ only in which migration matrix comes first in the .tpl. In model 8, matrix 0 (the present) has no migration and matrix 1 (older) has Ambatovy ↔︎ Ambavala migration; in model 9, it is the other way round. The switch happens at Tau_Change, whose range in the .est file runs from 0 to Tau_NodeS (paramInRange), so the switch is always more recent than node S. Model 10 does the same twice: Tau_Crash lies between 0 and Tau_NodeS, and Tau_Change between Tau_Crash and Tau_NodeS.

Q5: In words, what biological story does model 8 tell, and what story does model 9 tell? Why is the range of Tau_Change written as 0 to Tau_NodeS rather than as two fixed numbers? What could go wrong with a fixed range?

7. Running fastsimcoal2 (on Linux)

Read this section first and submit an analysis to the scheduler before beginning exercises from the beginning. This will allow some time for analyses to finish.

Each model is fit by maximum composite likelihood. fsc starts from random parameter values, simulates the expected spectrum, and climbs towards better parameters through a number of optimisation cycles. Because that climb can get stuck on a local peak, every model is run many times independently, and only the best run is kept.

A submission script is provided to help you run multiple replicates. It can be found at /home/fmt019/work/mwcgls/sfs/results/fsc-slurm.sh or where ever you might have copied the fsc directory on the cluster. In order to get all of the models done on time, we will have to work togehter! Go around the room and pick your model number in order. We need to make sure all models are run at least by one person. There are two things you will need to change: 1. The DATA_MODEL line. Change it to match your model. 2. Your DATA_DIR and RESULT_DIR paths. Probably safest to have it run in your own directly so you do not accidentally overwrite a colleague’s result.

#!/bin/bash
#SBATCH -p workq
#SBATCH --time=02:00:00 
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=4
#SBATCH --mem=4G
#SBATCH --job-name=fsc_mleh
#SBATCH --output=fsc_mleh.out
#SBATCH --error=fsc_mleh.err

DATA_DIR=/home/fmt019/work/mwcgls/sfs/data/northVSsouth
DATA_MODEL="2"

RESULT_DIR=/home/fmt019/work/mwcgls/sfs/results/runs

#How many independent replicates to run
#Each analysis could get stuck in a local optimum, so multiple independent starts are needed
N_REPS=10

#Provide the path to the fsc28 program
FSC=/home/fmt019/work/mwcgls/sfs/programs/fsc28_linux64/fsc28

#Generate a random number to ensure a different initial state
SEED=$SRANDOM

#Run fsc
for i in $(seq 1 ${N_REPS}); do
    mkdir -p "${RESULT_DIR}/${DATA_MODEL}_$i"
    MY_SEED=$SRANDOM
    RUN_DIR="${RESULT_DIR}/${DATA_MODEL}_$i"
    cp "${DATA_DIR}/${DATA_MODEL}/${DATA_MODEL}.tpl" "${DATA_DIR}/${DATA_MODEL}/${DATA_MODEL}.est" "${DATA_DIR}/${DATA_MODEL}/${DATA_MODEL}"_jointMAFpop*.obs "$RUN_DIR/"
    cd ${RUN_DIR}
    echo "Running command: $FSC -t "${DATA_MODEL}.tpl" -e "${DATA_MODEL}.est" -m -M -n 10000 -L 40 -c 4 -B 4 -r $SEED --foldedSFS"
    $FSC -t "${DATA_MODEL}.tpl" -e "${DATA_MODEL}.est" -m -M -n 10000 -L 40 -c 4 -B 4 -r $SEED --foldedSFS
done

The main program is run with this line $FSC -t "${DATA_MODEL}.tpl" -e "${DATA_MODEL}.est" -m -M -n 10000 -L 40 -c 4 -B 4 -r $SEED --foldedSFS. Below is a description of the options. We are doing a low number of simulations for this example, but a real analysis would want to increase that.

option	meaning
-t 0.tpl	the template file; the .obs files are found from its name
-e 0.est	the estimation file
-m	use the minor allele (folded) spectra, *_jointMAF*.obs (-d would mean derived/unfolded)
-M	estimate parameters by maximum composite likelihood
-n 10000	number of coalescent simulations used to approximate the expected SFS at each step
-L 10	number of optimisation (ECM) cycles
-c 2 -B 2	threads and batches
-q	quiet
--foldedSFS	fold the model’s spectra the same way the observed ones were folded (below)

If you started here and successfully launched your job, go back to the beginning to work through the content.

When your analysis is done, inspect the output files. Can you find the best composite lnL for your model? Record it here along with your name: workshop_results

Fold the model the way the data were folded

A folded spectrum counts each site by its minor allele, but “minor” has to be decided over some set of gene copies. ANGSD folded each pairwise spectrum over that pair’s own copies: 28 for Ambatovy × Tsinjoarivo. That is why the upper-right triangle of every observed spectrum is empty (section 3). By default, fsc’s -m decides the minor allele over all the copies in the model, 42 here. Its expected spectra therefore put sites into cells that are always empty in our data. For these models that is about 2% of the expected SNPs.

--foldedSFS (fsc 2.7 and later) makes fsc fold each pair the way ANGSD did. In a quick test with the same model and random seed, it shrank the gap between model 10’s best log10-likelihood and a perfect fit (MaxObsLhood, section 8) from about 25,900 to about 7,500 units. The paper used fsc version 2.6, before this option existed.

Many runs per model

The paper used -n 200000 -L 40 and 50–100 independent runs per model. A full run of every model is a job for a cluster: loop over models and replicates, each in its own folder so the outputs do not overwrite each other. The fsc-slurm.sh script helps manage this a little. It does not automate selecting the best run for you though. Here, we will rely on everyone working together. On your own, you would need to write a script to loop through each .bestlhoods file. Those files have the estimated parameters, then two likelihoods. This is the example from the fsc manual, for a one-population model with three parameters:

NPOP    NANC  TEXP  MaxEstLhood  MaxObsLhood
502642  542   5070  -90349.486   -90348.375

MaxEstLhood is the best likelihood the model found. MaxObsLhood is the likelihood you would get if the model predicted the observed spectrum perfectly, so the gap between them measures how well (or badly) the model fits.

8. Choosing a model

Each model has a maximum likelihood \(L\) and a number of free parameters \(k\). If one model contains another as a special case, its best possible likelihood can only be higher, so we need a penalty for complexity. The Akaike information criterion is

\[\text{AIC} = 2k - 2\ln L\]

and the model with the smallest AIC is preferred. To say how much it is preferred, we turn the differences \(\Delta_i = \text{AIC}_i - \text{AIC}_{min}\) into Akaike weights, which can be read as model probabilities:

\[w_i = \frac{\exp(-\Delta_i/2)}{\sum_j \exp(-\Delta_j/2)}\]

Table S12 of the paper gives the best log-likelihood and \(k\) for each of the eleven models:

s12 <- read.csv(file.path("data", "tableS12_model_selection_north_vs_south.csv"))

# our parameter counts from section 6 should match the paper's
stopifnot(all(s12$n_params == k))

s12$AIC       <- 2 * s12$n_params - 2 * s12$lnL
s12$delta_AIC <- s12$AIC - min(s12$AIC)
s12$weight    <- exp(-0.5 * s12$delta_AIC) / sum(exp(-0.5 * s12$delta_AIC))

# do we reproduce the published AIC?
all.equal(s12$AIC, s12$AIC_published)

s12[order(s12$AIC), c("model", "lnL", "n_params", "AIC", "delta_AIC", "weight")]

You should reproduce the published AIC, and find that model 10 gets a weight of 1 — every other model rounds to 0. The best model has gene flow between Ambavala and the southern ancestor, gene flow with Ambatovy that stopped some time before the present, and a recent size change in all three populations. Notice also the clear tiers: every model with ancestral gene flow (1, 2, 3, 4, 8, 9) beats every model without it (0, 5, 6, 7) by thousands of units.

You will reproduce this table with your own results. After a best composite lnL has been reported for each model, see if you can fill in the “Model Selection” sheet of the linked google sheet on workshop_results.

Q7: Model 3 contains model 2 as a special case: set the Tsinjoarivo–Ambatovy migration to (almost) zero and you have model 2. Yet model 3’s best log-likelihood is lower (−6,971,967 against −6,971,944). How can a model that contains another one fit worse? Think about how fsc gets its expected spectrum (-n) and how it searches.

Doing the same with your own fsc runs

When you have run fsc yourself (section 7), this chunk collects every .bestlhoods file, keeps the best replicate of each model, and builds the same kind of table. There is one trap:

Importantfastsimcoal2 reports log10 likelihoods

MaxEstLhood and MaxObsLhood are base-10 logarithms. AIC is defined with the natural logarithm, so multiply by \(\ln(10) \approx 2.303\) first.

Q9: Model 10 beats model 8 by 1,794 AIC units, a relative likelihood of about \(e^{-897}\) for model 8. Section 5 explained that this is a composite likelihood over ~30 million sites, many of them linked within RAD loci. Do you believe the evidence is really that strong? The paper estimated its confidence intervals with a block bootstrap: it resampled blocks of 10,000 sites rather than single sites. Why blocks?

9. Does the best model hold up? (optional)

Skip this section if you are short on time. You lose two checks that are separate questions from whether model 10 beat the other ten: does it actually reproduce the data, and do its estimates hold up when it is fitted again?

AIC tells you which model is least worst among the ones you tried. It does not tell you whether that model is any good. The way to find out is to simulate data from the fitted model and compare it with the real data.

Table S15 gives the maximum likelihood estimates for model 10. They are in fastsimcoal2’s own units: times in generations and sizes in gene copies (section 10 explains these units). With a generation time of 3.5 years, we can put them on a calendar:

s15 <- read.csv(file.path("data", "tableS15_parameters_north_vs_south.csv"))

GENERATION_TIME <- 3.5    # years

times <- s15[grepl("^tau", s15$parameter), ]
times$MLE_years     <- times$MLE * GENERATION_TIME
times$lower95_years <- times$lower95 * GENERATION_TIME
times$upper95_years <- times$upper95 * GENERATION_TIME
times

sizes <- s15[grepl("^theta", s15$parameter), ]
sizes$diploid_Ne <- sizes$MLE / 2        # gene copies -> diploid individuals
sizes[, c("parameter", "MLE", "diploid_Ne", "lower95", "upper95")]

# migration rates per generation: "a" = Ambatovy <-> Ambavala after node S,
# "b" = the southern ancestor S <-> Ambavala before node S
s15[grepl("^m_", s15$parameter), ]

Some numbers to think about:

• Root split (Ambavala vs. the south): about 79,000 generations, or ~277,000 years — long before the LGM.
• Tsinjoarivo–Ambatovy split (node S): about 2,200 generations, ~7,800 years, after the LGM.
• End of gene flow between Ambatovy and Ambavala: almost immediately after node S. And look at the rates for that short window: the estimates of m_a12 and m_a21 are about \(2\times10^{-9}\), effectively zero, although their confidence intervals reach up to \(10^{-5}\)–\(10^{-4}\). The gene flow the model is sure about is the older one, between Ambavala and the southern ancestor (m_b02, m_b20 ≈ 4–6 × 10⁻⁵). At the estimates, north–south gene flow ended about when Tsinjoarivo and Ambatovy split.
• Size change: 534 generations, ~1,900 years ago. Forward in time, Tsinjoarivo fell from 75,215 to 3,750 gene copies (20-fold) and Ambavala from 329,425 to 7,340 (45-fold), while Ambatovy stayed about the same.

Model 10 was simulated with these values with msprime. Here are its spectra next to the observed ones:

Observed spectra (top), spectra simulated under model 10 at the published estimates (middle), and their ratio (bottom). White means a perfect fit; red cells have more sites than the model expects, blue fewer.

You would find that model 10 reproduces the number of SNPs to within a few percent. The fit is good, not perfect, and the areas where the fit is poor is not random. Look at the bottom row: the clearest red cells sit along the Tsinjoarivo axis at high allele counts. The data have more sites where the minor allele is common in Tsinjoarivo but rare or absent in the other population than model 10 expects, and, more weakly, the same along the Ambavala axis. (Very sparse cells are noisy in a 10 Mb simulation, so do not read much into any single pale cell.)

Q10: What kind of history would produce more high-frequency private variants in Tsinjoarivo than model 10 predicts? Think about what drift in a small population does to variants that start out rare.

Q11: The size change is estimated at ~1,900 years ago. Its 95% confidence interval runs from 34 to 3,554 generations. Convert that to years. Does the interval exclude the LGM (19,000–26,500 years ago)? Does it exclude a decline before humans arrived (~2,000 years ago)? What would you write in a paper?

The same model, fitted again

Table S15 came from fastsimcoal 2.6 and the paper’s original files. We fitted model 10 again with fastsimcoal 2.8 and the files in data/. We used the paper’s settings (-n 200000 -L 40) plus --foldedSFS (section 7), with 10 independent replicates, which took about five hours on a laptop. Each replicate’s best estimates are in data/refit_model10_fsc28.csv, one row per replicate, with the same columns as a .bestlhoods file plus a rep column for the replicate number.

refit <- read.csv(file.path("data", "refit_model10_fsc28.csv"))
refit <- refit[order(refit$MaxEstLhood, decreasing = TRUE), ]   # best replicate first

# how far each replicate is behind the best one, and behind a perfect fit (log10 units)
refit$behind_best    <- round(refit$MaxEstLhood - refit$MaxEstLhood[1])
refit$gap_to_perfect <- round(refit$MaxEstLhood - refit$MaxObsLhood)
refit[, c("rep", "MaxEstLhood", "behind_best", "gap_to_perfect")]

You should see that replicate 8 found the highest likelihood. The other nine finished between 165 and 884 log10 units behind it, and the best one is 2,570 units short of a perfect fit. Only one of ten replicates reached the top, and that is exactly why the paper ran 50–100 per model.

The complete fastsimcoal2 output of that best replicate is in data/model10_best_fit/, under the names fsc gave it. Besides .bestlhoods, fsc writes its expected spectra at the estimates: same names as the .obs files, with a .txt extension. Here they are checked against the data, the same way we checked the msprime simulation above:

best_fit <- file.path("data", "model10_best_fit")     # or the folder of your own best replicate

fsc_check <- data.frame(pair = pair_names, FST_observed = NA, FST_fsc_expected = NA)
for (k in seq_along(observed)) {
  pair <- names(observed)[k]
  expected <- read_obs(file.path(best_fit, paste0("10_jointMAFpop", pair, ".txt")))

  # fsc rescales its expected spectrum so that the POLYMORPHIC cells sum to 1 (fsc manual),
  # so leave out the monomorphic cell and scale to the observed number of SNPs instead
  expected[1, 1] <- 0
  observed_snps  <- sum(observed[[pair]]) - observed[[pair]][1, 1]
  expected       <- expected * observed_snps / sum(expected)

  fsc_check$FST_observed[k]     <- round(hudson_fst(observed[[pair]]), 3)
  fsc_check$FST_fsc_expected[k] <- round(hudson_fst(expected), 3)
}
fsc_check

fastsimcoal2’s own expectation at the refit estimates reproduces the differentiation between all three pairs (\(F_{ST}\) of about 0.099, 0.143 and 0.083, against 0.088, 0.143 and 0.082 observed).

Now rescale every replicate to absolute time alongside Table S15:

GENERATION_TIME <- 3.5     # years

# the four event times: column in the refit file, and name in Table S15
time_params <- data.frame(quantity  = c("north-south split", "Tsinjoarivo-Ambatovy split",
                                        "end of gene flow", "size change"),
                          refit_col = c("Tau_NodeR", "Tau_NodeS", "Tau_Change", "Tau_Crash"),
                          s15_name  = c("tau_r", "tau_s", "tau_migration_end", "tau_crash"))

time_compare <- data.frame()
for (k in 1:nrow(time_params)) {
  years   <- refit[[time_params$refit_col[k]]] * GENERATION_TIME     # every replicate, in years
  s15_row <- s15[s15$parameter == time_params$s15_name[k], ]
  time_compare <- rbind(time_compare, data.frame(
    quantity       = time_params$quantity[k],
    refit_best     = round(years[1]),          # row 1 is the best replicate
    refit_lowest   = round(min(years)),
    refit_highest  = round(max(years)),
    paper_estimate = round(s15_row$MLE * GENERATION_TIME),
    paper_lower95  = round(s15_row$lower95 * GENERATION_TIME),
    paper_upper95  = round(s15_row$upper95 * GENERATION_TIME)))
}
time_compare

# size before / size after the recent size change: above 1 means a decline (forward in time)
fold <- data.frame(rep         = refit$rep,
                   Tsinjoarivo = refit$Theta_TsinChange  / refit$Theta_TsinCurr,
                   Ambatovy    = refit$Theta_AmbatChange / refit$Theta_AmbatCurr,
                   Ambavala    = refit$Theta_AmbavChange / refit$Theta_AmbavCurr)
s15_value  <- setNames(s15$MLE, s15$parameter)
paper_fold <- c(Tsinjoarivo = s15_value[["theta_Tsinjoarivo_before"]] / s15_value[["theta_Tsinjoarivo_after"]],
                Ambatovy    = s15_value[["theta_Ambatovy_before"]]    / s15_value[["theta_Ambatovy_after"]],
                Ambavala    = s15_value[["theta_Ambavala_before"]]    / s15_value[["theta_Ambavala_after"]])
round(fold, 1)
round(paper_fold, 1)

# long tables for the figure: one row per replicate and quantity
times_long <- data.frame()
for (k in 1:nrow(time_params)) {
  times_long <- rbind(times_long, data.frame(quantity = time_params$quantity[k],
                                             years    = refit[[time_params$refit_col[k]]] * GENERATION_TIME,
                                             best     = refit$rep == refit$rep[1]))
}
times_long$quantity   <- factor(times_long$quantity,   levels = time_params$quantity)
time_compare$quantity <- factor(time_compare$quantity, levels = time_params$quantity)

fold_long <- data.frame()
for (pop in c("Tsinjoarivo", "Ambatovy", "Ambavala")) {
  fold_long <- rbind(fold_long, data.frame(population = pop, fold = fold[[pop]],
                                           best = fold$rep == refit$rep[1]))
}
paper_fold_df <- data.frame(population = names(paper_fold), fold = paper_fold)
pop_order <- c("Tsinjoarivo", "Ambatovy", "Ambavala")       # same order as everywhere else
fold_long$population     <- factor(fold_long$population,     levels = pop_order)
paper_fold_df$population <- factor(paper_fold_df$population, levels = pop_order)

year_breaks <- c(100, 1000, 10000, 100000, 1000000)
year_labels <- c("100 yr", "1 ka", "10 ka", "100 ka", "1 Ma")

p_times <- ggplot() +
  geom_blank(data = time_compare, aes(x = quantity, y = paper_estimate))   +  # sets up the axes first
  # the two periods the paper asks about
  annotate("rect", xmin = -Inf, xmax = Inf, ymin = 19000, ymax = 26500, fill = "#4477AA", alpha = 0.15) +
  annotate("rect", xmin = -Inf, xmax = Inf, ymin = 1000, ymax = 2000, fill = "#DDAA33", alpha = 0.3) +
  annotate("text", x = 0.55, y = 22500, label = "LGM", hjust = 0, size = 3, colour = "grey30") +
  annotate("text", x = 0.55, y = 1400, label = "human arrival", hjust = 0, size = 3, colour = "grey30") +
  # Table S15: estimate (diamond) and 95% interval, drawn a little to the left
  geom_errorbar(data = time_compare, aes(x = quantity, ymin = paper_lower95, ymax = paper_upper95),
                width = 0.12, colour = "grey30", position = position_nudge(x = -0.22)) +
  geom_point(data = time_compare, aes(x = quantity, y = paper_estimate),
             shape = 23, size = 3.2, fill = "white", position = position_nudge(x = -0.22)) +
  # the refit: every replicate (grey) and the best one (red)
  geom_jitter(data = times_long[!times_long$best, ], aes(x = quantity, y = years),
              width = 0.08, height = 0, colour = "grey55", size = 1.8) +
  geom_point(data = times_long[times_long$best, ], aes(x = quantity, y = years), colour = "#CC3311", size = 3) +
  scale_y_log10(breaks = year_breaks, labels = year_labels) +
  labs(x = NULL, y = "years ago (log scale)") +
  theme_bw(base_size = 11) +
  theme(panel.grid.minor = element_blank())

p_fold <- ggplot() +
  geom_hline(yintercept = 1, linetype = "dashed", colour = "grey50") +
  geom_point(data = paper_fold_df, aes(x = population, y = fold),
             shape = 23, size = 3.2, fill = "white", position = position_nudge(x = -0.22)) +
  geom_jitter(data = fold_long[!fold_long$best, ], aes(x = population, y = fold),
              width = 0.08, height = 0, colour = "grey55", size = 1.8) +
  geom_point(data = fold_long[fold_long$best, ], aes(x = population, y = fold), colour = "#CC3311", size = 3) +
  scale_y_log10() +
  labs(x = NULL, y = "size before / size after (log scale)\nabove the dashed line = a decline") +
  theme_bw(base_size = 11) +
  theme(panel.grid.minor = element_blank())

p_times + p_fold + plot_layout(widths = c(2.2, 1))

Model 10 fitted again with fsc 2.8 (10 replicates) against the published estimates (Table S15). Red: the best replicate; grey: the other nine; white diamonds and bars: Table S15 estimates and 95% intervals.

Read the figure one quantity at a time:

• Robust: the Tsinjoarivo–Ambatovy split. Every replicate puts it at about 15,500–22,000 years ago. That is older than Table S15’s 7,800 years, though still inside its interval, and right at the end of the LGM.
• Robust: the ancestral gene flow. MIG_ANC_02 and MIG_ANC_20 come out at about 3–6 × 10⁻⁵ in every replicate, close to the published 5.5 and 3.7 × 10⁻⁵.
• Loosely constrained: the root. It lands anywhere between about 265,000 and 540,000 years ago.
• Fragile: the recent size change. The best replicate puts it at about 1,200 years ago, with 20- to 33-fold declines in Tsinjoarivo and Ambavala. The next four replicates, all within 300 log10 units, put it at 3,000–5,000 years ago with milder declines, and the weakest ones at 12,000–20,000 years with much smaller declines (at most about five-fold). Timing and size trade off against each other: a recent, sharp decline and an older, gentle one leave very similar spectra.
• Fragile: Ambatovy. It declines in the five best replicates (2- to 9-fold), but not consistently in the others, and not at all in Table S15.

Several things changed at once between the published fit and this one: the program version, the folding, the way the search is parameterized, and the number of replicates. Estimates that sit on a flat ridge of the likelihood, like the timing of the recent decline, move with any of them. Estimates on a sharp peak, like the southern split, stay put.

Q12: Which conclusions would you state confidently from this refit? The paper concluded that the recent decline probably happened after humans arrived. Does the refit support that, contradict it, or fail to settle it? What would you do next to find out?

10. From generations and gene copies to years and individuals

Every number fastsimcoal2 writes is in its own units. Sizes are in gene copies, so a population of \(N_e\) diploid individuals has \(N = 2N_e\) copies. Times are in generations. Migration rates (MIG_*) are the probability, per generation, that a lineage moves from one deme to another, backward in time. None of these is something you can put on a map or a calendar yet.

Why these units? The spectrum itself only pins down combinations of parameters with the mutation rate: \(\theta = 4N_e\mu = 2N\mu\), the products (time × μ), and the number of migrants \(N \times m\). Because the .tpl file fixes \(\mu = 1.52 \times 10^{-8}\) per site per generation, fsc can report gene copies and generations. A generation time (\(g = 3.5\) years for mouse lemurs) then turns generations into years. Each conversion is an assumption, and this section makes them explicit.

Reading model 10’s output

Each replicate writes its best estimates to <model>/<model>.bestlhoods: one header line of parameter names (the simple parameters, then the complex ones), then one line of values, ending with MaxEstLhood and MaxObsLhood. Our best fit is in data/model10_best_fit/, exactly as fsc wrote it, so you can follow along without running fsc:

est <- read.table(file.path("data", "model10_best_fit", "10.bestlhoods"), header = TRUE)
names(est)       # simple parameters, complex parameters, then the two log10 likelihoods

If you have run model 10 yourself, read every replicate and keep the best instead. This replaces est with your own fit, so your numbers below will differ somewhat from ours:

# every replicate of model 10 wrote runs/model10/rep<r>/10/10.bestlhoods
files <- Sys.glob(file.path("runs", "model10", "rep*", "10", "10.bestlhoods"))
fits  <- data.frame()
for (f in files) {
  fits <- rbind(fits, read.table(f, header = TRUE))   # one row per replicate
}
est <- fits[which.max(fits$MaxEstLhood), ]            # keep the best replicate

fsc also writes the model with these estimates substituted into the template, 10_maxL.par. Its historical events are the fitted history: times in generations and, with absoluteResize, sizes in gene copies. Those are exactly the numbers we are about to convert:

par_lines <- readLines(file.path("data", "model10_best_fit", "10_maxL.par"))
start <- grep("historical event", par_lines)[2]        # the line that says "6 historical event"
writeLines(par_lines[start:(start + 6)])

Rescaling

MU              <- 1.52e-8   # mutation rate per site per generation, as in the .tpl file
GENERATION_TIME <- 3.5       # years per generation

# times: fsc reports generations
times <- data.frame(event       = c("north-south split", "Tsinjoarivo-Ambatovy split",
                                    "end of gene flow", "size change"),
                    generations = c(est$Tau_NodeR, est$Tau_NodeS, est$Tau_Change, est$Tau_Crash))
times$years <- times$generations * GENERATION_TIME
times

# sizes: fsc reports gene copies; Ne diploid individuals carry 2 * Ne copies
sizes <- data.frame(population  = c("Tsinjoarivo now", "Ambatovy now", "Ambavala now",
                                    "Tsinjoarivo before the change", "Ambatovy before the change",
                                    "Ambavala before the change", "ancestor S", "root R"),
                    gene_copies = c(est$Theta_TsinCurr, est$Theta_AmbatCurr, est$Theta_AmbavCurr,
                                    est$Theta_TsinChange, est$Theta_AmbatChange, est$Theta_AmbavChange,
                                    est$Theta_NodeS, est$Theta_NodeR))
sizes$diploid_Ne     <- sizes$gene_copies / 2
sizes$theta_per_site <- signif(2 * sizes$gene_copies * MU, 3)   # = 4 * Ne * mu
sizes

# migration: MIG is the probability per generation that a lineage moves from the row deme to the
# column deme, backward in time. Forward in time, that is the fraction of the row population that
# arrives from the column population each generation, so the number of immigrant gene copies per
# generation is MIG times the size of the receiving (row) population.
migration <- data.frame(direction      = c("into ancestor S, from Ambavala", "into Ambavala, from ancestor S"),
                        MIG            = c(est$MIG_ANC_02, est$MIG_ANC_20),
                        recipient_size = c(est$Theta_NodeS, est$Theta_AmbavChange))
migration$immigrants_per_generation <- round(migration$MIG * migration$recipient_size, 2)
migration

# does the ancestral size predict the diversity we measured in section 3?
cat("theta per site in the ancestor S:", signif(2 * est$Theta_NodeS * MU, 3), "\n")
cat("observed pi per site (section 3):", one_d$pi, "\n")

You should get these dates for the best replicate: - the north–south split about 299,000 years ago; - the Tsinjoarivo–Ambatovy split about 18,600 years ago; - the end of gene flow about 6,200 years ago; - the size change about 1,160 years ago.

The ancestor S had about 147,000 gene copies, so \(N_e \approx 73{,}500\). That predicts \(\theta = 0.0045\) per site, right inside the observed \(\pi\) of 0.004–0.005. About 7 immigrant gene copies per generation (roughly 3.5 diploid migrants) moved in each direction between Ambavala and the southern ancestor. That is above the classic rule of thumb that about one migrant per generation is enough gene flow to keep populations from drifting apart.

Drift in coalescent units

Years are what we report, but drift is measured in coalescent units. Recall from lecture that one coalescent unit is \(2N_e\) generations, which is \(N\) generations when \(N\) counts gene copies. So the drift along a branch is its length in generations divided by its size, added up over epochs:

drift <- data.frame(
  branch      = c("Tsinjoarivo", "Tsinjoarivo", "Ambatovy", "Ambatovy",
                  "Ambavala", "Ambavala", "ancestor S"),
  epoch       = c("since the size change", "node S to the size change",
                  "since the size change", "node S to the size change",
                  "since the size change", "root to the size change", "root to node S"),
  generations = c(est$Tau_Crash, est$Tau_NodeS - est$Tau_Crash,
                  est$Tau_Crash, est$Tau_NodeS - est$Tau_Crash,
                  est$Tau_Crash, est$Tau_NodeR - est$Tau_Crash,
                  est$Tau_NodeR - est$Tau_NodeS),
  gene_copies = c(est$Theta_TsinCurr, est$Theta_TsinChange,
                  est$Theta_AmbatCurr, est$Theta_AmbatChange,
                  est$Theta_AmbavCurr, est$Theta_AmbavChange,
                  est$Theta_NodeS))
drift$coalescent_units <- round(drift$generations / drift$gene_copies, 3)
drift

# total drift along each branch
round(tapply(drift$coalescent_units, drift$branch, sum), 3)

In Tsinjoarivo, the 332 generations since the size change add 0.097 coalescent units of drift. That is more than the 4,988 generations before it (0.074), because the population is 20 times smaller. Ambatovy, still large, has drifted only 0.034 units since the split. This puts numbers on Q4: the short, recent crash adds more drift to Tsinjoarivo than the much longer stretch between the split and the crash, which is why the two southern populations differ so much for so young a split. Notice also that every branch is less than 1 coalescent unit long, even the one that is ~300,000 years old. Lecture showed that two lineages take 2 units, on average, to coalesce in a constant population. Most of these lineages had not finished coalescing when the populations split, which is why the joint spectra have many shared segragating sites.

Q13: Use the drift table to explain why a crash lasting only about 1,200 years can matter as much as the preceding 17,000 years.

Rescaling is not refitting

What if the mutation rate were wrong? The data fix \(N\mu\), \(T\mu\) and \(Nm\), not \(N\), \(T\) and \(m\) separately. A different \(\mu\) therefore simply rescales the estimates, with no need to fit again: sizes and times scale by \(\mu_{old}/\mu_{new}\), and migration rates by \(\mu_{new}/\mu_{old}\). (We checked this with fsc itself. Simulating model 10 at these estimates, and again with \(\mu\) doubled, sizes and times halved and migration rates doubled, gave the same SNP counts and \(F_{ST}\) to within simulation noise.)

# Re-express the estimates for a different mutation rate
rescale_mu <- function(est, mu_new) {
  factor <- MU / mu_new
  out <- est
  for (column in names(est)) {
    if (startsWith(column, "Theta_") || startsWith(column, "Tau_")) out[[column]] <- est[[column]] * factor
    if (startsWith(column, "MIG_"))                                  out[[column]] <- est[[column]] / factor
  }
  out
}
est_2x <- rescale_mu(est, 2 * MU)      # what if mu were twice as large?

# what changes ...
c(original = est$Tau_Crash, mu_doubled = est_2x$Tau_Crash)                              # generations
c(original = est$Theta_NodeS / 2, mu_doubled = est_2x$Theta_NodeS / 2)                  # diploid Ne
# ... and what does not
c(original = est$Theta_TsinChange / est$Theta_TsinCurr,
  mu_doubled = est_2x$Theta_TsinChange / est_2x$Theta_TsinCurr)                         # fold decline
c(original = est$Theta_NodeS * est$MIG_ANC_02,
  mu_doubled = est_2x$Theta_NodeS * est_2x$MIG_ANC_02)                                  # immigrants per generation
c(original = est$Tau_Crash / est$Theta_TsinCurr,
  mu_doubled = est_2x$Tau_Crash / est_2x$Theta_TsinCurr)                                # coalescent units

# Calendar dates under a range of mutation rates and generation times.
# These ranges are illustrative, chosen for this exercise -- not published uncertainty.
mu_values  <- c(1.0e-8, 1.52e-8, 2.0e-8)
gen_values <- seq(2.5, 4.5, by = 0.5)
grid <- data.frame()
for (mu_new in mu_values) {
  scaled <- rescale_mu(est, mu_new)
  for (g in gen_values) {
    grid <- rbind(grid, data.frame(mu = mu_new, generation_time = g, event = times$event,
                                   years = c(scaled$Tau_NodeR, scaled$Tau_NodeS,
                                             scaled$Tau_Change, scaled$Tau_Crash) * g))
  }
}
grid$event    <- factor(grid$event, levels = times$event)
mu_names      <- c("mu = 1.0e-8", "mu = 1.52e-8 (as in the paper)", "mu = 2.0e-8")
grid$mu_label <- factor(mu_names[match(grid$mu, mu_values)], levels = mu_names)

# for the figure: shade the LGM and human arrival only in panels whose dates reach them
periods <- data.frame(period = c("LGM", "human arrival"), from = c(19000, 1000), to = c(26500, 2000))
bands <- data.frame()
for (ev in levels(grid$event)) {
  lo <- min(grid$years[grid$event == ev])
  hi <- max(grid$years[grid$event == ev])
  for (k in 1:nrow(periods)) {
    if (periods$to[k] > lo && periods$from[k] < hi) {           # the period overlaps this panel
      bands <- rbind(bands, data.frame(event = ev, period = periods$period[k],
                                       ymin = max(periods$from[k], lo), ymax = min(periods$to[k], hi)))
    }
  }
}
bands$event <- factor(bands$event, levels = levels(grid$event))

# earliest and latest date of each event over the whole grid
data.frame(event    = levels(grid$event),
           earliest = round(tapply(grid$years, grid$event, min)),
           latest   = round(tapply(grid$years, grid$event, max)), row.names = NULL)

ggplot(grid, aes(x = generation_time, y = years, colour = mu_label)) +
  # the two periods the paper asks about, where they fall inside a panel
  geom_rect(data = bands, aes(xmin = -Inf, xmax = Inf, ymin = ymin, ymax = ymax, fill = period),
            inherit.aes = FALSE, alpha = 0.25) +
  scale_fill_manual(values = c("LGM" = "#4477AA", "human arrival" = "#DDAA33"), name = NULL) +
  geom_line(linewidth = 0.9) +
  geom_point(size = 1.4) +
  facet_wrap(~ event, nrow = 1, scales = "free_y") +
  scale_y_log10(labels = function(y) format(y, big.mark = ",", scientific = FALSE)) +
  scale_colour_manual(values = c("#4477AA", "black", "#CC6677"), name = NULL) +
  labs(x = "generation time (years)", y = "years ago (log scale)") +
  theme_bw(base_size = 11) +
  theme(legend.position = "bottom", panel.grid.minor = element_blank())

Absolute age of the four events in the best model-10 replicate, under illustrative mutation rates (lines) and generation times (x axis). Shaded where they fall inside a panel: the Last Glacial Maximum (19,000–26,500 years ago) and human arrival (1,000–2,000 years ago).

You should find that the size change stays between about 630 and 2,270 years ago across the whole grid, always long after the LGM. The Tsinjoarivo–Ambatovy split, by contrast, ranges from about 10,000 to 36,000 years ago. That puts it after the LGM, during it, or before it, depending on values of \(\mu\) and \(g\) that nobody knows that precisely.

NoteWhen fsc rescales for you

Our fits used the monomorphic sites, so the estimates are already in gene copies and generations for the \(\mu\) in the .tpl. If you fit only the polymorphic sites (option -0, or -l with a reference parameter), the data cannot fix the absolute scale. fsc then reports a scaling factor \(r = S_{obs}/(\mu T_{tot})\), by which you multiply sizes and times and divide migration rates (fsc manual, “Parameter rescaling”).

Q15: Using the figure, which conclusions survive the uncertainty in \(\mu\) and \(g\)? Can you say whether the decline happened after the LGM? After humans arrived? Can you say whether Tsinjoarivo and Ambatovy split before the LGM?

11. Takeaways

1. A joint SFS is a picture of shared history. The diagonal represents shared ancestry or gene flow, and the axes means private variation. The corners hold fixed differences. Each history leaves a different footprint, which is why spectra can distinguish models that a tree cannot.
2. fastsimcoal2 models are written backward in time, exactly like the coalescent from lecture: a split is a merge, sizes and samples are counted in gene copies, and the deme order in the .tpl sets the names of the .obs files.
3. Monomorphic sites and a mutation rate make the estimates absolute. They turn \(\theta\) into gene copies and generations; halving the copies gives diploid individuals, and a generation time turns generations into years. Dates and sizes scale with \(1/\mu\), and dates also with \(g\). Ratios, the order of events, \(Nm\) and drift in coalescent units do not depend on either, so those are the most robust results to report.
4. Know your likelihood’s units. fsc reports log10 likelihoods; AIC needs natural logs. And \(k\) is the number of free parameters in the .est file. Fold the model the way the data were folded: --foldedSFS for pairwise spectra from ANGSD.
5. Composite likelihoods are overconfident. A weight of 1 from ~30 million linked sites is not certainty. We did not discuss bootstrapping, but check that your conclusions are robust to plausible alternative hypotheses.
6. The best model can still be a poor fit to the data. Simulate from the best model and compare with the data. Here model 10 reproduces both the SNP counts and the \(F_{ST}\) results and that helps us have some confidence in the model, even if there is uncertainty regarding human impact versus paleoclimate (probably both).