Bio Ecological Genomics Conservation Genetics

Version Compatibility

Reference examples tested with: hierfstat 0.5+, adegenet 2.1+, detectRUNS 2.0+, poppr 2.9+, NeEstimator V2.1+, GONE2, SNeP 1.1+, SLiM 4+, msprime 1.3+, bcftools 1.19+

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('<pkg>') then ?function_name to verify parameters
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Conservation Genetics

"Assess the genetic health of my endangered population" -> Estimate Ne by time horizon (contemporary LD, recent trajectory, deep demographic history), measure inbreeding via F_ROH with length-class thresholds, decompose realized vs masked genetic load, and interpret against the modern 100/1000 Ne rule.

• R: hierfstat::basic.stats() for F-statistics and diversity
• R: detectRUNS::consecutiveRUNS.run() or bcftools roh -G30 for ROH detection
• CLI: Ne2.exe option_file.ne2 for NeEstimator V2 contemporary Ne (option-file driven, NOT CLI flags)
• CLI: gone2 for recent Ne trajectory from genome-wide LD (requires populated cM map)

The Single Most Important Modern Insight -- There is No Single Best Ne Estimator

The most common methodological failure in conservation genetics is treating "Ne" as one quantity. It is not. Each Ne estimator captures a different time horizon with different biases:

• LD-based (NeEstimator V2, Waples 2006): Last 1-few generations; biased downward at small N
• GONE / GONE2 (Santiago 2020): Recent 100-200 generations trajectory; requires populated cM genetic map
• Heterozygote-excess: Last 1 generation only; very low power
• Temporal (Jorde-Ryman): Between two time-sampled cohorts
• Stairway Plot 2 / dadi / fastsimcoal2: 10^3 to 10^5 generations from SFS
• PSMC: 10^4 to 10^6 generations from single high-coverage genome

A second cornerstone: Ne/Nc ratio is NOT 0.1 universally. Frankham 1995 reported a median of 0.1 in vertebrates, but Hauser & Carvalho 2008 Fish Fish 9:333-362 documented Ne/Nc spanning 2-6 orders of magnitude smaller than census (i.e., 10^-2 to 10^-6) in marine fish with sex-biased and sweepstakes-recruitment reproduction. Conservation papers that assume Ne/Nc = 0.1 are wrong for many taxa.

A third: the 50/500 rule was revised to 100/1000 by Frankham et al. 2014 Biol Conserv 170:56-63. The 50/500 numbers came from 1980 and underestimate the Ne needed for adaptive maintenance. The modern thresholds are Ne >= 100 for short-term inbreeding-fitness protection and Ne >= 1000 for long-term adaptive potential.

Algorithmic Taxonomy

Decision Tree by Scenario

Genetic Diversity and F-Statistics

Goal: Compute standard population-genetic metrics with Weir-Cockerham estimators and bootstrap confidence intervals.

Approach: Convert VCF to genind via adegenet, then to hierfstat format. Run basic.stats for Fis/Fst/Ho/He and pairwise.WCfst with bootstrapping for population differentiation.

library(hierfstat)
library(adegenet)
library(poppr)

# VCF/genepop/PLINK -> genind via adegenet
data_genind <- read.genepop('populations.gen')
data_hf <- genind2hierfstat(data_genind)

# F-statistics with Weir-Cockerham estimators
bstats <- basic.stats(data_hf)
cat('Overall Fis:', bstats$overall['Fis'], '\n')
cat('Overall Fst:', bstats$overall['Fst'], '\n')

# Pairwise FST with bootstrap CIs
pw_fst <- pairwise.WCfst(data_hf)
boot_fst <- boot.ppfst(data_hf, nboot = 1000)

# Rarefied allelic richness (corrects for unequal sample sizes)
ar <- allelic.richness(data_hf)

# Private alleles unique to each population
pa <- private_alleles(data_genind, count.alleles = TRUE)

Runs of Homozygosity with Length-Class Dating

Goal: Detect autozygous segments and date the inbreeding events by ROH length class.

Approach: Use bcftools roh -G30 (HMM-based, density-independent) for VCF input OR detectRUNS::consecutiveRUNS.run() for PLINK. Bin ROH by length to date inbreeding events: > 16 Mb means parents/grandparents; 4-16 Mb means ~5 generations; 1-4 Mb means ~10-20 generations; < 1 Mb is deep background.

library(detectRUNS)

# detectRUNS: SNP-by-SNP scanning from PLINK files
runs <- consecutiveRUNS.run(
    genotypeFile = 'genotypes.ped',
    mapFile = 'genotypes.map',
    minSNP = 20,        # Minimum SNPs per run; prevents short-stretch false positives
    minLengthBps = 1e6, # 1 Mb minimum; shorter ROH usually not IBD-derived
    maxGap = 1e6,
    maxOppRun = 1,
    maxMissRun = 2
)

# F_ROH per individual = sum(ROH length) / autosomal genome length
# > 16 Mb ROH: parents/grandparents inbreeding (very recent)
# 4-16 Mb: ~5 generations back
# 1-4 Mb: ~10-20 generations (historical bottleneck)
# < 1 Mb: deep background, ancient homozygosity
froh <- Froh_inbreeding(runs, mapFile = 'genotypes.map', genome_wide = TRUE)

# Bin ROH by length to estimate inbreeding timing
runs_summary <- summaryRuns(runs, genotypeFile = 'genotypes.ped',
                            mapFile = 'genotypes.map')
plot_DistributionRuns(runs)

# Alternative: bcftools roh with HMM (density-independent; better for low-density data)
# -G30: per-sample genotype likelihood threshold (30 = high-quality calls)
# --AF-tag AF: use AF tag from VCF for allele frequencies (else --AF-file)
bcftools roh -G30 --AF-tag AF -o roh_results.txt input.vcf.gz

NeEstimator V2 — Option-File API

Goal: Estimate contemporary Ne via LD method from a single sample (Do et al. 2014 Mol Ecol Resour 14:209-214).

Approach: NeEstimator V2 is option-file driven, NOT CLI-flag driven. Create a .ne2 text option file specifying input format, methods (LD/HetExcess/Coancestry/Temporal), Pcrit thresholds, and output. Run with Ne2.exe option_file.ne2.

# Example NeEstimator V2 option file (option_file.ne2)
# Line order matters; defaults can cause silent failures
# See NeEstimator V2 documentation: https://github.com/bunop/NeEstimator2.X

# Input format
1 0                                   # 1=GenePop, 2=FSTAT; second value reserved
input_genotypes.gen

# Methods (1 = run; 0 = skip)
1 0 0 0                               # LD only; others (HetExcess, Coancestry, Temporal) off

# Pcrit cutoffs (allele frequencies below threshold excluded)
3                                     # number of Pcrit values
0.05 0.02 0.01                        # 0.02 is standard; 0.05 conservative; 0.01 sensitive

# Mating system: 0 = random mating; 1 = monogamy
0

# Output
output_results.txt

# Run NeEstimator V2 (Java-based)
java -jar NeEstimator.jar option_file.ne2
# Verify INFO output line: confirms which methods ran and which were skipped

# Common failure: "Ne = infinity" reported
# Diagnosis: insufficient drift signal (population larger than method can detect)
# Try multiple Pcrit; with genomic data add SNeP physical-linkage correction

GONE2 — Recent Ne Trajectory with the cM-Column Trap

Goal: Estimate the Ne trajectory over the last ~100-200 generations from genome-wide LD across recombination rates.

Approach: GONE2 requires a genetic map with POPULATED cM positions (not just physical position). PLINK MAP files default to cM=0; if this is not corrected, GONE2 will silently produce nonsense. Verify the cM column is populated from a genetic-map file or use Hi-C-derived recombination map.

# GONE2 input: PLINK BED/BIM/FAM or VCF + populated MAP
# CRITICAL: BIM/MAP cM column must be populated; default cM=0 produces silent failure

# Check cM column populated
head genotypes.bim   # column 3 should NOT be all zeros

# Run GONE2
# -t 4: threads
# -u 0.05: upper recombination rate bound (exclude pairs with r > 0.05)
# Smaller -u focuses on more recent generations
./gone2 -t 4 -u 0.05 genotypes.vcf

# Output: OUTPUT_GONE2 with generation, Ne, CI_low, CI_high

# Parse GONE2 output and plot Ne trajectory
gone_out <- read.table('OUTPUT_GONE2', header = TRUE, sep = '\t')

pdf('gone2_ne_trajectory.pdf', width = 8, height = 5)
plot(gone_out$generation, gone_out$Ne,
     type = 'l', lwd = 2, col = 'blue',
     xlab = 'Generations ago', ylab = 'Effective population size (Ne)',
     main = 'Recent Ne Trajectory (GONE2)', log = 'y')
polygon(c(gone_out$generation, rev(gone_out$generation)),
        c(gone_out$CI_low, rev(gone_out$CI_high)),
        col = adjustcolor('blue', alpha = 0.2), border = NA)
dev.off()

SNeP — Physical-Linkage-Corrected LDNe for Genomic Data

Goal: Apply physical-linkage correction to LDNe when SNPs come from RAD-seq or WGS (where inter-locus r^2 reflects chromosomal linkage, not just demography).

Approach: SNeP (Barbato 2015 Front Genet 6:109) implements Waples & Do's physical-linkage correction. Use when SNPs are not pre-thinned to >= 1 cM apart.

# SNeP input: PLINK PED/MAP or VCF + map file
# Multi-threaded; supports several corrections (sample size, mutation, phasing, recombination)
./SNeP1.1 -ped genotypes.ped -map genotypes.map -threads 4 \
          -mutationrate 1.4e-8 -out snep_results.txt

# Output: Ne estimates per recombination-rate bin

Deep Demographic History

Goal: Reconstruct Ne trajectory over thousands of generations from genome-wide variation.

Approach: Pick the appropriate tool by data type. PSMC for a single high-coverage diploid genome. Stairway Plot 2 for SFS from population samples. dadi or fastsimcoal2 for multi-population history with composite likelihood.

# --- PSMC: single high-coverage WGS (>= 20x) ---
bcftools mpileup -C50 -Q 30 -q 30 -f reference.fa sample.bam | \
    bcftools call -c | vcfutils.pl vcf2fq -d 10 -D 100 > consensus.fq
fq2psmcfa -q20 consensus.fq > consensus.psmcfa
psmc -N25 -t15 -r5 -p '4+25*2+4+6' -o sample.psmc consensus.psmcfa
psmc_plot.pl -u 1.4e-8 -g 5 -p sample_psmc_plot sample.psmc

# --- Stairway Plot 2: SFS-based (works with RAD-seq or WGS) ---
# Build SFS from VCF (use easySFS or vcf2sfs)
# Create blueprint.txt specifying nseq, L, mu, generation time
java -cp stairway_plot_v2.jar Stairbuilder blueprint.txt
bash blueprint.sh

# --- fastsimcoal2: multi-population composite-likelihood SFS ---
# REQUIRES >= 50 independent optimization replicates; single-run is unreliable
for i in {1..50}; do
    fsc27 -t template.tpl -e template.est -m -L 50 -n 100000 -q
done
# Compare likelihoods across replicates; report best AND distribution

Genetic Load — Realized vs Masked Decomposition

Goal: Decompose total genetic load into realized (currently expressed) and masked (heterozygous, potential) components per Bertorelle 2022 framework.

Approach: Use SnpEff/VEP to annotate variants for predicted functional effect. Compute realized load from homozygous-derived-allele counts at deleterious positions; compute masked load from heterozygous counts. Hedrick & Garcia-Dorado 2016 distinguish purging vs drift dynamics; report both.

# Conceptual workflow (Bertorelle et al. 2022 NRG 23:492-503)
# Realized load: count_homozygous_deleterious / total_deleterious_sites
# Masked load: count_heterozygous_deleterious / total_deleterious_sites
# Total load = realized + 0.5 * masked (assuming partial dominance)

# Annotate VCF with SnpEff or VEP first to classify deleterious vs neutral
# Then per individual:
# realized_load <- sum(genotype == 2 & severity == 'HIGH') / sum(severity == 'HIGH')
# masked_load <- sum(genotype == 1 & severity == 'HIGH') / sum(severity == 'HIGH')

# Cite Hedrick & Garcia-Dorado 2016 TREE 31:940-952 for purging vs drift:
# Strong-s alleles purge under inbreeding (potentially good)
# Weak-s alleles fix by drift (definitely bad)
# Net effect depends on selection-coefficient distribution
# Empirical purging example: Robinson 2018 Curr Biol 28:3487-3494 (Channel Island foxes)
# Simulation-based load assessment: Kyriazis 2021 Evol Lett 5:33-47

Forward Simulation with Tree-Sequence Recording

Goal: Simulate selection and non-Wright-Fisher demography efficiently using SLiM 4 with tskit tree-sequence recording.

Approach: SLiM 4 with treeSeqOutput() records the genealogy; pyslim + tskit adds neutral mutations a posteriori. This is 5-100x faster than mutation-tracking (Haller, Galloway, Kelleher, Messer, Ralph 2019 Mol Ecol Resour 19:552-566). The Eidos scripting language is NOT Python.

# Reference: SLiM 4+, pyslim 1.0+, tskit 0.5+, msprime 1.3+
# Note: SLiM scripts (.slim files) use Eidos, not Python.
# See https://messerlab.org/slim/ for full Eidos syntax.

# After running a SLiM simulation with treeSeqOutput():
import tskit
import pyslim
import msprime

ts = tskit.load('simulation.trees')
ts = pyslim.update(ts)  # update to current SLiM/pyslim conventions

# Recapitate: add coalescence above the SLiM tree root (neutral burn-in)
ts_recap = pyslim.recapitate(ts, recombination_rate=1e-8,
                              ancestral_Ne=10000, random_seed=42)

# Add neutral mutations after-the-fact via msprime
ts_mut = msprime.sim_mutations(ts_recap, rate=1e-8, random_seed=42)

# Export VCF for downstream analyses
with open('simulation.vcf', 'w') as f:
    ts_mut.write_vcf(f)

Per-Method Failure Modes

LDNe applied to RAD-seq SNPs without physical-linkage correction

Trigger: Running NeEstimator V2 LD method on RAD-seq or WGS genotypes thinned only by MAF, with multiple SNPs per chromosome at close physical distance. For RAD-seq biases more broadly, see Andrews 2016 Nat Rev Genet 17:81-92.

Mechanism: LDNe assumes inter-locus r^2 reflects only demographic LD (drift, random mating). With physically linked SNPs, much r^2 is chromosomal, NOT demographic. The estimator interprets this as "more drift has happened" and reports a smaller Ne.

Symptom: LDNe estimate dramatically lower than ecologically reasonable; downward bias compared with temporal or other estimators on the same population.

Fix: Thin SNPs to >= 1 cM apart OR use SNeP (Barbato 2015) which corrects for physical linkage explicitly. NeEstimator V2 has a Chrom flag for chromosome-aware LDNe; document its use.

GONE2 silent failure with PLINK MAP cM=0 (the default)

Trigger: Running GONE/GONE2 on PLINK output where the MAP file's cM column was never populated (defaults to 0 because PLINK uses physical position).

Mechanism: GONE/GONE2 use the cM positions to interpret recombination distances; cM=0 for all loci means the algorithm cannot distinguish linked from unlinked SNPs and produces nonsense.

Symptom: GONE2 output shows extreme Ne values or completely flat trajectory; results inconsistent across runs or with other methods.

Fix: Populate the cM column from a species-specific genetic map, OR build one from Hi-C data, OR use the linkage-rate approximation from PLINK --cm-map. Verify with head genotypes.bim — column 3 should not be all zeros.