Performing Parametric Bootstrap in fastsimcoal2

What are limitations of parameter estimation?

The site frequency spectrum (SFS) describes the distribution of derived allele frequencies across a set of genomic sites and summarizes patterns of genetic variation within a population or among populations.While it captures many genetic signatures of demographic history, it represents a substantial reduction of the information contained in the underlying genomic data. Additionally, demographic events leave delayed genetic signatures because allele frequencies require time to respond to changes in population size, migration, or structure. This delay is particularly relevant in diploid organisms, where coalescent processes influence the expected distribution of variation. Combined with uncertainty in mutation rates, generation times, recombination rates, and the assumed demographic model, these factors introduce uncertainty into parameter estimates and limit their precision. fastsimcoal2 may converge on a single maximum likelihood estimate and identify parameter combinations that optimize the likelihood of the observed SFS. However, this estimate represents the most likely parameterization under the assumed model and does not necessarily correspond to the true historical value. This tutorial directly ties into the fastsimcoal2 protocol provided on my website under https://www.timlheller.com/fastsimcoal-introduction

What is Parametric Bootstrap?

Parametric bootstrap is a resampling approach used to estimate uncertainty around parameter estimates. It generates synthetic datasets under the fitted demographic model, re-estimates parameters from each simulated dataset using the same inference procedure, and evaluates the distribution of recovered parameter values. This technique works with the following assumptions

Assumptions

• The fitted demographic model is a sufficiently accurate approximation of the true evolutionary history.

• The maximum likelihood parameter estimates are sufficiently close to the true underlying parameters.

• The stochastic variation observed among simulated datasets reflects the uncertainty expected in the empirical dataset.
  

How to simulate datasets in fastsimcoal2?

fastsimcoal2 does have the options to simulate SFS based on a pre-defined model. This differs from the standard inference workflow, where fastsimcoal2 internally simulates SFS repeatedly to evaluate the likelihood of parameter combinations against the observed data. Here, simulated SFS are generated explicitly and treated as independent datasets for uncertainty estimation. UUsing the best-fitting model parameters from previous optimization runs, we can directly analyze the methodological stability of estimated parameters based on available and simulated data.

In fastsimcoal, this works by supplying a template file with all parameters estimated to generate SFS based on this evolutionary scenario. In my case, I used the *_maxL.par file with the parameters of highest likelihood from the last run to generate the required .par file. I replaced the number of blocks with the estimated number of independent genomic regions to simulate. This value was derived from the estimated genome size divided by the linkage disequilibrium cutoff: 1.4 Gb / 20 kb = approximately 70,000 independent blocks.

N_BLOCKS=70000
PREFIX=Test_FSC
PAR="Bootstrap_$REPLICATES"
sed "s/^1 0$/$N_BLOCKS 0/" ./"$PREFIX"/"$PREFIX"_maxL.par > $PAR.par

#With the blocksize set to the LD cutoff

BLOCKSIZE=2e4
sed -i "s/^FREQ 1 0 \(.*\) OUTEXP/DNA $BLOCKSIZE 0 \1 OUTEXP/" $PAR.par

I then generated 100 replicate SFS using this modified .par file

# Run fastsimcoal to generate bootstrap SFS
REPLICATES=100
FSC=./fsc28
$FSC -i $PAR.par -n $REPLICATES -j -d -s0 -x -c 24 -I -q

Each one of these 100 replicate folders required their identical .tpl and .est file, identical to the original analysis, to re-estimate parameters on synthetic datasets.

# Copy .tpl and .est files into each bootstrap replicate folder
for i in $(seq 1 $REPLICATES); do
replicate_dir=./"$PAR"/"$PAR"_"$i"
cp -rf "$PREFIX".est "$replicate_dir/$PAR.est"
cp -rf "$PREFIX".tpl "$replicate_dir/$PAR.tpl"
done
echo "Finished copying .est and .tpl files to $PAR"

How to re-analyze simulated dataset?

First, we create the folders for the summary files and define relevant parameters, such as the number of simulations ($SIMS) and the number of independent runs ($IND_RUNS) for each SFS estimation.

mkdir -p "$PAR"/Replicates
SIMS=10000
ITER=20
IND_RUNS=20

Then, we run the fastsimcoal2 analysis on each individual SFS, repeated by the number of $IND_RUNS

# Run Fastsimcoal on each bootstrap replicate
for tpl in ./"$PAR"/"$PAR"_*/"$PAR".tpl; do
run=$(dirname "$tpl")
replicate_id=$(basename "$run")
log_file="Replicates/$replicate_id.log"
# Move into the replicate directory
cd "$run" || exit
# Run fastsimcoal repetitions within the directory
### Clean old files
rm -f $output.runs.tsv
rm -rf $replicate_id.cumul.data.tsv
for i in $(seq 1 $IND_RUNS); do
$FSC -t "$PAR".tpl -e "$PAR".est -M -L "$ITER" -n "$SIMS" -d -I -x -y 3 -c 24 -q >> $replicate_id.cumul.data.tsv
done
cd - || exit
## Append to summary file
grep "Iter" "$PAR"/"$replicate_id"/"$replicate_id".cumul.data.tsv | sed "s;lhood=;;" | cut --complement -f1,2 | sort -k15 | head -n1 >> Summary.tsv
echo "Wrote Simulation $replicate_id Parameters to Summary File"
done

How to calculate Confidence Intervals from Parametric Bootstrap in fastsimcoal2?

After obtaining the final Summary.tsv file, containing the best fitting runs across these 20 repetitions for each of the 100 replicate SFS, we move to R to calculate Confidence Intervals (CI) and check the distribution of the replicate parameters.

## Load Summary.tsv into R and properly name the columns.
bootstrap = read.csv("Summary.tsv", header = F, sep = "\t")
boot_params <- bootstrap[, 1:13]

## Rename parameters for better visualization
param_names <- c("N_CAU", "N_BIT", "N_TRO", "N_VOE",
"T_CAUBIT", "T_BITVOE",
"m_BITCAU", "m_CAUBIT", "m_VOEBIT", "m_BITVOE",
"m_09CAUHYB", "m_19BITHYB", "m_19CAUHYB")
colnames(boot_params) <- param_names

Calculate 95% percentile confidence intervals across all estimated parameters.

CI <- t(apply(boot_params, 2, function(x) {
quantile(x, probs = c(0.025, 0.975), na.rm = TRUE)
}))

colnames(CI) <- c("CI_2.5", "CI_97.5")
CI_table <- cbind(
Median = apply(boot_params, 2, median),
Mean = apply(boot_params, 2, mean),
StdDev = apply(boot_params, 2, sd),
CI
)

print(CI_table)

Now, we want to plot all the replicate values to look for skews and properly set estimation boundaries in the .est file. Ideally, recovered parameter estimates should cluster around the original simulation values. However, distributions may be asymmetric or non normal, particularly for parameters with boundary constraints, low information content, or strong parameter correlations.

Parametric bootstrap confidence intervals should therefore not be interpreted as absolute uncertainty around the true demographic history. Instead, they quantify the expected variability of parameter estimates under the assumptions of the fitted model and sampling process. Narrow confidence intervals indicate parameter stability under the model, whereas broad intervals suggest limited information or parameter uncertainty.

## Check Histograms for the distribution and skew
par(mfrow = c(4, 4), mar = c(3, 3, 3, 1))
for (i in 1:13) {
hist(boot_params[, i],
main = param_names[i],
xlab = "",
col = "grey85",
border = "grey40")
# Add 95% CI lines
abline(v = CI[i, 1], col = "blue", lwd = 2, lty = 2)
abline(v = CI[i, 2], col = "blue", lwd = 2, lty = 2)
# Add median
abline(v = median(boot_params[, i]), col = "red", lwd = 2)
}

## Reset plotting layout
par(mfrow = c(1, 1))

And finally, we print our CI to include them in our report.

write.table(CI_table, "CI_bootstrap.csv", sep = ";", quote = F, col.names = NA)

How to cite this article?

Heller, T. L. (2026). Performing Parametric Bootstrap in fastsimcoal2. Tim L. Heller. https://www.timlheller.com/parametric-bootstrap-fastsimcoal

Impressum