Add pvalue cutoff option

R/colocboost.R

@@ -157,6 +157,7 @@ colocboost <- function(X = NULL, Y = NULL, # individual data
                        multicorrection_cut = 1,
                        ash_prior = "normal", # only applicable if func_multicorrection = lfsr
                        p.adjust.methods = "fdr",
+                       pv_cutoff = 1e-4, 
                        check_null = 0.1, # the cut off value for change conditional objective function
                        check_null_method = "profile",
                        check_null_max = 0.02,
@@ -565,6 +566,7 @@ colocboost <- function(X = NULL, Y = NULL, # individual data
                                      coverage = coverage,
                                      func_intw = func_intw,
                                      alpha = alpha,
+                                     pv_cutoff = pv_cutoff,
                                      check_null = check_null,
                                      check_null_method = check_null_method,
                                      check_null_max = check_null_max,

R/colocboost_assemble.R

@@ -26,6 +26,7 @@ colocboost_assemble <- function(cb_obj,
                                 check_null = 0.1,
                                 check_null_method = "profile",
                                 check_null_max=2e-5,
+                                pv_cutoff = 1e-4,
                                 dedup = TRUE,
                                 overlap = TRUE,
                                 n_purity = 100,
@@ -193,6 +194,9 @@ colocboost_assemble <- function(cb_obj,
         cb_obj$cb_model_para$alpha <- alpha
         cb_obj$cb_model_para$coverage <- coverage
         cos_results <- get_cos_details(cb_obj, coloc_out = past_out$cos$cos, data_info = data_info)
+        if (!is.null(pv_cutoff) & !is.null(cos_results$cos_results)){
+            cos_results <- cos_pvalue_filter(cos_results, data_info = data_info, pv_cutoff = pv_cutoff)
+        }
         cb_output <- list("vcp" = cos_results$vcp,
                           "cos_details" = cos_results$cos_results,
                           "data_info" = data_info,

R/colocboost_init.R

@@ -117,7 +117,11 @@ colocboost_init_data <- function(X, Y, dict_YX,
             }
             if (change_x){
                 dict_YX_final[i] == i
-                x_stand <- Rfast::standardise(x_tmp, center = intercept, scale = standardize)
+                if (!intercept&!standardize){
+                  x_stand <- x_tmp
+                } else {
+                  x_stand <- Rfast::standardise(x_tmp, center = intercept, scale = standardize)
+                }
                 x_stand[which(is.na(x_stand))] = 0
                 tmp$X <- x_stand
             }

R/colocboost_output.R

@@ -75,7 +75,7 @@ get_data_info <- function(cb_obj){
 
 #' @noRd
 #' @keywords cb_post_inference
-get_cos_details <- function(cb_obj, coloc_out, data_info = NULL){
+get_cos_details <- function(cb_obj, coloc_out, data_info = NULL, pv_cutoff = 1e-4){
 
     if (is.null(data_info))
         data_info <- get_data_info(cb_obj)
@@ -560,3 +560,86 @@ get_summary_table_fm <- function(cb_output, outcome_names = NULL, gene_name = NU
 
 }
 
+
+
+cos_pvalue_filter <- function(cos_results, data_info = NULL, pv_cutoff = 1e-4){
+
+  if (is.null(data_info))
+    data_info <- get_data_info(cb_obj)
+
+  n_cos <- length(cos_results$cos_results$cos$cos_index)
+  filtered_cos <- list("cos_index" = NULL,
+                       "cos_variables" = NULL)
+  filtered_traits <- list("outcome_index" = NULL,
+                          "outcome_name" = NULL)
+  filtered_cos_vcp <- list()
+  pp_remove <- c()
+  for (i in 1:n_cos){
+    cos_tmp <- cos_results$cos_results$cos$cos_index[[i]]
+    cos_trait <- cos_results$cos_results$cos_outcomes$outcome_index[[i]]
+    minPV <- sapply(cos_trait, function(tmp){
+      z <- data_info$z[[tmp]]
+      z <- z[cos_tmp]
+      pv <- pchisq(z^2, 1, lower.tail = FALSE)
+      min(pv)
+    })
+    pp <- which(minPV < 1e-4)
+    if (length(pp) == 0 | length(pp) == 1) {
+      pp_remove <- c(pp_remove, i)
+      next
+    }
+    filtered_cos$cos_index <- c(filtered_cos$cos_index, 
+                                cos_results$cos_results$cos$cos_index[i])
+    filtered_cos$cos_variables <- c(filtered_cos$cos_variables, 
+                                    cos_results$cos_results$cos$cos_variables[i])
+    filtered_cos_vcp <- c(filtered_cos_vcp, cos_results$cos_results$cos_vcp[i])
+    if (length(pp) == length(cos_trait)){
+      filtered_traits$outcome_index <- c(filtered_traits$outcome_index, 
+                                         cos_results$cos_results$cos_outcomes$outcome_index[i])
+      filtered_traits$outcome_name <- c(filtered_traits$outcome_name, 
+                                        cos_results$cos_results$cos_outcomes$outcome_name[i])
+
+    } else {
+      filtered_traits$outcome_index <- c(filtered_traits$outcome_index, 
+                                         list(cos_results$cos_results$cos_outcomes$outcome_index[[i]][pp]))
+      filtered_traits$outcome_name <- c(filtered_traits$outcome_name, 
+                                        list(cos_results$cos_results$cos_outcomes$outcome_name[[i]][pp]))
+    }
+  }
+  if (length(filtered_cos$cos_index)==0){
+    cos_results <- list("cos_results" = NULL, "vcp" = NULL)
+  } else {
+    cos_results$cos_results$cos <- filtered_cos
+    names(filtered_traits$outcome_index) <- names(cos_results$cos_results$cos$cos_index)
+    names(filtered_traits$outcome_name) <- names(cos_results$cos_results$cos$cos_index)
+    cos_results$cos_results$cos_outcomes <- filtered_traits
+    cos_results$cos_results$cos_vcp <- filtered_cos_vcp
+    coloc_hits <- coloc_hits_variablenames <- coloc_hits_names <- c()
+    for (i in 1:length(filtered_cos_vcp)){
+      inw <- filtered_cos_vcp[[i]]
+      pp1 <- which(inw == max(inw))
+      coloc_hits <- c(coloc_hits, pp1)
+      coloc_hits_variablenames <- c(coloc_hits_variablenames, data_info$variables[pp1])
+      if (length(pp1)==1){
+        coloc_hits_names <- c(coloc_hits_names, names(filtered_cos_vcp)[i])
+      } else {
+        coloc_hits_names <- c(coloc_hits_names, paste0(names(filtered_cos_vcp)[i], ".", 1:length(pp1)))
+      }
+    }
+    coloc_hits <- data.frame("top_index" = coloc_hits, "top_variables" = coloc_hits_variablenames) 
+    rownames(coloc_hits) <- coloc_hits_names
+    cos_results$cos_results$cos_top_variables <- coloc_hits
+
+    # - change purity
+    if (length(pp_remove)!=0){
+      cos_results$cos_results$cos_purity$min_abs_cor <- as.matrix(cos_results$cos_results$cos_purity$min_abs_cor)[-pp_remove, -pp_remove, drop=FALSE]
+      cos_results$cos_results$cos_purity$median_abs_cor <- as.matrix(cos_results$cos_results$cos_purity$median_abs_cor)[-pp_remove, -pp_remove, drop=FALSE]
+      cos_results$cos_results$cos_purity$max_abs_cor <- as.matrix(cos_results$cos_results$cos_purity$max_abs_cor)[-pp_remove, -pp_remove, drop=FALSE]
+    }
+
+  }
+  return(cos_results)
+
+}
+
+

R/colocboost_utils.R

@@ -59,7 +59,7 @@ check_null_post <- function(cb_obj,
             if (length(miss_idx)!=0){
                 xty <- XtY[-miss_idx] / scaling_factor
                 cs_beta <- cs_beta[-miss_idx]
-            } else { xty <- XtY  }
+            } else { xty <- XtY  / scaling_factor  }
 
             yty - 2*sum(cs_beta*xty) + sum( (xtx %*% as.matrix(cs_beta)) * cs_beta )
         }
@@ -174,8 +174,6 @@ check_null_post <- function(cb_obj,
     } 
     cs_change <- as.data.frame(cs_change)
     is_non_null <- which(rowSums( (check_cs_change >= check_null) * max_change ) != 0)
-    # is_non_null <- which(rowSums( (cs_change >= check_null) * (max_change >= check_null_max) ) != 0)
-    # is_non_null <- which(rowSums(cs_change >= check_null) != 0)
 
     ll = list("cs_change" = cs_change,
               "is_non_null" = is_non_null)