1.1 Mac Users

1.2 Windows Users

3.1 File formats

Let’s have a look at the file. While .fam and .bim are text file, the .bed is a binary file and cannot be opened with a text editor. These are Plink formats currently used in population genetics because they can effectively store a huge amount of genetic data. .fam file contains the sample information, .bim contains the genetic variant information (see the table below) and .bed codes the genotype calls. More information about the file formats in Plink can be found here.

3.2 Set the working directory

Before starting to work in R environment, we need to specify the working directory, which will be used by R to read and save files. In order to set our wirking directory we need to find the path leading to the directory called “popgen_tutorial”.

## For Mac users you path will be something like this:
#setwd("/Users/I/Desktop/Uni/forensic_genetics/popgen_tutorial")

## Fow Windows users the pathe will be something like this:
# setwd("c:/Documents/Desktop/Uni/forensic_genetics/popgen_tutorial")

# My path is:
setwd("/Users/bioserenina/Google Drive File Stream/Il mio Drive/lezioni_carlo/lezioni/esercitazione/popgen_tutorial")

Given that the genotype information are store in a binary file (.bed) which cannot be opened easily, we will make use of an R library specifically developed to work with genetic data: SNPRelate. You can find here a tutorial on its functions.

As a first step, you need to install the library, one of its dependency (the package “gdsfmt”) as well as other packages we will use in this tutorial using the command below.

# installing some packages
source("http://bioconductor.org/biocLite.R")
biocLite("gdsfmt")
biocLite("SNPRelate")
install.packages(c("ggplot2","RColorBrewer", "knitr", "kableExtra", "dplyr"))

Then, once you have installed a package, you need to “call” the library.

library(SNPRelate)
library(gdsfmt)

As reported in the SNPRelate tutorial, this library and the gdsfmt package can store the huge amount of genetic information into a different and more efficient file firmat: .gds.

To support efficient memory management for genome-wide numerical data, the gdsfmt package provides the genomic data structure (GDS) file format for array-oriented bioinformatic data, which is a container for storing annotation data and SNP genotypes. In this format each byte encodes up to four SNP genotypes thereby reducing file size and access time.

# create the "test.gds" file starting from our .bed, .bim, and .fam.
bed.fn <- "subsampled.bed"
fam.fn <- "subsampled.fam"
bim.fn <- "subsampled.bim"
snpgdsBED2GDS(bed.fn, fam.fn, bim.fn, "test.gds")

## Start snpgdsBED2GDS ...
##  BED file: "subsampled.bed" in the SNP-major mode (Sample X SNP)
##  FAM file: "subsampled.fam", DONE.
##  BIM file: "subsampled.bim", DONE.
## Wed May 20 09:51:52 2020     store sample id, snp id, position, and chromosome.
##  start writing: 786 samples, 29547 SNPs ...
##      Wed May 20 09:51:52 2020    0%
##      Wed May 20 09:51:52 2020    100%
## Wed May 20 09:51:52 2020     Done.
## Optimize the access efficiency ...
## Clean up the fragments of GDS file:
##     open the file 'test.gds' (5.7M)
##     # of fragments: 39
##     save to 'test.gds.tmp'
##     rename 'test.gds.tmp' (5.7M, reduced: 252B)
##     # of fragments: 18

snpgdsSummary("test.gds")

## The file name: /Volumes/GoogleDrive/Il mio Drive/lezioni_carlo/lezioni/esercitazione/popgen_tutorial/test.gds 
## The total number of samples: 786 
## The total number of SNPs: 29547 
## SNP genotypes are stored in SNP-major mode (Sample X SNP).

Once we have created the “test.gds” file, we need to open it and extract the information of interest, which are:

• samples information (sample.id)
• genetic variant information (snp.id)
• genotypes (g)

genofile <- snpgdsOpen("test.gds")
sample.id <- read.gdsn(index.gdsn(genofile, "sample.id"))
snp.id <- read.gdsn(index.gdsn(genofile, "snp.id"))
g <- read.gdsn(index.gdsn(genofile, "genotype"))  

4.1 Creation of the genotype matrix

Once we have finished the cleaning steps (already perfomed), we set up the genotypes matrix. Using the R commands dim(geno) we can see matrix size (rows and columns number), while in order to see part of the matrix we can use the command geno[1:5,1:5], which shows the first five rows and columns of the matrix.

geno = as.data.frame(g)
colnames(geno) = snp.id
rownames(geno) = sample.id

4.2 Missingness check

We know that the genotypes have already been cleaned by me, but the first rule in bioinformatics is “never trust the guy giving you some data”.

paste0("how many missing cells are there in our dataset? ", sum(geno==3))

## [1] "how many missing cells are there in our dataset? 0"

At this point, we need another information: the ancestry of the individuals in sample.id variable. You can find this information in the file “kgp3.pheno.tsv” and, with the commands below, you can add this information to the genotype matrix.

4.3 Adding population labels

# read the population labels file
pop_labels = read.table("kgp3.pheno.tsv", hea=T)

# add the population labels to geno
target = rownames(geno)
#head(pop_labels[match(target,pop_labels$SAMPLE),])
geno$POP = pop_labels[match(target,pop_labels$SAMPLE),"POP"]
geno$SAMPLE = pop_labels[match(target,pop_labels$SAMPLE),"SAMPLE"]
# check if the order of samples is correct: if the answer is TRUE can be happy
identical(rownames(geno), as.character(geno$SAMPLE))

## [1] TRUE

geno = geno[,c("POP", snp.id)]
#head(geno)[1:5]

How many populations do we have? As you can see with the R command table(geno$POP), there are 8 worldwide populations and 5 individuals of unknown ancestry.

# Divide the dataset in known populations and unknown individuals
geno_unk = geno[which(geno$POP %in% c("UNK1", "UNK2", "UNK3", "UNK4", "UNK5")),]
geno_pop = geno[!(geno$POP %in% c("UNK1", "UNK2", "UNK3", "UNK4", "UNK5")),]

In the end, our final dataset will be something like this, where each row is a sample, the first column is the population the sample comes from and the other columns are the SNPs:

5.1 Allele frequencies

At this point, the genotype matrix is ready and we can start the hard work: allele frequencies computation.

The frequency of an allele is defined as the number of the copies of the allelel in the population divided by the total number of gene copies in the population. In a diploid individual (in which all individuals carry two copies of each chromosome) with N individual, there are 2N gene copies. So the frequencies of alleles A1 and A2 are:

\(f_{A1}\) = \(\frac{N_{A1}}{2N}\)

and

\(f_{A2}\) = \(\frac{N_{A2}}{2N}\)

Let’s have a look at our genotypes, by taking the first five rows and five columns.

geno_pop[1:5,1:5]

##         POP rs116400033 rs115209712 rs112455420 rs375803869
## HG00096 GBR           2           0           0           0
## HG00097 GBR           0           0           0           0
## HG00099 GBR           0           1           0           0
## HG00100 GBR           0           0           1           0
## HG00101 GBR           0           0           0           0

There are four possible values stored in geno_pop: 0, 1, 2 and 3. For bi-allelic SNP sites, “0” indicates two A2 alleles, “1” indicates one A1 allele and one A2 allele, “2” indicates two A1 alleles, and “3” is a missing genotype.

Let’s see the second SNP in the geno_pop matrix: rs115209712.

Given that the geno_pop dataset harbours the count of one allele for each genotype (the cell), in order to compute the allele frequency, we simply need to sum all values in each column (the SNP) and divide it by 2N.

We can compute the frequency of the allele A across all samples or for each population separately.

# A1 frequency of the SNP rs115209712 in the entire dataset
sum(geno_pop$"rs115209712")/(2*dim(geno_pop)[1])

## [1] 0.1101152

# A2 frequency of the SNP rs115209712 in the entire dataset
1 - sum(geno_pop$"rs115209712")/(2*dim(geno_pop)[1])

## [1] 0.8898848

# Another way to compute the A2 frequency
(sum(geno_pop$"rs115209712" == 0)*2 + sum(geno_pop$"rs115209712" == 1))/(2*dim(geno_pop)[1])

## [1] 0.8898848

# A1 frequency of the SNP rs115209712 in the entire dataset
af_A1 = data.frame(matrix(NA,nrow=1, ncol=length(c("GBR", "FIN", "IBS", "PEL", "BEB", "YRI", "CHB", "TSI"))))
rownames(af_A1)= "rs115209712"
colnames(af_A1) = c("GBR", "FIN", "IBS", "PEL", "BEB", "YRI", "CHB", "TSI")

for (pop in c("GBR", "FIN", "IBS", "PEL", "BEB", "YRI", "CHB", "TSI")){
  pop.df = geno_pop[geno_pop$POP == pop,]
  n.alleles = 2*nrow(pop.df)
  snp ="rs115209712"
  af_A1[snp, pop] = sum(pop.df[, snp])/n.alleles
}

# A2 frequency of the SNP rs115209712 in the entire dataset
af_A2 = 1 - af_A1

The frequencies of the alleles A1 and A2 of SNP rs115209712 in each populations are:

If you compare the allele frequencies found in our dataset with the ones reported in Ensembl, you can see that they match.

5.2 Genotype frequencies

The netx step will be the computation of the genotype frequencies according to Hardy-Weinberg formulas.

The genotype frequencies under Hardy-Weinberg Equilibrium are:

gf_rs115209712 = data.frame(matrix(NA,nrow=3, ncol=length(c("GBR", "FIN", "IBS", "PEL", "BEB", "YRI", "CHB", "TSI"))))
rownames(gf_rs115209712)= c("GG", "GA", "AA")
colnames(gf_rs115209712) = c("GBR", "FIN", "IBS", "PEL", "BEB", "YRI", "CHB", "TSI")

for (pop in c("GBR", "FIN", "IBS", "PEL", "BEB", "YRI", "CHB", "TSI")){
  gf_rs115209712["GG",pop] = af_A1[pop]*af_A1[pop]
  gf_rs115209712["GA",pop] = 2*af_A1[pop]*af_A2[pop]
  gf_rs115209712["AA",pop] = af_A2[pop]*af_A2[pop]
}

As you can see the snp rs115209712 is not variable across different worldwide populations: in fact the homozygous genotype AA has an allele frequency higher that 0.71 everywhere, while the other two genotypes show little or no variation. From this rapid “visual inspection”, we can conclude that this SNP won’t surely be informative for ancestry inference.

Let’s try with another SNP, this time known in the literature for being informative for ancestry inference. The SNP is rs2814778 and has been reported byt this paper as a genetic marker informative in distinguishing European individuals. Let’s see if it works.

# allele frequencies
af_A1 = data.frame(matrix(NA,nrow=1, ncol=length(c("GBR", "FIN", "IBS", "PEL", "BEB", "YRI", "CHB", "TSI"))))
rownames(af_A1)= "rs2814778"
colnames(af_A1) = c("GBR", "FIN", "IBS", "PEL", "BEB", "YRI", "CHB", "TSI")

for (pop in c("GBR", "FIN", "IBS", "PEL", "BEB", "YRI", "CHB", "TSI")){
  pop.df = geno_pop[geno_pop$POP == pop,]
  n.alleles = 2*nrow(pop.df)
  snp ="rs2814778"
  af_A1[snp, pop] = sum(pop.df[, snp])/n.alleles
}
# A2 frequency of the SNP rs2814778 in the entire dataset
af_A2 = 1 - af_A1

# genotype frequencies
gf_rs2814778 = data.frame(matrix(NA,nrow=3, ncol=length(c("GBR", "FIN", "IBS", "PEL", "BEB", "YRI", "CHB", "TSI"))))
rownames(gf_rs2814778)= c("CC", "CT", "TT")
colnames(gf_rs2814778) = c("GBR", "FIN", "IBS", "PEL", "BEB", "YRI", "CHB", "TSI")

for (pop in c("GBR", "FIN", "IBS", "PEL", "BEB", "YRI", "CHB", "TSI")){
  gf_rs2814778["CC",pop] = af_A1[pop]*af_A1[pop]
  gf_rs2814778["CT",pop] = 2*af_A1[pop]*af_A2[pop]
  gf_rs2814778["TT",pop] = af_A2[pop]*af_A2[pop]
}

AS you can see from the table above and its Ensembl webpage, this SNP is informative in distinguishing Yoruba (African from Nigeria) individuals, because only in the Yoruba population the allele C and the genotype CC have almost reached the fixation.

6.1 Where does NA18520 come from?

unknown_sample = "NA18520"

random_markers = sample(snp.id,30)
snp_set = random_markers

# compute the allele frequency across all populations
af_A1 = data.frame(matrix(NA,nrow=length(snp_set), ncol=length(c("GBR", "FIN", "IBS", "PEL", "BEB", "YRI", "CHB", "TSI"))))
rownames(af_A1)= snp_set
colnames(af_A1) = c("GBR", "FIN", "IBS", "PEL", "BEB", "YRI", "CHB", "TSI")

for (pop in c("GBR", "FIN", "IBS", "PEL", "BEB", "YRI", "CHB", "TSI")){
  pop.df = geno_pop[geno_pop$POP == pop,]
  n.alleles = 2*nrow(pop.df)
  for (snp in snp_set){
    af_A1[snp, pop] = (1 + sum(pop.df[, snp]))/(2 + n.alleles)
  }
}

# compute the genotypic frequencies of the genotypes in the unknown sample
geno_prob = data.frame(matrix(NA,nrow=length(snp_set), ncol=length(c("GBR", "FIN", "IBS", "PEL", "BEB", "YRI", "CHB", "TSI"))))
rownames(geno_prob)= snp_set
colnames(geno_prob) = c("GBR", "FIN", "IBS", "PEL", "BEB", "YRI", "CHB", "TSI")

for (pop in c("GBR", "FIN", "IBS", "PEL", "BEB", "YRI", "CHB", "TSI")){
  #pop.df = geno_pop[geno_pop$POP == pop,]
  #n.alleles = 2*nrow(pop.df)
  for (snp in snp_set){
    p = af_A1[snp,pop]
    q = 1 - p
    if (geno_unk[unknown_sample,snp] == 0){
      geno_prob[snp,pop] = q*q
    } else if (geno_unk[unknown_sample,snp] == 1){
      geno_prob[snp,pop] = 2*p*q
    } else if (geno_unk[unknown_sample,snp] == 2){
      geno_prob[snp,pop] = p*p
    }
  }   
}
RMP = sort(apply(geno_prob, 2, prod), decreasing = T)
LR = RMP[1]/RMP[2]

## [1] "The YRI ancestry is 3627.38281716031 times more probable than BEB ancestry."

## [1] "The YRI ancestry is 134892248993230752 times more probable than BEB ancestry."

6.2 Where does NA18532 come from?

Now, we will compare the classification performances of the random and the Nano set on the individual NA18532.

The sample NA18532 is a Han Chinese individual from Beijing and while the random sets cannot infer his ancestry with high confidence (see the reduced LR with the random SNPs), the Global AIMs Nano set can accurately find his origins.

6.3 Where does HG03585 come from?

Now, we will compare the classification performances of tha random and the Nano set on the individual HG03585.

The sample HG03585 comes from Bangladesh. The Central-South Eurasia has been crossed for millennia by different populations and its populations show intermediate genetic frequencies between East Asia, Europe and Middle East. In the slide 11 of the pdf “Beyond identification (ancestry inference)” you can see that with K 4, 5 and 6 the Central/South Asia has the same colors of Europe (light blue) and East Asia (pink). However, the authors of the Global AIMs Nano set made an effort to distinguish this geographical area.

6.4 Where does NA20502 come from?

Now, we will compare the classification performances of tha random and the Nano set on the individual NA20502.

The sample NA20502 comes from Italy and particularly from Tuscany. None of the marker sets is able to correctly infer his ancestry. A random set of markers could never solve the difficult problem of distinguishing different countries from the same continents (in this case, Europe) and neither the Global AIMs Nano set can do it, because this set was developed with the aim of distinguishing the different continents. In fact, this set has correctly understood that this sample comes from Europe, since all the European populations have around the same probabilities of being the population of origin of NA20502. However, none of them prevails.

A question for you: where does the individual HG01620 come from?