Objectives of this practical session

  1. Summarize genetic diversity patterns.

  1. Summarize genetic differentiation with pairwise \(F_{ST}\)

  2. (OPTIONAL) Read Variant Call Format (VCF) files and create a file with genotypes for each SNP and individual.

To discuss:

0. Download the files and create a RStudio Project

Download the material from Moodle as a zip file. Uncompress the zip file.

Open RStudio, set the working directory as “Session>Set Working Directory>Choose Directory…” Select the directory of this TP (where all the files and folder downloaded from Moodle are).

To check that you are in the correct working directory, please copy and run the following command line. The output should be the folder where you have all the TP material.

getwd()

1. Install required packages

We will use the following packages

To install the required packages, copy paste the following command lines:

install.packages("vcfR")
install.packages("RColorBrewer")

Copy the next command lines to load the packages and check if they are installed.

library(vcfR)
library(RColorBrewer)

2. Analysis of Human exome data from Henn et al. (2015) PNAS

The folder “/Henn_et_al_data” contains exome data from seven human populations (file: Hennetal_genotypeMatrix.geno). This is a sub-set of 10,000 SNPs of the data published in Henn et al. (2015) PNAS (http://www.pnas.org/content/113/4/E440.abstract). This data was kindly provided by Stephan Peischl (http://www.bioinformatics.unibe.ch/about_us/staff/dr_peischl_stephan/index_eng.html), and cannot be used outside the scope of this course.

The seven human populations analysed are:
- Namibian San (SAN, https://en.wikipedia.org/wiki/San_people);
- Mbuti Pygmi from Democratic Republic of Congo (MBUTI, https://en.wikipedia.org/wiki/Mbuti_people);
- Algerian Mozabite (MOZABITE, https://en.wikipedia.org/wiki/Mozabite_people);
- Pakistani Pathan (PATHAN, https://en.wikipedia.org/wiki/Pathans_of_Punjab);
- Cambodian (CAMBODIAN);
- Siberian Yakut (YAKUT, https://en.wikipedia.org/wiki/Yakuts);
- Mexican Mayan (MAYA, https://en.wikipedia.org/wiki/Maya_peoples)

You can have a look at the input file “Hennetal_genotypeMatrix.geno” in the folder “Henn_et_al_data” by opening it with a text editor (for example, Notepad, Notepadd++, TextPad). For those of you not used to work with text editors, you can use Excel to open the file “Hennetal_genotypeMatrix_Excel.xlsx” in the folder “Henn_et_al_data”. This contains the genotype matrix data in Excel format.

Exercise 1
Have a look at the file with the data Hennetal_genotypeMatrix_Excel.xlsx, opening it with Excel.
- What are the values reported?
- What are rows and columns?

This is a genotypic data matrix where each entry codes for the genotypes as:

In this matrix each column corresponds to an individual, and each row corresponds to a site (position in the genome).

Note that in this case the mutations are said to be polarized, as we know what is the ancestral and derived state of each allele. This information was obtained by comparing, for each SNP, the alleles in modern humans with an outgroup (chimpanzee). It was considered that the allele found in chimpanzee reflect the ancestral state.

Note that in this case we do not really care what are the nucleotides. That is, imagine two SNPs where an individual is homozygote AA for SNP1 and heterozygote GT for SNP2. Those genotypes would be coded as 2 if the ancestral allele was C and the derived was A for SNP1, and coded as 1 for SNP2. When we consider the diversity we do not really gain information about the alleles. The information is in the patterns of diversity, and on sharing of alleles between individuals.

If the ancestral state is unknown, the genotypes would be coded as 0 for homozygote for the reference allele, 1 for heterozygote, and 2 for homozygote for alternative allele.

2.1 What is the most genetically diverse population?

2.1.1 Compute expected heterozygosity

We will now read the genotype matrix file with R and compute the expected heterozygosity for each sample/population.

The expected heterozygosity is computed for each SNP as the probability of taking two sequences from the population that are different, i.e. taking two sequences with two different alleles. This can then be averaged across SNPs to obtain a measure of genetic diversity for a given sample. In a sample of \(n_i\) sequences at a given site, the expected heterozygosity \(H_e\) can be computed as the mean across \(L\) sites:

\(H_e = \frac{1}{L} \sum_{i=1}^{L} 2 p_i (1-p_i)\)

where \(p_i\) is the frequency of one of the alleles at a particular site \(i\), and \((1-p_i)\) is the frequency of the alternative allele at site \(i\). Note that the sum \(\sum_{i=1}^{L}\) means that you will sum across all the values that \(i\) can take from \(1\) to \(L\). That is, the above formula is saying that we need to sum the heterozygosity for each SNP, \(2 p_1 (1-p_1) + 2 p_2 (1-p_2) + ... + 2p_L (1-p_L)\). Note that actually, the formula is a bit more complex as there is a correction factor because of the sample size, given by \(\frac{2n_i}{2n_i-1}\) for each site. NOTE: given that we have diploid individuals, each individual has two copies of a gene (one from the mother, one from the father), and hence the sample size is \(2 n_i\).

IMPORTANT: The name “expected heterozygosity” (He) is a source of confusion. The He is a summary statistic that is a measure of genetic diversity of a sample. Even though it is called expected heterozygosity this is something we observe! It is something we compute based on the allele frequencies in a sample, and hence it is observed. Thus, in some plots I called “observed expected heterozygosity” just to clarify that He is something we compute based on observed allele frequencies.

Copy the R commands below to:
- read the genotype table
- get the matrix with the genotypes for each sample
- compute the mean expected heterozygosity for each population.
- compute the mean observed heterozygosity for each population.

# load functions to compute the expected heterozygosity
source("utilfunctions.r")

# read the genotype matrix
geno <- as.matrix(read.table("./Henn_et_al_data/Hennetal_genotypeMatrix.geno", header=T, stringsAsFactors = F, na.strings = "NA"))

# get the label of the individuals
inds <- colnames(geno)

# vector with the population name
pop.names <- c("San","Mbuti","Mozabite","Pathan","Cambodian","Yakut","Maya")

# get the columns that correspond to each population
sample_inds <- get_index_ind(inds, pop.names)

# Get the expected heterozygosity for each sample
mean_exp_het <- ExpHet_mean(geno, sample_inds, pop.names)
mean_exp_het 
##        San      Mbuti   Mozabite     Pathan  Cambodian      Yakut       Maya 
## 0.13604545 0.13207806 0.11056417 0.10324971 0.09379727 0.09311416 0.08445250
# Compare the het for each pop in a barplot
# create a pallete with one color for each pop using a qualitative pallete
mycols <- brewer.pal("Paired", n=length(pop.names))
barplot(mean_exp_het, names.arg = pop.names, col=mycols, border=mycols, beside=T, ylab="mean observed Expected heterozygosity")

Exercise 2.1.1.
- What is the population with the higher expected heterozygosity?
- What is the population with the lower genetic diversity?
- What could explain differences between populations, i.e. why don’t all populations have the same expected heterozygosity?

2.1.2. Compute the observed heterozygosity for each sample

The observed heterozygosity is the actual proportion of heterozygote sites per individual. For a given population, the observed heterozygosity is the averaged across individuals.

Follow the script below to:
- Compute the observed heterozygosity for each individual
- Compute the mean observed heterozygosity for each sample

The script below gives you an example of how to compute the observed heterozygosity for one individual.

# compute the mean observed heterozygosity per population
mean_obs_het <- obs_het_mean(geno, sample_inds, pop.names)
mean_obs_het
##        San      Mbuti   Mozabite     Pathan  Cambodian      Yakut       Maya 
## 0.13450133 0.13032388 0.10911681 0.10095763 0.09328206 0.09182285 0.08418569

After computing the mean oberved heterozygosity for each population, you can compare it with the mean expected heterozygosity for each population using a barplot.

# Plot the mean_obs_het and mean_exp_het
barplot(rbind(mean_obs_het,mean_exp_het), col=mycols[1:2], beside=T, names.arg = pop.names)

The observed Heterozygosity (Ho) is light blue, and the observed expected heterozygosity (He) is dark blue.

Exercise 2.1.2.
- What is the location with higher observed heterozygosity?
- What is the location with lower observed heterozygosity?
- What is the best measure to quantify the genetic diversity of a population - the expected heterozygosity or the observed heterozygosity?

2.2. Estimate the inbreeding coefficient (\(F_{IS}\))

By comparing the observed heterozygosity with the expected heterozygosity, we can get an idea of whether the population from which we sampled the data is mating at random, i.e. according to Hardy-Weinberg equilibrium. Recall that if our samples came from populations under Hardy-Weinberg equilibrium, the observed heterozygosity (computed based on genotype frequencies) should be the same as the expected heterozygosity (computed based on allele frequencies).

We can obtain an estimate for the inbreeding coefficient \(F_{IS}\) for each individual based on the excess of homozygotes due to deviations from Hardy-Weinberg, as

\(F_{IS}=1-\frac{H_o}{H_e}\)

where \(H_o\) is the mean observed heterozygosity across sites for a given individual, and \(H_e\) is the observed mean expected heterozygosity of a given population.

Copy the command lines below to compute the estimator of \(F_{IS}\) for each population using the code provided below.

# Estimate the inbreeding coefficient per population
fis <- 1-(mean_obs_het/mean_exp_het)

# Plot the FIS per sample
barplot(fis, names.arg = pop.names, col=mycols, border=mycols, beside=T, ylab="mean FIS per pop", ylim=c(-0.5,0.5))

Exercise 2.2
- What are the range of values that \(F_{IS}\) can take, i.e. what is the minimum and maximum values?
- What is the population with higher inbreeding coefficient?
- Is there a population with evidence of random mating?
- What factors that lead to deviations from random mating can explain the observed values?

3. Estimate the pairwise \(F_{ST}\) among all pairs of populations

The \(F_{ST}\) can be seen as a measure of population structure, and is a function of the allele frequencies within each sub-population, and the overall allele frequencies across all sub-populations. \(F_{ST}\) varies from 0 to 1, where 0 indicates no genetic differentiation and 1 indicates that populations are fully differentiated and that fixed different alleles at all loci.

For a pair of populations, one common estimator of \(F_{ST}\) is based on comparing the probability of taking two different gene copies from within a sub-population, with the probability of taking two different gene copies from two different sub-populations. This is given by:

\(\hat{F}_{ST}=1-\frac{\hat{H}_w}{\hat{H}_b}\)

where \(H_w\) is the mean expected heterozygosity of the two sub-populations (this is the mean within-populations, hence the subscript \(w\)) with the total expected heterozygosity based on the probability of taking two different gene copies from each population \(H_b\) (this is the probability between populations). This is usually called the Hudson estimator (Hudson et al. 1992, Genetics http://www.genetics.org/content/132/2/583.long).

DETAIL: There are many estimators of \(F_{ST}\), and there is still research on finding the estimators with the best properties. That is, an estimator that will be close to the real population \(F_{ST}\), which is an unknown parameter.

Copy the R commands below to:
- read the genotype table
- get the matrix with the genotypes for each sample
- compute the pairwise \(F_{ST}\) for all pairs of populations.

# load functions to compute the expected heterozygosity
source("utilfunctions.r")

# read the genotype matrix
geno <- as.matrix(read.table("./Henn_et_al_data/Hennetal_genotypeMatrix.geno", header=T, stringsAsFactors = F, na.strings = "NA"))

# get the label of the individuals
inds <- colnames(geno)

# vector with the population name
pop.names <- c("San","Mbuti","Mozabite","Pathan","Cambodian","Yakut","Maya")

# get the columns that correspond to each population
sample_inds <- get_index_ind(inds, pop.names)
# Compute the pairwise FST
pairfst <- matrix(NA, ncol=length(pop.names), nrow=length(pop.names))
# go through each pair of populations
for(i in 1:(length(pop.names)-1)) {
  for(j in (i+1):length(pop.names)) {
    # call the function to compute the pairwise FST between each pair of populations
    pairfst[i,j] <-  getFst(geno[,sample_inds[[i]]],geno[,sample_inds[[j]]])       
  }
}

# print the matrix with the pairwise FST values
colnames(pairfst) <- pop.names
rownames(pairfst) <- pop.names
pairfst

# use the function to make the plot
plotFst(pairfst, pop.names)

Exercise 3
Answer the following questions, giving a justification:
- 3.1. What is the pair of populations with the higher genetic differentiation?
- 3.2. Give two examples of demographic history events or factors that lead to high genetic differentiation?
- 3.3. What is the pair of populations that likely experienced higher historical gene flow levels?
- 3.4. What is the pair of populations that likely diverged more recently?
- 3.5. What is the pair of populations that likely experienced less historical gene flow levels?
- 3.6. What is the pair of populations that likely had the older divergence?
- 3.7. Given the expected heterozygosity and \(F_{ST}\) values, what population likely have the smaller effective population size, i.e., experienced more genetic drift?
- 3.8. Are these results compatible with the Out of Africa hypothesis?

4. VCF (variant call format) files- OPTIONAL

Follow the command lines below to read a VCF file, and based on that create a matrix just with the genotypes, where each row is a site and each column an individual. In a genotypic data matrix each entry codes for the genotypes as 0 (homozygote for reference allele), 1 (heterozygote), 2 (homozygote for alternative allele).

In a VCF file the genotypes are coded using the following notation:

Note that since we are just interested in variation, usually we do not distinguish in the genotypes what is the reference and alternative allele. The reference and alternative allele are saved in a specific column of the VCF file.

While doing this section of the tutorial, you should try to understand the meaning of different sections paying special attention to the way genotypes are coded and to the depth of coverage (DP).

We will use as an example data from the diploid, oomycete pathogen Phytophthora infestans. This dataset contains data for 18 individuals genotyped at 22,031 SNPs. This data was statistically phased.

More info here: https://cran.r-project.org/web/packages/pinfsc50/index.html

Copy the lines of code below to the console.

# check that you are in the correct working directory
getwd()

# read VCF file
vcf <- read.vcfR("./VCF/pinf_sc50.vcf", verbose = FALSE )
# print info about the vcf
vcf

# Show the meta information
head(vcf@meta)

This shows that there are several statistics in this VCF file. For instance, in the FORMAT you have the field:

Then, let’s have a look at the fixed information. The fixed information contains information about:

# Fix region (info for each site and sample)
head(getFIX(vcf))

Now, we can have a look at the genotype region. The genotype region has information for the genotype of each individual, together with info for the depth of coverage and quality of the genotype for each individual at each site. Below, we will have a look at data for SNPs 1 to 6, looking at the columns with format and genotypes for individuals 1 to 3. We will also check the Depth of coverage field.

# Genotype region
# show data for sites 1 to 6, and show columns with format, and genotypes for individuals 1 to 3
vcf@gt[1:6, 1:3]

Exercise 4.1

After copying the command lines above to look at the fix information and the genotype region corresponding to sites 1 to 6 and individuals 1 to 3, answer the following questions.
Look at individual DDR7602 at position 136 of Supercontig_1.50.
- What is the most likely genotype for individual DDR7602 at position 136 of Supercontig_1.50?
- What is the depth of coverage for individual DDR7602 at position 136 of Supercontig_1.50?
- What is the genotype quality for individual DDR7602 at position 136 of Supercontig_1.50?
Look at individual BL2009P4_us23 at position 41 of Supercontig_1.50.
- What is the number of reads for the reference and alternative alleles?
- What is the most likely genotype?
- What is the genotype quality?

The next steps are done to convert the genotypes of VCF into a genotype matrix of 0, 1 and 2.

# Get a matrix with the genotypes
# 1st. discard indels
vcf <- extract.indels(vcf, return.indels = FALSE)

# 2nd. get only SNPs with 2 alleles
# alternative allele can be used to detect multi-allelic sites
# we just want SNPs, so we keep only sites with 1 letter in alternative
alt <- vcf@fix[,5]
snp_sites_i <- which(alt=="A" | alt=="T" | alt=="C" | alt=="G")
vcf <- vcf[snp_sites_i,]
vcf

# 3rd. get only the genotype field (GT)
gt <- extract.gt(vcf, element="GT", as.numeric = FALSE)
str(gt)

# check how many genotypes of different types we have
table(gt)

# Genotypes are coded as 0|0, 0|1, 1|0 and 1|1
# we want to transform this into the values 0, 1, 1, 2
# as these correspond to
# homozygote for reference allele (0|0), 
# heterozygote, heterozygote (0|1) or (1|0),
# homozygote for alternative allele (1|1)

# replace elements of matrix 
gt[gt=="0|0"] <- 0
gt[gt=="0|1"] <- 1
gt[gt=="1|0"] <- 1
gt[gt=="1|1"] <- 2

# 4th - tranform into numeric (this is optional)
gt2 <- apply(gt, c(1,2), as.numeric)
str(gt2)

# 5th - save the file (here just saving 10 lines to save time)
write.table(gt[1:10,], file="test.geno", quote = FALSE, col.names = names(vcf@gt[1,-1]), row.names = FALSE)

# 6th - plot with genotypes
# here showing just first 200 snps
nsites <- 200
image(1:nrow(gt2[1:nsites,]), 1:ncol(gt2[1:nsites,]), gt2[1:nsites,],
      xlab="sites", ylab="individuals", col=c("blue","red","yellow"))
legend("topright", c("GT=0","GT=1","GT=2", "NA"), fill=c("blue","red","yellow","white"))

Exercise 4.2
Based on the image of the genotypes across sites answer the following questions.
- Identify one individual with an unusually high proportion of missing data.
- Identify one individual with an unusually high proportion of heterozygote sites.
- Identify one SNP with unusually high proportion of missing data.