Bio Ecological Genomics Conservation Genetics

Version Compatibility

Reference examples tested with: 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 effective population size (Ne) trajectories, measure F-statistics and heterozygosity, detect inbreeding through runs of homozygosity, and reconstruct demographic history from SNP or microsatellite data.

• R: hierfstat::basic.stats() for F-statistics and diversity
• R: detectRUNS::consecutiveRUNS.run() for ROH detection
• CLI: gone2 for recent Ne trajectory from phased genotypes

Assesses genetic health of populations through diversity metrics, effective population size estimation, inbreeding detection, and demographic history reconstruction.

Genetic Diversity with hierfstat

Basic population genetics statistics using Weir & Cockerham estimators:

library(hierfstat)
library(adegenet)

# genind object from adegenet (VCF, PLINK, or genepop input)
# read.genepop() or vcfR::read.vcfR() + vcfR2genind()
data_genind <- read.genepop('populations.gen')

# Convert to hierfstat format
data_hf <- genind2hierfstat(data_genind)

# Basic F-statistics: Fis (inbreeding), Fst (differentiation), Ho, He
bstats <- basic.stats(data_hf)
cat('Overall Fis:', bstats$overall['Fis'], '\n')
cat('Overall Fst:', bstats$overall['Fst'], '\n')

# Per-population statistics
pop_stats <- data.frame(
    Ho = colMeans(bstats$Ho, na.rm = TRUE),
    He = colMeans(bstats$Hs, na.rm = TRUE),
    Fis = colMeans(bstats$Fis, na.rm = TRUE)
)
pop_stats

Pairwise Fst

# Bootstrap confidence intervals for pairwise Fst
# nboot=1000: standard for publication-quality CIs
pw_fst <- pairwise.WCfst(data_hf)
pw_fst

# Bootstrap p-values
boot_fst <- boot.ppfst(data_hf, nboot = 1000)
boot_fst$ll  # lower 95% CI
boot_fst$ul  # upper 95% CI

Allelic Richness

Rarefied allelic richness corrects for unequal sample sizes:

# min.n: rarefaction to smallest population sample size
ar <- allelic.richness(data_hf)
cat('Rarefied allelic richness per population:\n')
colMeans(ar$Ar, na.rm = TRUE)

Private Alleles

library(poppr)

# Private alleles unique to each population
pa <- private_alleles(data_genind, count.alleles = TRUE)
private_counts <- rowSums(pa > 0)
cat('Private alleles per population:\n')
print(private_counts)

Runs of Homozygosity (ROH)

Goal: Detect autozygous genomic segments and quantify individual-level inbreeding from SNP genotype data.

Approach: Scan genotypes with detectRUNS to identify consecutive homozygous stretches exceeding length and SNP thresholds, then compute F_ROH as the fraction of the genome covered by ROH.

Detects autozygous segments indicating recent inbreeding:

library(detectRUNS)

# Input: PLINK .ped/.map files
# consecutiveRUNS.run: file-based SNP-by-SNP scanning (preferred for SNP arrays/WGS)
runs <- consecutiveRUNS.run(
    genotypeFile = 'genotypes.ped',
    mapFile = 'genotypes.map',
    minSNP = 20,        # min SNPs in a run; 20 prevents false positives from short homozygous stretches
    minLengthBps = 1e6, # 1 Mb minimum; shorter ROH are often not due to IBD
    maxGap = 1e6,       # max gap between consecutive SNPs within a run
    maxOppRun = 1,      # max opposing homozygotes allowed (genotyping errors)
    maxMissRun = 2      # max missing genotypes allowed per run
)

# Summary statistics
summary_runs <- summaryRuns(runs, genotypeFile = 'genotypes.ped',
                            mapFile = 'genotypes.map')

# F_ROH: inbreeding coefficient from ROH
# F_ROH = total ROH length / autosomal genome length
# <0.0625: background; 0.0625-0.125: offspring of half-cousins; >0.125: moderate inbreeding
# Froh_inbreeding calculates genome length automatically from the map file
froh <- Froh_inbreeding(runs, mapFile = 'genotypes.map', genome_wide = TRUE)

cat('F_ROH per individual:\n')
print(froh)

ROH Length Classes

# Plot ROH length distribution by population
plot_DistributionRuns(runs)

Effective Population Size (Ne)

GONE2: Recent Ne Trajectory (Linkage Disequilibrium)

Estimates Ne changes over the last ~200 generations from LD patterns:

# GONE2 is a standalone CLI tool (esrud/GONE2), not an R package
# Input: PLINK bed/bim/fam or VCF
# Requires at least 10,000 SNPs and 50 individuals for reliable estimates

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

# Parse GONE2 output (tab-separated: generation, Ne, CI_low, CI_high)
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')
dev.off()

cat('Current Ne estimate:', gone_out$Ne[1], '\n')

NeEstimator: Contemporary Ne (LD Method)

# NeEstimator v2 uses an option file (.ne2), not command-line flags
# Create option file specifying input, method, and parameters:
#   Input: genepop or FSTAT format genotype file
#   Method: LD (linkage disequilibrium, single time point)
#   Pcrit: 0.02 (exclude alleles with frequency < 2%)
#     0.02 balances bias (lower pcrit) vs precision (higher pcrit)
#   Output: tab-separated with Ne point estimate and 95% CI (jackknife + parametric)
# Run: java -jar NeEstimator.jar option_file.ne2
# Build from GitHub: bunop/NeEstimator2.X

Stairway Plot 2: Demographic History from SFS

Reconstructs Ne over thousands of generations using the site frequency spectrum:

# Step 1: Generate folded SFS from VCF
# Blueprint file specifies SFS, mutation rate, generation time
# nseq: number of haploid sequences (2 * n_diploid_individuals)
# L: total number of monomorphic + polymorphic sites analyzed
# mutation_rate: per-generation per-site (e.g., 1.4e-8 for vertebrates)
# year_per_generation: species generation time

# Step 2: Run Stairway Plot 2
java -cp stairway_plot_v2.jar Stairbuilder blueprint.txt

# Step 3: Execute the generated bash script
bash blueprint.sh

# Output: Ne trajectory plot and table

PSMC: Whole-Genome Pairwise Coalescent

For a single diploid genome (whole-genome sequencing):

# Consensus FASTQ from BAM
bcftools mpileup -C50 -Q 30 -q 30 -f reference.fa sample.bam | \
    bcftools call -c | vcfutils.pl vcf2fq -d 10 -D 100 > consensus.fq

# PSMC input format
fq2psmcfa -q20 consensus.fq > consensus.psmcfa

# Run PSMC
# -N25: number of iterations (25 is standard)
# -t15: maximum 2N0 coalescent time
# -r5: initial theta/rho ratio
# -p: time interval pattern (4+25*2+4+6: standard for most organisms)
psmc -N25 -t15 -r5 -p '4+25*2+4+6' -o sample.psmc consensus.psmcfa

# Bootstrap (100 replicates for confidence intervals)
seq 100 | xargs -P 4 -I {} \
    psmc -N25 -t15 -r5 -p '4+25*2+4+6' -b \
    -o sample_bootstrap_{}.psmc consensus.psmcfa

# Plot with species-specific mutation rate and generation time
# -u: mutation rate per generation per site
# -g: generation time in years
psmc_plot.pl -u 1.4e-8 -g 5 -p sample_psmc_plot sample.psmc

Ne Interpretation Thresholds

Bottleneck Detection

Goal: Test whether a population has experienced a recent genetic bottleneck using heterozygosity excess.

Approach: Compare observed heterozygosity against the equilibrium heterozygosity expected from the observed number of alleles; after a bottleneck, alleles are lost faster than heterozygosity, creating a transient excess detectable by Wilcoxon test.

# Heterozygosity excess test (Cornuet & Luikart 1996)
# Under mutation-drift equilibrium, He ~ expected from allele number
# After bottleneck, alleles lost faster than He, creating excess
library(hierfstat)

# Compare observed He against equilibrium He expected from allele number
# After bottleneck, alleles lost faster than He, creating transient excess
# IAM (infinite alleles model) for SNPs; TPM preferred for microsatellites
bstats <- basic.stats(data_hf)
he_obs <- bstats$perloc$Hs

# Heq under IAM: expected He given observed allele count k and sample size
# For microsatellites with k>2 alleles, Heq is simulated (use BOTTLENECK software)
# For biallelic SNPs (k=2), use 2*p*(1-p) equilibrium or simulation-based approach
# Simple approximation for microsatellites (k alleles):
k_alleles <- sapply(data_hf[, -1], function(x) length(unique(na.omit(x))))
n_samples <- nrow(data_hf)
# Ewens sampling formula approximation (valid for microsatellites, not biallelic SNPs)
heq_iam <- 1 - 1 / k_alleles

# Wilcoxon signed-rank test: He > Heq indicates bottleneck
# NOTE: for biallelic SNPs, use BOTTLENECK software (simulates expected Heq given k and n)
wilcox.test(he_obs, heq_iam, paired = TRUE, alternative = 'greater')

Related Skills

• landscape-genomics - Adaptive variation and genotype-environment associations
• species-delimitation - Taxonomic unit definition for conservation
• population-genetics/population-structure - Population stratification
• population-genetics/selection-statistics - Selection signatures
• variant-calling/vcf-basics - VCF input from RADseq/WGS