Devel

.gitignore

@@ -18,6 +18,7 @@ methylTFR.BiocCheck/*
 .Trashes
 .fseventsd
 .TemporaryItems
+methylTFR_tmp/
 
 # R session artefacts
 .Rhistory

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)

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)

NEWS.md

@@ -1,4 +1,4 @@
-# methylTFR 0.99.4
+# methylTFR 0.99.6
 
 NEW FEATURES
 

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,

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,

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),

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(