diff --git a/.gitignore b/.gitignore index c691af8..b514a05 100644 --- a/.gitignore +++ b/.gitignore @@ -18,6 +18,7 @@ methylTFR.BiocCheck/* .Trashes .fseventsd .TemporaryItems +methylTFR_tmp/ # R session artefacts .Rhistory diff --git a/DESCRIPTION b/DESCRIPTION index 3cdb343..0e491d1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: methylTFR Title: Quantification of DNA Methylation Signatures in TFBS -Version: 0.99.4 +Version: 0.99.6 Date: 2026-08-19 Authors@R: c( person("Irem B.", "Gündüz", , "irembgunduz@gmail.com", role = c("aut", "cre"), @@ -27,6 +27,7 @@ Depends: SummarizedExperiment Imports: BiocGenerics, + BiocParallel, DelayedArray, GenomicRanges, ggplot2, @@ -35,7 +36,6 @@ Imports: logger, matrixStats, methods, - parallel, R.utils, S4Vectors, stats, @@ -49,6 +49,7 @@ Suggests: knitr, RefManageR, RnBeads, + RnBeads.hg38, rmarkdown, sessioninfo, testthat (>= 3.0.0) diff --git a/NAMESPACE b/NAMESPACE index 64ca8d4..1a5cdee 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -22,6 +22,10 @@ import(logger) importClassesFrom(SummarizedExperiment,SummarizedExperiment) importFrom(BiocGenerics,cbind) importFrom(BiocGenerics,rbind) +importFrom(BiocParallel,MulticoreParam) +importFrom(BiocParallel,SerialParam) +importFrom(BiocParallel,SnowParam) +importFrom(BiocParallel,bplapply) importFrom(DelayedArray,ArbitraryArrayGrid) importFrom(DelayedArray,DelayedArray) importFrom(DelayedArray,close) @@ -57,6 +61,7 @@ importFrom(ggplot2,xlim) importFrom(ggplot2,ylab) importFrom(logger,log_error) importFrom(logger,log_info) +importFrom(logger,log_success) importFrom(logger,log_warn) importFrom(matrixStats,colMads) importFrom(matrixStats,colMeans2) @@ -68,7 +73,6 @@ importFrom(methods,as) importFrom(methods,is) importFrom(methods,new) importFrom(methods,setMethod) -importFrom(parallel,mclapply) importFrom(stats,aggregate) importFrom(stats,aov) importFrom(stats,kruskal.test) diff --git a/NEWS.md b/NEWS.md index 8091345..b6fa53f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,4 @@ -# methylTFR 0.99.4 +# methylTFR 0.99.6 NEW FEATURES diff --git a/R/compute_deviations.R b/R/compute_deviations.R index a659636..c3a0370 100644 --- a/R/compute_deviations.R +++ b/R/compute_deviations.R @@ -1,3 +1,35 @@ +#' @title check_deviation_inputs +#' @description Validate the inputs of \code{computeDeviation}. +#' @param motif Motif name as a character string. +#' @param msites Methylation sites as a \code{GRanges} object. +#' @param tf_bindsites a \code{GRangesList} of TF binding site positions. +#' @param enhancer a \code{GRanges} restricting the analysis (optional). +#' @return Invisible \code{NULL}. Called for the errors it raises. +#' @importFrom methods is +#' @keywords internal +check_deviation_inputs <- function( + motif, msites, tf_bindsites, enhancer = NULL +) { + if (is.null(motif) || !is.character(motif)) { + stop("Please provide a valid motif name") + } + if (is.null(msites) || !is(msites, "GRanges")) { + stop( + "Please provide a valid methylation sites with ", + "read_methylome function" + ) + } + if (is.null(tf_bindsites) || + !any(c(!is(tf_bindsites, "GRangesList") || + !is.list(tf_bindsites)))) { + stop("Please provide a valid tf binding sites as GRangesList") + } + if (!is.null(enhancer) && !is(enhancer, "GRanges")) { + stop("Please provide a valid enhancer regions") + } + invisible(NULL) +} + #' @title computeDeviation #' @description computeDeviation is a function to calculate #' the deviation in transcription factor @@ -55,31 +87,18 @@ #' ) #' @export computeDeviation <- function( - motif, msites, tf_bindsites, gcfreqs, - enhancer = NULL, ignoreStrand = TRUE, - binMsites + motif, msites, tf_bindsites, gcfreqs, enhancer = NULL, + ignoreStrand = TRUE, binMsites ) { if (!is.logical(ignoreStrand)) { warning("Found invalid strand option, using the default") ignoreStrand <- TRUE } - if (is.null(motif) || !is.character(motif)) { - stop("Please provide a valid motif name") - } - if (is.null(msites) || !is(msites, "GRanges")) { - stop("Please provide a valid methylation - sites with read_methylome function") - } - if (is.null(tf_bindsites) || - !any(c(!is(tf_bindsites, "GRangesList") || - !is.list(tf_bindsites)))) { - stop("Please provide a valid tf binding sites as GRangesList") - } - if (!is.null(enhancer) && !is(enhancer, "GRanges")) { - stop("Please provide a valid enhancer regions") - } + check_deviation_inputs(motif, msites, tf_bindsites, enhancer) tfbs <- tf_bindsites[[motif]] - tfbs <- resize(tfbs, width(tfbs)[1] + 130, fix = "center") + tfbs <- resize(tfbs, width(tfbs)[1] + 130, + fix = "center" + ) gcfreq <- gcfreqs[[motif]] if (!is.null(enhancer)) { tfbs <- subsetByOverlaps(tfbs, enhancer, diff --git a/R/differential_analysis.R b/R/differential_analysis.R index 92ebab9..acdec53 100644 --- a/R/differential_analysis.R +++ b/R/differential_analysis.R @@ -1,3 +1,76 @@ +#' @title resolve_diff_groups +#' @description Derive and validate the group labels used by +#' \code{differential_deviation_test}. +#' @param deviations A matrix of deviation scores, motifs in rows. +#' @param groups Group labels, or NULL to take them from the column names. +#' @return The group labels as a \code{factor}. +#' @keywords internal +resolve_diff_groups <- function(deviations, groups) { + if (is.null(groups)) { + groups <- colnames(deviations) + } + if (is.null(groups)) { + stop( + "No group labels found. Provide 'groups', or supply deviations ", + "with column names identifying the groups." + ) + } + groups <- as.factor(groups) + if (length(groups) != ncol(deviations)) { + stop("'groups' must have one entry per column of 'deviations'") + } + if (nlevels(groups) < 2) { + stop("'groups' must contain at least two distinct groups") + } + return(groups) +} + +#' @title diff_pvalues +#' @description Test every motif for a difference between the groups, with +#' the test chosen from the number of groups and \code{parametric}. +#' @param deviations A matrix of deviation scores, motifs in rows. +#' @param groups The group labels as a \code{factor}. +#' @param parametric if TRUE, use a t-test or ANOVA, otherwise a Wilcoxon +#' or Kruskal-Wallis test. +#' @param alternative The alternative hypothesis of the two-group tests. +#' @return A numeric vector of p-values, one per motif. +#' @keywords internal +diff_pvalues <- function(deviations, groups, parametric, alternative) { + if (parametric) { + if (nlevels(groups) == 2) { + # t-test + return(apply(deviations, 1, t_helper, groups, alternative)) + } + # anova + return(apply(deviations, 1, anova_helper, groups)) + } + if (nlevels(groups) == 2) { + # wilcoxon + return(apply(deviations, 1, wilcoxon_helper, groups, alternative)) + } + # kruskal-wallis + return(apply(deviations, 1, kw_helper, groups)) +} + +#' @title group_mean_difference +#' @description Unsigned effect size: the difference of the group means for +#' two groups, and their range for more. +#' @param deviations A matrix of deviation scores, motifs in rows. +#' @param groups The group labels as a \code{factor}. +#' @return A numeric vector with one value per motif. +#' @keywords internal +group_mean_difference <- function(deviations, groups) { + group_means <- vapply( + levels(groups), + function(g) rowMeans(deviations[, groups == g, drop = FALSE]), + numeric(nrow(deviations)) + ) + if (nlevels(groups) == 2) { + return(abs(group_means[, 1] - group_means[, 2])) + } + return(apply(group_means, 1, function(x) max(x) - min(x))) +} + #' @title differential_deviation_test #' @description Differential analysis is to test which #' motifs are having significant @@ -53,12 +126,12 @@ #' ) #' @export differential_deviation_test <- function( - deviations, - groups = NULL, - motifs = rownames(deviations), - alternative = c("two.sided", "less", "greater"), - parametric = TRUE, - padjMethod = "BH" + deviations, + groups = NULL, + motifs = rownames(deviations), + alternative = c("two.sided", "less", "greater"), + parametric = TRUE, + padjMethod = "BH" ) { if (!any(class(deviations) %in% c("data.frame", "matrix", "methylTFRdeviations"))) { @@ -70,70 +143,15 @@ differential_deviation_test <- function( if (is(deviations, "methylTFRdeviations")) { deviations <- deviations(deviations) } - if (is.null(groups)) { - groups <- colnames(deviations) - } - if (is.null(groups)) { - stop( - "No group labels found. Provide 'groups', or supply deviations ", - "with column names identifying the groups." - ) - } - groups <- as.factor(groups) - if (length(groups) != ncol(deviations)) { - stop("'groups' must have one entry per column of 'deviations'") - } - if (nlevels(groups) < 2) { - stop("'groups' must contain at least two distinct groups") - } + groups <- resolve_diff_groups(deviations, groups) if (length(alternative) > 1) { stop( "Please indicate one of the alternatives only." ) } - if (parametric) { - if (nlevels(groups) == 2) { - # t-test - p_val <- apply( - deviations, 1, - t_helper, groups, alternative - ) - } else { - # anova - p_val <- apply( - deviations, 1, - anova_helper, groups - ) - } - } else { - if (nlevels(groups) == 2) { - # wilcoxon - p_val <- apply( - deviations, 1, - wilcoxon_helper, groups, alternative - ) - } else { - # kruskal-wallis - p_val <- apply( - deviations, 1, - kw_helper, groups - ) - } - } - p_adj <- p.adjust(p_val, - method = padjMethod - ) - # Compute group means - group_means <- vapply( - levels(groups), - function(g) rowMeans(deviations[, groups == g, drop = FALSE]), - numeric(nrow(deviations)) - ) - mean_diff <- if (nlevels(groups) == 2) { - abs(group_means[, 1] - group_means[, 2]) - } else { - apply(group_means, 1, function(x) max(x) - min(x)) - } + p_val <- diff_pvalues(deviations, groups, parametric, alternative) + p_adj <- p.adjust(p_val, method = padjMethod) + mean_diff <- group_mean_difference(deviations, groups) return(data.frame( motifs = motifs, diff --git a/R/expected_deviations.R b/R/expected_deviations.R index 48731b5..4a51462 100644 --- a/R/expected_deviations.R +++ b/R/expected_deviations.R @@ -31,18 +31,19 @@ #' bin_meth <- addGCBintoMethylome(msites, gcdist) #' @author Irem Gunduz addGCBintoMethylome <- function( - msites, - gcdist, - ignoreStrand = TRUE + msites, + gcdist, + ignoreStrand = TRUE ) { if (!is.logical(ignoreStrand)) { - warning("Found invalid strand option, - using the default") + warning("Found invalid strand option, using the default") ignoreStrand <- TRUE } if (is.null(msites) || !is(msites, "GRanges")) { - stop("Please provide a valid methylation - sites with read_methylome function") + stop( + "Please provide a valid methylation sites with ", + "read_methylome function" + ) } if (is.null(gcdist) || !is(gcdist, "GRanges")) { stop("Please provide a valid GC distribution") @@ -51,8 +52,7 @@ addGCBintoMethylome <- function( ignore.strand = ignoreStrand ) if (length(hits@from) == 0) { - stop("No methylation sites found - in the GC distribution") + stop("No methylation sites found in the GC distribution") } gcmap <- data.table( mscore = msites[hits@from]$score, @@ -81,12 +81,10 @@ addGCBintoMethylome <- function( #' @keywords internal computeExpectations <- function(binMsites, gcfreq) { if (!is.matrix(binMsites)) { - stop("Please provide a valid - GC bin frequency table as a matrix") + stop("Please provide a valid GC bin frequency table as a matrix") } if (!is.matrix(gcfreq)) { - stop("Please provide a valid - GC bin frequency table as a matrix") + stop("Please provide a valid GC bin frequency table as a matrix") } exp.data <- t(gcfreq) %*% binMsites[, 2] mpos <- round(seq(-floor(length(exp.data) / 2), diff --git a/R/memory_helpers.R b/R/memory_helpers.R index 007b43d..602f680 100644 --- a/R/memory_helpers.R +++ b/R/memory_helpers.R @@ -12,14 +12,18 @@ #' @importFrom logger log_info #' @keywords internal create_sink <- function( - files_list, motifs, temp_dir = "methylTFR_tmp", pattern = "methylTFR", - fileext = ".h5", verbose = TRUE + files_list, motifs, temp_dir = "methylTFR_tmp", pattern = "methylTFR", + fileext = ".h5", verbose = TRUE ) { # Create a temp sink if (!dir.exists(temp_dir)) { dir.create(temp_dir) } - tempfile <- tempfile(pattern = pattern, tmpdir = temp_dir, fileext = fileext) + tempfile <- tempfile( + pattern = pattern, + tmpdir = temp_dir, + fileext = fileext + ) # Create a sink for each region type sink <- HDF5Array::HDF5RealizationSink( diff --git a/R/methyltfr_core.R b/R/methyltfr_core.R index a551414..e7b989a 100644 --- a/R/methyltfr_core.R +++ b/R/methyltfr_core.R @@ -10,8 +10,8 @@ #' @importFrom methods is #' @keywords internal check_annotation_inputs <- function( - tf_bindsites, gcfreqs, gc_dist, - enhancer = NULL + tf_bindsites, gcfreqs, gc_dist, + enhancer = NULL ) { if (any(vapply( list(tf_bindsites, gcfreqs, gc_dist), is.null, logical(1) @@ -44,8 +44,8 @@ check_annotation_inputs <- function( #' @return A named list with the validated values. #' @keywords internal check_run_options <- function( - chunkSize = 20, threads = 1, - ignoreStrand = TRUE, cov_threshold = 1 + chunkSize = 20, threads = 1, + ignoreStrand = TRUE, cov_threshold = 1 ) { if (!is.logical(ignoreStrand)) { warning("Found invalid strand option, using the default") @@ -107,6 +107,178 @@ read_sample_annotation <- function(annfile, sampleColName) { } +#' @title check_core_inputs +#' @description Validate the arguments shared by both methylTFR entry points. +#' @param sample_ids A character vector of sample identifiers. +#' @param msites_fun A function of a single integer sample index. +#' @param samples A \code{data.frame} with one row per sample. +#' @return Invisible \code{NULL}. Called for the errors it raises. +#' @keywords internal +check_core_inputs <- function(sample_ids, msites_fun, samples) { + if (!is.character(sample_ids) || length(sample_ids) == 0) { + stop("No samples to process.") + } + if (!is.function(msites_fun)) { + stop("msites_fun must be a function of a single sample index.") + } + if (nrow(samples) != length(sample_ids)) { + stop("Sample annotation must have one row per sample.") + } + invisible(NULL) +} + +#' @title valid_core_motifs +#' @description Drop motifs whose binding sites are empty or whose GC bin +#' frequency matrix is missing. +#' @param tf_bindsites a \code{GRangesList} of TF binding site positions. +#' @param gcfreqs a \code{list} of GC bin frequency tables. +#' @return A character vector of the motif names that can be processed. +#' @importFrom logger log_info +#' @keywords internal +valid_core_motifs <- function(tf_bindsites, gcfreqs) { + motifs <- names(gcfreqs) + valid_motifs <- vapply(motifs, function(m) { + has_tfbs <- !is.null(tf_bindsites[[m]]) && + length(tf_bindsites[[m]]) > 0 + has_matrix <- !is.null(gcfreqs[[m]]) + return(has_tfbs && has_matrix) + }, logical(1)) + + if (any(!valid_motifs)) { + num_discarded <- sum(!valid_motifs) + log_info( + "Discarding ", num_discarded, + " motifs due to empty TFBS or missing matrix." + ) + motifs <- motifs[valid_motifs] + } + + if (length(motifs) == 0) { + stop("No valid motifs remaining after validation.") + } + return(motifs) +} + +#' @title bpparam_from_threads +#' @description Build the \pkg{BiocParallel} back-end used to spread the +#' motifs of one chunk over workers. +#' @details A forking back-end is used where the platform supports it and a +#' socket back-end on Windows, so that \code{threads} has the same meaning +#' on every platform. \code{threads = 1} runs serially in the current +#' process. +#' @param threads Thread count for parallel processing. +#' @return A \code{BiocParallelParam} object. +#' @importFrom BiocParallel SerialParam MulticoreParam SnowParam +#' @keywords internal +bpparam_from_threads <- function(threads) { + if (is.null(threads) || !is.numeric(threads) || threads <= 1) { + return(SerialParam()) + } + if (.Platform$OS.type == "windows") { + return(SnowParam(workers = threads)) + } + return(MulticoreParam(workers = threads)) +} + +#' @title process_core_sample +#' @description Compute and write the deviations of one sample, one motif +#' chunk at a time. +#' @param index Integer index of the sample within \code{sample_ids}. +#' @param sample_ids A character vector of sample identifiers. +#' @param msites_fun A function of a single integer sample index. +#' @param motif_chunks A \code{list} of character vectors of motif names. +#' @param tf_bindsites a \code{GRangesList} of TF binding site positions. +#' @param gcfreqs a \code{list} of GC bin frequency tables. +#' @param gc_dist a \code{GRanges} of the genome-wide GC distribution. +#' @param dev_grid,exp_grid The grids the blocks are written on. +#' @param dev_sink,exp_sink The sinks the blocks are written to. +#' @param BPPARAM A \code{BiocParallelParam} object. +#' @param enhancer a \code{GRanges} restricting the analysis (optional). +#' @param ignoreStrand if TRUE, strand information is ignored. +#' @return Invisible \code{NULL}. Called for its effect on the sinks. +#' @importFrom BiocParallel bplapply +#' @importFrom logger log_info +#' @importFrom methods is +#' @keywords internal +process_core_sample <- function( + index, sample_ids, msites_fun, motif_chunks, tf_bindsites, gcfreqs, + gc_dist, dev_grid, exp_grid, dev_sink, exp_sink, BPPARAM, enhancer, + ignoreStrand +) { + sample_name <- sample_ids[index] + msites <- msites_fun(index) + if (!is(msites, "GRanges")) { + stop( + "msites_fun did not return a GRanges object for sample ", + sample_name + ) + } + log_info("Processing ", sample_name) + bin_meth <- addGCBintoMethylome(msites, gc_dist, ignoreStrand) + + # Process motifs in chunks + for (j in seq_along(motif_chunks)) { + chunk_motifs <- motif_chunks[[j]] + + sample_deviations <- bplapply(chunk_motifs, + computeDeviation, + msites = msites, + tf_bindsites = tf_bindsites, + gcfreqs = gcfreqs, + binMsites = bin_meth, + enhancer = enhancer, + ignoreStrand = ignoreStrand, + BPPARAM = BPPARAM + ) + names(sample_deviations) <- chunk_motifs + + # Write the block to the sink + write_block_to_sink( + lapply(sample_deviations, function(x) x$dev), + dev_grid, index, j, dev_sink + ) + write_block_to_sink( + lapply(sample_deviations, function(x) x$exp_dev), + exp_grid, index, j, exp_sink + ) + rm(sample_deviations) + } + rm(msites) + cleanMem() + log_info("Finished processing ", sample_name) + invisible(NULL) +} + +#' @title assemble_core_result +#' @description Close the sinks and assemble the deviations, their row-wise +#' Z-scores and the expected deviations into a result object. +#' @param dev_sink The sink holding the bias-corrected deviations. +#' @param exp_sink The sink holding the expected deviations. +#' @param samples A \code{data.frame} with one row per sample. +#' @return a \code{methylTFRdeviations} object. +#' @importFrom SummarizedExperiment SummarizedExperiment +#' @importFrom S4Vectors DataFrame +#' @importFrom DelayedArray DelayedArray close +#' @importFrom methods as new +#' @keywords internal +assemble_core_result <- function(dev_sink, exp_sink, samples) { + DelayedArray::close(dev_sink) + DelayedArray::close(exp_sink) + deviation <- as.matrix(t(as(dev_sink, "DelayedArray"))) + exp_dev <- as.matrix(t(as(exp_sink, "DelayedArray"))) + + se <- SummarizedExperiment( + assays = list( + deviations = deviation, + z = computeRowZScore(deviation), + expected = exp_dev + ), + colData = samples, + rowData = DataFrame(motifs = row.names(deviation)) + ) + return(new("methylTFRdeviations", se)) +} + #' @title methyltfr_core #' @description Internal engine shared by \code{\link{run_methyltfr}} and #' \code{\link{run_methylTFR_RnBeads}}. It validates the motif set, allocates @@ -139,7 +311,6 @@ read_sample_annotation <- function(annfile, sampleColName) { #' row-wise Z-scores and expected deviations. #' @importFrom GenomicRanges GRanges findOverlaps width resize start end #' @importFrom IRanges subsetByOverlaps -#' @importFrom parallel mclapply #' @importFrom logger log_info #' @importFrom SummarizedExperiment SummarizedExperiment #' @importFrom S4Vectors DataFrame @@ -147,48 +318,12 @@ read_sample_annotation <- function(annfile, sampleColName) { #' @importFrom methods as new is #' @keywords internal methyltfr_core <- function( - sample_ids, - msites_fun, - samples, - tf_bindsites, - gcfreqs, - gc_dist, - chunkSize = 20, - threads = 1, - enhancer = NULL, - ignoreStrand = TRUE + sample_ids, msites_fun, samples, tf_bindsites, gcfreqs, gc_dist, + chunkSize = 20, threads = 1, enhancer = NULL, ignoreStrand = TRUE ) { - if (!is.character(sample_ids) || length(sample_ids) == 0) { - stop("No samples to process.") - } - if (!is.function(msites_fun)) { - stop("msites_fun must be a function of a single sample index.") - } - if (nrow(samples) != length(sample_ids)) { - stop("Sample annotation must have one row per sample.") - } - - motifs <- names(gcfreqs) - - # Validate motifs: discard if TFBS is empty or matrix is missing - valid_motifs <- vapply(motifs, function(m) { - has_tfbs <- !is.null(tf_bindsites[[m]]) && length(tf_bindsites[[m]]) > 0 - has_matrix <- !is.null(gcfreqs[[m]]) - return(has_tfbs && has_matrix) - }, logical(1)) - - if (any(!valid_motifs)) { - num_discarded <- sum(!valid_motifs) - log_info( - "Discarding ", num_discarded, - " motifs due to empty TFBS or missing matrix." - ) - motifs <- motifs[valid_motifs] - } - - if (length(motifs) == 0) { - stop("No valid motifs remaining after validation.") - } + check_core_inputs(sample_ids, msites_fun, samples) + motifs <- valid_core_motifs(tf_bindsites, gcfreqs) + BPPARAM <- bpparam_from_threads(threads) # Split the motifs into chunks numChunks <- ceiling(length(motifs) / chunkSize) @@ -213,67 +348,16 @@ methyltfr_core <- function( } for (i in seq_along(sample_ids)) { - sample_name <- sample_ids[i] - msites <- msites_fun(i) - if (!is(msites, "GRanges")) { - stop( - "msites_fun did not return a GRanges object for sample ", - sample_name - ) - } - log_info("Processing ", sample_name) - bin_meth <- addGCBintoMethylome( - msites, - gc_dist, ignoreStrand + process_core_sample( + index = i, sample_ids = sample_ids, msites_fun = msites_fun, + motif_chunks = motif_chunks, tf_bindsites = tf_bindsites, + gcfreqs = gcfreqs, gc_dist = gc_dist, dev_grid = dev_grid, + exp_grid = exp_grid, dev_sink = dev_sink, exp_sink = exp_sink, + BPPARAM = BPPARAM, enhancer = enhancer, + ignoreStrand = ignoreStrand ) - - # Process motifs in chunks - for (j in seq_along(motif_chunks)) { - chunk_motifs <- motif_chunks[[j]] - - sample_deviations <- mclapply(chunk_motifs, - computeDeviation, - msites = msites, - tf_bindsites = tf_bindsites, - gcfreqs = gcfreqs, - binMsites = bin_meth, - enhancer = enhancer, - mc.cores = threads, - ignoreStrand = ignoreStrand - ) - names(sample_deviations) <- chunk_motifs - - # Write the block to the sink - write_block_to_sink( - lapply(sample_deviations, function(x) x$dev), - dev_grid, i, j, dev_sink - ) - write_block_to_sink( - lapply(sample_deviations, function(x) x$exp_dev), - exp_grid, i, j, exp_sink - ) - rm(sample_deviations) - } - rm(msites) - cleanMem() - log_info("Finished processing ", sample_name) } log_success("Computed all deviations successfully") - # Close the sinks - DelayedArray::close(dev_sink) - DelayedArray::close(exp_sink) - deviation <- as.matrix(t(as(dev_sink, "DelayedArray"))) - exp_dev <- as.matrix(t(as(exp_sink, "DelayedArray"))) - - se <- SummarizedExperiment( - assays = list( - deviations = deviation, - z = computeRowZScore(deviation), - expected = exp_dev - ), - colData = samples, - rowData = DataFrame(motifs = row.names(deviation)) - ) - return(new("methylTFRdeviations", se)) + return(assemble_core_result(dev_sink, exp_sink, samples)) } diff --git a/R/plot_helpers.R b/R/plot_helpers.R index 7fdce33..ca571f1 100644 --- a/R/plot_helpers.R +++ b/R/plot_helpers.R @@ -14,9 +14,9 @@ #' @importFrom S4Vectors mcols #' @import data.table computeFootprint <- function( - motif_name, - tf_bindsites, - msites, enhancer = NULL + motif_name, + tf_bindsites, + msites, enhancer = NULL ) { tfbs <- tf_bindsites[[motif_name]] w <- width(tfbs)[1] @@ -76,8 +76,8 @@ computeFootprint <- function( #' @importFrom S4Vectors mcols #' @import data.table computeExpectedFootprint <- function( - motif, gcfreqs, gc_dist, - enhancer = NULL, msites + motif, gcfreqs, gc_dist, + enhancer = NULL, msites ) { gcfreq <- gcfreqs[[motif]] diff --git a/R/plots.R b/R/plots.R index f03d2a2..28da0bd 100644 --- a/R/plots.R +++ b/R/plots.R @@ -1,3 +1,82 @@ +#' @title check_footprint_inputs +#' @description Validate the inputs shared by the two footprint plotting +#' functions. +#' @param motif Motif name as a character string. +#' @param tf_bindsites a \code{GRangesList} of TF binding site positions. +#' @param msites Methylation sites as a \code{GRanges} object. +#' @param gc_dist a \code{GRanges} of the genome-wide GC distribution. +#' @param gcfreqs a \code{list} of GC bin frequency tables. +#' @param enhancer a \code{GRanges} restricting the analysis (optional). +#' @return Invisible \code{NULL}. Called for the errors it raises. +#' @importFrom methods is +#' @keywords internal +check_footprint_inputs <- function( + motif, tf_bindsites, msites, gc_dist, gcfreqs, enhancer = NULL +) { + if (is.null(msites)) { + stop( + "msites must be a data frame, ", + "please provide the methylation sites" + ) + } + if (is.null(motif) || !is.character(motif)) { + stop("Please provide a valid motif name") + } + if (is.null(tf_bindsites) || + !any(c(!is(tf_bindsites, "GRangesList") || + !is.list(tf_bindsites)))) { + stop("Please provide a valid tf binding sites as GRangesList") + } + if (!is.null(enhancer) && !is(enhancer, "GRanges")) { + stop("Please provide a valid enhancer regions") + } + if (is.null(gc_dist) || !is(gc_dist, "GRanges")) { + stop("Please provide a valid gc_dist as GRanges") + } + if (is.null(gcfreqs) || !is(gcfreqs, "list")) { + stop("Please provide a valid gcfreqs as list") + } + if (!is(msites, "GRanges")) { + stop( + "Please provide a valid methylation ", + "sites with read_methylome function" + ) + } + invisible(NULL) +} + +#' @title expected_footprint_plot +#' @description Draw the observed and expected methylation profiles. +#' @param combined_data A \code{data.table} with the columns \code{x}, +#' \code{avg_methyl} and \code{type}. +#' @param motif Motif name as a character string. +#' @param sample_name Sample label used in the title. +#' @return A \code{ggplot} object. +#' @importFrom ggplot2 ggplot aes geom_line geom_point xlab ylab ggtitle +#' @importFrom ggplot2 theme_classic scale_color_manual theme +#' @keywords internal +expected_footprint_plot <- function(combined_data, motif, sample_name) { + ggplot(combined_data, aes( + x = x, + y = avg_methyl, + color = type + )) + + geom_line() + + geom_point() + + xlab("Distance from motif center") + + ylab("Methylation level") + + theme_classic() + + ggtitle(paste( + "TF footprint for", motif, "in", + sample_name + )) + + scale_color_manual(values = c( + "Expected" = "blue", + "Observed" = "red" + )) + + theme(legend.position = "bottom") +} + #' @title plotExpectedFootprint #' @description Creates a footprint plot of expected #' vs observed methylation @@ -58,52 +137,18 @@ #' #' @importFrom ggplot2 ggplot geom_point geom_line ggtitle theme_classic plotExpectedFootprint <- function( - motif, tf_bindsites, msites, - sample_name = NULL, gc_dist, gcfreqs, - enhancer = NULL, returnPlotData = FALSE + motif, tf_bindsites, msites, + sample_name = NULL, gc_dist, gcfreqs, + enhancer = NULL, returnPlotData = FALSE ) { - if (is.null(msites)) { - stop( - "msites must be a data frame, ", - "please provide the methylation sites" - ) - } # If sample label is not provided, use file name as label if (is.null(sample_name)) { sample_name <- "sample" } - # Check if motif is valid - if (is.null(motif) || !is.character(motif)) { - stop("Please provide a valid motif name") - } - - # Check if tf_bindsites is a GRangesList - if (is.null(tf_bindsites) || - !any(c(!is(tf_bindsites, "GRangesList") || - !is.list(tf_bindsites)))) { - stop("Please provide a valid tf binding sites as GRangesList") - } - # Check if enhancer is a GRanges - if (!is.null(enhancer) && !is(enhancer, "GRanges")) { - stop("Please provide a valid enhancer regions") - } - - # Check if gc_dist is a GRanges - if (is.null(gc_dist) || !is(gc_dist, "GRanges")) { - stop("Please provide a valid gc_dist as GRanges") - } + check_footprint_inputs( + motif, tf_bindsites, msites, gc_dist, gcfreqs, enhancer + ) - # Check if gcfreqs is a list - if (is.null(gcfreqs) || !is(gcfreqs, "list")) { - stop("Please provide a valid gcfreqs as list") - } - # Check if msites is a GRanges - if (is.null(msites) || !is(msites, "GRanges")) { - stop( - "Please provide a valid methylation ", - "sites with read_methylome function" - ) - } # Compute footprint for the motif plot_data <- computeFootprint( motif, @@ -124,26 +169,7 @@ plotExpectedFootprint <- function( plot_data[, .(x, avg_methyl, type)] )) - # Generate the footprint plot - p1 <- ggplot(combined_data, aes( - x = x, - y = avg_methyl, - color = type - )) + - geom_line() + - geom_point() + - xlab("Distance from motif center") + - ylab("Methylation level") + - theme_classic() + - ggtitle(paste( - "TF footprint for", motif, "in", - sample_name - )) + - scale_color_manual(values = c( - "Expected" = "blue", - "Observed" = "red" - )) + - theme(legend.position = "bottom") + p1 <- expected_footprint_plot(combined_data, motif, sample_name) if (returnPlotData) { return(list(plot = p1, plotDF = combined_data)) } else { @@ -151,6 +177,73 @@ plotExpectedFootprint <- function( } } +#' @title footprint_difference +#' @description Combine the observed and expected profiles into a single +#' corrected curve. +#' @param plot_data A \code{data.table} with the columns \code{x}, +#' \code{avg_methyl} and \code{type}. +#' @param method Either \code{"substraction"} or \code{"division"}. +#' @param sample_name Sample label used in the curve label. +#' @return A \code{list} with the corrected \code{data} and the axis +#' \code{lab}. +#' @keywords internal +footprint_difference <- function(plot_data, method, sample_name) { + if (method == "substraction") { + # Calculate observed - expected methylation + difference_data <- plot_data[, .( + avg_methyl = + avg_methyl[type == "Observed"] - avg_methyl[type == "Expected"] + ), + by = x + ] + difference_data[, type := paste( + "Obs. sub. Exp.", + sample_name + )] + lab <- "(Observed - Expected)" + } else { + # Calculate observed / expected methylation + difference_data <- plot_data[, .( + avg_methyl = + avg_methyl[type == "Observed"] / avg_methyl[type == "Expected"] + ), + by = x + ] + difference_data[, type := paste( + "Obs. div. Exp.", + sample_name + )] + lab <- "(Observed / Expected)" + } + return(list(data = difference_data, lab = lab)) +} + +#' @title normalise_footprint_flank +#' @description Normalise a corrected footprint against its outer flanking +#' windows, on the same scale as the statistic itself. +#' @details Dividing a difference by its flank mean would rescale it by an +#' arbitrary factor, because that mean is near zero: for a typical footprint +#' it inflates the curve by an order of magnitude and can invert its sign. +#' @param difference_data The corrected footprint as a \code{data.table}. +#' @param method Either \code{"substraction"} or \code{"division"}. +#' @param flankNorm Width of the flanking window used for normalisation. +#' @return The normalised \code{data.table}. +#' @keywords internal +normalise_footprint_flank <- function(difference_data, method, flankNorm) { + if (is.null(flankNorm) || flankNorm <= 0) { + return(difference_data) + } + flank <- max(abs(difference_data$x), na.rm = TRUE) + idx <- abs(difference_data$x) >= flank - flankNorm + norm_factor <- mean(difference_data$avg_methyl[idx], na.rm = TRUE) + if (method == "substraction") { + difference_data[, avg_methyl := avg_methyl - norm_factor] + } else { + difference_data[, avg_methyl := avg_methyl / norm_factor] + } + return(difference_data) +} + #' @title plotMotifFootprint #' @description Creates a footprint plot of bias #' corrected methylation for a given motif and sample @@ -211,11 +304,11 @@ plotExpectedFootprint <- function( #' #' @importFrom ggplot2 ggplot geom_point geom_line ggtitle theme_classic plotMotifFootprint <- function( - motif, - tf_bindsites, msites, - sample_name = NULL, gc_dist, gcfreqs, - enhancer = NULL, method = "division", - flankNorm = 50 + motif, + tf_bindsites, msites, + sample_name = NULL, gc_dist, gcfreqs, + enhancer = NULL, method = "division", + flankNorm = 50 ) { if (is.null(method) || !method %in% c("substraction", "division")) { @@ -230,48 +323,12 @@ plotMotifFootprint <- function( enhancer = enhancer, returnPlotData = TRUE )$plotDF - if (method == "substraction") { - # Calculate observed - expected methylation - difference_data <- plot_data[, .( - avg_methyl = - avg_methyl[type == "Observed"] - avg_methyl[type == "Expected"] - ), - by = x - ] - difference_data[, type := paste( - "Obs. sub. Exp.", - sample_name - )] - lab <- "(Observed - Expected)" - } else if (method == "division") { - # Calculate observed / expected methylation - difference_data <- plot_data[, .( - avg_methyl = - avg_methyl[type == "Observed"] / avg_methyl[type == "Expected"] - ), - by = x - ] - difference_data[, type := paste( - "Obs. div. Exp.", - sample_name - )] - lab <- "(Observed / Expected)" - } - if (!is.null(flankNorm) && flankNorm > 0) { - # Normalise against the outer flanking windows, on the same scale - # as the statistic itself. Dividing a difference by its flank mean - # would rescale it by an arbitrary factor, because that mean is - # near zero: for a typical footprint it inflates the curve by an - # order of magnitude and can invert its sign. - flank <- max(abs(difference_data$x), na.rm = TRUE) - idx <- abs(difference_data$x) >= flank - flankNorm - norm_factor <- mean(difference_data$avg_methyl[idx], na.rm = TRUE) - if (method == "substraction") { - difference_data[, avg_methyl := avg_methyl - norm_factor] - } else { - difference_data[, avg_methyl := avg_methyl / norm_factor] - } - } + diff <- footprint_difference(plot_data, method, sample_name) + difference_data <- diff$data + lab <- diff$lab + difference_data <- normalise_footprint_flank( + difference_data, method, flankNorm + ) p_combined <- ggplot(difference_data, aes( x = x, diff --git a/R/rnbeads_interface.R b/R/rnbeads_interface.R index ee8df49..e453608 100644 --- a/R/rnbeads_interface.R +++ b/R/rnbeads_interface.R @@ -1,3 +1,47 @@ +#' @title check_rnb_inputs +#' @description Check that \pkg{RnBeads} is available and that the object +#' handed in is an \code{RnBSet}. +#' @param rnb_set The object to validate. +#' @return Invisible \code{NULL}. Called for the errors it raises. +#' @importFrom methods is +#' @keywords internal +check_rnb_inputs <- function(rnb_set) { + if (!requireNamespace("RnBeads", quietly = TRUE)) { + stop( + "The RnBeads package is required for run_methylTFR_RnBeads(). ", + "Install it with BiocManager::install('RnBeads'), or export ", + "your samples and use run_methyltfr() instead." + ) + } + if (is.null(rnb_set) || !is(rnb_set, "RnBSet")) { + stop("Please provide a valid RnBSet object") + } + invisible(NULL) +} + +#' @title resolve_rnb_sample_ann +#' @description Default the sample annotation to the phenotype table of the +#' RnBeads object and check its shape. +#' @param rnb_set A preprocessed \code{RnBSet} object. +#' @param sample_ann Sample annotation supplied by the caller, or NULL. +#' @param sample_ids Character vector of sample identifiers. +#' @return The sample annotation as a \code{data.frame}. +#' @keywords internal +resolve_rnb_sample_ann <- function(rnb_set, sample_ann, sample_ids) { + if (is.null(sample_ann)) { + sample_ann <- as.data.frame(RnBeads::pheno(rnb_set), + stringsAsFactors = FALSE + ) + } + if (!is.data.frame(sample_ann)) { + stop("sample_ann must be a data.frame") + } + if (nrow(sample_ann) != length(sample_ids)) { + stop("sample_ann must have one row per sample in the RnBSet object") + } + return(sample_ann) +} + #' @title run_methylTFR_RnBeads #' @description Run the methylTFR workflow directly on a preprocessed #' \pkg{RnBeads} object, without exporting per-sample BED files first. @@ -60,59 +104,59 @@ #' @seealso \code{\link{run_methyltfr}} for the file-based entry point. #' @author Irem Gunduz #' @examples -#' # Not run: requires the RnBeads package, an hg38 annotation package and a -#' # preprocessed RnBeads set. -#' \donttest{ -#' if (requireNamespace("RnBeads", quietly = TRUE)) { -#' # rnb_set <- RnBeads::load.rnb.set("reports/rnb.set_preprocessed") -#' # gcfreqs <- getGCfreq(motifSet = "jaspar2020") -#' # gc_dist <- getGenomeGC("hg38") -#' # tf_bindsites <- getTFbindsites(motifSet = "jaspar2020") -#' # -#' # deviations <- run_methylTFR_RnBeads( -#' # rnb_set = rnb_set, -#' # tf_bindsites = tf_bindsites, -#' # gcfreqs = gcfreqs, -#' # gc_dist = gc_dist, -#' # threads = 8, -#' # chunkSize = 15 -#' # ) -#' } +#' # A minimal end-to-end run on the BATF example data bundled with the +#' # package. RnBeads and its hg38 annotation build the input object; both +#' # are optional dependencies. +#' if (requireNamespace("RnBeads", quietly = TRUE) && +#' requireNamespace("RnBeads.hg38", quietly = TRUE)) { +#' load(system.file("extdata", "example_data.rda", package = "methylTFR")) +#' load(system.file( +#' "extdata", "BATF_tf_bindsites.rda", +#' package = "methylTFR" +#' )) +#' load(system.file("extdata", "BATF_gcfreqs.rda", package = "methylTFR")) +#' load(system.file("extdata", "gcdist_subset.rda", package = "methylTFR")) +#' +#' # RnBiseqSet() takes methylation as a fraction and coverage as counts, +#' # with one column per sample. +#' sites <- data.frame( +#' chromosome = as.character(GenomicRanges::seqnames(msites)), +#' position = GenomicRanges::start(msites), +#' strand = "*", +#' stringsAsFactors = FALSE +#' ) +#' rnb_set <- RnBeads::RnBiseqSet( +#' pheno = data.frame( +#' sampleName = "sample_1", stringsAsFactors = FALSE +#' ), +#' sites = sites, +#' meth = matrix(msites$score, ncol = 1), +#' covg = matrix(msites$coverage, ncol = 1), +#' assembly = "hg38", +#' summarize.regions = FALSE +#' ) +#' +#' devs <- run_methylTFR_RnBeads( +#' rnb_set = rnb_set, +#' tf_bindsites = tf_bindsites, +#' gcfreqs = gcfreqs, +#' gc_dist = gcdist +#' ) +#' deviations(devs) #' } #' @export run_methylTFR_RnBeads <- function( - rnb_set, tf_bindsites = NULL, - gcfreqs = NULL, gc_dist = NULL, - chunkSize = 20, threads = 1, - enhancer = NULL, ignoreStrand = TRUE, - cov_threshold = 1, sample_ann = NULL + rnb_set, tf_bindsites = NULL, gcfreqs = NULL, gc_dist = NULL, + chunkSize = 20, threads = 1, enhancer = NULL, ignoreStrand = TRUE, + cov_threshold = 1, sample_ann = NULL ) { - if (!requireNamespace("RnBeads", quietly = TRUE)) { - stop( - "The RnBeads package is required for run_methylTFR_RnBeads(). ", - "Install it with BiocManager::install('RnBeads'), or export ", - "your samples and use run_methyltfr() instead." - ) - } - if (is.null(rnb_set) || !is(rnb_set, "RnBSet")) { - stop("Please provide a valid RnBSet object") - } + check_rnb_inputs(rnb_set) check_annotation_inputs(tf_bindsites, gcfreqs, gc_dist, enhancer) opts <- check_run_options(chunkSize, threads, ignoreStrand, cov_threshold) sample_ids <- rnb_sample_ids(rnb_set) - if (is.null(sample_ann)) { - sample_ann <- as.data.frame(RnBeads::pheno(rnb_set), - stringsAsFactors = FALSE - ) - } - if (!is.data.frame(sample_ann)) { - stop("sample_ann must be a data.frame") - } - if (nrow(sample_ann) != length(sample_ids)) { - stop("sample_ann must have one row per sample in the RnBSet object") - } + sample_ann <- resolve_rnb_sample_ann(rnb_set, sample_ann, sample_ids) sites_gr <- rnb_sites_to_granges(rnb_set, opts$ignoreStrand) has_covg <- rnb_has_coverage(rnb_set) @@ -268,8 +312,8 @@ rnb_has_coverage <- function(rnb_set) { #' @importFrom logger log_warn #' @keywords internal rnb_sample_msites <- function( - rnb_set, sites_gr, index, - cov_threshold = 1, has_covg = TRUE + rnb_set, sites_gr, index, + cov_threshold = 1, has_covg = TRUE ) { index <- as.integer(index) mvals <- rnb_column(RnBeads::meth, rnb_set, index) diff --git a/R/run_methyltfr.R b/R/run_methyltfr.R index 54a0177..9a957dc 100644 --- a/R/run_methyltfr.R +++ b/R/run_methyltfr.R @@ -1,3 +1,83 @@ +#' @title check_file_run_inputs +#' @description Validate the file-specific arguments of +#' \code{run_methyltfr}. +#' @param filetype File type of the methylation call files. +#' @param sampleColName Column name holding the file names. +#' @param full_path if TRUE, the annotation file holds full paths. +#' @return Invisible \code{NULL}. Called for the errors it raises. +#' @keywords internal +check_file_run_inputs <- function(filetype, sampleColName, full_path) { + if (!tolower(filetype) %in% c( + "bissnp", "epp", "allc", "bismarkcytosine", + "bismarkcov", "encode" + )) { + stop("Please provide a valid file type") + } + if (is.null(sampleColName) || !is.character(sampleColName)) { + stop("Please provide a valid sample column name") + } + if (!is.logical(full_path)) { + stop( + "Invalid full path flag detected, ", + "please provide a valid logical value" + ) + } + invisible(NULL) +} + +#' @title resolve_annotation_file +#' @description Work out the path of the sample annotation file. +#' @param annfile Explicit path to the annotation file, or NULL. +#' @param sample_ann Name of the annotation file inside \code{sample_dir}. +#' @param sample_dir Directory holding the methylation call files. +#' @return The path of the annotation file as a character string. +#' @keywords internal +resolve_annotation_file <- function(annfile, sample_ann, sample_dir) { + if (!is.null(annfile) && is.character(annfile)) { + return(annfile) + } + if (is.null(sample_ann) || !is.character(sample_ann)) { + stop("Please provide a valid sample annotation file") + } + if (is.null(sample_dir) || !is.character(sample_dir)) { + stop("Please provide a valid sample directory") + } + if (!dir.exists(sample_dir)) { + stop( + "Sample directory does not exist, ", + "please check the directory path" + ) + } + return(file.path(sample_dir, sample_ann)) +} + +#' @title locate_sample_files +#' @description Build and check the list of per-sample methylation files. +#' @param samples The sample annotation as a \code{data.frame}. +#' @param sample_dir Directory holding the methylation call files. +#' @param sampleColName Column name holding the file names. +#' @param full_path if TRUE, the annotation file holds full paths. +#' @return A character vector of existing file paths. +#' @importFrom logger log_success +#' @keywords internal +locate_sample_files <- function( + samples, sample_dir, sampleColName, full_path +) { + if (full_path) { + files_list <- samples[, sampleColName] + } else { + files_list <- file.path(sample_dir, samples[, sampleColName]) + } + if (!all(file.exists(files_list))) { + stop( + "Some of the files does not exist, ", + "please check the file path!" + ) + } + log_success("The samples are successfully located") + return(files_list) +} + #' @title run_methyltfr #' @description This function is a wrapper function to #' calculate the deviation @@ -34,7 +114,6 @@ #' @importFrom GenomicRanges GRanges findOverlaps width resize start end #' @importFrom IRanges subsetByOverlaps #' @importFrom data.table data.table -#' @importFrom parallel mclapply #' @importFrom logger log_info log_error #' @importFrom SummarizedExperiment SummarizedExperiment #' @importFrom S4Vectors DataFrame @@ -47,65 +126,69 @@ #' bias-corrected deviation and Z-scores #' @seealso \code{\link{run_methylTFR_RnBeads}} for running methylTFR #' directly on a preprocessed RnBeads object. +#' @examples +#' # A minimal end-to-end run on the BATF example data bundled with the +#' # package. The annotation objects cover a single motif, so the result +#' # has one row. +#' load(system.file("extdata", "example_data.rda", package = "methylTFR")) +#' load(system.file("extdata", "BATF_tf_bindsites.rda", package = "methylTFR")) +#' load(system.file("extdata", "BATF_gcfreqs.rda", package = "methylTFR")) +#' load(system.file("extdata", "gcdist_subset.rda", package = "methylTFR")) +#' +#' # run_methyltfr() reads per-sample calls from disk, so the bundled sites +#' # are written out as a bismarkCov file first. +#' sample_dir <- tempfile("methylTFR_example") +#' dir.create(sample_dir) +#' n_meth <- round(msites$score * msites$coverage) +#' write.table( +#' data.frame( +#' chr = as.character(GenomicRanges::seqnames(msites)), +#' start = GenomicRanges::start(msites), +#' end = GenomicRanges::end(msites), +#' percent = msites$score * 100, +#' meth = n_meth, +#' unmeth = msites$coverage - n_meth +#' ), +#' file.path(sample_dir, "sample_1.cov"), +#' sep = "\t", row.names = FALSE, col.names = FALSE, quote = FALSE +#' ) +#' write.table( +#' data.frame(sampleName = "sample_1", bedFile = "sample_1.cov"), +#' file.path(sample_dir, "samples.tsv"), +#' sep = "\t", row.names = FALSE, quote = FALSE +#' ) +#' +#' devs <- run_methyltfr( +#' sample_ann = "samples.tsv", +#' sample_dir = sample_dir, +#' tf_bindsites = tf_bindsites, +#' gcfreqs = gcfreqs, +#' gc_dist = gcdist, +#' filetype = "bismarkcov" +#' ) +#' deviations(devs) +#' +#' unlink(sample_dir, recursive = TRUE) #' @export run_methyltfr <- function( - sample_ann, sample_dir, tf_bindsites = NULL, - gcfreqs = NULL, gc_dist = NULL, - sampleColName = "bedFile", chunkSize = 20, - full_path = FALSE, annfile = NULL, threads = 1, - enhancer = NULL, filetype = NULL, - ignoreStrand = TRUE, cov_threshold = 1 + sample_ann, sample_dir, tf_bindsites = NULL, + gcfreqs = NULL, gc_dist = NULL, + sampleColName = "bedFile", chunkSize = 20, + full_path = FALSE, annfile = NULL, threads = 1, + enhancer = NULL, filetype = NULL, + ignoreStrand = TRUE, cov_threshold = 1 ) { - if (!tolower(filetype) %in% c( - "bissnp", "epp", "allc", "bismarkcytosine", - "bismarkcov", "encode" - )) { - stop("Please provide a valid file type") - } - if (is.null(sampleColName) || !is.character(sampleColName)) { - stop("Please provide a valid sample column name") - } - if (!is.logical(full_path)) { - stop( - "Invalid full path flag detected, ", - "please provide a valid logical value" - ) - } + check_file_run_inputs(filetype, sampleColName, full_path) check_annotation_inputs(tf_bindsites, gcfreqs, gc_dist, enhancer) opts <- check_run_options( chunkSize, threads, ignoreStrand, cov_threshold ) - if (is.null(annfile) || !is.character(annfile)) { - if (is.null(sample_ann) || !is.character(sample_ann)) { - stop("Please provide a valid sample annotation file") - } - if (is.null(sample_dir) || !is.character(sample_dir)) { - stop("Please provide a valid sample directory") - } - if (!dir.exists(sample_dir)) { - stop( - "Sample directory does not exist, ", - "please check the directory path" - ) - } - annfile <- file.path(sample_dir, sample_ann) - } - + annfile <- resolve_annotation_file(annfile, sample_ann, sample_dir) samples <- read_sample_annotation(annfile, sampleColName) - - if (full_path) { - files_list <- samples[, sampleColName] - } else { - files_list <- file.path(sample_dir, samples[, sampleColName]) - } - if (!all(file.exists(files_list))) { - stop( - "Some of the files does not exist, ", - "please check the file path!" - ) - } - log_success("The samples are successfully located") + files_list <- locate_sample_files( + samples, sample_dir, sampleColName, full_path + ) # Per-sample reader handed to the shared engine msites_fun <- function(i) { diff --git a/R/variability.R b/R/variability.R index dc8a531..d524ac8 100644 --- a/R/variability.R +++ b/R/variability.R @@ -52,6 +52,111 @@ calibrateDeviations <- function(devs, method = c("robust", "gaussian")) { } +#' @title check_variability_inputs +#' @description Validate the inputs of \code{computeZScoreVariability} and +#' return the deviation scores as a matrix. +#' @param object A \code{methylTFRdeviations} object, matrix or data.frame. +#' @param bootstrap Logical, whether bootstrap bounds were requested. +#' @param conf_level Confidence level of the bootstrap bounds. +#' @return The deviation scores as a numeric matrix. +#' @importFrom methods is +#' @keywords internal +check_variability_inputs <- function(object, bootstrap, conf_level) { + if (is(object, "methylTFRdeviations")) { + devs <- deviations(object) + } else if (is.matrix(object) || is.data.frame(object)) { + devs <- as.matrix(object) + } else { + stop( + "object must be a methylTFRdeviations object, ", + "a matrix or a data.frame" + ) + } + if (!is.numeric(devs)) { + stop("The deviation scores must be numeric") + } + if (ncol(devs) < 3) { + stop( + "At least three samples are required to estimate variability; ", + "found ", ncol(devs), "." + ) + } + if (!is.logical(bootstrap) || length(bootstrap) != 1) { + stop("bootstrap must be a single logical value") + } + if (!is.numeric(conf_level) || conf_level <= 0 || conf_level >= 1) { + stop("conf_level must be a number between 0 and 1") + } + if (nrow(devs) < 50) { + message( + "Only ", nrow(devs), " motifs supplied. The within-sample null is ", + "estimated across motifs, so variability estimates from small ", + "motif sets should be treated as indicative only." + ) + } + return(devs) +} + +#' @title variability_test +#' @description Summarise the calibrated Z-scores into a variability per +#' motif and test each against the chi-squared null. +#' @param z The calibrated Z-score matrix. +#' @param motif_names Character vector of motif names. +#' @param padjMethod Multiple testing correction passed to +#' \code{stats::p.adjust}. +#' @return A \code{data.frame} with one row per motif. +#' @keywords internal +variability_test <- function(z, motif_names, padjMethod) { + variability <- matrixStats::rowSds(z, na.rm = TRUE) + n_obs <- rowSums(!is.na(z)) + + p_value <- rep(NA_real_, length(variability)) + testable <- is.finite(variability) & n_obs > 1 + p_value[testable] <- stats::pchisq( + variability[testable]^2 * (n_obs[testable] - 1), + df = n_obs[testable] - 1, + lower.tail = FALSE + ) + + return(data.frame( + motifs = motif_names, + variability = variability, + p_value = p_value, + p_value_adjusted = stats::p.adjust(p_value, method = padjMethod), + stringsAsFactors = FALSE, + row.names = NULL + )) +} + +#' @title bootstrap_variability +#' @description Add bootstrap confidence bounds to a variability table. +#' @param z The calibrated Z-score matrix. +#' @param res The variability table to add the bounds to. +#' @param niterations Number of bootstrap iterations. +#' @param conf_level Confidence level of the bounds. +#' @return \code{res} with the two bound columns added. +#' @keywords internal +bootstrap_variability <- function(z, res, niterations, conf_level) { + if (!is.numeric(niterations) || niterations < 2) { + stop("niterations must be a number greater than 1") + } + niterations <- as.integer(niterations) + boot <- vapply(seq_len(niterations), function(k) { + idx <- sample.int(ncol(z), ncol(z), replace = TRUE) + matrixStats::rowSds(z[, idx, drop = FALSE], na.rm = TRUE) + }, numeric(nrow(z))) + alpha <- (1 - conf_level) / 2 + res$bootstrap_lower_bound <- apply( + boot, 1, stats::quantile, + probs = alpha, na.rm = TRUE + ) + res$bootstrap_upper_bound <- apply( + boot, 1, stats::quantile, + probs = 1 - alpha, na.rm = TRUE + ) + return(res) +} + #' @title computeZScoreVariability #' @description Identify transcription factor motifs whose methylTFR activity #' varies across samples more than expected by chance. @@ -119,47 +224,11 @@ calibrateDeviations <- function(devs, method = c("robust", "gaussian")) { #' @author Irem Gunduz #' @export computeZScoreVariability <- function( - object, - method = c("robust", "gaussian"), - bootstrap = FALSE, - niterations = 1000L, - conf_level = 0.95, - padjMethod = "BH" + object, method = c("robust", "gaussian"), bootstrap = FALSE, + niterations = 1000L, conf_level = 0.95, padjMethod = "BH" ) { method <- match.arg(method) - if (is(object, "methylTFRdeviations")) { - devs <- deviations(object) - } else if (is.matrix(object) || is.data.frame(object)) { - devs <- as.matrix(object) - } else { - stop( - "object must be a methylTFRdeviations object, ", - "a matrix or a data.frame" - ) - } - if (!is.numeric(devs)) { - stop("The deviation scores must be numeric") - } - if (ncol(devs) < 3) { - stop( - "At least three samples are required to estimate variability; ", - "found ", ncol(devs), "." - ) - } - if (!is.logical(bootstrap) || length(bootstrap) != 1) { - stop("bootstrap must be a single logical value") - } - if (!is.numeric(conf_level) || conf_level <= 0 || conf_level >= 1) { - stop("conf_level must be a number between 0 and 1") - } - - if (nrow(devs) < 50) { - message( - "Only ", nrow(devs), " motifs supplied. The within-sample null is ", - "estimated across motifs, so variability estimates from small ", - "motif sets should be treated as indicative only." - ) - } + devs <- check_variability_inputs(object, bootstrap, conf_level) motif_names <- rownames(devs) if (is.null(motif_names)) { @@ -167,46 +236,10 @@ computeZScoreVariability <- function( } z <- calibrateDeviations(devs, method = method) - variability <- matrixStats::rowSds(z, na.rm = TRUE) - n_obs <- rowSums(!is.na(z)) - - p_value <- rep(NA_real_, length(variability)) - testable <- is.finite(variability) & n_obs > 1 - p_value[testable] <- stats::pchisq( - variability[testable]^2 * (n_obs[testable] - 1), - df = n_obs[testable] - 1, - lower.tail = FALSE - ) - p_adj <- stats::p.adjust(p_value, method = padjMethod) - - res <- data.frame( - motifs = motif_names, - variability = variability, - p_value = p_value, - p_value_adjusted = p_adj, - stringsAsFactors = FALSE, - row.names = NULL - ) + res <- variability_test(z, motif_names, padjMethod) if (bootstrap) { - if (!is.numeric(niterations) || niterations < 2) { - stop("niterations must be a number greater than 1") - } - niterations <- as.integer(niterations) - boot <- vapply(seq_len(niterations), function(k) { - idx <- sample.int(ncol(z), ncol(z), replace = TRUE) - matrixStats::rowSds(z[, idx, drop = FALSE], na.rm = TRUE) - }, numeric(nrow(z))) - alpha <- (1 - conf_level) / 2 - res$bootstrap_lower_bound <- apply( - boot, 1, stats::quantile, - probs = alpha, na.rm = TRUE - ) - res$bootstrap_upper_bound <- apply( - boot, 1, stats::quantile, - probs = 1 - alpha, na.rm = TRUE - ) + res <- bootstrap_variability(z, res, niterations, conf_level) } - return(res) } diff --git a/docs/404.html b/docs/404.html index 8eef24d..939fca4 100644 --- a/docs/404.html +++ b/docs/404.html @@ -20,7 +20,7 @@ methylTFR - 0.99.4 + 0.99.5 + + + + + +
+
+
+ +

assemble_core_result

+ Source: R/methyltfr_core.R +
assemble_core_result.Rd
+
+ +
+

Close the sinks and assemble the deviations, their row-wise +Z-scores and the expected deviations into a result object.

+
+ +
+

Usage

+ 
assemble_core_result(dev_sink, exp_sink, samples)
+
+ +
+

Arguments

+ + +
dev_sink
+

The sink holding the bias-corrected deviations.

+ + +
exp_sink
+

The sink holding the expected deviations.

+ + +
samples
+

A data.frame with one row per sample.

+ +
+
+

Value

+

a methylTFRdeviations object.

+
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/bootstrap_variability.html b/docs/reference/bootstrap_variability.html new file mode 100644 index 0000000..8e9a31c --- /dev/null +++ b/docs/reference/bootstrap_variability.html @@ -0,0 +1,95 @@ + +bootstrap_variability — bootstrap_variability • methylTFR + Skip to contents + + +
+
+
+ +

bootstrap_variability

+ Source: R/variability.R +
bootstrap_variability.Rd
+
+ +
+

Add bootstrap confidence bounds to a variability table.

+
+ +
+

Usage

+ 
bootstrap_variability(z, res, niterations, conf_level)
+
+ +
+

Arguments

+ + +
z
+

The calibrated Z-score matrix.

+ + +
res
+

The variability table to add the bounds to.

+ + +
niterations
+

Number of bootstrap iterations.

+ + +
conf_level
+

Confidence level of the bounds.

+ +
+
+

Value

+

res with the two bound columns added.

+
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/bpparam_from_threads.html b/docs/reference/bpparam_from_threads.html new file mode 100644 index 0000000..5268870 --- /dev/null +++ b/docs/reference/bpparam_from_threads.html @@ -0,0 +1,93 @@ + +bpparam_from_threads — bpparam_from_threads • methylTFR + Skip to contents + + +
+
+
+ +

bpparam_from_threads

+ Source: R/methyltfr_core.R +
bpparam_from_threads.Rd
+
+ +
+

Build the BiocParallel back-end used to spread the +motifs of one chunk over workers.

+
+ +
+

Usage

+ 
bpparam_from_threads(threads)
+
+ +
+

Arguments

+ + +
threads
+

Thread count for parallel processing.

+ +
+
+

Value

+

A BiocParallelParam object.

+
+
+

Details

+

A forking back-end is used where the platform supports it and a +socket back-end on Windows, so that threads has the same meaning +on every platform. threads = 1 runs serially in the current +process.

+
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/calibrateDeviations.html b/docs/reference/calibrateDeviations.html index aa723a4..9c22631 100644 --- a/docs/reference/calibrateDeviations.html +++ b/docs/reference/calibrateDeviations.html @@ -13,7 +13,7 @@ methylTFR - 0.99.4 + 0.99.5 + + + + + +
+
+
+ +

check_core_inputs

+ Source: R/methyltfr_core.R +
check_core_inputs.Rd
+
+ +
+

Validate the arguments shared by both methylTFR entry points.

+
+ +
+

Usage

+ 
check_core_inputs(sample_ids, msites_fun, samples)
+
+ +
+

Arguments

+ + +
sample_ids
+

A character vector of sample identifiers.

+ + +
msites_fun
+

A function of a single integer sample index.

+ + +
samples
+

A data.frame with one row per sample.

+ +
+
+

Value

+

Invisible NULL. Called for the errors it raises.

+
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/check_deviation_inputs.html b/docs/reference/check_deviation_inputs.html new file mode 100644 index 0000000..afe30ba --- /dev/null +++ b/docs/reference/check_deviation_inputs.html @@ -0,0 +1,95 @@ + +check_deviation_inputs — check_deviation_inputs • methylTFR + Skip to contents + + +
+
+
+ +

check_deviation_inputs

+ Source: R/compute_deviations.R +
check_deviation_inputs.Rd
+
+ +
+

Validate the inputs of computeDeviation.

+
+ +
+

Usage

+ 
check_deviation_inputs(motif, msites, tf_bindsites, enhancer = NULL)
+
+ +
+

Arguments

+ + +
motif
+

Motif name as a character string.

+ + +
msites
+

Methylation sites as a GRanges object.

+ + +
tf_bindsites
+

a GRangesList of TF binding site positions.

+ + +
enhancer
+

a GRanges restricting the analysis (optional).

+ +
+
+

Value

+

Invisible NULL. Called for the errors it raises.

+
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/check_file_run_inputs.html b/docs/reference/check_file_run_inputs.html new file mode 100644 index 0000000..680a013 --- /dev/null +++ b/docs/reference/check_file_run_inputs.html @@ -0,0 +1,94 @@ + +check_file_run_inputs — check_file_run_inputs • methylTFR + Skip to contents + + +
+
+
+ +

check_file_run_inputs

+ Source: R/run_methyltfr.R +
check_file_run_inputs.Rd
+
+ +
+

Validate the file-specific arguments of +run_methyltfr.

+
+ +
+

Usage

+ 
check_file_run_inputs(filetype, sampleColName, full_path)
+
+ +
+

Arguments

+ + +
filetype
+

File type of the methylation call files.

+ + +
sampleColName
+

Column name holding the file names.

+ + +
full_path
+

if TRUE, the annotation file holds full paths.

+ +
+
+

Value

+

Invisible NULL. Called for the errors it raises.

+
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/check_footprint_inputs.html b/docs/reference/check_footprint_inputs.html new file mode 100644 index 0000000..ebd282b --- /dev/null +++ b/docs/reference/check_footprint_inputs.html @@ -0,0 +1,113 @@ + +check_footprint_inputs — check_footprint_inputs • methylTFR + Skip to contents + + +
+
+
+ +

check_footprint_inputs

+ Source: R/plots.R +
check_footprint_inputs.Rd
+
+ +
+

Validate the inputs shared by the two footprint plotting +functions.

+
+ +
+

Usage

+ 
check_footprint_inputs(
+  motif,
+  tf_bindsites,
+  msites,
+  gc_dist,
+  gcfreqs,
+  enhancer = NULL
+)
+
+ +
+

Arguments

+ + +
motif
+

Motif name as a character string.

+ + +
tf_bindsites
+

a GRangesList of TF binding site positions.

+ + +
msites
+

Methylation sites as a GRanges object.

+ + +
gc_dist
+

a GRanges of the genome-wide GC distribution.

+ + +
gcfreqs
+

a list of GC bin frequency tables.

+ + +
enhancer
+

a GRanges restricting the analysis (optional).

+ +
+
+

Value

+

Invisible NULL. Called for the errors it raises.

+
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/check_rnb_inputs.html b/docs/reference/check_rnb_inputs.html new file mode 100644 index 0000000..08527ed --- /dev/null +++ b/docs/reference/check_rnb_inputs.html @@ -0,0 +1,86 @@ + +check_rnb_inputs — check_rnb_inputs • methylTFR + Skip to contents + + +
+
+
+ +

check_rnb_inputs

+ Source: R/rnbeads_interface.R +
check_rnb_inputs.Rd
+
+ +
+

Check that RnBeads is available and that the object +handed in is an RnBSet.

+
+ +
+

Usage

+ 
check_rnb_inputs(rnb_set)
+
+ +
+

Arguments

+ + +
rnb_set
+

The object to validate.

+ +
+
+

Value

+

Invisible NULL. Called for the errors it raises.

+
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/check_run_options.html b/docs/reference/check_run_options.html index dfd641f..cf1cc90 100644 --- a/docs/reference/check_run_options.html +++ b/docs/reference/check_run_options.html @@ -9,7 +9,7 @@ methylTFR - 0.99.4 + 0.99.5 + + + + + +
+
+
+ +

check_variability_inputs

+ Source: R/variability.R +
check_variability_inputs.Rd
+
+ +
+

Validate the inputs of computeZScoreVariability and +return the deviation scores as a matrix.

+
+ +
+

Usage

+ 
check_variability_inputs(object, bootstrap, conf_level)
+
+ +
+

Arguments

+ + +
object
+

A methylTFRdeviations object, matrix or data.frame.

+ + +
bootstrap
+

Logical, whether bootstrap bounds were requested.

+ + +
conf_level
+

Confidence level of the bootstrap bounds.

+ +
+
+

Value

+

The deviation scores as a numeric matrix.

+
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/cleanMem.html b/docs/reference/cleanMem.html index 77fe125..3cf2045 100644 --- a/docs/reference/cleanMem.html +++ b/docs/reference/cleanMem.html @@ -7,7 +7,7 @@ methylTFR - 0.99.4 + 0.99.5 + + + + + +
+
+
+ +

diff_pvalues

+ Source: R/differential_analysis.R +
diff_pvalues.Rd
+
+ +
+

Test every motif for a difference between the groups, with +the test chosen from the number of groups and parametric.

+
+ +
+

Usage

+ 
diff_pvalues(deviations, groups, parametric, alternative)
+
+ +
+

Arguments

+ + +
deviations
+

A matrix of deviation scores, motifs in rows.

+ + +
groups
+

The group labels as a factor.

+ + +
parametric
+

if TRUE, use a t-test or ANOVA, otherwise a Wilcoxon +or Kruskal-Wallis test.

+ + +
alternative
+

The alternative hypothesis of the two-group tests.

+ +
+
+

Value

+

A numeric vector of p-values, one per motif.

+
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/differential_deviation_test.html b/docs/reference/differential_deviation_test.html index 7b7153e..ac4b22a 100644 --- a/docs/reference/differential_deviation_test.html +++ b/docs/reference/differential_deviation_test.html @@ -13,7 +13,7 @@ methylTFR - 0.99.4 + 0.99.5 + + + + + +
+
+
+ +

expected_footprint_plot

+ Source: R/plots.R +
expected_footprint_plot.Rd
+
+ +
+

Draw the observed and expected methylation profiles.

+
+ +
+

Usage

+ 
expected_footprint_plot(combined_data, motif, sample_name)
+
+ +
+

Arguments

+ + +
combined_data
+

A data.table with the columns x, +avg_methyl and type.

+ + +
motif
+

Motif name as a character string.

+ + +
sample_name
+

Sample label used in the title.

+ +
+
+

Value

+

A ggplot object.

+
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/footprint_difference.html b/docs/reference/footprint_difference.html new file mode 100644 index 0000000..d69dc47 --- /dev/null +++ b/docs/reference/footprint_difference.html @@ -0,0 +1,96 @@ + +footprint_difference — footprint_difference • methylTFR + Skip to contents + + +
+
+
+ +

footprint_difference

+ Source: R/plots.R +
footprint_difference.Rd
+
+ +
+

Combine the observed and expected profiles into a single +corrected curve.

+
+ +
+

Usage

+ 
footprint_difference(plot_data, method, sample_name)
+
+ +
+

Arguments

+ + +
plot_data
+

A data.table with the columns x, +avg_methyl and type.

+ + +
method
+

Either "substraction" or "division".

+ + +
sample_name
+

Sample label used in the curve label.

+ +
+
+

Value

+

A list with the corrected data and the axis +lab.

+
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/granges_helper.html b/docs/reference/granges_helper.html index 5990e22..2658db7 100644 --- a/docs/reference/granges_helper.html +++ b/docs/reference/granges_helper.html @@ -7,7 +7,7 @@ methylTFR - 0.99.4 + 0.99.5 + + + + + +
+
+
+ +

group_mean_difference

+ Source: R/differential_analysis.R +
group_mean_difference.Rd
+
+ +
+

Unsigned effect size: the difference of the group means for +two groups, and their range for more.

+
+ +
+

Usage

+ 
group_mean_difference(deviations, groups)
+
+ +
+

Arguments

+ + +
deviations
+

A matrix of deviation scores, motifs in rows.

+ + +
groups
+

The group labels as a factor.

+ +
+
+

Value

+

A numeric vector with one value per motif.

+
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/index.html b/docs/reference/index.html index 77e268e..8689181 100644 --- a/docs/reference/index.html +++ b/docs/reference/index.html @@ -7,7 +7,7 @@ methylTFR - 0.99.4 + 0.99.5 + + + + + +
+
+
+ +

locate_sample_files

+ Source: R/run_methyltfr.R +
locate_sample_files.Rd
+
+ +
+

Build and check the list of per-sample methylation files.

+
+ +
+

Usage

+ 
locate_sample_files(samples, sample_dir, sampleColName, full_path)
+
+ +
+

Arguments

+ + +
samples
+

The sample annotation as a data.frame.

+ + +
sample_dir
+

Directory holding the methylation call files.

+ + +
sampleColName
+

Column name holding the file names.

+ + +
full_path
+

if TRUE, the annotation file holds full paths.

+ +
+
+

Value

+

A character vector of existing file paths.

+
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/methylTFRdeviations-class.html b/docs/reference/methylTFRdeviations-class.html index 2ba6fab..1450538 100644 --- a/docs/reference/methylTFRdeviations-class.html +++ b/docs/reference/methylTFRdeviations-class.html @@ -9,7 +9,7 @@ methylTFR - 0.99.4 + 0.99.5 + + + + + +
+
+
+ +

normalise_footprint_flank

+ Source: R/plots.R +
normalise_footprint_flank.Rd
+
+ +
+

Normalise a corrected footprint against its outer flanking +windows, on the same scale as the statistic itself.

+
+ +
+

Usage

+ 
normalise_footprint_flank(difference_data, method, flankNorm)
+
+ +
+

Arguments

+ + +
difference_data
+

The corrected footprint as a data.table.

+ + +
method
+

Either "substraction" or "division".

+ + +
flankNorm
+

Width of the flanking window used for normalisation.

+ +
+
+

Value

+

The normalised data.table.

+
+
+

Details

+

Dividing a difference by its flank mean would rescale it by an +arbitrary factor, because that mean is near zero: for a typical footprint +it inflates the curve by an order of magnitude and can invert its sign.

+
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/plotExpectedFootprint.html b/docs/reference/plotExpectedFootprint.html index 02575b0..749ddb6 100644 --- a/docs/reference/plotExpectedFootprint.html +++ b/docs/reference/plotExpectedFootprint.html @@ -11,7 +11,7 @@ methylTFR - 0.99.4 + 0.99.5 + + + + + +
+
+
+ +

process_core_sample

+ Source: R/methyltfr_core.R +
process_core_sample.Rd
+
+ +
+

Compute and write the deviations of one sample, one motif +chunk at a time.

+
+ +
+

Usage

+ 
process_core_sample(
+  index,
+  sample_ids,
+  msites_fun,
+  motif_chunks,
+  tf_bindsites,
+  gcfreqs,
+  gc_dist,
+  dev_grid,
+  exp_grid,
+  dev_sink,
+  exp_sink,
+  BPPARAM,
+  enhancer,
+  ignoreStrand
+)
+
+ +
+

Arguments

+ + +
index
+

Integer index of the sample within sample_ids.

+ + +
sample_ids
+

A character vector of sample identifiers.

+ + +
msites_fun
+

A function of a single integer sample index.

+ + +
motif_chunks
+

A list of character vectors of motif names.

+ + +
tf_bindsites
+

a GRangesList of TF binding site positions.

+ + +
gcfreqs
+

a list of GC bin frequency tables.

+ + +
gc_dist
+

a GRanges of the genome-wide GC distribution.

+ + +
dev_grid, exp_grid
+

The grids the blocks are written on.

+ + +
dev_sink, exp_sink
+

The sinks the blocks are written to.

+ + +
BPPARAM
+

A BiocParallelParam object.

+ + +
enhancer
+

a GRanges restricting the analysis (optional).

+ + +
ignoreStrand
+

if TRUE, strand information is ignored.

+ +
+
+

Value

+

Invisible NULL. Called for its effect on the sinks.

+
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/rbind-methylTFRdeviations-method.html b/docs/reference/rbind-methylTFRdeviations-method.html index 573edba..77c93d5 100644 --- a/docs/reference/rbind-methylTFRdeviations-method.html +++ b/docs/reference/rbind-methylTFRdeviations-method.html @@ -7,7 +7,7 @@ methylTFR - 0.99.4 + 0.99.5 + + + + + +
+
+
+ +

resolve_annotation_file

+ Source: R/run_methyltfr.R +
resolve_annotation_file.Rd
+
+ +
+

Work out the path of the sample annotation file.

+
+ +
+

Usage

+ 
resolve_annotation_file(annfile, sample_ann, sample_dir)
+
+ +
+

Arguments

+ + +
annfile
+

Explicit path to the annotation file, or NULL.

+ + +
sample_ann
+

Name of the annotation file inside sample_dir.

+ + +
sample_dir
+

Directory holding the methylation call files.

+ +
+
+

Value

+

The path of the annotation file as a character string.

+
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/resolve_diff_groups.html b/docs/reference/resolve_diff_groups.html new file mode 100644 index 0000000..36eb8a0 --- /dev/null +++ b/docs/reference/resolve_diff_groups.html @@ -0,0 +1,90 @@ + +resolve_diff_groups — resolve_diff_groups • methylTFR + Skip to contents + + +
+
+
+ +

resolve_diff_groups

+ Source: R/differential_analysis.R +
resolve_diff_groups.Rd
+
+ +
+

Derive and validate the group labels used by +differential_deviation_test.

+
+ +
+

Usage

+ 
resolve_diff_groups(deviations, groups)
+
+ +
+

Arguments

+ + +
deviations
+

A matrix of deviation scores, motifs in rows.

+ + +
groups
+

Group labels, or NULL to take them from the column names.

+ +
+
+

Value

+

The group labels as a factor.

+
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/resolve_rnb_sample_ann.html b/docs/reference/resolve_rnb_sample_ann.html new file mode 100644 index 0000000..a88ba2c --- /dev/null +++ b/docs/reference/resolve_rnb_sample_ann.html @@ -0,0 +1,94 @@ + +resolve_rnb_sample_ann — resolve_rnb_sample_ann • methylTFR + Skip to contents + + +
+
+
+ +

resolve_rnb_sample_ann

+ Source: R/rnbeads_interface.R +
resolve_rnb_sample_ann.Rd
+
+ +
+

Default the sample annotation to the phenotype table of the +RnBeads object and check its shape.

+
+ +
+

Usage

+ 
resolve_rnb_sample_ann(rnb_set, sample_ann, sample_ids)
+
+ +
+

Arguments

+ + +
rnb_set
+

A preprocessed RnBSet object.

+ + +
sample_ann
+

Sample annotation supplied by the caller, or NULL.

+ + +
sample_ids
+

Character vector of sample identifiers.

+ +
+
+

Value

+

The sample annotation as a data.frame.

+
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/rnb_column.html b/docs/reference/rnb_column.html index 522fa57..099754f 100644 --- a/docs/reference/rnb_column.html +++ b/docs/reference/rnb_column.html @@ -11,7 +11,7 @@ methylTFR - 0.99.4 + 0.99.5 + + + + + +
+
+
+ +

valid_core_motifs

+ Source: R/methyltfr_core.R +
valid_core_motifs.Rd
+
+ +
+

Drop motifs whose binding sites are empty or whose GC bin +frequency matrix is missing.

+
+ +
+

Usage

+ 
valid_core_motifs(tf_bindsites, gcfreqs)
+
+ +
+

Arguments

+ + +
tf_bindsites
+

a GRangesList of TF binding site positions.

+ + +
gcfreqs
+

a list of GC bin frequency tables.

+ +
+
+

Value

+

A character vector of the motif names that can be processed.

+
+ +
+ + +
+

Developed by Irem B. Gündüz, Sarath Kumar Murugan, Fabian Muller.

+
+ +
+

Site built with pkgdown 2.1.3.

+
+ +
+ + + + + + + diff --git a/docs/reference/variability_test.html b/docs/reference/variability_test.html new file mode 100644 index 0000000..833176c --- /dev/null +++ b/docs/reference/variability_test.html @@ -0,0 +1,95 @@ + +variability_test — variability_test • methylTFR + Skip to contents + + +
+
+
+ +

variability_test

+ Source: R/variability.R +
variability_test.Rd
+
+ +
+

Summarise the calibrated Z-scores into a variability per +motif and test each against the chi-squared null.

+
+ +
+

Usage

+ 
variability_test(z, motif_names, padjMethod)
+
+ +
+

Arguments

+ + +
z
+

The calibrated Z-score matrix.

+ + +
motif_names
+

Character vector of motif names.

+ + +
padjMethod
+

Multiple testing correction passed to +stats::p.adjust.

+ +
+
+

Value

+

A data.frame with one row per motif.

+
+ +
+ + +
+

Developed by Irem B. Gündüz, Sarath Kumar Murugan, Fabian Muller.

+
+ +
+

Site built with pkgdown 2.1.3.

+
+ +
+ + + + + + + diff --git a/docs/reference/wilcoxon_helper.html b/docs/reference/wilcoxon_helper.html index 631d838..58a0a4c 100644 --- a/docs/reference/wilcoxon_helper.html +++ b/docs/reference/wilcoxon_helper.html @@ -7,7 +7,7 @@ methylTFR - 0.99.4 + 0.99.5