Last updated: 2023-11-29

Checks: 7 0

Knit directory: multigroup_ctwas_analysis/

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(20231112)

The command set.seed(20231112) 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: 8883dec

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 8883dec. 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:    .Rhistory
    Ignored:    .Rproj.user/

Unstaged changes:
    Modified:   analysis/detect_LD_mismatch_gwas_UKBBref_susie_rss.Rmd

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/detect_LD_mismatch_gwas_UKBBref_DENTIST_susie_rss.Rmd) and HTML (docs/detect_LD_mismatch_gwas_UKBBref_DENTIST_susie_rss.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 8883dec kevinlkx 2023-11-29 update susie_rss LDL result with chr6 result
html c9fee4b kevinlkx 2023-11-29 Build site.
Rmd c195252 kevinlkx 2023-11-29 compare DENTIST vs. susie_rss results

Load packages and functions

library(ctwas)
library(susieR)
library(foreach)
library(data.table)
library(tidyverse)

LD Regions (ldetect blocks)

regions <- system.file("extdata/ldetect", "EUR.b38.bed", package = "ctwas")
regions_df <- read.table(regions, header = T)
regions_df <- regions_df %>% dplyr::arrange(chr, start, stop) %>% dplyr::mutate(locus = 1:nrow(regions_df))

LDL

trait <- "LDL"

Compare DENTIST and SuSiE RSS pvalues in chr22

load DENTIST results

CHR=22
dentist.dir <- paste0("/project2/xinhe/shared_data/multigroup_ctwas/DENTIST/", trait)
dentist.res.file <- file.path(dentist.dir, paste0("LDL-ukb-d-30780_irnt.ukb_chr", CHR, ".b38.DENTIST.full.txt"))
dentist.chr.df <- data.table::fread(dentist.res.file)
colnames(dentist.chr.df) <- c("rsID", "chisq", "LP", "ifDup")

dentist.chr.snps <- dentist.chr.df$rsID
dentist_detected_snps <- dentist.chr.df$rsID[which(dentist.chr.df$LP > -log10(5e-8))]
cat(sprintf("%d variants with DENTIST result in chr%s \n", length(dentist.chr.snps), CHR))
cat(sprintf("%d detected variants with DENTIST pvalue < 5e-8.\n", length(dentist_detected_snps)))
# 118796 variants with DENTIST result in chr22 
# 0 detected variants with DENTIST pvalue < 5e-8.

load SuSiE RSS result

select_loci <- paste0("chr", CHR)
susie_rss_dir <- paste0("/project2/xinhe/shared_data/multigroup_ctwas/ld_mismatch_susie_rss/", trait)
condz_dist_chr <- readRDS(file.path(susie_rss_dir, paste0(trait, ".condz.dist.", select_loci, "loci.rds")))
condz_dist_chr.df <- do.call(rbind.data.frame, condz_dist_chr)

susierss.chr.snps <- condz_dist_chr.df$id
susierss_detected_snps <- condz_dist_chr.df$id[which(condz_dist_chr.df$p_diff < 5e-8)]
susierss_detected_flipped_snps <- condz_dist_chr.df$id[which(condz_dist_chr.df$logLR > 2 & abs(condz_dist_chr.df$z) > 2)]
cat(sprintf("%d variants with susie_rss result in chr%s \n", length(susierss.chr.snps), CHR))
cat(sprintf("%d detected variants with susie_rss pvalue < 5e-8.\n", length(susierss_detected_snps)))
cat(sprintf("%d detected variants with susie_rss allele flipped.\n", length(susierss_detected_flipped_snps)))
# 120845 variants with susie_rss result in chr22 
# 93 detected variants with susie_rss pvalue < 5e-8.
# 3 detected variants with susie_rss allele flipped.
cat(sprintf("%d variants with DENTIST result \n", length(dentist.chr.snps)))
cat(sprintf("%d variants with susie_rss result \n", length(susierss.chr.snps)))
cat(sprintf("%d variants with both DENTIST and susie_rss result \n", length(intersect(dentist.chr.snps, susierss.chr.snps))))
cat(sprintf("%d variants detected by both DENTIST and susie_rss (pvalue < 5e-8).\n", length(intersect(dentist_detected_snps, susierss_detected_snps))))

common.snps <- intersect(dentist.chr.snps, susierss.chr.snps)
df <- data.frame(rsID = common.snps, dentist.pval = NA, susierss.pval = NA)
df$dentist.pval <- dentist.chr.df$LP[match(common.snps, dentist.chr.df$rsID)]
df$susierss.pval <- -log10(condz_dist_chr.df$p_diff[match(common.snps, condz_dist_chr.df$id)])

ggplot(df, aes(x = dentist.pval, y = susierss.pval)) +
  geom_point(alpha=0.3) +
  xlim(0, 100) + ylim(0, 100) +
  labs(x = "DENTIST -log10P", y = "SuSiE RSS -log10P", title = paste0(trait, " chr", CHR)) +
  geom_abline(intercept = 0, slope = 1) + 
  geom_vline(xintercept = -log10(5e-8), col = "red") +
  geom_hline(yintercept = -log10(5e-8), col = "red") +
  theme_bw()
# Warning: Removed 6 rows containing missing values (`geom_point()`).

Version Author Date
c9fee4b kevinlkx 2023-11-29 
# 118796 variants with DENTIST result 
# 120845 variants with susie_rss result 
# 118796 variants with both DENTIST and susie_rss result 
# 0 variants detected by both DENTIST and susie_rss (pvalue < 5e-8).

Genome-wide results

DENTIST results genome-wide

dentist.dir <- paste0("/project2/xinhe/shared_data/multigroup_ctwas/DENTIST/", trait)
dentist.res <- foreach(CHR=1:22) %do% {
  dentist.res.file <- file.path(dentist.dir, paste0("LDL-ukb-d-30780_irnt.ukb_chr", CHR, ".b38.DENTIST.full.txt"))
  if(file.exists(dentist.res.file)){
    dentist.chr.res <- data.table::fread(dentist.res.file)
    colnames(dentist.chr.res) <- c("rsID", "chisq", "LP", "ifDup")
    cbind(chr = CHR, dentist.chr.res)
  }else{
    NULL
  }
}
dentist.res.df <- do.call(rbind.data.frame, dentist.res)
data.table::fwrite(dentist.res.df, file.path(dentist.dir, paste0("LDL-ukb-d-30780_irnt.ukb_chrs.b38.DENTIST.full.txt.gz")), sep = "\t", col.names = TRUE)
dentist.dir <- paste0("/project2/xinhe/shared_data/multigroup_ctwas/DENTIST/", trait)
dentist.res.df <- data.table::fread(file.path(dentist.dir, paste0("LDL-ukb-d-30780_irnt.ukb_chrs.b38.DENTIST.full.txt.gz")))
cat(sprintf("%d variants with DENTIST result in total \n", length(dentist.res.df$rsID)))
cat(length(which(dentist.res.df$LP > -log10(5e-8))), "variants genome-wide with DENTIST pval < 5e-8")
# 8687900 variants with DENTIST result in total 
# 292 variants genome-wide with DENTIST pval < 5e-8

SuSiE RSS results genome-wide

susie_rss_dir <- paste0("/project2/xinhe/shared_data/multigroup_ctwas/ld_mismatch_susie_rss/", trait)
condz_dist.res <- foreach(CHR=1:22) %do% {
  select_loci <- paste0("chr", CHR)
  condz_dist.file <- file.path(susie_rss_dir, paste0(trait, ".condz.dist.", select_loci, "loci.rds"))
  if(file.exists(condz_dist.file)){
    condz_dist_chr <- readRDS(condz_dist.file)
    condz_dist_chr.df <- do.call(rbind.data.frame, condz_dist_chr)
    cbind(chr = CHR, condz_dist_chr.df)
  }else{
    NULL
  }
}
# summary(condz_dist.res)
condz_dist_all.df <- do.call(rbind.data.frame, condz_dist.res)
saveRDS(condz_dist_all.df, file.path(susie_rss_dir, paste0(trait, ".susie_rss.condz.dist.rds")))
susie_rss_dir <- paste0("/project2/xinhe/shared_data/multigroup_ctwas/ld_mismatch_susie_rss/", trait)
condz_dist_all.df <- readRDS(file.path(susie_rss_dir, paste0(trait, ".susie_rss.condz.dist.rds")))
susierss_detected_snps <- condz_dist_all.df$id[which(condz_dist_all.df$p_diff < 5e-8)]
susierss_detected_flipped_snps <- condz_dist_all.df$id[which(condz_dist_all.df$logLR > 2 & abs(condz_dist_all.df$z) > 2)]
cat(sprintf("%d variants with susie_rss result in total \n", length(condz_dist_all.df$id)))
cat(sprintf("%d detected variants with susie_rss pvalue < 5e-8.\n", length(susierss_detected_snps)))
cat(sprintf("%d detected variants with susie_rss allele flipped.\n", length(susierss_detected_flipped_snps)))
# 8841628 variants with susie_rss result in total 
# 4607 detected variants with susie_rss pvalue < 5e-8.
# 8 detected variants with susie_rss allele flipped.

aFib

trait <- "aFib"

Compare DENTIST and SuSiE RSS pvalues in chr22

load DENTIST results

CHR=22
dentist.dir <- paste0("/project2/xinhe/shared_data/multigroup_ctwas/DENTIST/", trait)
dentist.res.file <- file.path(dentist.dir, paste0("aFib-ebi-a-GCST006414.ukb_chr", CHR, ".b38.DENTIST.full.txt"))
dentist.chr.df <- data.table::fread(dentist.res.file)
colnames(dentist.chr.df) <- c("rsID", "chisq", "LP", "ifDup")

dentist.chr.snps <- dentist.chr.df$rsID
dentist_detected_snps <- dentist.chr.df$rsID[which(dentist.chr.df$LP > -log10(5e-8))]
cat(sprintf("%d variants with DENTIST result in chr%s \n", length(dentist.chr.snps), CHR))
cat(sprintf("%d detected variants with DENTIST pvalue < 5e-8.\n", length(dentist_detected_snps)))
# 109507 variants with DENTIST result in chr22 
# 2802 detected variants with DENTIST pvalue < 5e-8.

load SuSiE RSS result

select_loci <- paste0("chr", CHR)
susie_rss_dir <- paste0("/project2/xinhe/shared_data/multigroup_ctwas/ld_mismatch_susie_rss/", trait)
condz_dist_chr <- readRDS(file.path(susie_rss_dir, paste0(trait, ".condz.dist.", select_loci, "loci.rds")))
condz_dist_chr.df <- do.call(rbind.data.frame, condz_dist_chr)

susierss.chr.snps <- condz_dist_chr.df$id
susierss_detected_snps <- condz_dist_chr.df$id[which(condz_dist_chr.df$p_diff < 5e-8)]
susierss_detected_flipped_snps <- condz_dist_chr.df$id[which(condz_dist_chr.df$logLR > 2 & abs(condz_dist_chr.df$z) > 2)]
cat(sprintf("%d variants with susie_rss result in chr%s \n", length(susierss.chr.snps), CHR))
cat(sprintf("%d detected variants with susie_rss pvalue < 5e-8.\n", length(susierss_detected_snps)))
cat(sprintf("%d detected variants with susie_rss allele flipped.\n", length(susierss_detected_flipped_snps)))
# 110716 variants with susie_rss result in chr22 
# 449 detected variants with susie_rss pvalue < 5e-8.
# 12 detected variants with susie_rss allele flipped.
cat(sprintf("%d variants with DENTIST result \n", length(dentist.chr.snps)))
cat(sprintf("%d variants with susie_rss result \n", length(susierss.chr.snps)))
cat(sprintf("%d variants with both DENTIST and susie_rss result \n", length(intersect(dentist.chr.snps, susierss.chr.snps))))
cat(sprintf("%d variants detected by both DENTIST and susie_rss (pvalue < 5e-8).\n", length(intersect(dentist_detected_snps, susierss_detected_snps))))

common.snps <- intersect(dentist.chr.snps, susierss.chr.snps)
df <- data.frame(rsID = common.snps, dentist.pval = NA, susierss.pval = NA)
df$dentist.pval <- dentist.chr.df$LP[match(common.snps, dentist.chr.df$rsID)]
df$susierss.pval <- -log10(condz_dist_chr.df$p_diff[match(common.snps, condz_dist_chr.df$id)])

ggplot(df, aes(x = dentist.pval, y = susierss.pval)) +
  geom_point(alpha=0.3) +
  xlim(0, 100) + ylim(0, 100) +
  labs(x = "DENTIST -log10P", y = "SuSiE RSS -log10P", title = paste0(trait, " chr", CHR)) +
  geom_abline(intercept = 0, slope = 1) + 
  geom_vline(xintercept = -log10(5e-8), col = "red") +
  geom_hline(yintercept = -log10(5e-8), col = "red") +
  theme_bw()
# Warning: Removed 23 rows containing missing values (`geom_point()`).

Version Author Date
c9fee4b kevinlkx 2023-11-29 
# 109507 variants with DENTIST result 
# 110716 variants with susie_rss result 
# 109504 variants with both DENTIST and susie_rss result 
# 366 variants detected by both DENTIST and susie_rss (pvalue < 5e-8).

Genome-wide results

DENTIST results genome-wide

dentist.dir <- paste0("/project2/xinhe/shared_data/multigroup_ctwas/DENTIST/", trait)
dentist.res <- foreach(CHR=1:22) %do% {
  dentist.res.file <- file.path(dentist.dir, paste0("aFib-ebi-a-GCST006414.ukb_chr", CHR, ".b38.DENTIST.full.txt"))
  if(file.exists(dentist.res.file)){
    dentist.chr.res <- data.table::fread(dentist.res.file)
    colnames(dentist.chr.res) <- c("rsID", "chisq", "LP", "ifDup")
    cbind(chr = CHR, dentist.chr.res)
  }else{
    NULL
  }
}
dentist.res.df <- do.call(rbind.data.frame, dentist.res)
data.table::fwrite(dentist.res.df, file.path(dentist.dir, paste0("aFib-ebi-a-GCST006414.ukb_chrs.b38.DENTIST.full.txt.gz")), sep = "\t", col.names = TRUE)
dentist.dir <- paste0("/project2/xinhe/shared_data/multigroup_ctwas/DENTIST/", trait)
dentist.res.df <- data.table::fread(file.path(dentist.dir, paste0("aFib-ebi-a-GCST006414.ukb_chrs.b38.DENTIST.full.txt.gz")))
cat(sprintf("%d variants with DENTIST result in total \n", length(dentist.res.df$rsID)))
cat(length(which(dentist.res.df$LP > -log10(5e-8))), "variants genome-wide with DENTIST pval < 5e-8")
# 8230059 variants with DENTIST result in total 
# 321808 variants genome-wide with DENTIST pval < 5e-8

SuSiE RSS results genome-wide

susie_rss_dir <- paste0("/project2/xinhe/shared_data/multigroup_ctwas/ld_mismatch_susie_rss/", trait)
condz_dist.res <- foreach(CHR=1:22) %do% {
  select_loci <- paste0("chr", CHR)
  condz_dist.file <- file.path(susie_rss_dir, paste0(trait, ".condz.dist.", select_loci, "loci.rds"))
  if(file.exists(condz_dist.file)){
    condz_dist_chr <- readRDS(condz_dist.file)
    condz_dist_chr.df <- do.call(rbind.data.frame, condz_dist_chr)
    cbind(chr = CHR, condz_dist_chr.df)
  }else{
    NULL
  }
}
# summary(condz_dist.res)
condz_dist_all.df <- do.call(rbind.data.frame, condz_dist.res)
saveRDS(condz_dist_all.df, file.path(susie_rss_dir, paste0(trait, ".susie_rss.condz.dist.rds")))
susie_rss_dir <- paste0("/project2/xinhe/shared_data/multigroup_ctwas/ld_mismatch_susie_rss/", trait)
condz_dist_all.df <- readRDS(file.path(susie_rss_dir, paste0(trait, ".susie_rss.condz.dist.rds")))
susierss_detected_snps <- condz_dist_all.df$id[which(condz_dist_all.df$p_diff < 5e-8)]
susierss_detected_flipped_snps <- condz_dist_all.df$id[which(condz_dist_all.df$logLR > 2 & abs(condz_dist_all.df$z) > 2)]
cat(sprintf("%d variants with susie_rss result in total \n", length(condz_dist_all.df$id)))
cat(sprintf("%d detected variants with susie_rss pvalue < 5e-8.\n", length(susierss_detected_snps)))
cat(sprintf("%d detected variants with susie_rss allele flipped.\n", length(susierss_detected_flipped_snps)))
# 8258166 variants with susie_rss result in total 
# 38330 detected variants with susie_rss pvalue < 5e-8.
# 439 detected variants with susie_rss allele flipped.

sessionInfo()
# R version 4.2.0 (2022-04-22)
# Platform: x86_64-pc-linux-gnu (64-bit)
# Running under: CentOS Linux 7 (Core)
# 
# Matrix products: default
# BLAS/LAPACK: /software/openblas-0.3.13-el7-x86_64/lib/libopenblas_haswellp-r0.3.13.so
# 
# locale:
#  [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C         LC_TIME=C           
#  [4] LC_COLLATE=C         LC_MONETARY=C        LC_MESSAGES=C       
#  [7] LC_PAPER=C           LC_NAME=C            LC_ADDRESS=C        
# [10] LC_TELEPHONE=C       LC_MEASUREMENT=C     LC_IDENTIFICATION=C 
# 
# attached base packages:
# [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#  [1] forcats_1.0.0     stringr_1.5.0     dplyr_1.1.0       purrr_1.0.1      
#  [5] readr_2.1.4       tidyr_1.3.0       tibble_3.1.8      ggplot2_3.4.1    
#  [9] tidyverse_1.3.2   data.table_1.14.6 foreach_1.5.2     susieR_0.12.35   
# [13] ctwas_0.1.35      workflowr_1.7.0  
# 
# loaded via a namespace (and not attached):
#  [1] matrixStats_0.63.0  fs_1.6.1            lubridate_1.9.2    
#  [4] httr_1.4.4          rprojroot_2.0.3     tools_4.2.0        
#  [7] backports_1.4.1     bslib_0.4.2         utf8_1.2.3         
# [10] R6_2.5.1            irlba_2.3.5         DBI_1.1.3          
# [13] colorspace_2.1-0    withr_2.5.0         tidyselect_1.2.0   
# [16] processx_3.8.0      compiler_4.2.0      git2r_0.30.1       
# [19] cli_3.6.0           rvest_1.0.3         logging_0.10-108   
# [22] xml2_1.3.3          labeling_0.4.2      sass_0.4.5         
# [25] scales_1.2.1        callr_3.7.3         mixsqp_0.3-43      
# [28] digest_0.6.31       R.utils_2.12.2      rmarkdown_2.20     
# [31] pkgconfig_2.0.3     htmltools_0.5.4     highr_0.10         
# [34] dbplyr_2.3.0        fastmap_1.1.0       rlang_1.0.6        
# [37] readxl_1.4.2        rstudioapi_0.14     jquerylib_0.1.4    
# [40] generics_0.1.3      farver_2.1.1        jsonlite_1.8.4     
# [43] R.oo_1.25.0         googlesheets4_1.0.1 magrittr_2.0.3     
# [46] Matrix_1.5-3        Rcpp_1.0.10         munsell_0.5.0      
# [49] fansi_1.0.4         R.methodsS3_1.8.2   lifecycle_1.0.3    
# [52] stringi_1.7.12      whisker_0.4         yaml_2.3.7         
# [55] plyr_1.8.7          grid_4.2.0          promises_1.2.0.1   
# [58] crayon_1.5.2        lattice_0.20-45     haven_2.5.1        
# [61] hms_1.1.2           knitr_1.42          ps_1.7.2           
# [64] pillar_1.8.1        codetools_0.2-18    reprex_2.0.2       
# [67] glue_1.6.2          evaluate_0.20       getPass_0.2-2      
# [70] modelr_0.1.10       vctrs_0.5.2         tzdb_0.3.0         
# [73] httpuv_1.6.5        cellranger_1.1.0    gtable_0.3.1       
# [76] pgenlibr_0.3.3      reshape_0.8.9       assertthat_0.2.1   
# [79] cachem_1.0.6        xfun_0.37           broom_1.0.3        
# [82] later_1.3.0         googledrive_2.0.0   gargle_1.3.0       
# [85] iterators_1.0.14    timechange_0.2.0    ellipsis_0.3.2