cTWAS Simulation Power: N=226K

Simulation Info

#number of samples (original)
print(n.ori)

[1] 2e+05

#number of samples after filtering
print(n)

[1] 225582

#number of SNPs
print(p)

[1] 6228138

#number of genes
print(J)

[1] 8021

Parameter Estimation

configtag <- 1
runtag = "ukb-s400.226-adi"

simutags <- paste(1, 1:5, sep = "-")
plot_par(configtag, runtag, simutags)

simulations  1-1 1-2 1-3 1-4 1-5 : mean gene PVE: 0.04556472 , mean SNP PVE: 0.4569035 

Past versions of unnamed-chunk-6-1.png

PIP Calibration

plot_PIP(configtag, runtag, simutags)

Past versions of unnamed-chunk-7-1.png

Number of Causal and Non-Causal Genes in PIP Bins

phenofs <- paste0(outputdir, "ukb-s400.226-adi", "_simu", simutags, "-pheno.Rd")
susieIfs <- paste0(outputdir, runtag, "_simu",simutags, "_config", configtag,".susieIrss.txt")

cau <- lapply(phenofs, function(x){load(x);get_causal_id(phenores)})

pipfs <- susieIfs

df <- NULL
for (i in 1:length(pipfs)) {
  res <- fread(pipfs[i], header = T)
  res <- data.frame(res[res$type  =="gene", ])
  res$ifcausal <- ifelse(res$id %in% cau[[i]], 1, 0)
  res$runtag <- i
  res <- res[complete.cases(res),]
  df <- rbind(df, res)
}

bins_start <- 0:9/10
bins_end <- bins_start + 0.1
n_causal <- rep(NA, length(bins_start))
n_noncausal <- rep(NA, length(bins_start))

for (i in 1:length(bins_start)){
  n_causal[i] <- sum(df$susie_pip >= bins_start[i] & df$susie_pip <= bins_end[i] & df$ifcausal==1)
  n_noncausal[i] <- sum(df$susie_pip >= bins_start[i] & df$susie_pip <= bins_end[i] & df$ifcausal==0)
}

n_causal <- n_causal / length(simutags)
n_noncausal <- n_noncausal / length(simutags)

#number of causal and non-causal genes in PIP bins
rbind(bins_start, bins_end, n_causal, n_noncausal)

              [,1]  [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
bins_start     0.0   0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8   0.9
bins_end       0.1   0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9   1.0
n_causal      38.6   7.2  3.8  2.2  3.2  3.0  4.8  7.2  6.8  53.0
n_noncausal 7686.2 237.4 39.2  9.6  5.8  3.8  3.0  3.2  1.2   1.8

PIP Calibration

nca_plot <- function(pips, ifcausal, runtag = NULL, mode = c("PIP", "FDR"), xmin =0, main = mode[1], ...){
  # ifcausal:0,1, runtag: for adding std.
  if (is.null(runtag)){
    runtag <- rep(1, length(pips))
  }

  if (mode == "PIP"){
    bins <- seq(0, 1, by=0.1)[1:10]
  } else if (mode == "FDR"){
    bins <- c(0, 0.01, 0.05, 0.1, seq(0.2, 1, by=0.1))[1:12]
  }

  calist <- list()
  nonlist <- list()
  for (rt in unique(runtag)){
    pips.rt <- pips[runtag == rt]
    ifcausal.rt <- ifcausal[runtag == rt]
    res <- .obn(pips.rt, ifcausal.rt, mode = mode)
    calist[[rt]] <- cbind(res[["ncausal"]], "causal", bins)
    nonlist[[rt]] <- cbind(res[["nnoncausal"]], "noncausal", bins)
  }

  df <- rbind(do.call(rbind, calist),
          do.call(rbind, nonlist))
  
  #df[df[,2]=="noncausal",2] <- "Non-causal"
  
  
  df <- data.frame("count"= as.numeric(df[,1]),
                   "ifcausal" = factor(df[,2], levels = c("noncausal", "causal")),
                   "bins" = as.numeric(df[,3]))

  if (mode == "PIP"){
    ymax <- 1.1* (max(df[df$bins > xmin & df$ifcausal == "causal", "count"], na.rm = T) + max(df[df$bins > xmin & df$ifcausal == "noncausal", "count"], na.rm = T))
  } else {
    ymax <- 1.1* (max(df[df$bins < xmin & df$ifcausal == "causal", "count"], na.rm = T) + max(df[df$bins < xmin & df$ifcausal == "noncausal", "count"], na.rm = T))
  }
  
  df <- df[df$bins>=0.5,]
  
  levels(df$ifcausal) <- c("Non-causal", "Causal")

  fig <- ggbarplot(df, x = "bins", y = "count", add = "mean_se", fill = "ifcausal", 
                   palette=rev(c("#F8766D","#00BFC4")),
                   ylim=c(0,ymax), 
                   main="",
                   ylab="Gene counts",
                   xlab="PIP bins"
                   )
  fig <- ggpar(fig, legend.title="")
  return(fig)
}

ncausal_plot <- function(phenofs, pipfs, main = "PIP"){

  cau <- lapply(phenofs, function(x) {load(x);get_causal_id(phenores)})

  df <- NULL
  for (i in 1:length(pipfs)) {
    res <- fread(pipfs[i], header = T)
    res <- data.frame(res[res$type  =="gene", ])
    res$ifcausal <- ifelse(res$id %in% cau[[i]], 1, 0)
    res$runtag <- i
    res <- res[complete.cases(res),]
    df <- rbind(df, res)
  }

  fig <- nca_plot(df$susie_pip, df$ifcausal, df$runtag, mode ="PIP", xmin = 0.5, main = main)
  
  return(fig)
}

phenofs <- paste0(outputdir, "ukb-s400.226-adi", "_simu", simutags, "-pheno.Rd")
susieIfs <- paste0(outputdir, runtag, "_simu",simutags, "_config", configtag,".susieIrss.txt")

ncausal_plot(phenofs, susieIfs) 

Past versions of unnamed-chunk-10-1.png

# pdf(file = "./panel_e.pdf", width=3.5, height=3.5)
# 
# ncausal_plot(phenofs, susieIfs) 
# 
# dev.off()

Session information

sessionInfo()

R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
 [1] ggpubr_0.4.0      plotrix_3.7-6     cowplot_1.0.0    
 [4] stringr_1.4.0     plyr_1.8.4        tidyr_1.1.0      
 [7] plotly_4.9.0      ggplot2_3.3.3     data.table_1.14.0
[10] ctwas_0.1.29     

loaded via a namespace (and not attached):
 [1] httr_1.4.1        jsonlite_1.6      viridisLite_0.3.0
 [4] foreach_1.5.1     pgenlibr_0.3.1    carData_3.0-2    
 [7] logging_0.10-108  cellranger_1.1.0  yaml_2.2.0       
[10] pillar_1.6.1      backports_1.1.4   lattice_0.20-38  
[13] glue_1.4.2        digest_0.6.20     promises_1.0.1   
[16] ggsignif_0.5.0    colorspace_1.4-1  htmltools_0.3.6  
[19] httpuv_1.5.1      Matrix_1.2-18     pkgconfig_2.0.3  
[22] broom_0.7.9       haven_2.3.1       purrr_0.3.4      
[25] scales_1.1.0      whisker_0.3-2     openxlsx_4.1.0.1 
[28] later_0.8.0       rio_0.5.16        git2r_0.26.1     
[31] tibble_3.1.2      farver_2.1.0      generics_0.0.2   
[34] car_3.0-5         ellipsis_0.3.2    withr_2.4.1      
[37] lazyeval_0.2.2    magrittr_2.0.1    crayon_1.4.1     
[40] readxl_1.3.1      evaluate_0.14     fs_1.3.1         
[43] fansi_0.5.0       rstatix_0.7.0     forcats_0.4.0    
[46] foreign_0.8-71    tools_3.6.1       hms_1.1.0        
[49] lifecycle_1.0.0   munsell_0.5.0     ggsci_2.9        
[52] zip_2.0.3         compiler_3.6.1    rlang_0.4.11     
[55] grid_3.6.1        iterators_1.0.13  htmlwidgets_1.3  
[58] labeling_0.3      rmarkdown_1.13    gtable_0.3.0     
[61] codetools_0.2-16  abind_1.4-5       DBI_1.1.1        
[64] curl_3.3          R6_2.5.0          gridExtra_2.3    
[67] knitr_1.23        dplyr_1.0.7       utf8_1.2.1       
[70] workflowr_1.6.2   rprojroot_2.0.2   stringi_1.4.3    
[73] Rcpp_1.0.6        vctrs_0.3.8       tidyselect_1.1.0 
[76] xfun_0.8