bio-ecological-genomics-conservation-genetics

What this skill does

## 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

| Method | Time horizon | Strength | Bias mode |
|--------|--------------|----------|-----------|
| LDNe (NeEstimator V2 LD method) | Last 1-few generations | Single sample; SNP-friendly | Biased downward at small N; sensitive to physical linkage |
| SNeP (LD with physical-linkage correction) | Last 1-few generations | Corrects LDNe for chromosomal linkage in genomic data | Requires phased genotypes |
| GONE / GONE2 | Recent 100-200 generations | Detects bottleneck/expansion trajectory | Requires populated cM genetic map column (NOT physical position) |
| Heterozygote-excess | Last 1 generation | Detects recent bottleneck | Very low power; requires high heterozygosity |
| Temporal (Jorde-Ryman, Pollak) | Between two cohorts | Direct drift estimate | Requires temporal samples |
| Sibship (Wang COLONY) | Current generation | Pedigree-based from kinship | Computationally expensive |
| dadi (Gutenkunst 2009) | 10^3 - 10^5 generations | Diffusion SFS-based; analytical | Local-optima trap; <=3 populations practical |
| fastsimcoal2 (Excoffier 2013) | 10^3 - 10^5 generations | Simulation-based composite-likelihood | Optimization needs >= 50 replicates |
| moments (Jouganous 2017) | 10^3 - 10^5 generations | Moment-based ODE; faster than dadi | Same dimensionality limits |
| momi2 (Kamm 2020) | 10^3 - 10^5 generations | Analytical SFS via moments | Stiff at very recent times |
| Stairway Plot 2 | 10^3 - 10^5 generations | Folded SFS; no parametric model | Requires explicit mutation rate and generation time |
| PSMC | 10^4 - 10^6 generations | Whole-genome pairwise coalescent | Requires >= 20x WGS coverage |
| MSMC2 | Multi-individual deep coalescent | Better recent resolution than PSMC | Requires phased data; computational |
| msprime (Kelleher 2016) | Simulation (not inference) | Gold-standard neutral simulator | Used to generate data under a fitted model |
| SLiM 3/4 (Haller & Messer 2019) | Forward simulation | Non-WF, age structure, spatial | Eidos scripting; tree-sequence recording essential for speed |

## Decision Tree by Scenario

| Scenario | Recommended approach | Why |
|----------|---------------------|-----|
| Contemporary Ne from a single SNP sample | NeEstimator V2 LD method with Pcrit = 0.02 | Standard; widely accepted by conservation reviewers |
| LDNe applied to RAD-seq / WGS | Use SNeP (Barbato 2015) OR thin SNPs to >= 1 cM apart | LDNe assumes inter-locus r^2 reflects demography only; physical linkage inflates it |
| Recent Ne trajectory (bottleneck detection) | GONE2 with populated cM genetic-map column | GONE2 fails silently if cM column is zero (PLINK default uses physical position) |
| Deep demographic history (mammal/vertebrate) | PSMC with >= 20x WGS or Stairway Plot 2 from SFS | PSMC requires high coverage; Stairway Plot 2 needs only the SFS |
| Multi-population demographic inference | dadi / fastsimcoal2 with >= 50 independent optimization replicates | Likelihood surfaces have local optima; single-replicate inference is unreliable |
| Inbreeding from genomic data | F_ROH with length-class thresholds (bcftools roh OR detectRUNS) | F_IS is genome-averaged; F_ROH partitions inbreeding by time |
| Genetic-load assessment | Decompose realized vs masked via Bertorelle 2022 framework | Single "load" estimate hides the dynamics |
| Define management unit (MU) or ESU | Moritz 1994 reciprocal monophyly + nuclear divergence | Do NOT use BPP/BFD* species-delimitation methods (Sukumaran-Knowles 2017 oversplitting) |
| Forward simulation of selection + demography | SLiM 4 with tree-sequence recording + pyslim/tskit | Tree-sequence recording is 5-100x faster than mutation tracking |
| Comparing Ne across species | Cite Ne/Nc taxonomic variation (Hauser & Carvalho 2008) | The Ne/Nc = 0.1 default is wrong for many taxa |
| Genetic rescue decision | Frankham 2015 meta-analysis; usually favors rescue | Outbreeding-depression fear is overstated for most cases |

## 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.

```r
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