Last updated: 2022-07-17

Checks: 7 0

Knit directory: SISG2022_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(20220530)

The command set.seed(20220530) 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: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: ddfa683

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 ddfa683. 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:    data/.DS_Store

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 ddfa683 Joelle Mbatchou 2022-07-17 fix plink version
html 68b8d58 Joelle Mbatchou 2022-07-17 Build site.
Rmd 4b95c5b Joelle Mbatchou 2022-07-17 session 1-2 praticals

Before you begin:

  • Make sure that R is installed on your computer
  • For this lab, we will use the following R libraries:
require(data.table)
require(dplyr)
require(tidyr)
require(bigsnpr)
require(ggplot2)

Population Structure Inference

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 genotype data at autosomal SNPs 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)

Exercises

Here are some things to look at:

  1. Examine the dataset:
  • How many samples are present?
  • How many SNPs?
  • What is the number of samples in each population?
  1. Get the first 10 principal components (PCs) in PLINK using all SNPs. The basic command would look like
plink2 --bfile <plink_bed_prefix> --pca 10 --out <output_prefix>

This generates a file <output_prefix>.eigenvec containing the PCs (eigenvectors) as well as another file <output_prefix>.eigenval containing the top eigenvalues.

  • Make a scatterplot of the first two PCs with each point colored by population membership.
  • Interpret the first two PCs, what ancestries are they reflecting?
  • Make a scree plot of the eigenvalues for the first 10 PCs. Approximate the proportion of variance explained by the first two PCs.
  1. 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 would look like
# run PCA
obj.bed <- bed(bedfile = <plink_bed_file>)
pc.out <- bed_autoSVD(
  obj.bed, 
  thr.r2 = <r2_threshold>, 
  k = <number_of_PCs>, 
  min.mac = <min_MAC>
)
# plot PC2 vs PC1
plot(pc.out, type = "scores", scores = 1:2)
# scree plot
plot(pc.out) 
# plot SNP loadings (should be centered at 0)
plot(pc.out, type = "loadings", scores = 1:<number_of_PCs>, coeff = 0.4)
  • Run PCA and 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?
  • Check the SNP loadings for the first 10 PCs.

(Hint: This tutorial document from bigsnpr might be helpful)

  1. Predict the proportional Native American and European Ancestry for the HapMap MXL from the PCA output in Question 3 using one of the principal components. (Which PC is most appropriate for this analysis?) Assume that the HapMap MXL have negligible African Ancestry.

  2. Make a barplot of the proportional ancestry estimates from question 4.

Extra: 6. Check if there are samples related 2nd degree or closer. If so, run PCA as in Question 3 removing these samples then project the remaining samples onto the PC space. The basic command would look like

# check for 3rd degree relateds or closer
snp_plinkKINGQC(
  plink2.path = "/usr/bin/plink2", 
  bedfile.in = <plink_bed_prefix>, 
  thr.king = 2^-3.5,
  make.bed = FALSE
)

(Hint: This returns a data frame which contains all pairs of individuals related 3rd degree or closer. We can then remove them when calling bed_autoSVD() using the ind.row argument. Finally, you can use bed_projectSelfPCA() to project related samples on the PC space.)


sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Catalina 10.15.7

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

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

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

other attached packages:
[1] workflowr_1.7.0

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.8.3     bslib_0.3.1      jquerylib_0.1.4  compiler_3.6.1  
 [5] pillar_1.7.0     later_1.3.0      git2r_0.30.1     tools_3.6.1     
 [9] getPass_0.2-2    digest_0.6.25    jsonlite_1.7.2   evaluate_0.15   
[13] tibble_2.1.3     lifecycle_1.0.1  pkgconfig_2.0.3  rlang_1.0.1     
[17] cli_3.1.1        rstudioapi_0.13  yaml_2.2.1       xfun_0.31       
[21] fastmap_1.1.0    httr_1.4.3       stringr_1.4.0    knitr_1.39      
[25] sass_0.4.0       fs_1.5.2         vctrs_0.3.8      rprojroot_2.0.3 
[29] glue_1.6.1       R6_2.4.1         processx_3.5.3   fansi_0.4.1     
[33] rmarkdown_2.14   callr_3.7.0      magrittr_1.5     whisker_0.4     
[37] ps_1.7.0         promises_1.2.0.1 htmltools_0.5.2  ellipsis_0.3.2  
[41] httpuv_1.6.5     utf8_1.1.4       stringi_1.4.6    crayon_1.3.4