diff --git a/NAMESPACE b/NAMESPACE index 0ae59be..0d3d4dd 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -5,7 +5,7 @@ export(colocboost_plot) export(get_cormat) export(get_cos) export(get_cos_summary) -export(get_strong_colocalization) +export(get_robust_colocalization) importFrom(grDevices,adjustcolor) importFrom(graphics,abline) importFrom(graphics,axis) diff --git a/R/colocboost.R b/R/colocboost.R index e88ba5d..d2b0313 100644 --- a/R/colocboost.R +++ b/R/colocboost.R @@ -179,7 +179,7 @@ colocboost <- function(X = NULL, Y = NULL, # individual data LD_free = FALSE, output_level = 1) { ###################### ---- one module for data object - message("Starting checking the input data.") + message("Validating input data.") # - check if all missing check_individual <- (is.null(X) & is.null(Y)) check_sumstat <- (is.null(sumstat) & (is.null(effect_est) & is.null(effect_se))) @@ -604,7 +604,7 @@ colocboost <- function(X = NULL, Y = NULL, # individual data ) # --- post-processing of the colocboost updates - message("Starting assemble analyses and results summary.") + message("Performing inference on colocalization events.") cb_output <- colocboost_assemble(cb_obj, coverage = coverage, weight_fudge_factor = weight_fudge_factor, diff --git a/R/colocboost_output.R b/R/colocboost_output.R index eda9e2a..2073949 100644 --- a/R/colocboost_output.R +++ b/R/colocboost_output.R @@ -129,11 +129,11 @@ get_cos_summary <- function(cb_output, } -#' @rdname get_strong_colocalization +#' @rdname get_robust_colocalization #' -#' @title Get colocalization summary table from a ColocBoost output. +#' @title Recalibrate and summarize robust colocalization events. #' -#' @description `get_strong_colocalization` get the colocalization by discarding the weaker colocalization events or colocalized outcomes +#' @description `get_robust_colocalization` get the colocalization by discarding the weaker colocalization events or colocalized outcomes #' #' @param cb_output Output object from `colocboost` analysis #' @param cos_npc_cutoff Minimum threshold of normalized probability of colocalization (NPC) for CoS. @@ -171,12 +171,12 @@ get_cos_summary <- function(cb_output, #' } #' res <- colocboost(X = X, Y = Y) #' res$cos_details$cos$cos_index -#' filter_res <- get_strong_colocalization(res, cos_npc_cutoff = 0.5, npc_outcome_cutoff = 0.2) +#' filter_res <- get_robust_colocalization(res, cos_npc_cutoff = 0.5, npc_outcome_cutoff = 0.2) #' filter_res$cos_details$cos$cos_index #' #' @family colocboost_inference #' @export -get_strong_colocalization <- function(cb_output, +get_robust_colocalization <- function(cb_output, cos_npc_cutoff = 0.5, npc_outcome_cutoff = 0.2, pvalue_cutoff = NULL, @@ -458,7 +458,7 @@ get_ucos_summary <- function(cb_output, outcome_names = NULL, region_name = NULL return(summary_table) } -#' Extract CoS simply change the coverage without checking purity +#' Extract CoS at different coverages, without filtering by purity #' #' @description `get_cos` get the colocalization confidence sets (CoS) with different coverage. #' @@ -793,6 +793,10 @@ get_model_info <- function(cb_obj, outcome_names = NULL) { n_updates <- cb_obj$cb_model_para$num_updates model_coveraged <- cb_obj$cb_model_para$coveraged jk_update <- cb_obj$cb_model_para$real_update_jk + if (!is.null(jk_update)){ + rownames(jk_update) <- paste0("jk_star_", 1:nrow(jk_update)) + colnames(jk_update) <- outcome_names + } outcome_proximity_obj <- lapply(cb_obj$cb_model, function(cb) cb$obj_path) outcome_coupled_best_update_obj <- lapply(cb_obj$cb_model, function(cb) cb$obj_single) outcome_profile_loglik <- lapply(cb_obj$cb_model, function(cb) cb$profile_loglike_each) @@ -805,7 +809,7 @@ get_model_info <- function(cb_obj, outcome_names = NULL) { "outcome_profile_loglik" = outcome_profile_loglik, "outcome_proximity_obj" = outcome_proximity_obj, "outcome_coupled_best_update_obj" = outcome_coupled_best_update_obj, - "jk_update" = jk_update + "jk_star" = jk_update ) return(ll) } diff --git a/R/colocboost_workhorse.R b/R/colocboost_workhorse.R index 8be177e..6303024 100644 --- a/R/colocboost_workhorse.R +++ b/R/colocboost_workhorse.R @@ -80,7 +80,7 @@ colocboost_workhorse <- function(cb_data, # - if all outcomes do not have signals, STOP message(paste0( "Using multiple testing correction method: ", func_multi_test, - ". Stop ColocBoost since no outcomes ", focal_outcome_idx, "have association signals." + ". Stop ColocBoost since no outcomes have association signals." )) } else { message(paste0( @@ -92,7 +92,7 @@ colocboost_workhorse <- function(cb_data, } if (!is.null(focal_outcome_idx) & (M != 1)) { if (sum(cb_model_para$true_stop == focal_outcome_idx) != 0) { - message(paste("Stop ColocBoost since the focal outcome", focal_outcome_idx, "do not have association signals.")) + warning(paste("Stop ColocBoost since the focal outcome", focal_outcome_idx, "do not have association signals.")) cb_model_para$update_y <- 0 } } @@ -106,7 +106,7 @@ colocboost_workhorse <- function(cb_data, } if (M == 1) { # single effect with or without LD matrix - message("Running colocboost with assumption of one causal per outcome!") + message("Running ColocBoost with assumption of one causal per outcome per region!") cb_obj <- colocboost_one_causal(cb_model, cb_model_para, cb_data) cb_obj$cb_model_para$coveraged <- "one_causal" } else { @@ -133,10 +133,10 @@ colocboost_workhorse <- function(cb_data, if (length(pos_rtr_stop) != 0) { cb_model_para$update_y[pos.update[pos_rtr_stop]] <- 0 message(paste( - "Boosting iterations for outcome", paste(pos.update[pos_rtr_stop], collapse = ", "), + "Gradient boosting for outcome", paste(pos.update[pos_rtr_stop], collapse = ", "), "stop since rtr < 0 or max(correlation) > 1 after", m, "iterations!", "Results for this locus are not stable, please check if mismatch between sumstat and LD!", - "See details in website." + "See details in tutorial website." )) } @@ -166,14 +166,11 @@ colocboost_workhorse <- function(cb_data, } } } - # if (!is.null(focal_outcome_idx)){ - # if (!is.na(stop[focal_outcome_idx]) & stop[focal_outcome_idx]){ stop = TRUE } - # } if (all(length(stop) == 1 & stop)) { cb_model_para$update_y <- 0 if (cb_model_para$L == 1) { - message(paste("Boosting iterations for outcome 1 converge after", m, "iterations!")) + message(paste("Gradient boosting for outcome 1 converge after", m, "iterations!")) } } else { if (all(!stop[pos.update])) { @@ -194,25 +191,25 @@ colocboost_workhorse <- function(cb_data, if (!is.null(focal_outcome_idx)) { if (focal_outcome_idx %in% cb_model_para$true_stop) { message(paste( - "Boosting iterations for focal outcome", focal_outcome_idx, - "converge after", m, "iterations!" + "Gradient boosting for focal outcome", focal_outcome_idx, + "converged after", m, "iterations!" )) if (length(setdiff(cb_model_para$true_stop, focal_outcome_idx)) != 0) { message(paste( - "Boosting iterations for outcome", paste(setdiff(cb_model_para$true_stop, focal_outcome_idx), collapse = ", "), - "converge after", m, "iterations!" + "Gradient boosting for outcome", paste(setdiff(cb_model_para$true_stop, focal_outcome_idx), collapse = ", "), + "converged after", m, "iterations!" )) } } else { message(paste( - "Boosting iterations for outcome", paste(cb_model_para$true_stop, collapse = ", "), - "converge after", m, "iterations!" + "Gradient boosting for outcome", paste(cb_model_para$true_stop, collapse = ", "), + "converged after", m, "iterations!" )) } } else { message(paste( - "Boosting iterations for outcome", paste(cb_model_para$true_stop, collapse = ", "), - "converge after", m, "iterations!" + "Gradient boosting for outcome", paste(cb_model_para$true_stop, collapse = ", "), + "converged after", m, "iterations!" )) } } @@ -220,7 +217,7 @@ colocboost_workhorse <- function(cb_data, } } if (m %% 1000 == 0) { - message(paste("Boosting at", m, "iterations, still updating.")) + message(paste("Gradient boosting at", m, "iterations, still updating.")) } } @@ -239,7 +236,7 @@ colocboost_workhorse <- function(cb_data, ####### --------------------------------------------- # calculate objective function of Y for the last iteration. cb_model <- boost_obj_last(cb_data, cb_model, cb_model_para) - warning(paste("COLOC-BOOST updates did not converge in", M, "iterations; checkpoint at last iteration")) + warning(paste("ColocBoost updates did not converge in", M, "iterations; checkpoint at last iteration")) cb_model_para$coveraged <- FALSE } diff --git a/man/get_cos.Rd b/man/get_cos.Rd index b4a94d0..6b53e23 100644 --- a/man/get_cos.Rd +++ b/man/get_cos.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/colocboost_output.R \name{get_cos} \alias{get_cos} -\title{Extract CoS simply change the coverage without checking purity} +\title{Extract CoS at different coverages, without filtering by purity} \usage{ get_cos(cb_output, coverage = 0.95) } diff --git a/man/get_cos_summary.Rd b/man/get_cos_summary.Rd index 73e4cf9..88e9525 100644 --- a/man/get_cos_summary.Rd +++ b/man/get_cos_summary.Rd @@ -63,6 +63,6 @@ get_cos_summary(res) } \seealso{ Other colocboost_inference: -\code{\link{get_strong_colocalization}()} +\code{\link{get_robust_colocalization}()} } \concept{colocboost_inference} diff --git a/man/get_strong_colocalization.Rd b/man/get_robust_colocalization.Rd similarity index 88% rename from man/get_strong_colocalization.Rd rename to man/get_robust_colocalization.Rd index 5bc2c93..d5d9b68 100644 --- a/man/get_strong_colocalization.Rd +++ b/man/get_robust_colocalization.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/colocboost_output.R -\name{get_strong_colocalization} -\alias{get_strong_colocalization} -\title{Get colocalization summary table from a ColocBoost output.} +\name{get_robust_colocalization} +\alias{get_robust_colocalization} +\title{Recalibrate and summarize robust colocalization events.} \usage{ -get_strong_colocalization( +get_robust_colocalization( cb_output, cos_npc_cutoff = 0.5, npc_outcome_cutoff = 0.2, @@ -36,7 +36,7 @@ A \code{"colocboost"} object with some or all of the following elements: \item{model_info}{A object with detailed information for colocboost model} } \description{ -\code{get_strong_colocalization} get the colocalization by discarding the weaker colocalization events or colocalized outcomes +\code{get_robust_colocalization} get the colocalization by discarding the weaker colocalization events or colocalized outcomes } \examples{ # colocboost example @@ -59,7 +59,7 @@ for (l in 1:L) { } res <- colocboost(X = X, Y = Y) res$cos_details$cos$cos_index -filter_res <- get_strong_colocalization(res, cos_npc_cutoff = 0.5, npc_outcome_cutoff = 0.2) +filter_res <- get_robust_colocalization(res, cos_npc_cutoff = 0.5, npc_outcome_cutoff = 0.2) filter_res$cos_details$cos$cos_index } diff --git a/tests/testthat/test_inference.R b/tests/testthat/test_inference.R index 8c92f1d..72bb0f1 100644 --- a/tests/testthat/test_inference.R +++ b/tests/testthat/test_inference.R @@ -45,7 +45,6 @@ generate_test_result <- function(n = 100, p = 20, L = 2, seed = 42) { # Test colocboost_plot function test_that("colocboost_plot handles different plot options", { - skip_on_cran() # Generate a test colocboost results cb_res <- generate_test_result() @@ -62,8 +61,7 @@ test_that("colocboost_plot handles different plot options", { # Test get_cos_summary function test_that("get_cos_summary handles different parameters", { - skip_on_cran() - + # Generate a test colocboost results cb_res <- generate_test_result() @@ -85,8 +83,7 @@ test_that("get_cos_summary handles different parameters", { # Test for get_strong_colocalization test_that("get_strong_colocalization filters results correctly", { - skip_on_cran() - + # Generate a test colocboost results cb_res <- generate_test_result() diff --git a/vignettes/Interpret_ColocBoost_Output.Rmd b/vignettes/Interpret_ColocBoost_Output.Rmd index 9c7f2cf..fdf29b2 100644 --- a/vignettes/Interpret_ColocBoost_Output.Rmd +++ b/vignettes/Interpret_ColocBoost_Output.Rmd @@ -76,14 +76,16 @@ In `cos_summary`, for each 95% CoS, the `cos_npc` column provides a normalized p `min_npc_outcome` column provides the minimum normalized probability among colocalized traits. Those two metrices are measured as an empirical evidence of colocalization both in CoS-level and in trait-level. To obtain the best minimal colocalization configuration can be defined by using both `cos_npc` and `npc_outcome` -See detailed usage of this function in [link](https://statfungen.github.io/colocboost/reference/get_strong_colocalization.html). +See detailed usage of this function in [link](https://statfungen.github.io/colocboost/reference/get_robust_colocalization.html). ```{r run-strong-colocalization} -filter_res <- get_strong_colocalization(res, cos_npc_cutoff = 0.5, npc_outcome_cutoff = 0.2) +filter_res <- get_robust_colocalization(res, cos_npc_cutoff = 0.5, npc_outcome_cutoff = 0.2) ``` -The output from `get_strong_colocalization` is the same as output from `colocboost`, which can be directly used in any post inference and visualization. +- The output from `get_robust_colocalization` is the same as output from `colocboost`, which can be directly used in any post inference and visualization. +- `npc=0.5` or `npc_outcome = 0.2` maintains robust colocalization signals for cases when many traits are evaluated. +Higher thresholds can be specified if users only want to focus on strong colocalization events. ## 3. More details on ColocBoost output diff --git a/vignettes/Visualization_ColocBoost_Output.Rmd b/vignettes/Visualization_ColocBoost_Output.Rmd index 779e534..2551c71 100644 --- a/vignettes/Visualization_ColocBoost_Output.Rmd +++ b/vignettes/Visualization_ColocBoost_Output.Rmd @@ -83,6 +83,7 @@ colocboost_plot(res, show_top_variables = TRUE) ### 2.3. Plot CoS variants to uncolocalized traits to diagnostic the colocalization. There are three options available for plotting the CoS variants to uncolocalized traits: + - `show_cos_to_uncoloc = FALSE` (default), if `TRUE` will plot all CoS variants to all uncolocalized traits. - `show_cos_to_uncoloc_idx = NULL` (default), if specified, will plot the specified CoS variants to all uncolocalized traits. - `show_cos_to_uncoloc_outcome = NULL` (default), if specified, will plot the all CoS variants to the specified uncolocalized traits.