diff --git a/.Rbuildignore b/.Rbuildignore index a296df3..bd1d928 100755 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -10,6 +10,7 @@ ^docs$ ^tools$ ^doc$ +^SplineOmics\.BiocCheck$ ^Meta$ ^docker$ ^dev$ diff --git a/R/cluster_hits.R b/R/cluster_hits.R index 5a198ae..d4e2904 100755 --- a/R/cluster_hits.R +++ b/R/cluster_hits.R @@ -193,6 +193,9 @@ cluster_hits <- function( genes = genes, spline_params = spline_params, adj_pthresholds = adj_pthresholds, + adj_pthresh_avrg_diff_conditions = adj_pthresh_avrg_diff_conditions, + adj_pthresh_interaction_condition_time = + adj_pthresh_interaction_condition_time, report_dir = report_dir, mode = mode, report_info = report_info, @@ -438,6 +441,8 @@ perform_clustering <- function( #' @param adj_pthresholds Numeric vector, containing a float < 1 > 0 as each #' value. There is one float for every level, and this is #' the adj. p-value threshold. +#' @param adj_pthresh_avrg_diff_conditions Float +#' @param adj_pthresh_interaction_condition_time Float #' @param report_dir A character string specifying the report directory. #' @param mode A character string specifying the mode #' ('isolated' or 'integrated'). @@ -489,6 +494,8 @@ make_clustering_report <- function( genes, spline_params, adj_pthresholds, + adj_pthresh_avrg_diff_conditions, + adj_pthresh_interaction_condition_time, report_dir, mode, report_info, @@ -727,6 +734,9 @@ make_clustering_report <- function( topTables = topTables, enrichr_format = enrichr_format, adj_pthresholds = adj_pthresholds, + adj_pthresh_avrg_diff_conditions = adj_pthresh_avrg_diff_conditions, + adj_pthresh_interaction_condition_time = + adj_pthresh_interaction_condition_time, report_type = "cluster_hits", feature_name_columns = feature_name_columns, mode = mode, @@ -1194,7 +1204,7 @@ remove_batch_effect_cluster_hits <- function( level <- unique_levels[level_index] level_columns <- which(meta[[condition]] == level) data_copy <- data[, level_columns] - meta_copy <- subset(meta, meta[[condition]] == level) + meta_copy <- meta[meta[[condition]] == level, , drop = FALSE] } else { data_copy <- data # Take the full data for mode == "integrated" meta_copy <- meta @@ -1229,7 +1239,7 @@ remove_batch_effect_cluster_hits <- function( if (mode == "isolated") { level <- unique_levels[level_index] - meta_copy <- subset(meta, meta[[condition]] == level) + meta_copy <- meta[meta[[condition]] == level, , drop = FALSE] if (!is.null(meta_batch2_column) && length(unique(meta_copy[[meta_batch2_column]])) > 1) { @@ -1243,10 +1253,16 @@ remove_batch_effect_cluster_hits <- function( } } - data_copy <- do.call(limma::removeBatchEffect, args) + data_copy <- do.call( + limma::removeBatchEffect, + args + ) # For mode == "integrated", all elements are identical - datas <- c(datas, list(data_copy)) + datas <- c( + datas, + list(data_copy) + ) } } else { # no meta batch column specified, just return right data @@ -2071,9 +2087,13 @@ plot_spline_comparisons <- function( adj_pthresh_interaction ) { - # Sort and prepare data + # Sort and prepare data (sorting based on feature name for easy navigation) time_effect_1 <- time_effect_1 |> dplyr::arrange(.data$feature_names) time_effect_2 <- time_effect_2 |> dplyr::arrange(.data$feature_names) + avrg_diff_conditions <- + avrg_diff_conditions |> dplyr::arrange(.data$feature_names) + interaction_condition_time <- + interaction_condition_time |> dplyr::arrange(.data$feature_names) # Get relevant parameters DoF <- which(names(time_effect_1) == "AveExpr") - 1 @@ -2350,6 +2370,8 @@ prepare_gene_lists_for_enrichr <- function( #' @param level_headers_info A list of header information for each level. #' @param spline_params A list of spline parameters. #' @param adj_pthresholds Float vector with values for any level for adj.p.tresh +#' @param adj_pthresh_avrg_diff_conditions Float +#' @param adj_pthresh_interaction_condition_time Float #' @param mode A character string specifying the mode #' ('isolated' or 'integrated'). #' @param report_info A named list containg the report info fields. Here used @@ -2370,6 +2392,8 @@ build_cluster_hits_report <- function( level_headers_info, spline_params, adj_pthresholds, + adj_pthresh_avrg_diff_conditions, + adj_pthresh_interaction_condition_time, mode, report_info, output_file_path @@ -2578,11 +2602,9 @@ build_cluster_hits_report <- function( pb$tick() } - + # Add sections for limma_result_2_and_3_plots - if (!is.null(limma_result_2_and_3_plots)) { - adj_pthresh_avrg_diff_conditions <- 0.05 - adj_pthresh_interaction_condition_time <- 0.05 + if (length(limma_result_2_and_3_plots) > 0) { # Create a new main header for the limma result plots header_index <- header_index + 1 diff --git a/R/run_limma_splines.R b/R/run_limma_splines.R index 66481c9..77ccb40 100755 --- a/R/run_limma_splines.R +++ b/R/run_limma_splines.R @@ -247,6 +247,7 @@ between_level <- function( condition, compared_levels[2] ) + condition_only <- limma::topTable( fit, coef = factor_only_contrast_coeff, @@ -258,6 +259,7 @@ between_level <- function( top_table = condition_only, fit = fit ) + top_table_condition_only <- process_top_table( condition_only_resuls, feature_names @@ -277,7 +279,7 @@ between_level <- function( ":X", seq_len(num_matching_columns) ) - + condition_time <- limma::topTable( fit, coef = factor_time_contrast_coeffs, @@ -342,11 +344,11 @@ within_level <- function( padjust_method, mode ) { - + if (mode == "isolated") { samples <- which(meta[[condition]] == level) data_copy <- data[, samples] - meta_copy <- subset(meta, meta[[condition]] == level) + meta_copy <- meta[meta[[condition]] == level, , drop = FALSE] } else { # mode == "integrated" data_copy <- data meta_copy <- meta diff --git a/R/utils_input_validation.R b/R/utils_input_validation.R index 21e942c..b4984ca 100755 --- a/R/utils_input_validation.R +++ b/R/utils_input_validation.R @@ -422,10 +422,17 @@ InputControl <- R6::R6Class("InputControl", check_design_formula = function() { formula <- self$args[["design"]] - meta <- self$args$meta - meta_index <- self$args$meta_index + meta <- self$args[["meta"]] + meta_index <- self$args[["meta_index"]] + + # Not strictly required + meta_batch_column <- self$args[["meta_batch_column"]] + meta_batch2_column <- self$args[["meta_batch2_column"]] - required_args <- list(formula, meta) + required_args <- list( + formula, + meta + ) if (any(sapply(required_args, is.null))) { return(NULL) @@ -475,16 +482,43 @@ InputControl <- R6::R6Class("InputControl", missing_columns <- setdiff(formula_terms, names(meta)) if (length(missing_columns) > 0) { if (!is.null(meta_index)) { - stop(sprintf("%s (data/meta pair index: %s): %s", + stop_call_false(sprintf("%s (data/meta pair index: %s): %s", "The following design columns are missing in meta", meta_index, - paste(missing_columns, collapse = ", ")), - call. = FALSE) + paste(missing_columns, collapse = ", "))) } else { - stop(paste("The following design columns are missing in meta:", - paste(missing_columns, collapse = ", ")), - call. = FALSE) + stop_call_false(paste("The following design columns are missing in meta:", + paste(missing_columns, collapse = ", "))) + } + } + + # Convert formula to string for regex checking + formula_str <- as.character(formula) + + # Check if batch column is provided and validate its presence in the formula + if (!is.null(meta_batch_column)) { + if (!grepl(meta_batch_column, formula_str)) { + stop_call_false( + paste("The batch effect column", meta_batch_column, + "is provided but not present in the design formula. ", + "Please ensure that if you specify a batch column, ", + "it is included in the design formula to", + "remove batch effects.") + ) + } + } + + # Check if the second batch column is provided and validate its presence + if (!is.null(meta_batch2_column)) { + if (!grepl(meta_batch2_column, formula_str)) { + stop_call_false( + paste("The second batch effect column", meta_batch2_column, + "is provided but not present in the design formula.", + "Please ensure that if you specify a second batch column,", + "it is included in the design formula", + "to remove batch effects.") + ) } } diff --git a/R/utils_report_generation.R b/R/utils_report_generation.R index 1332ee8..48c4f10 100755 --- a/R/utils_report_generation.R +++ b/R/utils_report_generation.R @@ -30,6 +30,8 @@ #' @param spline_params A list of spline parameters, such as dof and type. #' @param adj_pthresholds Numeric vector with the values for the adj.p.tresholds #' for each level. +#' @param adj_pthresh_avrg_diff_conditions Float, only for cluster_hits() +#' @param adj_pthresh_interaction_condition_time Float, only for cluster_hits() #' @param report_type A character string specifying the report type #' ('screen_limma_hyperparams' or 'cluster_hits'). #' @param feature_name_columns Character vector with the column names of the @@ -69,6 +71,8 @@ generate_report_html <- function( level_headers_info = NA, spline_params = NA, adj_pthresholds = NA, + adj_pthresh_avrg_diff_conditions = NA, # only for cluster_hits() + adj_pthresh_interaction_condition_time = NA, # only for cluster_hits() report_type = "explore_data", feature_name_columns = NA, # only for cluster_hits() mode = NA, @@ -79,7 +83,7 @@ generate_report_html <- function( ) { feature_names_formula <- "" - + if (report_type == "explore_data") { if (filename == "explore_data") { title <- "explore data" @@ -142,6 +146,7 @@ generate_report_html <- function( "meta_condition", "meta_batch", "limma_design", + "mode", "analyst_name", "contact_info", "project_name", @@ -187,7 +192,9 @@ generate_report_html <- function( '