Merge pull request #123 from xueweic/main

Update about post-filtering of uCoS

Author: gaow

NAMESPACE

@@ -10,6 +10,7 @@ export(get_cos_purity)
 export(get_cos_summary)
 export(get_hierarchical_clusters)
 export(get_robust_colocalization)
+export(get_robust_ucos)
 export(get_ucos_summary)
 importFrom(grDevices,adjustcolor)
 importFrom(graphics,abline)

R/colocboost.R

@@ -210,7 +210,12 @@ colocboost <- function(X = NULL, Y = NULL, # individual data
     dict_YX = dict_YX, dict_sumstatLD = dict_sumstatLD,
     effect_est = effect_est, effect_se = effect_se, effect_n = effect_n,
     overlap_variables = overlap_variables,
-    M = M, min_abs_corr = min_abs_corr
+    M = M, min_abs_corr = min_abs_corr,
+    jk_equiv_corr = jk_equiv_corr, 
+    jk_equiv_loglik = jk_equiv_loglik,
+    func_simplex = func_simplex,
+    cos_npc_cutoff = cos_npc_cutoff,
+    npc_outcome_cutoff = npc_outcome_cutoff
   )
   if (is.null(validated_data)) {
     return(NULL)
@@ -235,6 +240,12 @@ colocboost <- function(X = NULL, Y = NULL, # individual data
   jk_equiv_corr <- validated_data$jk_equiv_corr
   jk_equiv_loglik <- validated_data$jk_equiv_loglik
   func_simplex <- validated_data$func_simplex
+  cos_npc_cutoff <- validated_data$cos_npc_cutoff
+  npc_outcome_cutoff <- validated_data$npc_outcome_cutoff
+  if (M == 1){
+      cos_npc_cutoff <- 0
+      npc_outcome_cutoff <- 0
+  }
 
   # - initial colocboost object
   keep_variables <- c(keep_variable_individual, keep_variable_sumstat)
@@ -338,6 +349,20 @@ colocboost <- function(X = NULL, Y = NULL, # individual data
     weight_fudge_factor = weight_fudge_factor,
     coverage = coverage
   )
+  # ---- post filtering of the colocboost results (get robust trait-specific events)
+  if ("ucos_details" %in% names(cb_output)) {
+    if (is.null(pvalue_cutoff)){
+        pvalue_cutoff_ucos <- NULL
+    } else {
+        # only keep the suggestive significant results
+        pvalue_cutoff_ucos <- ifelse(pvalue_cutoff > 1e-5, 1e-5, pvalue_cutoff)
+    }
+    cb_output <- get_robust_ucos(
+        cb_output, 
+        npc_outcome_cutoff = npc_outcome_cutoff, 
+        pvalue_cutoff = pvalue_cutoff_ucos
+    )
+  }
 
   return(cb_output)
 }
@@ -387,7 +412,11 @@ colocboost_validate_input_data <- function(X = NULL, Y = NULL,
                                            dict_YX = NULL, dict_sumstatLD = NULL,
                                            effect_est = NULL, effect_se = NULL, effect_n = NULL,
                                            overlap_variables = FALSE,
-                                           M = 500, min_abs_corr = 0.5) {
+                                           M = 500, min_abs_corr = 0.5,
+                                           jk_equiv_corr = 0.8, jk_equiv_loglik = 1,
+                                           func_simplex = "LD_z2z", 
+                                           cos_npc_cutoff = 0.2,
+                                           npc_outcome_cutoff = 0.2) {
 
   # - check individual level data
   if (!is.null(X) & !is.null(Y)) {
@@ -624,9 +653,11 @@ colocboost_validate_input_data <- function(X = NULL, Y = NULL,
     # --- check input of LD
     M_updated <- M
     min_abs_corr_updated <- min_abs_corr
-    jk_equiv_corr_updated <- 0.8
-    jk_equiv_loglik_updated <- 1
-    func_simplex_updated <- "LD_z2z"
+    jk_equiv_corr_updated <- jk_equiv_corr
+    jk_equiv_loglik_updated <- jk_equiv_loglik
+    func_simplex_updated <- func_simplex
+    cos_npc_cutoff_updated <- cos_npc_cutoff
+    npc_outcome_cutoff_updated <- npc_outcome_cutoff
 
     if (is.null(LD)) {
       # if no LD input, set diagonal matrix to LD
@@ -644,6 +675,8 @@ colocboost_validate_input_data <- function(X = NULL, Y = NULL,
       jk_equiv_corr_updated <- 0
       jk_equiv_loglik_updated <- 0.1
       func_simplex_updated <- "only_z2z"
+      cos_npc_cutoff_updated <- 0
+      npc_outcome_cutoff_updated <- 0
 
     } else {
 
@@ -804,9 +837,11 @@ colocboost_validate_input_data <- function(X = NULL, Y = NULL,
     Z <- N_sumstat <- Var_y <- SeBhat <- sumstatLD_dict <- keep_variable_sumstat <- NULL
     M_updated <- M
     min_abs_corr_updated <- min_abs_corr
-    jk_equiv_corr_updated <- 0.8
-    jk_equiv_loglik_updated <- 1
-    func_simplex_updated <- "LD_z2z"
+    jk_equiv_corr_updated <- jk_equiv_corr
+    jk_equiv_loglik_updated <- jk_equiv_loglik
+    func_simplex_updated <- func_simplex
+    cos_npc_cutoff_updated = cos_npc_cutoff
+    npc_outcome_cutoff_updated = npc_outcome_cutoff
   }
 
   return(list(
@@ -826,7 +861,9 @@ colocboost_validate_input_data <- function(X = NULL, Y = NULL,
     min_abs_corr = min_abs_corr_updated,
     jk_equiv_corr = jk_equiv_corr_updated,
     jk_equiv_loglik = jk_equiv_loglik_updated,
-    func_simplex = func_simplex_updated
+    func_simplex = func_simplex_updated,
+    cos_npc_cutoff = cos_npc_cutoff_updated,
+    npc_outcome_cutoff = npc_outcome_cutoff_updated
   ))
 }
 

R/colocboost_inference.R

@@ -269,7 +269,7 @@ check_null_post <- function(cb_obj,
     ld_feature <- sqrt(abs(ld_jk))
     # - calculate delta
     delta <- boost_KL_delta(
-      z = z, ld_feature = ld_feature, adj_dep = adj_dep,
+      z = z, ld_feature = ld_feature, adj_dep = adj_dep, 
       func_simplex = func_simplex, lambda = lambda
     )
     scaling_factor <- if (!is.null(N)) (N - 1) else 1
@@ -575,7 +575,6 @@ get_cos_evidence <- function(cb_obj, coloc_out, data_info) {
   }
 
   avWeight <- coloc_out$avWeight
-  cs_change <- coloc_out$cs_change
   check_null_max <- sapply(cb_obj$cb_model, function(cb) cb$check_null_max)
   outcome_names <- data_info$outcome_info$outcome_names
   n_cos <- length(avWeight)
@@ -603,3 +602,99 @@ get_cos_evidence <- function(cb_obj, coloc_out, data_info) {
   return(list(normalization_evidence = normalization_evidence, npc = npc))
 }
 
+
+
+#' Function to get the evidence of trait-specific ucos
+#' @keywords cb_post_inference
+#' @noRd
+#' @importFrom stats var
+#' @importFrom utils tail
+get_ucos_evidence <- function(cb_obj, ucoloc_info) {
+    
+  get_ucos_profile <- function(cs_beta, outcome_idx, X = NULL, Y = NULL, N = NULL,
+                              XtX = NULL, YtY = NULL, XtY = NULL, miss_idx = NULL, adj_dep = 1) {
+    if (!is.null(X)) {
+      cos_profile <- mean((Y - X %*% as.matrix(cs_beta))^2) * N / (N - 1)
+      yty <- var(Y)
+    } else if (!is.null(XtY)) {
+      scaling_factor <- if (!is.null(N)) (N - 1) else 1
+      beta_scaling <- if (!is.null(N)) 1 else 100
+      cs_beta <- cs_beta / beta_scaling
+      yty <- YtY / scaling_factor
+      xtx <- XtX
+      if (length(miss_idx) != 0) {
+        xty <- XtY[-miss_idx] / scaling_factor
+        cs_beta <- cs_beta[-miss_idx]
+      } else {
+        xty <- XtY / scaling_factor
+      }
+      if (length(xtx) == 1){
+        cos_profile <- (yty - 2 * sum(cs_beta * xty) + sum(cs_beta^2)) * adj_dep
+      } else {
+        cos_profile <- (yty - 2 * sum(cs_beta * xty) + sum((xtx %*% as.matrix(cs_beta)) * cs_beta)) * adj_dep
+      }
+    }
+    delta <- yty - cos_profile
+    if (delta <= 0) {
+      warning(paste(
+        "Warning message: potential sumstat & LD mismatch may happens for outcome", outcome_idx,
+        ". Using logLR  = uCoS(profile) - max(profile). Please check our website https://statfungen.github.io/colocboost/articles/."
+      ))
+    }
+    cos_profile
+  }
+
+  get_outcome_profile_change <- function(outcome_idx, ucos, cb_obj) {
+    extract_last <- function(lst) {
+      tail(lst, n = 1)
+    }
+    cb_data <- cb_obj$cb_data
+    cs_beta <- rep(0, cb_obj$cb_model_para$P)
+    cs_beta[ucos] <- cb_obj$cb_model[[outcome_idx]]$beta[ucos]
+    X_dict <- cb_data$dict[outcome_idx]
+    cos_profile <- get_ucos_profile(cs_beta, outcome_idx,
+      X = cb_data$data[[X_dict]]$X, Y = cb_data$data[[outcome_idx]]$Y,
+      XtX = cb_data$data[[X_dict]]$XtX, XtY = cb_data$data[[outcome_idx]]$XtY,
+      YtY = cb_data$data[[outcome_idx]]$YtY, N = cb_data$data[[outcome_idx]]$N,
+      miss_idx = cb_data$data[[outcome_idx]]$variable_miss,
+      adj_dep = cb_data$data[[outcome_idx]]$dependency
+    )
+    max_profile <- max(cb_obj$cb_model[[outcome_idx]]$profile_loglike_each)
+    ifelse(max_profile < cos_profile, 0, max_profile - cos_profile)
+  }
+
+  # - Calculate best configuration likelihood explained by minimal configuration
+  get_normalization_evidence <- function(profile_change, null_max, outcomes, outcome_names) {
+    # Define the exponential likelihood ratio normalization (ELRN)
+    logLR_normalization <- function(ratio) {
+      1 - exp(-2 * ratio)
+    }
+
+    ratio <- profile_change / null_max
+    prob <- logLR_normalization(ratio)
+    df <- data.frame(outcome = outcome_names, outcomes_index = outcomes, relative_logLR = ratio, npc_outcome = prob)
+    return(df)
+  }
+
+  check_null_max_ucos <- sapply(cb_obj$cb_model, function(cb) cb$check_null_max_ucos)
+  n_ucos <- length(ucoloc_info$ucos)
+  normalization_evidence <- list()
+  for (i in 1:n_ucos) {
+    outcome_idx <- ucoloc_info$outcome[[i]]
+    outcome_name <- ucoloc_info$outcome_name[[i]]
+    # most likely cos
+    ucos <- ucoloc_info$ucos[[i]]
+    profile_change_outcome <- get_outcome_profile_change(outcome_idx, ucos, cb_obj)
+    normalization_evidence[[i]] <- get_normalization_evidence(
+      profile_change = profile_change_outcome,
+      null_max = check_null_max_ucos[outcome_idx],
+      outcome_idx, outcome_name
+    )
+  }
+  normalization_evidence <- do.call(rbind, normalization_evidence)
+  rownames(normalization_evidence) <- names(ucoloc_info$ucos)
+  return(normalization_evidence)
+}
+
+
+