Introduce you to how genetic data can inform us about genetic diversity, evolutionary processes, and demographic histories of populations
Work with real Next Generation Sequencing (NGS) data without needing bioinformatics expertise
Learn to analyze Single Nucleotide Polymorphisms (SNPs) from a prepared Variant Call Format (VCF) file
Convert a VCF file into a genlight object (the format used by R’s adegenet package for population genetics)
Perform basic population genetic analyses including:
Principal Component Analysis (PCA) to visualize population structure
Observed heterozygosity per population
Inbreeding coefficients (FIS) per population
Pairwise FST values (genetic distance between populations)
Understand how these metrics can be used in conservation to:
Identify which populations are genetically similar or distinct
Detect populations with low genetic diversity
Pinpoint populations that may need increased conservation efforts
Identify suitable candidate populations for translocation
The Fennoscandian Arctic fox (Vulpes lagopus) population has been persecuted extensively for its fur and is now in a fragmented state due to climate change
We will compare it to the larger, continuous population in Russia (see Figure 1 for subpopulation locations)
As you work through the exercise, consider:
What evolutionary processes (genetic drift, gene flow, selection) may be shaping the genetic structure
What demographic events (bottlenecks, expansion, isolation) could have caused these processes
What you would expect from a large outbreeding population versus a small fragmented one
The data has been subset to one chromosome and randomly sampled for informative SNPs to make it lighter on your computer
This subset should represent the larger dataset, but the reduced resolution may affect the patterns observed in the final output
We need to install some packages to complete the exercise, take some time to make sure you have them loaded and that they are working.
# STEP 1: Install regular packages one by one
# Just copy and paste each line, one at a time
install.packages("vcfR")
install.packages("tidyverse")
install.packages("RColorBrewer")
install.packages("ggrepel")
install.packages("readxl")
# STEP 2: Install BiocManager (special installer)
install.packages("BiocManager")
# Load BiocManager
library(BiocManager)
# STEP 3: Install genetics packages using BiocManager
BiocManager::install("dartR")
BiocManager::install("adegenet")
# Load libraries
library(vcfR)
library(tidyverse)
library(RColorBrewer)
library(ggrepel);
library(readxl)
library(dartR)
library(adegenet)
Now we will import the data in vcf format and the accompanying metadata, see Figure 1 for how a vcf file is formatted below.
# STEP 1: Set your folder
# CHANGE THIS to your folder path!
setwd("/home/christopheralancockerill/Documents/Work/PhD_work/Teaching/Conservation of Populations/") # YOUR folder path
# STEP 2: Read your genetic data
vcf <- read.vcfR("popgen_practical_scaff22.biallelic.0missing_subset.vcf")
# STEP 3: Convert to genlight format
gl_data <- vcfR2genlight(vcf)
# STEP 4: Add population info
metadata <- read_excel("popgen_metadata.xlsx")
pop(gl_data) <- metadata$pop
ploidy(gl_data) <- 2
#check the file
vcf.gl
Now we will plot the data using ggplot and look to see if there is any populations structure.
# STEP 1: Run PCA (takes 2-3 minutes)
pca_result <- glPca(gl_data, nf = 5)
# STEP 2: Create a simple data frame for plotting
pca_data <- data.frame(PC1 = pca_result$scores[,1],
PC2 = pca_result$scores[,2],
Population = pop(gl_data))
# STEP 3: Calculate variance explained
var_exp <- round(pca_result$eig[1:2]/sum(pca_result$eig)*100, 1)
# STEP 4: Create the plot (this is the main result!)
ggplot(pca_data, aes(x = PC1, y = PC2, color = Population)) +
geom_point(size = 3) +
labs(x = paste0("PC1 (", var_exp[1], "%)"),
y = paste0("PC2 (", var_exp[2], "%)")) +
theme_bw() + scale_color_brewer(palette = "Set1")
# STEP 5: Save your plot
ggsave("my_pca_plot.pdf", width = 8, height = 6)
Now we are going to calculate some traditional populations genetics estimates from our data. You are probably familiar with Hardy-Weinberg equilibrium.
Briefly, Hardy-Weinberg equilibrium (HWE) is a fundamental principle in population genetics that describes the expected genotype frequencies in a population under certain conditions. It serves as a null hypothesis against which observed genotype frequencies can be compared to assess departures from expected values, which may indicate evolutionary processes such as genetic drift, migration, mutation, or natural selection.
Remember the principles?
#Now we will explore genetic variation within the populations with the package DartR
# STEP 1: Clean up data (remove monomorphic loci)
gl_clean <- gl.filter.monomorphs(gl_data)
gl_clean <- gl.recalc.metrics(gl_clean) # Calculate required metrics
gl_clean@other$loc.metrics <- as.data.frame(gl_clean@other$loc.metrics) # Fix format
# STEP 2: Calculate heterozygosity
gl.report.heterozygosity(gl_clean, method = 'pop')
Now lastly, It seems that these populations of arctic foxes vary in genetic variation and show substantial populations structure. Lets take a quick look at pairwise fst values to see if we can find anything interesting.
#Now lets take a look at the genetic distance of the populations by calculating Fst
#this might take a few minutes.
# STEP 3: Calculate Fst (genetic distance between populations)
gl.fst.pop(gl_clean)
In this exercise, we have shown you how to:
Assignment questions