Sensitivity Analysis - estimation of genomic inflation

1 Paths

dir.single.cohort <- "../DATASETS"
dir.gasparoni <- file.path(dir.single.cohort, "GASPARONI/step9_single_cpg_pval")
dir.london <- file.path(dir.single.cohort, "LONDON/step7_single_cpg_pval")
dir.mtsinai <- file.path(dir.single.cohort, "MtSinai/step7_single_cpg_pval")
dir.rosmap <- file.path(dir.single.cohort, "ROSMAP/step9_single_cpg_pval")
dir.meta.single.cpg <- "../meta_analysis_single_cpg_results/"
dir.data <- "./DATASETS"
dir.fig <- "./FIGURES"

2 Estimation of inflation

library(dplyr)
library(bacon)
library(GWASTools)

2.1 Auxiliary functions

estimation_of_inflation <- function(data){
  ### 1. Compute genomic inflation factor before bacon adjustment
  data$zvalue <- data$Estimate / data$StdErr
  data$chisq <- (data$zvalue) ^ 2

  # inflation factor - last term is median from chisq distrn with 1 df  
  inflationFactor <- median(data$chisq,na.rm = TRUE) / qchisq(0.5, 1)
  print("lambda")
  print(inflationFactor)

  # genome-wide sig cpgs
  sig <- ifelse(data$pValue < 2.4e-7, 1, 0)
  # table(sig)  # 1 sig

  ### 2. bacon analysis
  bc <- bacon(
    teststatistics = NULL,
    effectsizes =  data$Estimate,
    standarderrors = data$StdErr,
    na.exclude = TRUE
  )

  # inflation factor
  print("lambda.bacon")
  print(inflation(bc))
  
  ### 3. Create final dataset
  data.with.inflation <- data.frame(
    data,
    Estimate.bacon = bacon::es(bc),
    StdErr.bacon = bacon::se(bc),
    pValue.bacon = pval(bc),
    fdr.bacon = p.adjust(pval(bc), method = "fdr"),
    stringsAsFactors = FALSE
  )
  data.with.inflation <- data.with.inflation %>% select(-c(zvalue, chisq))
  return(
    list("data.with.inflation" = data.with.inflation, 
          "inflationFactor" =   inflationFactor,
          "estimatedInflation" = inflation(bc)
         )
    )
}

2.2 GASPARONI

est.inflation <- estimation_of_inflation(res)

## [1] "lambda"
## [1] 1.002218
## [1] "lambda.bacon"
##   sigma.0 
## 0.9852269

res.with.inflation <- est.inflation$data.with.inflation


cohort <- "GASPARONI"
write.csv(
    res.with.inflation,
    file.path(dir.data, cohort, "single_cpg_pVal_df.csv"),
    row.names = FALSE
)

# genome-wide sig cpgs
sig.fdr <- ifelse(res.with.inflation$fdr.bacon < 0.05, 1, 0)
# table(sig.fdr)   # 2 sig

2.3 LONDON

est.inflation <- estimation_of_inflation(res)

## [1] "lambda"
## [1] 1.208283
## [1] "lambda.bacon"
##  sigma.0 
## 1.082232

res.with.inflation <- est.inflation$data.with.inflation

cohort <- "LONDON"

write.csv(
    res.with.inflation,
    file.path(dir.data, cohort, "single_cpg_pVal_df.csv"),
    row.names = FALSE
)

# genome-wide sig cpgs
sig.fdr <- ifelse(res.with.inflation$fdr.bacon < 0.05, 1, 0)
# table(sig.fdr)   # 11 sig

2.4 MTSINAI

est.inflation <- estimation_of_inflation(res)

## [1] "lambda"
## [1] 1.248826
## [1] "lambda.bacon"
##  sigma.0 
## 1.077693

res.with.inflation <- est.inflation$data.with.inflation

cohort <- "MTSINAI"

write.csv(
    res.with.inflation,
    file.path(dir.data, cohort, "single_cpg_pVal_df.csv"),
    row.names = FALSE
)

# genome-wide sig cpgs
sig.fdr <- ifelse(res.with.inflation$fdr.bacon < 0.05, 1, 0)
# table(sig.fdr)   # 0 sig

2.5 ROSMAP

est.inflation <- estimation_of_inflation(res)

## [1] "lambda"
## [1] 1.015642
## [1] "lambda.bacon"
##  sigma.0 
## 1.001666

res.with.inflation <- est.inflation$data.with.inflation

cohort <- "ROSMAP"

write.csv(
    res.with.inflation,
    file.path(dir.data, cohort, "single_cpg_pVal_df.csv"),
    row.names = FALSE
)

# genome-wide sig cpgs
sig.fdr <- ifelse(res.with.inflation$fdr.bacon < 0.05, 1, 0)
# table(sig.fdr)   # 34 sig

3 Session information

devtools::session_info()

## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 4.0.2 (2020-06-22)
##  os       macOS Catalina 10.15.6      
##  system   x86_64, darwin17.0          
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  ctype    en_US.UTF-8                 
##  tz       America/New_York            
##  date     2020-08-11                  
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package      * version  date       lib source        
##  assertthat     0.2.1    2019-03-21 [1] CRAN (R 4.0.2)
##  backports      1.1.8    2020-06-17 [1] CRAN (R 4.0.2)
##  bacon        * 1.16.0   2020-04-27 [1] Bioconductor  
##  base64enc      0.1-3    2015-07-28 [1] CRAN (R 4.0.2)
##  Biobase      * 2.48.0   2020-04-27 [1] Bioconductor  
##  BiocGenerics * 0.34.0   2020-04-27 [1] Bioconductor  
##  BiocParallel * 1.22.0   2020-04-27 [1] Bioconductor  
##  bit            4.0.3    2020-07-30 [1] CRAN (R 4.0.2)
##  bit64          4.0.2    2020-07-30 [1] CRAN (R 4.0.2)
##  blob           1.2.1    2020-01-20 [1] CRAN (R 4.0.2)
##  broom          0.7.0    2020-07-09 [1] CRAN (R 4.0.2)
##  callr          3.4.3    2020-03-28 [1] CRAN (R 4.0.2)
##  cli            2.0.2    2020-02-28 [1] CRAN (R 4.0.2)
##  colorspace     1.4-1    2019-03-18 [1] CRAN (R 4.0.2)
##  conquer        1.0.1    2020-05-06 [1] CRAN (R 4.0.2)
##  crayon         1.3.4    2017-09-16 [1] CRAN (R 4.0.2)
##  DBI            1.1.0    2019-12-15 [1] CRAN (R 4.0.2)
##  desc           1.2.0    2018-05-01 [1] CRAN (R 4.0.2)
##  devtools       2.3.1    2020-07-21 [1] CRAN (R 4.0.2)
##  digest         0.6.25   2020-02-23 [1] CRAN (R 4.0.2)
##  DNAcopy        1.62.0   2020-04-27 [1] Bioconductor  
##  dplyr        * 1.0.1    2020-07-31 [1] CRAN (R 4.0.2)
##  ellipse      * 0.4.2    2020-05-27 [1] CRAN (R 4.0.2)
##  ellipsis       0.3.1    2020-05-15 [1] CRAN (R 4.0.2)
##  evaluate       0.14     2019-05-28 [1] CRAN (R 4.0.1)
##  fansi          0.4.1    2020-01-08 [1] CRAN (R 4.0.2)
##  fs             1.5.0    2020-07-31 [1] CRAN (R 4.0.2)
##  gdsfmt         1.24.1   2020-06-16 [1] Bioconductor  
##  generics       0.0.2    2018-11-29 [1] CRAN (R 4.0.2)
##  ggplot2      * 3.3.2    2020-06-19 [1] CRAN (R 4.0.2)
##  glue           1.4.1    2020-05-13 [1] CRAN (R 4.0.2)
##  gtable         0.3.0    2019-03-25 [1] CRAN (R 4.0.2)
##  GWASExactHW    1.01     2013-01-05 [1] CRAN (R 4.0.2)
##  GWASTools    * 1.34.0   2020-04-27 [1] Bioconductor  
##  htmltools      0.5.0    2020-06-16 [1] CRAN (R 4.0.2)
##  knitr          1.29     2020-06-23 [1] CRAN (R 4.0.2)
##  lattice        0.20-41  2020-04-02 [1] CRAN (R 4.0.2)
##  lifecycle      0.2.0    2020-03-06 [1] CRAN (R 4.0.2)
##  lmtest         0.9-37   2019-04-30 [1] CRAN (R 4.0.2)
##  logistf        1.23     2018-07-19 [1] CRAN (R 4.0.2)
##  magrittr       1.5      2014-11-22 [1] CRAN (R 4.0.2)
##  Matrix         1.2-18   2019-11-27 [1] CRAN (R 4.0.2)
##  MatrixModels   0.4-1    2015-08-22 [1] CRAN (R 4.0.2)
##  matrixStats    0.56.0   2020-03-13 [1] CRAN (R 4.0.2)
##  memoise        1.1.0    2017-04-21 [1] CRAN (R 4.0.2)
##  mgcv           1.8-31   2019-11-09 [1] CRAN (R 4.0.2)
##  mice           3.10.0.1 2020-08-02 [1] CRAN (R 4.0.2)
##  munsell        0.5.0    2018-06-12 [1] CRAN (R 4.0.2)
##  nlme           3.1-148  2020-05-24 [1] CRAN (R 4.0.2)
##  pillar         1.4.6    2020-07-10 [1] CRAN (R 4.0.2)
##  pkgbuild       1.1.0    2020-07-13 [1] CRAN (R 4.0.2)
##  pkgconfig      2.0.3    2019-09-22 [1] CRAN (R 4.0.2)
##  pkgload        1.1.0    2020-05-29 [1] CRAN (R 4.0.2)
##  prettyunits    1.1.1    2020-01-24 [1] CRAN (R 4.0.2)
##  processx       3.4.3    2020-07-05 [1] CRAN (R 4.0.2)
##  ps             1.3.3    2020-05-08 [1] CRAN (R 4.0.2)
##  purrr          0.3.4    2020-04-17 [1] CRAN (R 4.0.2)
##  quantreg       5.61     2020-07-09 [1] CRAN (R 4.0.2)
##  quantsmooth    1.54.0   2020-04-27 [1] Bioconductor  
##  R6             2.4.1    2019-11-12 [1] CRAN (R 4.0.2)
##  Rcpp           1.0.5    2020-07-06 [1] CRAN (R 4.0.2)
##  remotes        2.2.0    2020-07-21 [1] CRAN (R 4.0.2)
##  rlang          0.4.7    2020-07-09 [1] CRAN (R 4.0.2)
##  rmarkdown      2.3      2020-06-18 [1] CRAN (R 4.0.2)
##  rprojroot      1.3-2    2018-01-03 [1] CRAN (R 4.0.2)
##  RSQLite        2.2.0    2020-01-07 [1] CRAN (R 4.0.2)
##  sandwich       2.5-1    2019-04-06 [1] CRAN (R 4.0.2)
##  scales         1.1.1    2020-05-11 [1] CRAN (R 4.0.2)
##  sessioninfo    1.1.1    2018-11-05 [1] CRAN (R 4.0.2)
##  SparseM        1.78     2019-12-13 [1] CRAN (R 4.0.2)
##  stringi        1.4.6    2020-02-17 [1] CRAN (R 4.0.2)
##  stringr        1.4.0    2019-02-10 [1] CRAN (R 4.0.2)
##  survival       3.2-3    2020-06-13 [1] CRAN (R 4.0.2)
##  testthat       2.3.2    2020-03-02 [1] CRAN (R 4.0.2)
##  tibble         3.0.3    2020-07-10 [1] CRAN (R 4.0.2)
##  tidyr          1.1.1    2020-07-31 [1] CRAN (R 4.0.2)
##  tidyselect     1.1.0    2020-05-11 [1] CRAN (R 4.0.2)
##  usethis        1.6.1    2020-04-29 [1] CRAN (R 4.0.2)
##  vctrs          0.3.2    2020-07-15 [1] CRAN (R 4.0.2)
##  withr          2.2.0    2020-04-20 [1] CRAN (R 4.0.2)
##  xfun           0.16     2020-07-24 [1] CRAN (R 4.0.2)
##  yaml           2.2.1    2020-02-01 [1] CRAN (R 4.0.2)
##  zoo            1.8-8    2020-05-02 [1] CRAN (R 4.0.2)
## 
## [1] /Library/Frameworks/R.framework/Versions/4.0/Resources/library