diff --git a/R/colocboost.R b/R/colocboost.R index 91890f1..19904ae 100644 --- a/R/colocboost.R +++ b/R/colocboost.R @@ -157,7 +157,6 @@ colocboost <- function(X = NULL, Y = NULL, # individual data multicorrection_cut = 1, ash_prior = "normal", # only applicable if func_multicorrection = lfsr p.adjust.methods = "fdr", - npc_cutoff = 0.7, check_null = 0.1, # the cut off value for change conditional objective function check_null_method = "profile", check_null_max = 0.02, @@ -568,7 +567,6 @@ colocboost <- function(X = NULL, Y = NULL, # individual data coverage = coverage, func_intw = func_intw, alpha = alpha, - npc_cutoff = npc_cutoff, check_null = check_null, check_null_method = check_null_method, check_null_max = check_null_max, diff --git a/R/colocboost_assemble.R b/R/colocboost_assemble.R index 9de7f39..4677859 100644 --- a/R/colocboost_assemble.R +++ b/R/colocboost_assemble.R @@ -26,7 +26,6 @@ colocboost_assemble <- function(cb_obj, check_null = 0.1, check_null_method = "profile", check_null_max=2e-5, - npc_cutoff = 0.7, dedup = TRUE, overlap = TRUE, n_purity = 100, @@ -193,7 +192,7 @@ colocboost_assemble <- function(cb_obj, # - colocalization results 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, npc_cutoff = npc_cutoff) + cos_results <- get_cos_details(cb_obj, coloc_out = past_out$cos$cos, data_info = data_info) cb_output <- list("vcp" = cos_results$vcp, "cos_details" = cos_results$cos_results, "data_info" = data_info, diff --git a/R/colocboost_output.R b/R/colocboost_output.R index 7f2899f..8a84e9b 100644 --- a/R/colocboost_output.R +++ b/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, npc_cutoff = 0.7){ +get_cos_details <- function(cb_obj, coloc_out, data_info = NULL){ if (is.null(data_info)) data_info <- get_data_info(cb_obj) @@ -86,7 +86,7 @@ get_cos_details <- function(cb_obj, coloc_out, data_info = NULL, npc_cutoff = 0. if (length(coloc_sets)!=0){ # - colocalization outcome configurations - tmp <- get_cos_evidence(cb_obj, coloc_out, data_info, npc_cutoff = npc_cutoff) + tmp <- get_cos_evidence(cb_obj, coloc_out, data_info) normalization_evidence <- tmp$normalization_evidence npc <- tmp$npc diff --git a/R/colocboost_utils.R b/R/colocboost_utils.R index 000bea0..4b5731b 100644 --- a/R/colocboost_utils.R +++ b/R/colocboost_utils.R @@ -339,7 +339,7 @@ get_between_purity = function (pos1, pos2, X=NULL, Xcorr=NULL, N = NULL, miss_id } -get_cos_evidence <- function(cb_obj, coloc_out, data_info, npc_cutoff = 0.7){ +get_cos_evidence <- function(cb_obj, coloc_out, data_info){ get_cos_config <- function(w, config_idx, alpha = 1.5, coverage = 0.95){ @@ -395,7 +395,7 @@ get_cos_evidence <- function(cb_obj, coloc_out, data_info, npc_cutoff = 0.7){ } # - calculate best configuration likelihood explained by minimal configuration - get_normalization_evidence <- function(profile_change, null_max, outcomes, outcome_names, npc_cutoff = 0.7) { + 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 ) } @@ -405,10 +405,16 @@ get_cos_evidence <- function(cb_obj, coloc_out, data_info, npc_cutoff = 0.7){ rownames(df) <- outcome_names[outcomes] sorted_df <- df[order(-df$relative_logLR), ] sorted_df$strong_effect <- sorted_df$relative_logLR >= 1 - sorted_df$suggestive_effect <- sorted_df$npc_outcome >= npc_cutoff return(sorted_df) } + get_npuc <- function(npc_outcome){ + # npuc_outcome <- sapply(1:length(npc_outcome), function(i) npc_outcome[i]*prod(1-npc_outcome[-i]) ) + # 1 - prod(1-npuc_outcome) + max_idx <- which.max(npc_outcome) + npc_outcome[max_idx] * prod(1 - npc_outcome[-max_idx]) + } + avWeight <- coloc_out$avWeight cs_change <- coloc_out$cs_change check_null_max <- sapply(cb_obj$cb_model, function(cb) cb$check_null_max) @@ -426,12 +432,10 @@ get_cos_evidence <- function(cb_obj, coloc_out, data_info, npc_cutoff = 0.7){ }) normalization_evidence[[i]] <- get_normalization_evidence(profile_change = profile_change_outcome, null_max = check_null_max[outcomes], - outcomes, outcome_names, npc_cutoff = npc_cutoff) + outcomes, outcome_names) # - calcualte CoS colocalization probability (CCP) and CoS uncolocalization probability (UCCP) - npc_outcome <- normalization_evidence[[i]]$npc_outcome - max_idx <- which.max(npc_outcome) - npuc <- npc_outcome[max_idx] * prod(1 - npc_outcome[-max_idx]) + npuc <- get_npuc(normalization_evidence[[i]]$npc_outcome) npc[i] <- 1 - npuc } names(normalization_evidence) <- names(npc) <- names(coloc_out$cos)