Last updated: 2023-08-10
Checks: 7 0
Knit directory: fitnessGWAS/
This reproducible R Markdown analysis was created with workflowr (version 1.7.0.4). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
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.
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(20180914) 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.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
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 6b63c18. 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: .Rapp.history
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: .httr-oauth
Ignored: .pversion
Ignored: analysis/.DS_Store
Ignored: code/.DS_Store
Ignored: code/Drosophila_GWAS.Rmd
Ignored: data/.DS_Store
Ignored: data/derived/
Ignored: data/input/.DS_Store
Ignored: data/input/.pversion
Ignored: data/input/dgrp.fb557.annot.txt
Ignored: data/input/dgrp2.bed
Ignored: data/input/dgrp2.bim
Ignored: data/input/dgrp2.fam
Ignored: data/input/huang_transcriptome/
Ignored: figures/.DS_Store
Untracked files:
Untracked: JAGS-4.3.0/
Untracked: biv_mod.jags
Untracked: old_analyses/
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/gwas_adaptive_shrinkage.Rmd) and HTML
(docs/gwas_adaptive_shrinkage.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 | 6b63c18 | lukeholman | 2023-08-10 | wflow_publish("analysis/gwas_adaptive_shrinkage.Rmd") |
| html | 16ba947 | lukeholman | 2023-07-31 | Build site. |
| Rmd | 926d76f | lukeholman | 2023-07-31 | wflow_publish("analysis/gwas_adaptive_shrinkage.Rmd") |
| html | 3a5750d | lukeholman | 2023-03-21 | Build site. |
| Rmd | a5026e7 | lukeholman | 2023-03-21 | wflow_publish("analysis/gwas_adaptive_shrinkage.Rmd") |
| html | 5078e01 | lukeholman | 2023-03-21 | Build site. |
| Rmd | 0ab0260 | lukeholman | 2023-03-21 | wflow_publish("analysis/gwas_adaptive_shrinkage.Rmd") |
| html | 7449a90 | lukeholman | 2021-10-01 | Build site. |
| html | 8d14298 | lukeholman | 2021-09-26 | Build site. |
| Rmd | af15dd6 | lukeholman | 2021-09-26 | Commit Sept 2021 |
| Rmd | 3855d33 | lukeholman | 2021-03-04 | big fist commit 2021 |
| html | 3855d33 | lukeholman | 2021-03-04 | big fist commit 2021 |
| Rmd | 8d54ea5 | Luke Holman | 2018-12-23 | Initial commit |
| html | 8d54ea5 | Luke Holman | 2018-12-23 | Initial commit |
library(tidyverse)
library(ashr)
library(mashr)
library(glue)
library(kableExtra)
library(gridExtra)
# Get the list of SNPs (or SNP clumps in 100% LD) which are in approx. LD with one another
get_mashr_data <- function(){
SNPs_in_LD <- read.table("data/derived/SNPs_in_LD", header = FALSE)
data_for_mashr <- read_csv("data/derived/all_univariate_GEMMA_results.csv") %>%
filter(str_detect(SNPs, ",") | SNPs %in% SNPs_in_LD$V1) %>%
select(SNPs, starts_with("beta"), starts_with("SE"))
single_snps <- data_for_mashr %>% filter(!str_detect(SNPs, ","))
multiple_snps <- data_for_mashr %>% filter(str_detect(SNPs, ","))
multiple_snps <- multiple_snps[enframe(strsplit(multiple_snps$SNPs, split = ", ")) %>%
unnest(value) %>%
filter(value %in% SNPs_in_LD$V1) %>%
pull(name), ]
bind_rows(single_snps, multiple_snps) %>% arrange(SNPs)
}mashrThe following code runs the package mashr, which
attempts to infer more accurate estimates of the true SNP effect sizes
by applying shrinkage based on the multivariate structure of the data.
mashr accepts a matrix of four effect sizes for each locus
(one fore each fitness trait) as well as the four associated standard
errors. mashr then uses Bayesian multivariate mixture
models to attempt to improve the estimates of the true effect sizes by
applying shrinkage. It can also calculate the local false sign rate and
other useful things. For more information on mashr, see the
package
website, and the associated paper by Urbut, Wang, and Stephens (doi:10.1101/096552).
There are two ways to run mashr: fitting a mixture model
where the covariance matrices describing the true relationships between
the effects sizes are estimated from the data by a deconvolution
algorithm (termed the ‘data-driven’ approach), and an alternative
approach where these covariance matrices are defined a priori
by the user (the ‘canonical’ approach). The following code runs both
approaches, though only the analysis using the “data-driven” approach is
presented in the paper (we felt that the canonical approach did not add
any new insights, and the likelihood of the mashr model was
far better when using the data-driven approach).
The mashr analysis use the default ‘null-biased’ prior,
meaning that our prior is that loci with no effect on any of the fitness
components are 10-fold more common than any of the other possibilities
(reflecting our expectation that genetic polymorphisms with strong
effects on fitness are probably rare).
Because mashr models are computationally intensive, and
to avoid issues of pseudoreplication in downstream analyses, we ran
mashr on a subset of the 1,207,357 loci that were analysed
using GEMMA. We pruned this list to a subset of 209,053 loci that were
in approximate linkage disequilibrium with one another.
run_mashr <- function(beta_and_se, mashr_mode, ED_p_cutoff = NULL){
print(str(beta_and_se))
mashr_setup <- function(beta_and_se){
betas <- beta_and_se %>% select(starts_with("beta")) %>% as.matrix()
SEs <- beta_and_se %>% select(starts_with("SE")) %>% as.matrix()
rownames(betas) <- beta_and_se$SNP
rownames(SEs) <- beta_and_se$SNP
mash_set_data(betas, SEs)
}
mash_data <- mashr_setup(beta_and_se)
# Setting mashr_mode == "ED" makes mashr select the covariance matrices for us, using the
# software Extreme Deconvolution. This software "reconstructs the error-deconvolved or 'underlying'
# distribution function common to all samples, even when the individual data points are samples from different distributions"
# Following the mashr vignette, we initialise the algorithm in ED using the principal components of the strongest effects in the dataset
# Reference for ED: https://arxiv.org/abs/0905.2979
if(mashr_mode == "ED"){
# Find the strongest effects in the data
m.1by1 <- mash_1by1(mash_data)
strong <- get_significant_results(m.1by1, thresh = ED_p_cutoff)
# Obtain data-driven covariance matrices by running Extreme Deconvolution
U.pca <- cov_pca(mash_data, npc = 4, subset = strong)
U <- cov_ed(mash_data, U.pca, subset = strong)
}
# Otherwise, we define the covariance matrices ourselves (a long list of a priori interesting matrices are checked)
if(mashr_mode == "canonical"){
make_SA_matrix <- function(r) matrix(c(1,1,r,r,1,1,r,r,r,r,1,1,r,r,1,1), ncol=4)
make_age_antag_matrix <- function(r) matrix(c(1,r,1,r,r,1,r,1,1,r,1,r,r,1,r,1), ncol=4)
make_sex_specific <- function(mat, sex){
if(sex == "F") {mat[, 3:4] <- 0; mat[3:4, ] <- 0}
if(sex == "M") {mat[, 1:2] <- 0; mat[1:2, ] <- 0}
mat
}
make_age_specific <- function(mat, age){
if(age == "early") {mat[, c(2,4)] <- 0; mat[c(2,4), ] <- 0}
if(age == "late") {mat[, c(1,3)] <- 0; mat[c(1,3), ] <- 0}
mat
}
add_matrix <- function(mat, mat_list, name){
mat_list[[length(mat_list) + 1]] <- mat
names(mat_list)[length(mat_list)] <- name
mat_list
}
id_matrix <- matrix(1, ncol=4, nrow=4)
# Get the mashr default canonical covariance matrices: this includes the ones
# called "null", "uniform", and "same sign" in the list that precedes this code chunk
U <- cov_canonical(mash_data)
# And now our custom covariance matrices:
# Identical across ages, but sex-antagonistic
U <- make_SA_matrix(-0.25) %>% add_matrix(U, "Sex_antag_0.25")
U <- make_SA_matrix(-0.5) %>% add_matrix(U, "Sex_antag_0.5")
U <- make_SA_matrix(-0.75) %>% add_matrix(U, "Sex_antag_0.75")
U <- make_SA_matrix(-1) %>% add_matrix(U, "Sex_antag_1.0")
# Identical across sexes, but age-antagonistic
U <- make_age_antag_matrix(-0.25) %>% add_matrix(U, "Age_antag_0.25")
U <- make_age_antag_matrix(-0.5) %>% add_matrix(U, "Age_antag_0.5")
U <- make_age_antag_matrix(-0.75) %>% add_matrix(U, "Age_antag_0.75")
U <- make_age_antag_matrix(-1) %>% add_matrix(U, "Age_antag_1.0")
# Sex-specific, identical effect in young and old
U <- id_matrix %>% make_sex_specific("F") %>% add_matrix(U, "Female_specific_1")
U <- id_matrix %>% make_sex_specific("M") %>% add_matrix(U, "Male_specific_1")
# Age-specific, identical effect in males and females
U <- id_matrix %>% make_age_specific("early") %>% add_matrix(U, "Early_life_specific_1")
U <- id_matrix %>% make_age_specific("late") %>% add_matrix(U, "Late_life_specific_1")
# Positively correlated but variable effect across ages, and also sex-specific
U <- make_age_antag_matrix(0.25) %>% make_sex_specific("F") %>% add_matrix(U, "Female_specific_0.25")
U <- make_age_antag_matrix(0.5) %>% make_sex_specific("F") %>% add_matrix(U, "Female_specific_0.5")
U <- make_age_antag_matrix(0.75) %>% make_sex_specific("F") %>% add_matrix(U, "Female_specific_0.75")
U <- make_age_antag_matrix(0.25) %>% make_sex_specific("M") %>% add_matrix(U, "Male_specific_0.25")
U <- make_age_antag_matrix(0.5) %>% make_sex_specific("M") %>% add_matrix(U, "Male_specific_0.5")
U <- make_age_antag_matrix(0.75) %>% make_sex_specific("M") %>% add_matrix(U, "Male_specific_0.75")
# Negatively correlated across ages, and also sex-specific
U <- make_age_antag_matrix(-0.25) %>% make_sex_specific("F") %>% add_matrix(U, "Female_specific_age_antag_0.25")
U <- make_age_antag_matrix(-0.5) %>% make_sex_specific("F") %>% add_matrix(U, "Female_specific_age_antag_0.5")
U <- make_age_antag_matrix(-0.75) %>% make_sex_specific("F") %>% add_matrix(U, "Female_specific_age_antag_0.75")
U <- make_age_antag_matrix(-0.25) %>% make_sex_specific("M") %>% add_matrix(U, "Male_specific_age_antag_0.25")
U <- make_age_antag_matrix(-0.5) %>% make_sex_specific("M") %>% add_matrix(U, "Male_specific_age_antag_0.5")
U <- make_age_antag_matrix(-0.75) %>% make_sex_specific("M") %>% add_matrix(U, "Male_specific_age_antag_0.75")
# Positively correlated but variable effect across sexes, and also age-specific
U <- make_SA_matrix(0.25) %>% make_age_specific("early") %>% add_matrix(U, "Early_life_specific_0.25")
U <- make_SA_matrix(0.5) %>% make_age_specific("early") %>% add_matrix(U, "Early_life_specific_0.5")
U <- make_SA_matrix(0.75) %>% make_age_specific("early") %>% add_matrix(U, "Early_life_specific_0.75")
U <- make_SA_matrix(0.25) %>% make_age_specific("late") %>% add_matrix(U, "Late_life_specific_0.25")
U <- make_SA_matrix(0.5) %>% make_age_specific("late") %>% add_matrix(U, "Late_life_specific_0.5")
U <- make_SA_matrix(0.75) %>% make_age_specific("late") %>% add_matrix(U, "Late_life_specific_0.75")
# Negatively correlated but variable effect across sexes, and also age-specific
U <- make_SA_matrix(-0.25) %>% make_age_specific("early") %>% add_matrix(U, "Early_life_antag_0.25")
U <- make_SA_matrix(-0.5) %>% make_age_specific("early") %>% add_matrix(U, "Early_life_antag_0.5")
U <- make_SA_matrix(-0.75) %>% make_age_specific("early") %>% add_matrix(U, "Early_life_antag_0.75")
U <- make_SA_matrix(-0.25) %>% make_age_specific("late") %>% add_matrix(U, "Late_life_antag_0.25")
U <- make_SA_matrix(-0.5) %>% make_age_specific("late") %>% add_matrix(U, "Late_life_antag_0.5")
U <- make_SA_matrix(-0.75) %>% make_age_specific("late") %>% add_matrix(U, "Late_life_antag_0.75")
}
return(mash(data = mash_data, Ulist = U)) # Run mashr
}
if(!file.exists("data/derived/mashr_results_canonical.rds")){
data_for_mashr <- get_mashr_data()
print("Starting the data-driven analysis")
run_mashr(data_for_mashr, mashr_mode = "ED", ED_p_cutoff = 0.2) %>%
saveRDS(file = "data/derived/mashr_results_ED.rds")
print("Starting the canonical analysis")
run_mashr(data_for_mashr, mashr_mode = "canonical") %>%
saveRDS(file = "data/derived/mashr_results_canonical.rds")
mashr_results_ED <- read_rds("data/derived/mashr_results_ED.rds")
mashr_results_canonical <- read_rds("data/derived/mashr_results_canonical.rds")
} else {
data_for_mashr <- get_mashr_data()
mashr_results_ED <- read_rds("data/derived/mashr_results_ED.rds")
mashr_results_canonical <- read_rds("data/derived/mashr_results_canonical.rds")
}mashr results to the databaseHere, we make a single large dataframe holding all of the ‘raw’
results from the univariate GEMMA analysis, and the corresponding
“shrinked” results from mashr (for both the canonical and
data-driven mashr analyses). Because it is so large, we add
this sheet of results to the database, allowing memory-efficient
searching, joins, etc.
all_univariate_lmm_results <- read_csv("data/derived/all_univariate_GEMMA_results.csv") %>%
rename_at(vars(-SNPs), ~ str_c(., "_raw"))
mashr_snps <- data_for_mashr$SNPs
# canonical_estimates <- get_pm(mashr_results_canonical) %>%
# as_tibble() %>%
# rename_all(~str_c(., "_mashr_canonical"))
ED_estimates <- get_pm(mashr_results_ED) %>%
as_tibble() %>%
rename_all(~str_c(., "_mashr_ED"))
# lfsr_canonical <- get_lfsr(mashr_results_canonical) %>%
# as_tibble() %>%
# rename_all(~str_replace_all(., "beta", "LFSR")) %>%
# rename_all(~str_c(., "_mashr_canonical"))
lfsr_ED <- get_lfsr(mashr_results_ED) %>%
as_tibble() %>%
rename_all(~str_replace_all(., "beta", "LFSR")) %>%
rename_all(~str_c(., "_mashr_ED"))
all_mashr_results <- bind_cols(
tibble(SNPs = mashr_snps),
# canonical_estimates,
ED_estimates,
# lfsr_canonical,
lfsr_ED)
all_univariate_lmm_results <- left_join(
all_univariate_lmm_results,
all_mashr_results, by = "SNPs")
nested <- all_univariate_lmm_results %>%
filter(str_detect(SNPs, ", "))
split_snps <- strsplit(nested$SNPs, split = ", ")
nested <- lapply(1:nrow(nested),
function(i) {
data.frame(SNP = split_snps[[i]],
SNP_clump = nested$SNPs[i],
nested[i,] %>% select(-SNPs), stringsAsFactors = FALSE)
}) %>%
do.call("rbind", .) %>% as_tibble()
rm(split_snps)
all_univariate_lmm_results <- all_univariate_lmm_results %>%
filter(!str_detect(SNPs, ", ")) %>%
rename(SNP_clump = SNPs) %>% mutate(SNP = SNP_clump) %>%
select(SNP, SNP_clump, everything()) %>%
bind_rows(nested) %>%
arrange(SNP)
# Merge in the mixture proportions
# all_univariate_lmm_results <-
# all_univariate_lmm_results %>%
# left_join(mixture_assignment_probabilities, by = "SNP_clump")
db <- DBI::dbConnect(RSQLite::SQLite(),
"data/derived/annotations.sqlite3", create = FALSE)
DBI::dbRemoveTable(db, "univariate_lmm_results")
db %>% copy_to(all_univariate_lmm_results,
"univariate_lmm_results", temporary = FALSE)
sessionInfo()R version 4.2.2 (2022-10-31)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur ... 10.16
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2/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] gridExtra_2.3 kableExtra_1.3.4 glue_1.6.2 mashr_0.2.69
[5] ashr_2.2-54 lubridate_1.9.2 forcats_1.0.0 stringr_1.5.0
[9] dplyr_1.1.0 purrr_1.0.1 readr_2.1.4 tidyr_1.3.0
[13] tibble_3.1.8 ggplot2_3.4.1 tidyverse_2.0.0 workflowr_1.7.0.4
loaded via a namespace (and not attached):
[1] httr_1.4.5 sass_0.4.5 bit64_4.0.5 vroom_1.6.1
[5] jsonlite_1.8.4 viridisLite_0.4.1 bslib_0.4.2 rmeta_3.0
[9] assertthat_0.2.1 getPass_0.2-2 mixsqp_0.3-48 yaml_2.3.7
[13] pillar_1.8.1 lattice_0.20-45 digest_0.6.31 promises_1.2.0.1
[17] rvest_1.0.3 colorspace_2.1-0 htmltools_0.5.4 httpuv_1.6.9
[21] Matrix_1.5-1 plyr_1.8.8 pkgconfig_2.0.3 invgamma_1.1
[25] mvtnorm_1.1-3 scales_1.2.1 webshot_0.5.4 processx_3.8.0
[29] svglite_2.1.1 whisker_0.4.1 later_1.3.0 tzdb_0.3.0
[33] timechange_0.2.0 git2r_0.31.0 generics_0.1.3 ellipsis_0.3.2
[37] cachem_1.0.7 withr_2.5.0 cli_3.6.0 crayon_1.5.2
[41] magrittr_2.0.3 evaluate_0.20 ps_1.7.2 fs_1.6.1
[45] fansi_1.0.4 xml2_1.3.3 truncnorm_1.0-8 tools_4.2.2
[49] hms_1.1.2 lifecycle_1.0.3 munsell_0.5.0 irlba_2.3.5.1
[53] callr_3.7.3 compiler_4.2.2 jquerylib_0.1.4 systemfonts_1.0.4
[57] rlang_1.0.6 grid_4.2.2 rstudioapi_0.14 rmarkdown_2.20
[61] gtable_0.3.1 abind_1.4-5 R6_2.5.1 knitr_1.42
[65] bit_4.0.5 fastmap_1.1.1 utf8_1.2.2 rprojroot_2.0.3
[69] stringi_1.7.12 parallel_4.2.2 SQUAREM_2021.1 Rcpp_1.0.10
[73] vctrs_0.5.2 tidyselect_1.2.0 xfun_0.37