{
"cells": [
{
"cell_type": "markdown",
"id": "876323d6-5e7b-40c0-9597-5f8bedff28c7",
"metadata": {
"tags": []
},
"source": [
"# Population structure\n",
"\n",
"With the advent of SNP data it is possible to precisely infer the genetic distance across individuals or populations. As written in the book, one way of doing it is by comparing each SNP from each individual against every other individual. This comparison produces the so called: covariance matrix, which in genetic terms means the number of shared polymorphisms across individuals.\n",
"There are many ways to visualize this data, in this tutorial you will be exposed to `Principal Component Analysis` and `Admixture` software.\n",
"\n",
"We will use the R package `SNPRelate`, which can easily handle vcf files and do the PCA. If you want to explore a bit more on the functionality of the package access [here](https://www.rdocumentation.org/packages/SNPRelate/versions/1.6.4).\n",
"\n",
"## How to make this notebook work\n",
"\n",
"* In this notebook we will use both the `command line bash` commands and `R` to setup the file folders.\n",
"* Having to shift between two languages, you need to choose a kernel every time we shift from one language to another. A kernel contains a programming language and the necessary packages to run the course material. To choose a kernel, go on the menu on the top of the page and select `Kernel --> Change Kernel`, and then select the preferred one. We will shift between two kernels, and along the code in this notebook you will see a picture telling you to change kernel. The two pictures are below:\n",
"\n",
"<img src=\"Images/bash.png\" alt=\"Bash\" width=\"80\"> Shift to the `Bash` kernel\n",
"\n",
"<img src=\"Images/R.png\" alt=\"R\" width=\"80\"> Shift to the `R` kernel\n",
"\n",
"\n",
"## Learning outcomes\n",
"\n",
"At the end of this tutorial you will be able to\n",
"\n",
"- **Extract** information from a `vcf` file and create a PCA projection\n",
"- **Look at the effect of** LD pruning to reveal population structure\n",
"\n",
"## Setting up folders\n",
"\n",
"Here we setup a link to the `Data` folder and create the `Results` folder.\n",
"\n",
"<img src=\"Images/bash.png\" alt=\"Bash\" width=\"80\"> Choose the `Bash` kernel"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c31784ba-e790-455c-b772-4c12dd43fdc8",
"metadata": {
"tags": [],
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"ln -sf ../Data"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ef9e23c2-307e-4f43-ac2a-3e68025acf61",
"metadata": {
"tags": [],
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"unzip Data/hapmap.zip -d Data"
]
},
{
"cell_type": "markdown",
"id": "2c078d72-c9fd-4045-9c0c-373dabf0da8a",
"metadata": {},
"source": [
"# PCA from VCF data file"
]
},
{
"cell_type": "markdown",
"id": "4e0d238b-d21e-41e3-a23a-9ee3f5b8b4e8",
"metadata": {},
"source": [
"<img src=\"Images/R.png\" alt=\"R\" width=\"80\"> Shift to the `R` kernel"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d724bcd3-f6f2-42a7-863f-a6c43946f575",
"metadata": {
"tags": [],
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"library(SNPRelate)\n",
"library(ggplot2)"
]
},
{
"cell_type": "markdown",
"id": "100a9384-b0e2-4af5-9fc2-55581b0e2266",
"metadata": {},
"source": [
"## Import data and calculate PCA\n",
"\n",
"We read the metadata about the samples (geographic locations) and transform the `vcf` file into `gds` format, from which the package `SNPRelate` can calculate the PCA projection of the data."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1d7bd102-f4e9-4f4a-b9bd-0211d3210a22",
"metadata": {
"tags": [],
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"# Reading the metadata information \n",
"info = read.csv(\"Data/sample_infos_accessionnb.csv\", header = T, sep = ';')\n",
"\n",
"# Setting the directory of the VCF file \n",
"vcf.fn <- \"Data/Allvariants_135_145_chr2.vcf\"\n",
"\n",
"# Transforming the vcf file to gds format\n",
"snpgdsVCF2GDS(vcf.fn, \"Data/Allvariants_135_145_chr2_2.gds\", method=\"biallelic.only\")\n",
"\n",
"#Read the file and calculate the PCA\n",
"genofile <- snpgdsOpen(\"Data/Allvariants_135_145_chr2_2.gds\", FALSE, TRUE, TRUE)\n",
"pca <- snpgdsPCA(genofile)\n",
"summary(pca)"
]
},
{
"cell_type": "markdown",
"id": "1bfe101e-8197-4b77-bba9-b64682655cc1",
"metadata": {},
"source": [
"**Q.1** How many individuals and snps does this dataset have? What is an eigenvector and an eigenvalue?"
]
},
{
"cell_type": "markdown",
"id": "0f033932-59e8-4a6b-904c-71e295bc4e00",
"metadata": {},
"source": [
"The `pca` object just created is a list containing various elements."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "24333347-934f-4331-bddc-e0fc3bdc4b0d",
"metadata": {
"tags": [],
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"ls(pca)"
]
},
{
"cell_type": "markdown",
"id": "651c485f-39d7-4f99-9d70-8b2235da3790",
"metadata": {},
"source": [
"We use `pca$eigenvect` to plot the PCA. We extract also `pca$sample.id` to match the geographic locations in the metadata with the samples in `pca`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6efb680e-9b95-4fe8-9e6d-89d2b76f0c03",
"metadata": {
"tags": [],
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"eigenvectors = as.data.frame(pca$eigenvect[,1:5])\n",
"colnames(eigenvectors) = as.vector(sprintf(\"PC%s\", seq(1:ncol(eigenvectors))))\n",
"pca$sample.id = sub(\"_chr2_piece_dedup\", \"\", pca$sample.id)\n",
"\n",
"# Matching the sample names with their origin and population\n",
"rownames(info) <- info[,\"ENA.RUN\"]\n",
"eigenvectors <- cbind(eigenvectors, info[pca$sample.id, c(\"population\",\"region\")])\n"
]
},
{
"cell_type": "markdown",
"id": "38b4baf6-bb56-4a09-9250-10dbe707ee06",
"metadata": {},
"source": [
"In the end, we have created a table called `eigenvectors` containing the PCA coordinates and some metadata"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a83b5b02-53cd-4e48-b0c0-57d288c8f74f",
"metadata": {
"tags": [],
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"head(eigenvectors)"
]
},
{
"cell_type": "markdown",
"id": "439cea05-4081-44dd-b8cb-c02fbb68acb7",
"metadata": {},
"source": [
"Let's first look at how much of the variance of the data is explained by each eigenvector:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e5b69906-0939-4610-a541-ca204fe51a3a",
"metadata": {
"tags": [],
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"# Variance proportion:\n",
"pca$pca_percent <- pca$varprop*100\n",
"ggplot( NULL, aes(x=seq(1, length(pca$eigenval)), y=pca$pca_percent, label=sprintf(\"%0.2f\", round(pca$pca_percent, digits = 2))) ) +\n",
" geom_line() + geom_point() + \n",
" geom_text(nudge_y = .3, nudge_x = 1.5, check_overlap = T)"
]
},
{
"cell_type": "markdown",
"id": "b51bda09-95f8-473c-83fd-57cdaed85301",
"metadata": {},
"source": [
"**Q.2** How many PC's do we need in order to explain 50% of the variance of the data? Can you make an accumulative plot of the variance explained PC?"
]
},
{
"cell_type": "markdown",
"id": "546fb504-4a4d-4f66-9e5a-06440baaacfd",
"metadata": {},
"source": [
"## Visualization\n",
"\n",
"We plot now the first two PCA coordinates and label them by Population, with color by region. We can see how only africans are separated from the rest, but the PCA is quite confused and cannot distinguish EastAsia and WestEurasia."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5aff2315-9eae-4d81-b023-ad0be1d9b5a1",
"metadata": {
"tags": [],
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"ggplot(data = eigenvectors, aes(x = PC1, y = PC2, col = region)) + \n",
" geom_point(size=3,alpha=0.5) + geom_text( aes(label=population), col=\"black\") +\n",
" scale_color_manual(values = c(\"#FF1BB3\",\"#A7FF5B\",\"#99554D\")) +\n",
" theme_bw()"
]
},
{
"cell_type": "markdown",
"id": "4028f5df-3122-4195-a8b8-ad53f591455f",
"metadata": {},
"source": [
"**Q.3** Try to plot PC2 and PC3. Do you see the same patterns? What is the correlation between PC2 and PC3 (hint use the function cor())?\n",
"\n",
"**Q.4** Try also to color the graph based on population. What do you observe?"
]
},
{
"cell_type": "markdown",
"id": "79a54b0d-e337-4e6e-aa58-13d4402aaea4",
"metadata": {},
"source": [
"## LD pruning"
]
},
{
"cell_type": "markdown",
"id": "2c760205-a58a-4c10-8481-07567744fd1f",
"metadata": {},
"source": [
"We implement LD pruning to eliminate those SNPs that are in high linkage disequilibrium. What happens is that the structure of the populations will change. According to ([Bergovich et al, 2024, Biorxiv](https://www.biorxiv.org/content/10.1101/2024.05.02.592187v1)), this is not good practice and removes a lot of the population information."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6c3b9585-034e-4d31-9e09-29dfe9ed5848",
"metadata": {
"tags": [],
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"set.seed(1000)\n",
"\n",
"# This function prune the snps with a thrshold of maximum 0.3 of LD\n",
"snpset <- snpgdsLDpruning(genofile, ld.threshold=0.3)\n",
"\n",
"# Get all selected snp's ids\n",
"snpset.id <- unlist(snpset)\n",
"\n",
"pca_pruned <- snpgdsPCA(genofile, snp.id=snpset.id, num.thread=2)\n",
"\n",
"#add metadata\n",
"eigenvectors = as.data.frame(pca_pruned$eigenvect)\n",
"colnames(eigenvectors) = as.vector(sprintf(\"PC%s\", seq(1:nrow(pca$eigenvect))))\n",
"pca_pruned$sample.id = sub(\"_chr2_piece_dedup\", \"\", pca$sample.id)\n",
"eigenvectors <- cbind(eigenvectors, info[pca$sample.id, c(\"population\",\"region\")])\n",
"\n",
"#plot\n",
"ggplot(data = eigenvectors, aes(x = PC3, y = PC2, col = region, label=population)) + \n",
" geom_text(hjust=1, vjust=0, angle=45) +\n",
" geom_point(size=3,alpha=0.5) +\n",
" scale_color_manual(values = c(\"#FF1BB3\",\"#A7FF5B\",\"#99554D\")) +\n",
" theme_bw() + coord_flip()"
]
},
{
"cell_type": "markdown",
"id": "cb3b87ec-184f-4fa5-a8bb-d6b93aab775e",
"metadata": {},
"source": [
"**Q.5** Implement different LD thresholds (0.1, 0.2, 0.3, 0.4, 0.5). How many SNPs are left after each filtering threshold? Are these SNPs linked?"
]
},
{
"cell_type": "markdown",
"id": "5f6556c8-3e0d-4bdc-b0d6-ba3f24985059",
"metadata": {},
"source": [
"Now we are going to convert this GDS file into a plink format, to be later used in the admixture exercise:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cb187178-546a-4b9a-96fc-02d38fe8b870",
"metadata": {
"tags": [],
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"snpgdsGDS2BED(genofile, \"Data/chr2_135_145_flt_prunned.gds\", sample.id=NULL, snp.id=snpset.id, snpfirstdim=NULL, verbose=TRUE)"
]
},
{
"cell_type": "markdown",
"id": "e349bc4a-c1c8-4d2c-902b-fe73e24bbfbb",
"metadata": {},
"source": [
"Save the data for later"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3378d6c7-266d-4cc2-8e46-706649887d5c",
"metadata": {
"tags": [],
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"save(pca, pca_pruned, info, genofile, file = \"Results/data.Rdata\")"
]
},
{
"cell_type": "markdown",
"id": "87d54719-b44f-433c-b149-48d5685bd6b6",
"metadata": {
"tags": []
},
"source": [
"# Admixture Estimation\n",
"\n",
"<img src=\"Images/bash.png\" alt=\"Bash\" width=\"80\"> Choose the `Bash` kernel\n",
"\n",
"`Admixture` is a program for estimating ancestry in a model based manner from SNP genotype datasets, where individuals are unrelated. The input format required by the software is in binary PLINK (`.bed`) file. That is why we converted our vcf file into `.bed`.\n",
"\n",
"Now with adjusted format and pruned snps, we are ready to run the admixture analysis. We believe that our individuals are derived from three ancestral populations:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "586b3d34-77e0-4d69-9929-69b5f58a0893",
"metadata": {
"tags": [],
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"admixture Data/chr2_135_145_flt_prunned.gds.bed 3\n",
"#move output files with other results\n",
"mv chr2_135_145_flt_prunned.gds.3.P chr2_135_145_flt_prunned.gds.3.Q Results/"
]
},
{
"cell_type": "markdown",
"id": "477da2b3-53e3-4a76-acc7-90efc465c251",
"metadata": {},
"source": [
"**Q.6** Have a look at the Fst across populations, that is printed in the terminal. Would you guess which populations are Pop0, Pop1 and Pop2 referring to?\n",
"\n",
"After running admixture, 2 outputs are generated:\n",
"\n",
"- `Q`: the ancestry fractions\n",
"\n",
"- `P`: the allele frequencies of the inferred ancestral populations\n",
"\n",
"Sometimes we may have no priori about `K`, one good way of choosing the best `K` is by doing a cross-validation procedure impletemented in admixture as follow:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "73273dc7-80c9-4a65-8657-ef3fadef0513",
"metadata": {
"scrolled": true,
"tags": [],
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"for K in 1 2 3 4 5\n",
"do \n",
" admixture --cv Data/chr2_135_145_flt_prunned.gds.bed $K | tee log${K}.out\n",
" mv chr2_135_145_flt_prunned.gds.$K.* log$K.out Results/\n",
"done"
]
},
{
"cell_type": "markdown",
"id": "7d4163b1-4da9-41a8-a9b4-6e28a1331a01",
"metadata": {},
"source": [
"Have a look at the Cross Validation error of each `K`:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4741e4db-e6ce-4496-9034-f2cce303dace",
"metadata": {
"tags": [],
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"grep -h CV Results/log*.out"
]
},
{
"cell_type": "markdown",
"id": "ff6a401e-8f6e-45a0-8fcd-3fd79b3ace78",
"metadata": {},
"source": [
"Save it in a text file:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "10933cd3-cab7-4f41-8454-30c745276c71",
"metadata": {
"tags": [],
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"grep -h CV Results/log*.out > Results/CV_logs.txt"
]
},
{
"cell_type": "markdown",
"id": "ba39de62-ae36-4451-8ed6-00d5c7fa7557",
"metadata": {},
"source": [
"<img src=\"Images/R.png\" alt=\"R\" width=\"80\"> Shift to the `R` kernel\n",
"\n",
"Look at the distribution of CV error."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d018c6d2-ab0d-4398-ad55-019b8ef8e7db",
"metadata": {
"tags": [],
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"library(ggplot2)\n",
"\n",
"CV = read.table('Results/CV_logs.txt')\n",
"\n",
"p <- ggplot(data = CV, aes(x = V3, y = V4, group = 1)) + \n",
" geom_line() + geom_point() + theme_bw() + \n",
" labs(x = 'Number of clusters', y = 'Cross validation error')\n",
"\n",
"p"
]
},
{
"cell_type": "markdown",
"id": "73a0c22d-8489-4739-88d0-ef56bfa82b4f",
"metadata": {},
"source": [
"**Q.7** What do you understand of Cross validation error? Based on this graph, what is the best `K`?\n",
"\n",
"Plotting the `Q` estimates. Choose the `K` that makes more sense to you and substitute it in the first line of code (right now it is `K=3`)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "09b21282-b5a3-469d-9e5f-6a6a8944dc7e",
"metadata": {
"tags": [],
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"tbl = read.table(\"Results/chr2_135_145_flt_prunned.gds.3.Q\")\n",
"ord = tbl[order(tbl$V1,tbl$V2,tbl$V3),]\n",
"bp = barplot(t(as.matrix(ord)), legend.text = c(\"African\", \"2\", \"3\"),\n",
" space = c(0.2),\n",
" col=rainbow(3),\n",
" xlab=\"Population #\", \n",
" ylab=\"Ancestry\",\n",
" border=NA,\n",
" las=2)"
]
},
{
"cell_type": "markdown",
"id": "ec624904-90d3-4809-90cf-2ac391f57e60",
"metadata": {},
"source": [
"Note: Here we order the X-axis based on proportions for the first population component. However, you will see that in the HapMap data all the individuals show some portion of this component and the different individuals are more admixed in general, i.e they are no longer explained by mostly one component, it’s not useful to use that kind of ordering anymore to interpret the plot. Instead, we should keep the original order, since the files are originally ordered by population, and we should plot each population on the X axis to be able to interpret the plot. This can be achieved with something of the type:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "84b6ac8a-01eb-4d85-994e-1b22511cf0f4",
"metadata": {
"tags": [],
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"library(dplyr)\n",
"library(ggplot2)\n",
"load(\"Results/data.Rdata\")\n",
"\n",
"#resize plot\n",
"options(repr.plot.width = 12, repr.plot.height = 8)\n",
"\n",
"K=3\n",
"\n",
"tbl = read.table( paste(\"Results/chr2_135_145_flt_prunned.gds.\",K,\".Q\", sep=\"\") )\n",
"\n",
"origin <- rep( paste( rep(\"Pop\", K), as.character(c(0:(K-1))), sep=\"\" ) , each=dim(tbl)[1] )\n",
"population <- info[ pca$sample.id, ]$population\n",
"region <- info[ pca$sample.id, ]$region\n",
"regpop <- make.unique( paste(region, population) )\n",
"\n",
"tbl <- as.data.frame(unlist(tbl))\n",
"colnames(tbl) <- 'Admixture_fraction'\n",
"tbl['origin'] = origin\n",
"tbl['population'] = rep(population,K)\n",
"tbl['region'] = region\n",
"tbl['region_population'] = regpop\n",
"\n",
"ggplot(tbl, aes(fill=origin, y=Admixture_fraction, x=region_population)) +\n",
" geom_bar(position=\"stack\", stat=\"identity\") + \n",
" theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))"
]
},
{
"cell_type": "markdown",
"id": "c4d44c92-f57b-4512-96eb-e4a72c922c19",
"metadata": {
"tags": []
},
"source": [
"**Q.8** How many clusters do you identify in this plot? Does that agree with what was found using PCA?\n",
"\n",
"In the following part of this exercise you will do both analysis (PCA and Admixture) using a different dataset. The data comes from the HAPMAP Consortium, to learn more about the populations studied in this project access here. The vcf file `hapmap.vcf`, an information file `relationships_w_pops_121708.txt`, as well as `.bim, .bed, .fam` files (only to be used if you get stuck during the exercise) are available for the admixture analysis. This dataset is placed here:\n",
"\n",
"`Data/assignment`"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "R",
"language": "R",
"name": "ir"
},
"language_info": {
"codemirror_mode": "r",
"file_extension": ".r",
"mimetype": "text/x-r-source",
"name": "R",
"pygments_lexer": "r",
"version": "4.2.2"
}
},
"nbformat": 4,
"nbformat_minor": 5
}