Merge pull request #18 from leaffur/working

pre-filtering improvement

Author: leaffur

CASE_R_package.Rproj

@@ -1,4 +1,5 @@
 Version: 1.0
+ProjectId: 2241ace3-378b-4e48-bb0a-ef8f1b97ac60
 
 RestoreWorkspace: Default
 SaveWorkspace: Default

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: CASE
 Title: Cell-type-specific And Shared EQTL fine-mapping
-Version: 0.3.0
+Version: 0.3.1
 Authors@R: c(
     person("Chen", "Lin", , "c.lin@yale.edu", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0001-9821-2578")),
@@ -18,7 +18,7 @@ Imports:
     mvtnorm,
     stats
 Encoding: UTF-8
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 Suggests: 
     knitr,
     rmarkdown,

R/CASE.R

@@ -81,14 +81,30 @@ CASE <- function(Z = NULL, R, hatB = NULL, hatS = NULL, N, V = NULL, cs = TRUE,
     hatS = hatBS$hatS
 
     # V = estimate_null_correlation_simple(mash_set_data(do.call(rbind, raw.data$hatB), do.call(rbind, raw.data$hatS)))
-
-    m1 <- CASE_train(hatB = hatB, hatS = hatS, V = V, R = R, N = N, Z = NULL, verbose = verbose, ...)
+    sump = 2 - 2 * pnorm(abs(hatB / hatS))
+    sumfdr = apply(sump, 2, function(x) p.adjust(x, method = "fdr"))
+    C = ncol(sumfdr)
+    
+    # no strong signals in all cell types, or low SNP numbers 
+    if (sum(sumfdr <= 0.2) == 0){
+      if (verbose){
+        cat("No FDR-significant variants in the inputs.", "\n")
+      }
+      m1 = list(pi = 1, n.iter = 0)
+    } else if (nrow(sump) <= 20){
+      if (verbose){
+        cat("Too few SNPs in the data (<=20).", "\n")
+      }
+      m1 = list(pi = 1, n.iter = 0)
+    } else{
+      m1 <- CASE_train(hatB = hatB, hatS = hatS, V = V, R = R, N = N, Z = NULL, verbose = verbose, ...)
+    }
 
     res <- CASE_test(hatB = hatB, hatS = hatS, R = R, N = N, CASE_training = m1, Z = NULL, verbose = verbose, ...)
     t2 = Sys.time()
     res$time = difftime(t2, t1, units = "secs")
     if (cs){
-      res$sets <- get_credible_sets(res$pip, R = R, verbose = verbose)
+      res$sets <- get_credible_sets(pips = res$pip, R = R, verbose = verbose)
     }
 
     return(res)

R/CASE_models.R

@@ -80,10 +80,17 @@ CASE_train <- function(Z = NULL, R, hatB = NULL, hatS = NULL, N, V = NULL, verbo
   ## Train with only marginally significant SNPs
   M0 = nrow(R)
   ME_p <- 2 - 2 * pnorm(abs(hatB / hatS), 0, 1)
-  significant_thres = ifelse("significant_thres" %in% names(args), args$significant_thres, 1e-1)
-  idx <- which(ME_p <= significant_thres, arr.ind = TRUE)[, 1] %>% unique
+  significant_thres = ifelse("significant_thres" %in% names(args), 
+                             args$significant_thres, 
+                             ifelse(sqrt(6) / sqrt(diag(N)) >= 0.05, sqrt(6) / sqrt(diag(N)), 0.05))
+  if (C == 1){
+    idx <- which(ME_p <= significant_thres)
+  } else{
+    idx <- which(t(apply(ME_p, 1, function(x) x <= significant_thres)), arr.ind = TRUE)[, 1] %>% unique
+  }
+  
   if (length(idx) == 1){
-    idx = c(idx - 1, idx, idx + 1)
+    idx = idx + (-2):2
     idx = idx[idx >= 1 & idx <= nrow(hatB)]
   }
   hatB = hatB[idx, ]
@@ -96,7 +103,7 @@ CASE_train <- function(Z = NULL, R, hatB = NULL, hatS = NULL, N, V = NULL, verbo
     if (verbose){
       cat("No marginally significant variants in the inputs.", "\n")
     }
-    return(list(pi = 1, U = list(matrix(0, C, C)), V = V, n.iter = 0))
+    return(list(pi = 1, n.iter = 0))
   }
 
   J <- 0
@@ -310,10 +317,6 @@ CASE_test <- function(Z = NULL, R, hatB = NULL, hatS = NULL, N, CASE_training, v
     cat("Start Posterior Analysis.", "\n")
   }
 
-  U = CASE_training$U
-  V = CASE_training$V
-  pi = CASE_training$pi
-  
   if (is.null(Z)){
     Z = hatB / hatS
   }
@@ -322,15 +325,19 @@ CASE_test <- function(Z = NULL, R, hatB = NULL, hatS = NULL, N, CASE_training, v
   hatBS = transform_Z(Z, N)
   hatB = hatBS$hatB
   hatS = hatBS$hatS
-
+  
   C <- ncol(hatB)
   M <- nrow(R)
 
+  pi = CASE_training$pi
   if (length(pi) <= 1){
     post_mean = pip = matrix(0, M, C)
-    return(list(pi = pi, U = U, V = V, pip = pip, post_mean = post_mean))
+    return(list(pip = pip, post_mean = post_mean))
   }
 
+  U = CASE_training$U
+  V = CASE_training$V
+
   L = length(U)
   ## MC step
   gBc = gB_coef(U, V)