Last updated: 2024-06-10
Checks: 5 2
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.
The R Markdown is untracked by Git. To know which version of the R
Markdown file created these results, you’ll want to first commit it to
the Git repo. If you’re still working on the analysis, you can ignore
this warning. When you’re finished, you can run
wflow_publish to commit the R Markdown file and build the
HTML.
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.
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.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
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 |
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 a58f9a0. 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: analysis/.DS_Store
Ignored: data/sim_rels_geno.bed
Ignored: exe/
Ignored: lectures/
Ignored: mk_website.R
Ignored: tmp/
Untracked files:
Untracked: .Rhistory
Untracked: analysis/SISGM15_prac4.Rmd
Untracked: analysis/SISGM15_prac4Solution.Rmd
Untracked: analysis/SISGM15_prac5.Rmd
Untracked: analysis/SISGM15_prac5Solution.Rmd
Untracked: analysis/SISGM15_prac6.Rmd
Untracked: analysis/SISGM15_prac6Solution.Rmd
Untracked: analysis/SISGM15_prac9.Rmd
Untracked: analysis/SISGM15_prac9Solution.Rmd
Untracked: analysis/Session01_practical_Key.Rmd
Untracked: analysis/Session02_practical.Rmd
Untracked: analysis/Session02_practical_Key.Rmd
Untracked: analysis/Session03_practical.Rmd
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: code/Session01_exercises.R
Untracked: code/Session02_exercises.R
Untracked: code/Session03_exercises.R
Untracked: code/Session07_exercises.R
Untracked: code/Session08_exercises.R
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 |
|---|---|---|---|---|
| html | 3499240 | Joelle Mbatchou | 2024-06-05 | update session pages |
Before you begin:
library(data.table)
library(dplyr)
library(bigsnpr)
library(ggplot2)The R template to do the exercises is here.
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:
File with ancestry labels assignment for each sample: Population_Sample_Info.txt
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 2HGDP_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 -9When 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 NAMHere are some things to look at:
str(HGDP_fam)str(HGDP_bim)table(HGDP_ancestry_df$Population)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.
PC_df <- left_join(fam_pop_info, fread("pca_plink.eigenvec"), by = c("V1" = "#FID", "V2" = "IID"))
with(PC_df, plot(x=PC1, y=PC2, color = Population))Interpret the first two PCs, what ancestries are they reflecting?
Read in the eigenvalues and a scree plot corresponding to the first 10 PCs. Estimate the proportion of variance explained by the first two PCs.
eigenvalues_df <- fread("pca_plink.eigenval", header = FALSE)
with(eigenvalues_df, plot(x=1:10, y=V1, xlab = "PC", ylab = "Eigenvalue"))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 likeobj.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
)# plot PC2 vs PC1
plot(pca.bigsnpr, type = "scores", scores = 1:2)
# scree plot
plot(pca.bigsnpr)
# plot SNP loadings for the first 5 PCs
plot(pca.bigsnpr, type = "loadings", scores = 1:5, coeff = 0.4)plot(pca.bigsnpr, type = "scores", scores = 1:2) +
aes(color = fam_pop_info$Population) +
labs(color = "Population")Predict the proportional Native American and European Ancestry for the HapMap MXL from the PCA output in Question 6 using one of the principal components. (Which PC is most appropriate for this analysis?) Assume that the HapMap MXL have negligible African Ancestry.
Make a barplot of the proportional ancestry estimates from question 4.
# check for 3rd 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 for 3rd degree relateds
make.bed = FALSE
)This returns a data frame which contains all pairs of individuals related 3rd degree or closer.
bed_autoSVD()
using the ind.row argument.ind.rel <- match(c(relatedness_info$IID1, relatedness_info$IID2), obj.bed$fam$sample.ID)
ind.norel <- rows_along(obj.bed)[-ind.rel]
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
)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/Chicago
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.2 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