Practical 2 - Population Structure Inference

Last updated: 2024-06-13

Checks: 6 1

Knit directory: SISG2024_Association_Mapping/

This reproducible R Markdown analysis was created with workflowr (version 1.7.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20230530)

The command set.seed(20230530) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: absolute

Using absolute paths to the files within your workflowr project makes it difficult for you and others to run your code on a different machine. Change the absolute path(s) below to the suggested relative path(s) to make your code more reproducible.

absolute	relative
/Users/joelle.mbatchou/SISG/2024/SISG2024_Association_Mapping/data/	data

Repository version: 5558558

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 5558558. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .DS_Store
    Ignored:    analysis/.DS_Store
    Ignored:    data/sim_rels_geno.bed
    Ignored:    exe/
    Ignored:    gwas_plink.log
    Ignored:    gwas_regenie.log
    Ignored:    lectures/
    Ignored:    mk_website.R
    Ignored:    step2_gwas_regenie.log
    Ignored:    tmp/

Untracked files:
    Untracked:  .Rhistory
    Untracked:  analysis/SISGM15_prac4Solution.Rmd
    Untracked:  analysis/SISGM15_prac5Solution.Rmd
    Untracked:  analysis/SISGM15_prac6Solution.Rmd
    Untracked:  analysis/SISGM15_prac9Solution.Rmd
    Untracked:  analysis/Session02_practical_Key_cache/
    Untracked:  analysis/Session03_practical_Key.Rmd
    Untracked:  analysis/Session07_practical.Rmd
    Untracked:  analysis/Session07_practical_Key.Rmd
    Untracked:  analysis/Session08_practical.Rmd
    Untracked:  analysis/Session08_practical_Key.Rmd
    Untracked:  data/run_regenie.r
    Untracked:  notes.txt

Unstaged changes:
    Modified:   .gitignore

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/Session02_practical.Rmd) and HTML (docs/Session02_practical.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
Rmd	81864cc	Joelle Mbatchou	2024-06-12	update practicals
html	81864cc	Joelle Mbatchou	2024-06-12	update practicals
Rmd	6e92dae	Joelle Mbatchou	2024-06-10	update practical session 2
html	6e92dae	Joelle Mbatchou	2024-06-10	update practical session 2
Rmd	364d0cf	Joelle Mbatchou	2024-06-10	update practical session 2
html	364d0cf	Joelle Mbatchou	2024-06-10	update practical session 2
html	3499240	Joelle Mbatchou	2024-06-05	update session pages

Before you begin:

Make sure that R is installed on your computer
For this lab, we will use the following R libraries:

library(data.table)
library(dplyr)
library(bigsnpr)
library(ggplot2)

Introduction

We will be working with a subset of the genotype data from the Human Genome Diversity Panel (HGDP) and HapMap.

The file “YRI_CEU_ASW_MEX_NAM.bed” is a binary file in PLINK BED format with accompanying BIM and FAM files. It contains the genotype data at autosomal SNPs (i.e. chromosomes 1-22) for:

Native American samples from HGDP
Four population samples from HapMap:
- Yoruba in Ibadan, Nigeria (YRI)
- Utah residents with ancestry from Northern and Western Europe (CEU)
- Mexican Americans in Los Angeles, California (MXL)
- African Americans from the south-western United States (ASW)

File with ancestry labels assignment for each sample: Population_Sample_Info.txt

Data preparation

Let’s first load the HGDP data into the R session. We need to define the path to the directory containing the PLINK BED and the ancestry label files (change the path to the file location).

# change this to the directory on your machine
HGDP_dir <- "/SISGM19/data/"

Also specify the path to the PLINK2 binary

plink2_binary <- "/SISGM19/bin/plink2"

We can now read the PLINK BED and FAM files (recall the BED file is a binary file):

HGDP_bim <- fread(sprintf("%s/YRI_CEU_ASW_MEX_NAM.bim", HGDP_dir), header = FALSE)
head(HGDP_bim, 3)

   V1        V2 V3      V4 V5 V6
1:  1 rs9442372  0 1008567  1  2
2:  1 rs2887286  0 1145994  1  2
3:  1 rs3813199  0 1148140  1  2

HGDP_fam <- fread(sprintf("%s/YRI_CEU_ASW_MEX_NAM.fam", HGDP_dir), header = FALSE)
head(HGDP_fam, 3)

     V1        V2 V3 V4 V5 V6
1: 1432 HGDP00702  0  0  2 -9
2: 1433 HGDP00703  0  0  1 -9
3: 1434 HGDP00704  0  0  2 -9

When reading the ancestry label file, we need to make sure the order of samples matches that in the PLINK data:

HGDP_ancestry_df <- fread(sprintf("%s/Population_Sample_Info.txt", HGDP_dir))
HGDP_ancestry_df <- left_join(HGDP_fam[,c("V1","V2")], HGDP_ancestry_df, by = c("V1" = "FID", "V2" = "IID"))
head(HGDP_ancestry_df, 3)

     V1        V2 Population
1: 1432 HGDP00702        NAM
2: 1433 HGDP00703        NAM
3: 1434 HGDP00704        NAM

Exercises

Here are some things to look at:

Examine the dataset:
- How many samples are present? str(HGDP_fam)
- How many SNPs? str(HGDP_bim)
- What is the number of samples in each population? table(HGDP_ancestry_df$Population)
Get the first 10 principal components (PCs) in PLINK using all SNPs.

cmd <- sprintf("%s --bfile %s/YRI_CEU_ASW_MEX_NAM --pca 10 --out pca_plink", plink2_binary, HGDP_dir)
system(cmd)

This generates two files pca_plink.eigenvec containing the PCs (eigenvectors), and pca_plink.eigenval containing the top eigenvalues.

Read in the PCs in R and make a scatterplot of the first two PCs with each point colored by population membership.

PC_df <- left_join(HGDP_ancestry_df, fread("pca_plink.eigenvec"), by = c("V1" = "#FID", "V2" = "IID"))
ggplot(PC_df, aes(x=PC1, y=PC2, color = Population)) +
  geom_point()

Interpret the first two PCs, what ancestries are they reflecting?
Read in the eigenvalues and make a scree plot corresponding for these first 10 PCs. Estimate the proportion of variance explained by the first two PCs.

eigenvalues_df <- fread("pca_plink.eigenval", header = FALSE)
# Make scree plot
ggplot(eigenvalues_df, aes(x = 1:10, y = V1)) +
  geom_point() +
  geom_line() +
  scale_x_continuous(breaks = 1:10) +
  labs(x = "PC", y = "Eigenvalue")

Now redo Question 2 above using the bigsnpr R package specifying a $r^2$ threshold of 0.2 (i.e. LD pruning) as well as a minimum minor allele count (MAC) of 20. The basic command is

obj.bed <- bed(bedfile = sprintf("%s/YRI_CEU_ASW_MEX_NAM.bed", HGDP_dir))
pca.bigsnpr <- bed_autoSVD(
  obj.bed, 
  thr.r2 = 0.1, # R^2 threshold
  k = 10, # number of PCs
  min.mac = 10 # minimum minor allele count (MAC) filter
)

You can evaluate the PCA results using

# plot PC2 vs PC1
plot(pca.bigsnpr, type = "scores", scores = 1:2)
# scree plot
plot(pca.bigsnpr) 
# plot SNP loadings for the first two PCs
plot(pca.bigsnpr, type = "loadings", loadings = 1:2, coeff = 0.4)

Make a scatter plot of the first two principal components (PCs) with each point colored according to population membership. Does the plot change from the one in Question 2?

plot(pca.bigsnpr, type = "scores", scores = 1:2) +
  aes(color = HGDP_ancestry_df$Population) +
  labs(color = "Population")

Check the SNP loadings for the first 5 PCs.

Predict the proportion of Native American ancestry for the HapMap MXL based on the PCA output in Question 3 using one of the principal components and assuming that the HapMap MXL have negligible African Ancestry. Which PC is most appropriate for this analysis?

Hint: To compute the average PC2 value for individuals of CEU ancestry:

ceu.mean <- with(PC_df, mean(PC2[Population == "CEU"]))

Do the same for individuals of NAM ancestry. How can you express distances of the MXL individuals relative to those means based on the chosen PC?

Make a barplot of the proportional Native American/European ancestry estimates based on your output of question 7.

# this line assumes the proportion of NAM ancestry is stored in a variable called "mxl.prop.nam"
sorted.props.nam <- sort(mxl.prop.nam, decreasing = TRUE)

df_mxl_props <- data.frame(
  anc.props = c(sorted.props.nam, 1 - sorted.props.nam),
  x = rep(1:length(sorted.props.nam), times = 2),
  population.labels = rep(c("NAM", "CEU"), each = length(sorted.props.nam))
)
ggplot(df_mxl_props, aes(x = x, y = anc.props, fill = population.labels)) +
  geom_bar(position="stack", stat="identity") +
  labs(x="Sample", y = "Ancestry Proportion", fill = "Population")

Extra

Check if there are samples related 2nd degree or closer. If so, run PCA as in Question 6 removing these samples then project the remaining samples onto the PC space.

The command to check relatedness is:

# check for 2nd degree relateds or closer
relatedness_info <- snp_plinkKINGQC(
  plink2.path = plink2_binary, 
  bedfile.in = sprintf("%s/YRI_CEU_ASW_MEX_NAM.bed", HGDP_dir), 
  thr.king = 2^-3.5, # threshold to identify 2nd degree relateds
  make.bed = FALSE
)

This returns a data frame which contains all pairs of individuals related 2nd degree or closer.

We can then remove them when calling bed_autoSVD() using the ind.row argument.

ind.rel <- match(c(relatedness_info$IID1, relatedness_info$IID2), obj.bed$fam$sample.ID) # relateds
ind.norel <- rows_along(obj.bed)[-ind.rel] # unrelateds
# Run PCA on unrelateds
obj.svd2 <- bed_autoSVD(
  obj.bed, 
  thr.r2 = 0.1, # R^2 threshold
  k = 10, # number of PCs
  min.mac = 10, # minimum minor allele count (MAC) filter
  ind.row = ind.norel
)

Use bed_projectSelfPCA() to project related samples on the PC space. (Hint: This tutorial document from bigsnpr will be helpful – see the last section ‘Project remaining individuals’)

sessionInfo()

R version 4.3.0 (2023-04-21)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS 14.5

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] ggplot2_3.4.2     bigsnpr_1.12.9    bigstatsr_1.5.12  dplyr_1.1.2      
[5] data.table_1.14.8

loaded via a namespace (and not attached):
 [1] sass_0.4.6         utf8_1.2.3         generics_0.1.3     bigsparser_0.6.1  
 [5] lattice_0.21-8     stringi_1.7.12     digest_0.6.31      magrittr_2.0.3    
 [9] evaluate_0.21      grid_4.3.0         iterators_1.0.14   fastmap_1.1.1     
[13] Matrix_1.5-4       doParallel_1.0.17  foreach_1.5.2      rprojroot_2.0.3   
[17] workflowr_1.7.0    jsonlite_1.8.5     whisker_0.4.1      promises_1.2.0.1  
[21] doRNG_1.8.6        fansi_1.0.4        scales_1.2.1       codetools_0.2-19  
[25] jquerylib_0.1.4    cli_3.6.1          rlang_1.1.1        cowplot_1.1.1     
[29] munsell_0.5.0      withr_2.5.0        cachem_1.0.8       yaml_2.3.7        
[33] flock_0.7          parallel_4.3.0     tools_4.3.0        bigassertr_0.1.6  
[37] colorspace_2.1-0   httpuv_1.6.11      rngtools_1.5.2     vctrs_0.6.2       
[41] R6_2.5.1           lifecycle_1.0.3    bigparallelr_0.3.2 git2r_0.32.0      
[45] stringr_1.5.0      fs_1.6.2           pkgconfig_2.0.3    pillar_1.9.0      
[49] bslib_0.5.0        later_1.3.1        gtable_0.3.3       glue_1.6.2        
[53] Rcpp_1.0.10        xfun_0.39          tibble_3.2.1       tidyselect_1.2.0  
[57] highr_0.10         rstudioapi_0.14    knitr_1.43         htmltools_0.5.5   
[61] rmarkdown_2.22     compiler_4.3.0