Important update about post-filtering

R/colocboost.R

@@ -89,6 +89,9 @@
 #' @param output_level When \code{output_level = 1}, return basic cos details for colocalization results
 #'                     When \code{output_level = 2}, return the ucos details for the single specific effects.
 #'                     When \code{output_level = 3}, return the entire Colocboost model to diagnostic results (more space).
+#' @param cos_npc_cutoff Minimum threshold of normalized probability of colocalization (NPC) for CoS.
+#' @param npc_outcome_cutoff Minimum threshold of normalized probability of colocalized traits in each CoS.
+#' @param pvalue_cutoff Maximum threshold of marginal p-values of colocalized variants on colocalized traits in each CoS.
 #'
 #' @return A \code{"colocboost"} object with some or all of the following elements:
 #'
@@ -184,7 +187,12 @@ colocboost <- function(X = NULL, Y = NULL, # individual data
                        check_null_max_ucos = 0.015, # the smallest value of change of profile loglikelihood for each outcome in uCoS.
                        weaker_effect = TRUE,
                        LD_free = FALSE,
-                       output_level = 1) {
+                       output_level = 1,
+                       ###### - Post filtering parameters
+                       cos_npc_cutoff = 0.2, # remove the CoS with cos_npc less than this cutoff
+                       npc_outcome_cutoff = 0.2, # remove the colocalized outcome in CoS if npc_outcome less than this cutoff
+                       pvalue_cutoff = 1e-3 # remove the colocalized outcome in CoS if pvalue greater than this cutoff
+                       ) {
   ###################### ---- one module for data object
   message("Validating input data.")
   # - check if all missing
@@ -320,6 +328,17 @@ colocboost <- function(X = NULL, Y = NULL, # individual data
     output_level = output_level
   )
   class(cb_output) <- "colocboost"
+
+  # ---- post filtering of the colocboost results (get robust colocalization events)
+  cb_output <- get_robust_colocalization(
+    cb_output = cb_output,
+    cos_npc_cutoff = cos_npc_cutoff,
+    npc_outcome_cutoff = npc_outcome_cutoff,
+    pvalue_cutoff = pvalue_cutoff,
+    weight_fudge_factor = weight_fudge_factor,
+    coverage = coverage
+  )
+
   return(cb_output)
 }
 
@@ -646,20 +665,26 @@ colocboost_validate_input_data <- function(X = NULL, Y = NULL,
           warning('Error: Please provide dict_sumstatLD: you have ', length(sumstat), 
                   ' sumstats but only ', length(LD), ' LD matrices')
           return(NULL)
+        } else {
+          # - dict for sumstat to LD mapping
+          sumstatLD_dict <- rep(NA, length(sumstat))
+          for (i in 1:length(sumstat)) {
+            tmp <- unique(dict_sumstatLD[dict_sumstatLD[, 1] == i, 2])
+            if (length(tmp) == 0) {
+              warning(paste("Error: You don't provide matched LD for sumstat", i))
+              return(NULL)
+            } else if (length(tmp) != 1) {
+              warning(paste("Error: You provide multiple matched LD for sumstat", i))
+              return(NULL)
+            } else {
+              sumstatLD_dict[i] <- tmp
+            }
+          }
+          if (max(sumstatLD_dict) > length(LD)) {
+            warning("Error: You don't provide enough LD matrices!")
+            return(NULL)
+          }
         }
-        if (length(dict_sumstatLD) != length(sumstat)) {
-          warning('Error: dict_sumstatLD must have length ', length(sumstat))
-          return(NULL)
-        }
-        if (any(is.na(dict_sumstatLD))) {
-          warning('Error: dict_sumstatLD contains NA values')
-          return(NULL)
-        }
-        if (any(dict_sumstatLD < 1) || any(dict_sumstatLD > length(LD))) {
-          warning('Error: dict_sumstatLD values must be between 1 and ', length(LD))
-          return(NULL)
-        }
-        sumstatLD_dict <- as.integer(dict_sumstatLD)
       }
 
       # === Filter variants for each sumstat ===

tests/testthat/test_utils.R

@@ -19,12 +19,12 @@ generate_test_result <- function(n = 100, p = 20, L = 2, seed = 42) {
 
   if (L == 1) {
     # Single trait case
-    true_beta[5, 1] <- 0.7  # SNP5 affects the trait
-    true_beta[10, 1] <- 0.6 # SNP10 also affects the trait
+    true_beta[5, 1] <- 1  # SNP5 affects the trait
+    true_beta[10, 1] <- 1 # SNP10 also affects the trait
   } else {
     # Multi-trait case
-    true_beta[5, 1] <- 0.7  # SNP5 affects trait 1
-    true_beta[5, 2] <- 0.6  # SNP5 also affects trait 2 (colocalized)
+    true_beta[5, 1] <- 1  # SNP5 affects trait 1
+    true_beta[5, 2] <- 1  # SNP5 also affects trait 2 (colocalized)
     true_beta[10, 2] <- 0.5 # SNP10 only affects trait 2
   }
 

vignettes/announcements.Rmd

@@ -14,7 +14,8 @@ vignette: >
 - *May 2, 2025*: `colocboost` R package is available on [CRAN](https://CRAN.R-project.org/package=colocboost).
 
 ## Software updates
-- `v1.0.7` Improvements to ColocBoost (check out the full details in [PR](https://github.com/StatFunGen/colocboost/pull/116)). 
+- `v1.0.7` (**Important update**) Improvements to ColocBoost (check out the full details in [PR](https://github.com/StatFunGen/colocboost/pull/116) and [PR](https://github.com/StatFunGen/colocboost/pull/121)). 
+  - Enhanced `colocboost` main function with post-filtering and only keep the robust colocalization events.
   - Enhanced `colocboost_plot` function with flexible highlighting options and new visualization styles.
   - Optimized performance and computational efficiency
   - Improved documentation and examples for the wrapper pipeline