From 079f79fabde0fadd86bc6b05e1d5694169f806d9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=E2=80=9Cigunduz=E2=80=9D?= <irembgunduz@gmail.com>
Date: Sun, 30 Aug 2026 15:13:04 +0200
Subject: [PATCH 1/4] array support for rnbeads wrapper

---
 R/rnbeads_interface.R        | 394 +++++++++++++++++------------------
 man/rnb_annotation_table.Rd  |  27 +++
 man/rnb_annotation_target.Rd |  19 ++
 man/rnb_column.Rd            |   7 +-
 man/rnb_has_coverage.Rd      |   3 +-
 man/rnb_has_dpval.Rd         |  19 ++
 man/rnb_quality_mode.Rd      |  22 ++
 man/rnb_sample_ids.Rd        |   8 -
 man/rnb_sample_msites.Rd     |  25 ++-
 man/rnb_sites_to_granges.Rd  |  10 +-
 man/run_methylTFR_RnBeads.Rd | 117 ++---------
 11 files changed, 326 insertions(+), 325 deletions(-)
 create mode 100644 man/rnb_annotation_table.Rd
 create mode 100644 man/rnb_annotation_target.Rd
 create mode 100644 man/rnb_has_dpval.Rd
 create mode 100644 man/rnb_quality_mode.Rd

diff --git a/R/rnbeads_interface.R b/R/rnbeads_interface.R
index e453608..29ab37c 100644
--- a/R/rnbeads_interface.R
+++ b/R/rnbeads_interface.R
@@ -19,6 +19,7 @@ check_rnb_inputs <- function(rnb_set) {
     invisible(NULL)
 }
 
+
 #' @title resolve_rnb_sample_ann
 #' @description Default the sample annotation to the phenotype table of the
 #' RnBeads object and check its shape.
@@ -42,130 +43,47 @@ resolve_rnb_sample_ann <- function(rnb_set, sample_ann, sample_ids) {
     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.
 #'
-#' This is the RnBeads-based counterpart to \code{\link{run_methyltfr}}. Both
-#' functions share the same engine and produce numerically identical results
-#' for the same underlying methylation calls; they differ only in where the
-#' per-sample methylation levels come from.
-#'
 #' @details
-#' Methylation calls are always read at single-cytosine resolution
-#' (\code{type = "sites"}). Region-level summaries such as \code{tiling1kb} or
-#' \code{distal} cannot be used, because methylTFR needs base-resolution calls
-#' to build the footprint around each motif centre. To restrict the analysis to
-#' a set of regulatory regions, pass those regions through the \code{enhancer}
-#' argument instead.
-#'
-#' Samples are processed one at a time and methylation levels are pulled from
-#' the RnBeads object column by column, so disk-backed (\code{ff}-managed)
-#' RnBeads sets are never loaded into memory in full.
+#' Methylation calls are read at single-cytosine resolution. Sequencing
+#' sets (\code{RnBiseqSet}) are filtered by coverage, array sets
+#' (\code{RnBeadSet}) by detection p-value. The number of sites retained
+#' per sample is reported through \pkg{logger}.
 #'
-#' Coverage filtering is applied only when the object carries coverage
-#' information, which is the case for sequencing-based sets
-#' (\code{RnBiseqSet}). For array-based sets \code{cov_threshold} is ignored
-#' and a message is emitted.
-#'
-#' Note that RnBeads site annotation is 1-based while
-#' \code{\link{read_methylome}} reads 0-based BED coordinates as-is. The
-#' resulting one-base offset is not corrected here, since deviation scores
-#' aggregate methylation over windows of tens to hundreds of bases and are
-#' insensitive to a uniform single-base shift.
-#'
-#' @param rnb_set A preprocessed \code{RnBSet} object, for example the output
-#' of \code{rnb.run.preprocessing} or a set loaded with
-#' \code{RnBeads::load.rnb.set}.
-#' @param tf_bindsites a \code{GRangesList} object contains
-#'  tf binding sites positions
-#' @param gcfreqs a \code{list} of GC bin frequency tables
-#'  (matrices for multiple motif)
-#' @param gc_dist a \code{GRanges} object contains
-#' Genome wide GC distribution
-#' @param chunkSize Chunk size for parallel processing
-#'  of motifs (default: 20)
-#' @param threads Thread count for parallel processing
-#' @param enhancer a \code{GRanges} object specifying
-#' regions such as distal regulatory elements (optional)
-#' @param ignoreStrand if TRUE, it ignores strand info from annotation
-#' @param cov_threshold numeric, coverage threshold used to filter out low
-#' coverage sites, default is 1. Ignored for objects without coverage
-#' information.
-#' @param sample_ann Optional \code{data.frame} of sample annotation with one
-#' row per sample, used as \code{colData}. Defaults to
-#' \code{RnBeads::pheno(rnb_set)}.
-#' @return a \code{methylTFRdeviations} object with
-#' bias-corrected deviation and Z-scores
-#' @importFrom GenomicRanges GRanges
-#' @importFrom IRanges IRanges
-#' @importFrom logger log_info log_warn
+#' @param rnb_set A preprocessed \code{RnBSet} object.
+#' @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 chunkSize Chunk size for parallel processing of motifs.
+#' @param threads Thread count for parallel processing.
+#' @param enhancer an optional \code{GRanges} of regions to restrict to.
+#' @param ignoreStrand if TRUE, strand information is ignored.
+#' @param cov_threshold numeric, minimum coverage of a retained site.
+#' @param dpval_threshold numeric, maximum detection p-value of a
+#' retained probe.
+#' @param sample_ann Optional \code{data.frame} of sample annotation.
+#' @return a \code{methylTFRdeviations} object with bias-corrected
+#' deviations and Z-scores.
+#' @importFrom logger log_info
 #' @importFrom methods is
-#' @seealso \code{\link{run_methyltfr}} for the file-based entry point.
-#' @author Irem Gunduz
-#' @examples
-#' # 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
+    cov_threshold = 1, dpval_threshold = 0.05, sample_ann = NULL
 ) {
     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)
-
     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)
-    if (!has_covg) {
-        log_warn(
-            "The RnBSet object does not carry coverage information; ",
-            "cov_threshold is ignored."
-        )
-    }
+    mode <- rnb_quality_mode(rnb_set, dpval_threshold)
     log_info(
         "Found ", length(sites_gr), " sites across ",
         length(sample_ids), " samples"
@@ -173,24 +91,17 @@ run_methylTFR_RnBeads <- function(
 
     msites_fun <- function(i) {
         rnb_sample_msites(
-            rnb_set = rnb_set,
-            sites_gr = sites_gr,
-            index = i,
-            cov_threshold = opts$cov_threshold,
-            has_covg = has_covg
+            rnb_set = rnb_set, sites_gr = sites_gr, index = i,
+            cov_threshold = opts$cov_threshold, has_covg = mode$has_covg,
+            dpval_threshold = dpval_threshold, has_dpval = mode$has_dpval
         )
     }
 
     methyltfr_core(
-        sample_ids = sample_ids,
-        msites_fun = msites_fun,
-        samples = sample_ann,
-        tf_bindsites = tf_bindsites,
-        gcfreqs = gcfreqs,
-        gc_dist = gc_dist,
-        chunkSize = opts$chunkSize,
-        threads = opts$threads,
-        enhancer = enhancer,
+        sample_ids = sample_ids, msites_fun = msites_fun,
+        samples = sample_ann, tf_bindsites = tf_bindsites,
+        gcfreqs = gcfreqs, gc_dist = gc_dist, chunkSize = opts$chunkSize,
+        threads = opts$threads, enhancer = enhancer,
         ignoreStrand = opts$ignoreStrand
     )
 }
@@ -198,12 +109,6 @@ run_methylTFR_RnBeads <- function(
 
 #' @title rnb_sample_ids
 #' @description Determine the sample identifiers of an RnBeads object.
-#' @details RnBeads exports \code{samples} with \code{exportMethods} rather
-#' than \code{export}, so the generic is not reachable as
-#' \code{RnBeads::samples} and referring to it that way fails
-#' \code{R CMD check}. The identifiers are therefore taken from the column
-#' names of the methylation matrix, falling back to the row names of the
-#' sample annotation.
 #' @param rnb_set An \code{RnBSet} object.
 #' @return A character vector of sample identifiers.
 #' @keywords internal
@@ -212,63 +117,110 @@ rnb_sample_ids <- function(rnb_set) {
         colnames(RnBeads::meth(rnb_set, type = "sites", i = 1L)),
         error = function(e) NULL
     )
-    if (is.null(ids) || length(ids) == 0) {
-        ids <- tryCatch(rownames(RnBeads::pheno(rnb_set)),
+    if (length(ids) == 0) {
+        ids <- tryCatch(
+            rownames(RnBeads::pheno(rnb_set)),
             error = function(e) NULL
         )
     }
-    if (is.null(ids) || length(ids) == 0) {
-        nsamples <- tryCatch(nrow(RnBeads::pheno(rnb_set)),
+    if (length(ids) == 0) {
+        nsamples <- tryCatch(
+            nrow(RnBeads::pheno(rnb_set)),
             error = function(e) 0L
         )
-        if (!is.null(nsamples) && nsamples > 0) {
+        if (length(nsamples) == 1 && nsamples > 0) {
             ids <- paste0("sample_", seq_len(nsamples))
         }
     }
-    if (is.null(ids) || length(ids) == 0) {
+    if (length(ids) == 0) {
+        stop("Could not determine sample identifiers from the RnBSet object.")
+    }
+    return(as.character(ids))
+}
+
+
+#' @title rnb_annotation_target
+#' @description Resolve the annotation target of an RnBeads object.
+#' @param rnb_set An \code{RnBSet} object.
+#' @return A character scalar naming the annotation target, \code{"sites"}
+#' for sequencing sets and the array platform for array sets.
+#' @importFrom methods is
+#' @keywords internal
+rnb_annotation_target <- function(rnb_set) {
+    if (is(rnb_set, "RnBeadSet")) rnb_set@target else "sites"
+}
+
+
+#' @title rnb_annotation_table
+#' @description Look up the site or probe annotation of an RnBeads object.
+#' @details The annotation stored in the object is used when available,
+#' otherwise the genome-wide track registered for \code{target} is
+#' subset to the sites of the object.
+#' @param rnb_set An \code{RnBSet} object.
+#' @param target Character scalar naming the annotation target.
+#' @param assembly Character scalar naming the genome assembly.
+#' @return A \code{data.frame} with one row per site or probe.
+#' @keywords internal
+rnb_annotation_table <- function(rnb_set, target, assembly) {
+    pull <- function(expr) tryCatch(expr, error = function(e) NULL)
+
+    ann <- pull(RnBeads::annotation(rnb_set, type = target))
+    if (is.null(ann)) {
+        ann <- pull(RnBeads::annotation(rnb_set, type = "sites"))
+    }
+    if (is.null(ann) && target != "sites") {
+        ann <- pull({
+            track <- RnBeads::rnb.get.annotation(target, assembly)
+            RnBeads::rnb.annotation2data.frame(track)[rnb_set@sites, ]
+        })
+    }
+    if (is.null(ann) || nrow(ann) == 0) {
         stop(
-            "Could not determine sample identifiers from the RnBSet object; ",
-            "pass them explicitly via sample_ann."
+            "Could not extract coordinates for target '", target,
+            "' and assembly '", assembly, "'. Load the matching RnBeads ",
+            "annotation, or a custom annotation for this array, before ",
+            "calling methylTFR."
         )
     }
-    return(as.character(ids))
+    return(ann)
 }
 
 
 #' @title rnb_sites_to_granges
-#' @description Build a \code{GRanges} object of the site annotation of an
-#' RnBeads object. The order of the ranges matches the row order of the
-#' methylation matrix returned by \code{RnBeads::meth}.
+#' @description Build a \code{GRanges} object of the site or probe
+#' annotation of an RnBeads object.
 #' @param rnb_set An \code{RnBSet} object.
-#' @param ignoreStrand if TRUE, all ranges are returned with strand \code{"*"}.
-#' @return A \code{GRanges} object with one range per site.
+#' @param ignoreStrand if TRUE, all ranges are returned with strand
+#' \code{"*"}.
+#' @return A \code{GRanges} object with one range per site or probe.
 #' @importFrom GenomicRanges GRanges
 #' @importFrom IRanges IRanges
+#' @importFrom logger log_info
+#' @importFrom methods is
 #' @keywords internal
 rnb_sites_to_granges <- function(rnb_set, ignoreStrand = TRUE) {
-    ann <- RnBeads::annotation(rnb_set, type = "sites")
-    if (is.null(ann) || nrow(ann) == 0) {
-        stop("The RnBSet object does not contain any site annotation")
-    }
-    required <- c("Chromosome", "Start")
-    missing_cols <- setdiff(required, colnames(ann))
+    assembly <- rnb_set@assembly
+    target <- rnb_annotation_target(rnb_set)
+    log_info("Annotation target: ", target, " | assembly: ", assembly)
+    invisible(requireNamespace(paste0("RnBeads.", assembly), quietly = TRUE))
+
+    ann <- rnb_annotation_table(rnb_set, target, assembly)
+    missing_cols <- setdiff(c("Chromosome", "Start"), colnames(ann))
     if (length(missing_cols) > 0) {
-        stop(
-            "Unexpected RnBeads site annotation, missing column(s): ",
-            paste(missing_cols, collapse = ", ")
-        )
+        stop("RnBeads annotation misses column(s): ", toString(missing_cols))
     }
+
     ends <- if ("End" %in% colnames(ann)) ann$End else ann$Start
     strands <- "*"
     if (!ignoreStrand && "Strand" %in% colnames(ann)) {
         strands <- as.character(ann$Strand)
         strands[is.na(strands) | !strands %in% c("+", "-")] <- "*"
     }
+
     GenomicRanges::GRanges(
         seqnames = as.character(ann$Chromosome),
         ranges = IRanges::IRanges(
-            start = as.integer(ann$Start),
-            end = as.integer(ends)
+            start = as.integer(ann$Start), end = as.integer(ends)
         ),
         strand = strands
     )
@@ -276,13 +228,12 @@ rnb_sites_to_granges <- function(rnb_set, ignoreStrand = TRUE) {
 
 
 #' @title rnb_has_coverage
-#' @description Test whether an RnBeads object carries coverage information.
+#' @description Test whether an RnBeads object carries coverage
+#' information.
 #' @param rnb_set An \code{RnBSet} object.
 #' @return A logical scalar.
 #' @keywords internal
 rnb_has_coverage <- function(rnb_set) {
-    # Try the subsetting form first so that large disk-backed sets are not
-    # materialised, then fall back for RnBeads versions without the j argument.
     res <- tryCatch(
         !is.null(RnBeads::covg(rnb_set, type = "sites", j = 1L)),
         error = function(e) NULL
@@ -297,63 +248,54 @@ rnb_has_coverage <- function(rnb_set) {
 }
 
 
-#' @title rnb_sample_msites
-#' @description Extract the methylation calls of a single sample from an
-#' RnBeads object as a \code{GRanges} object in the layout expected by
-#' \code{\link{computeDeviation}}.
+#' @title rnb_has_dpval
+#' @description Test whether an RnBeads object carries detection
+#' p-values.
 #' @param rnb_set An \code{RnBSet} object.
-#' @param sites_gr A \code{GRanges} object of site positions, as returned by
-#' \code{rnb_sites_to_granges}.
-#' @param index Integer index of the sample to extract.
-#' @param cov_threshold numeric coverage threshold.
-#' @param has_covg logical, whether the object carries coverage information.
-#' @return A \code{GRanges} object with \code{score} and \code{coverage}
-#' metadata columns, restricted to sites with a non-missing methylation call.
-#' @importFrom logger log_warn
+#' @return A logical scalar.
 #' @keywords internal
-rnb_sample_msites <- function(
-    rnb_set, sites_gr, index,
-    cov_threshold = 1, has_covg = TRUE
-) {
-    index <- as.integer(index)
-    mvals <- rnb_column(RnBeads::meth, rnb_set, index)
-    if (length(mvals) != length(sites_gr)) {
-        stop(
-            "The number of methylation values does not match the number of ",
-            "annotated sites; the RnBSet object appears to be inconsistent."
+rnb_has_dpval <- function(rnb_set) {
+    res <- tryCatch(
+        !is.null(RnBeads::dpval(rnb_set, type = "sites", j = 1L)),
+        error = function(e) NULL
+    )
+    if (is.null(res)) {
+        res <- tryCatch(
+            !is.null(RnBeads::dpval(rnb_set, type = "sites")),
+            error = function(e) FALSE
         )
     }
-    keep <- !is.na(mvals)
-    if (has_covg) {
-        cvals <- rnb_column(RnBeads::covg, rnb_set, index)
-        if (length(cvals) != length(sites_gr)) {
-            stop(
-                "The number of coverage values does not match the number of ",
-                "annotated sites."
-            )
-        }
-        keep <- keep & !is.na(cvals) & cvals >= cov_threshold
-    } else {
-        cvals <- rep(NA_real_, length(sites_gr))
-    }
-    if (!any(keep)) {
-        stop(
-            "No sites passed the coverage threshold for sample index ", index
+    return(isTRUE(res))
+}
+
+
+#' @title rnb_quality_mode
+#' @description Decide which per-site quality filter applies to an
+#' RnBeads object and report it.
+#' @param rnb_set An \code{RnBSet} object.
+#' @param dpval_threshold numeric detection p-value threshold.
+#' @return A list with the logical flags \code{has_covg} and
+#' \code{has_dpval}.
+#' @importFrom logger log_info log_warn
+#' @keywords internal
+rnb_quality_mode <- function(rnb_set, dpval_threshold) {
+    has_covg <- rnb_has_coverage(rnb_set)
+    has_dpval <- rnb_has_dpval(rnb_set)
+    if (has_dpval) {
+        log_info("Filtering probes at detection p-value <= ", dpval_threshold)
+    } else if (!has_covg) {
+        log_warn(
+            "RnBSet carries neither coverage nor detection p-values, ",
+            "no quality filtering is applied"
         )
     }
-    gr <- sites_gr[keep]
-    gr$score <- as.numeric(mvals[keep])
-    gr$coverage <- as.numeric(cvals[keep])
-    return(gr)
+    list(has_covg = has_covg, has_dpval = has_dpval)
 }
 
 
 #' @title rnb_column
-#' @description Extract a single sample column from an RnBeads accessor,
-#' falling back to full extraction on RnBeads versions that do not support
-#' column subsetting.
-#' @param accessor An RnBeads accessor function, either \code{RnBeads::meth} or
-#' \code{RnBeads::covg}.
+#' @description Extract a single sample column from an RnBeads accessor.
+#' @param accessor An RnBeads accessor function.
 #' @param rnb_set An \code{RnBSet} object.
 #' @param index Integer index of the sample to extract.
 #' @return A numeric vector with one value per site.
@@ -371,3 +313,53 @@ rnb_column <- function(accessor, rnb_set, index) {
     }
     return(as.numeric(as.vector(vals)))
 }
+
+
+#' @title rnb_sample_msites
+#' @description Extract the methylation calls of a single sample from an
+#' RnBeads object as a \code{GRanges} object.
+#' @param rnb_set An \code{RnBSet} object.
+#' @param sites_gr A \code{GRanges} object of site positions.
+#' @param index Integer index of the sample to extract.
+#' @param cov_threshold numeric coverage threshold.
+#' @param has_covg logical, whether coverage filtering applies.
+#' @param dpval_threshold numeric detection p-value threshold.
+#' @param has_dpval logical, whether detection p-value filtering applies.
+#' @return A \code{GRanges} object restricted to valid methylation calls.
+#' @importFrom logger log_info
+#' @keywords internal
+rnb_sample_msites <- function(
+    rnb_set, sites_gr, index, cov_threshold = 1, has_covg = TRUE,
+    dpval_threshold = 0.05, has_dpval = FALSE
+) {
+    index <- as.integer(index)
+    mvals <- rnb_column(RnBeads::meth, rnb_set, index)
+    if (length(mvals) != length(sites_gr)) {
+        stop("Methylation values do not match annotated sites.")
+    }
+
+    keep <- !is.na(mvals)
+    if (has_covg) {
+        cvals <- rnb_column(RnBeads::covg, rnb_set, index)
+        keep <- keep & !is.na(cvals) & cvals >= cov_threshold
+    } else {
+        cvals <- rep(NA_real_, length(sites_gr))
+    }
+    if (has_dpval) {
+        dvals <- rnb_column(RnBeads::dpval, rnb_set, index)
+        keep <- keep & !is.na(dvals) & dvals <= dpval_threshold
+    }
+
+    log_info(
+        "Sample ", index, ": ", sum(keep), " of ", length(keep),
+        " sites retained (", round(100 * mean(keep), 2), "%)"
+    )
+    if (!any(keep)) {
+        stop("No sites passed the quality thresholds for sample index ", index)
+    }
+
+    gr <- sites_gr[keep]
+    gr$score <- as.numeric(mvals[keep])
+    gr$coverage <- as.numeric(cvals[keep])
+    return(gr)
+}
diff --git a/man/rnb_annotation_table.Rd b/man/rnb_annotation_table.Rd
new file mode 100644
index 0000000..ccb9b89
--- /dev/null
+++ b/man/rnb_annotation_table.Rd
@@ -0,0 +1,27 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/rnbeads_interface.R
+\name{rnb_annotation_table}
+\alias{rnb_annotation_table}
+\title{rnb_annotation_table}
+\usage{
+rnb_annotation_table(rnb_set, target, assembly)
+}
+\arguments{
+\item{rnb_set}{An \code{RnBSet} object.}
+
+\item{target}{Character scalar naming the annotation target.}
+
+\item{assembly}{Character scalar naming the genome assembly.}
+}
+\value{
+A \code{data.frame} with one row per site or probe.
+}
+\description{
+Look up the site or probe annotation of an RnBeads object.
+}
+\details{
+The annotation stored in the object is used when available,
+otherwise the genome-wide track registered for \code{target} is
+subset to the sites of the object.
+}
+\keyword{internal}
diff --git a/man/rnb_annotation_target.Rd b/man/rnb_annotation_target.Rd
new file mode 100644
index 0000000..c78696f
--- /dev/null
+++ b/man/rnb_annotation_target.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/rnbeads_interface.R
+\name{rnb_annotation_target}
+\alias{rnb_annotation_target}
+\title{rnb_annotation_target}
+\usage{
+rnb_annotation_target(rnb_set)
+}
+\arguments{
+\item{rnb_set}{An \code{RnBSet} object.}
+}
+\value{
+A character scalar naming the annotation target, \code{"sites"}
+for sequencing sets and the array platform for array sets.
+}
+\description{
+Resolve the annotation target of an RnBeads object.
+}
+\keyword{internal}
diff --git a/man/rnb_column.Rd b/man/rnb_column.Rd
index 06692b9..b1dcd18 100644
--- a/man/rnb_column.Rd
+++ b/man/rnb_column.Rd
@@ -7,8 +7,7 @@
 rnb_column(accessor, rnb_set, index)
 }
 \arguments{
-\item{accessor}{An RnBeads accessor function, either \code{RnBeads::meth} or
-\code{RnBeads::covg}.}
+\item{accessor}{An RnBeads accessor function.}
 
 \item{rnb_set}{An \code{RnBSet} object.}
 
@@ -18,8 +17,6 @@ rnb_column(accessor, rnb_set, index)
 A numeric vector with one value per site.
 }
 \description{
-Extract a single sample column from an RnBeads accessor,
-falling back to full extraction on RnBeads versions that do not support
-column subsetting.
+Extract a single sample column from an RnBeads accessor.
 }
 \keyword{internal}
diff --git a/man/rnb_has_coverage.Rd b/man/rnb_has_coverage.Rd
index 15ee028..81f6a73 100644
--- a/man/rnb_has_coverage.Rd
+++ b/man/rnb_has_coverage.Rd
@@ -13,6 +13,7 @@ rnb_has_coverage(rnb_set)
 A logical scalar.
 }
 \description{
-Test whether an RnBeads object carries coverage information.
+Test whether an RnBeads object carries coverage
+information.
 }
 \keyword{internal}
diff --git a/man/rnb_has_dpval.Rd b/man/rnb_has_dpval.Rd
new file mode 100644
index 0000000..49d3afd
--- /dev/null
+++ b/man/rnb_has_dpval.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/rnbeads_interface.R
+\name{rnb_has_dpval}
+\alias{rnb_has_dpval}
+\title{rnb_has_dpval}
+\usage{
+rnb_has_dpval(rnb_set)
+}
+\arguments{
+\item{rnb_set}{An \code{RnBSet} object.}
+}
+\value{
+A logical scalar.
+}
+\description{
+Test whether an RnBeads object carries detection
+p-values.
+}
+\keyword{internal}
diff --git a/man/rnb_quality_mode.Rd b/man/rnb_quality_mode.Rd
new file mode 100644
index 0000000..1f4ca2a
--- /dev/null
+++ b/man/rnb_quality_mode.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/rnbeads_interface.R
+\name{rnb_quality_mode}
+\alias{rnb_quality_mode}
+\title{rnb_quality_mode}
+\usage{
+rnb_quality_mode(rnb_set, dpval_threshold)
+}
+\arguments{
+\item{rnb_set}{An \code{RnBSet} object.}
+
+\item{dpval_threshold}{numeric detection p-value threshold.}
+}
+\value{
+A list with the logical flags \code{has_covg} and
+\code{has_dpval}.
+}
+\description{
+Decide which per-site quality filter applies to an
+RnBeads object and report it.
+}
+\keyword{internal}
diff --git a/man/rnb_sample_ids.Rd b/man/rnb_sample_ids.Rd
index 97046aa..00d2551 100644
--- a/man/rnb_sample_ids.Rd
+++ b/man/rnb_sample_ids.Rd
@@ -15,12 +15,4 @@ A character vector of sample identifiers.
 \description{
 Determine the sample identifiers of an RnBeads object.
 }
-\details{
-RnBeads exports \code{samples} with \code{exportMethods} rather
-than \code{export}, so the generic is not reachable as
-\code{RnBeads::samples} and referring to it that way fails
-\code{R CMD check}. The identifiers are therefore taken from the column
-names of the methylation matrix, falling back to the row names of the
-sample annotation.
-}
 \keyword{internal}
diff --git a/man/rnb_sample_msites.Rd b/man/rnb_sample_msites.Rd
index 6331c3b..b595f12 100644
--- a/man/rnb_sample_msites.Rd
+++ b/man/rnb_sample_msites.Rd
@@ -4,27 +4,36 @@
 \alias{rnb_sample_msites}
 \title{rnb_sample_msites}
 \usage{
-rnb_sample_msites(rnb_set, sites_gr, index, cov_threshold = 1, has_covg = TRUE)
+rnb_sample_msites(
+  rnb_set,
+  sites_gr,
+  index,
+  cov_threshold = 1,
+  has_covg = TRUE,
+  dpval_threshold = 0.05,
+  has_dpval = FALSE
+)
 }
 \arguments{
 \item{rnb_set}{An \code{RnBSet} object.}
 
-\item{sites_gr}{A \code{GRanges} object of site positions, as returned by
-\code{rnb_sites_to_granges}.}
+\item{sites_gr}{A \code{GRanges} object of site positions.}
 
 \item{index}{Integer index of the sample to extract.}
 
 \item{cov_threshold}{numeric coverage threshold.}
 
-\item{has_covg}{logical, whether the object carries coverage information.}
+\item{has_covg}{logical, whether coverage filtering applies.}
+
+\item{dpval_threshold}{numeric detection p-value threshold.}
+
+\item{has_dpval}{logical, whether detection p-value filtering applies.}
 }
 \value{
-A \code{GRanges} object with \code{score} and \code{coverage}
-metadata columns, restricted to sites with a non-missing methylation call.
+A \code{GRanges} object restricted to valid methylation calls.
 }
 \description{
 Extract the methylation calls of a single sample from an
-RnBeads object as a \code{GRanges} object in the layout expected by
-\code{\link{computeDeviation}}.
+RnBeads object as a \code{GRanges} object.
 }
 \keyword{internal}
diff --git a/man/rnb_sites_to_granges.Rd b/man/rnb_sites_to_granges.Rd
index 6dfd084..85c036d 100644
--- a/man/rnb_sites_to_granges.Rd
+++ b/man/rnb_sites_to_granges.Rd
@@ -9,14 +9,14 @@ rnb_sites_to_granges(rnb_set, ignoreStrand = TRUE)
 \arguments{
 \item{rnb_set}{An \code{RnBSet} object.}
 
-\item{ignoreStrand}{if TRUE, all ranges are returned with strand \code{"*"}.}
+\item{ignoreStrand}{if TRUE, all ranges are returned with strand
+\code{"*"}.}
 }
 \value{
-A \code{GRanges} object with one range per site.
+A \code{GRanges} object with one range per site or probe.
 }
 \description{
-Build a \code{GRanges} object of the site annotation of an
-RnBeads object. The order of the ranges matches the row order of the
-methylation matrix returned by \code{RnBeads::meth}.
+Build a \code{GRanges} object of the site or probe
+annotation of an RnBeads object.
 }
 \keyword{internal}
diff --git a/man/run_methylTFR_RnBeads.Rd b/man/run_methylTFR_RnBeads.Rd
index 659e2b5..a5a8b13 100644
--- a/man/run_methylTFR_RnBeads.Rd
+++ b/man/run_methylTFR_RnBeads.Rd
@@ -14,122 +14,45 @@ run_methylTFR_RnBeads(
   enhancer = NULL,
   ignoreStrand = TRUE,
   cov_threshold = 1,
+  dpval_threshold = 0.05,
   sample_ann = NULL
 )
 }
 \arguments{
-\item{rnb_set}{A preprocessed \code{RnBSet} object, for example the output
-of \code{rnb.run.preprocessing} or a set loaded with
-\code{RnBeads::load.rnb.set}.}
+\item{rnb_set}{A preprocessed \code{RnBSet} object.}
 
-\item{tf_bindsites}{a \code{GRangesList} object contains
-tf binding sites positions}
+\item{tf_bindsites}{a \code{GRangesList} of TF binding site positions.}
 
-\item{gcfreqs}{a \code{list} of GC bin frequency tables
-(matrices for multiple motif)}
+\item{gcfreqs}{a \code{list} of GC bin frequency tables.}
 
-\item{gc_dist}{a \code{GRanges} object contains
-Genome wide GC distribution}
+\item{gc_dist}{a \code{GRanges} of the genome-wide GC distribution.}
 
-\item{chunkSize}{Chunk size for parallel processing
-of motifs (default: 20)}
+\item{chunkSize}{Chunk size for parallel processing of motifs.}
 
-\item{threads}{Thread count for parallel processing}
+\item{threads}{Thread count for parallel processing.}
 
-\item{enhancer}{a \code{GRanges} object specifying
-regions such as distal regulatory elements (optional)}
+\item{enhancer}{an optional \code{GRanges} of regions to restrict to.}
 
-\item{ignoreStrand}{if TRUE, it ignores strand info from annotation}
+\item{ignoreStrand}{if TRUE, strand information is ignored.}
 
-\item{cov_threshold}{numeric, coverage threshold used to filter out low
-coverage sites, default is 1. Ignored for objects without coverage
-information.}
+\item{cov_threshold}{numeric, minimum coverage of a retained site.}
 
-\item{sample_ann}{Optional \code{data.frame} of sample annotation with one
-row per sample, used as \code{colData}. Defaults to
-\code{RnBeads::pheno(rnb_set)}.}
+\item{dpval_threshold}{numeric, maximum detection p-value of a
+retained probe.}
+
+\item{sample_ann}{Optional \code{data.frame} of sample annotation.}
 }
 \value{
-a \code{methylTFRdeviations} object with
-bias-corrected deviation and Z-scores
+a \code{methylTFRdeviations} object with bias-corrected
+deviations and Z-scores.
 }
 \description{
 Run the methylTFR workflow directly on a preprocessed
 \pkg{RnBeads} object, without exporting per-sample BED files first.
-
-This is the RnBeads-based counterpart to \code{\link{run_methyltfr}}. Both
-functions share the same engine and produce numerically identical results
-for the same underlying methylation calls; they differ only in where the
-per-sample methylation levels come from.
 }
 \details{
-Methylation calls are always read at single-cytosine resolution
-(\code{type = "sites"}). Region-level summaries such as \code{tiling1kb} or
-\code{distal} cannot be used, because methylTFR needs base-resolution calls
-to build the footprint around each motif centre. To restrict the analysis to
-a set of regulatory regions, pass those regions through the \code{enhancer}
-argument instead.
-
-Samples are processed one at a time and methylation levels are pulled from
-the RnBeads object column by column, so disk-backed (\code{ff}-managed)
-RnBeads sets are never loaded into memory in full.
-
-Coverage filtering is applied only when the object carries coverage
-information, which is the case for sequencing-based sets
-(\code{RnBiseqSet}). For array-based sets \code{cov_threshold} is ignored
-and a message is emitted.
-
-Note that RnBeads site annotation is 1-based while
-\code{\link{read_methylome}} reads 0-based BED coordinates as-is. The
-resulting one-base offset is not corrected here, since deviation scores
-aggregate methylation over windows of tens to hundreds of bases and are
-insensitive to a uniform single-base shift.
-}
-\examples{
-# 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)
-}
-}
-\seealso{
-\code{\link{run_methyltfr}} for the file-based entry point.
-}
-\author{
-Irem Gunduz
+Methylation calls are read at single-cytosine resolution. Sequencing
+sets (\code{RnBiseqSet}) are filtered by coverage, array sets
+(\code{RnBeadSet}) by detection p-value. The number of sites retained
+per sample is reported through \pkg{logger}.
 }

From 6cf3cca3aac1e50d84a9419ed100f9a0ffd5f2b5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=E2=80=9Cigunduz=E2=80=9D?= <irembgunduz@gmail.com>
Date: Sun, 30 Aug 2026 15:29:55 +0200
Subject: [PATCH 2/4] Fix small bug

---
 R/rnbeads_interface.R       | 89 ++++++++++++++++++++++++++++++++++---
 man/rnb_annotation_by_id.Rd | 29 ++++++++++++
 man/rnb_annotation_table.Rd |  2 +-
 man/rnb_probe_ids.Rd        | 22 +++++++++
 4 files changed, 136 insertions(+), 6 deletions(-)
 create mode 100644 man/rnb_annotation_by_id.Rd
 create mode 100644 man/rnb_probe_ids.Rd

diff --git a/R/rnbeads_interface.R b/R/rnbeads_interface.R
index 29ab37c..a8c4cfd 100644
--- a/R/rnbeads_interface.R
+++ b/R/rnbeads_interface.R
@@ -68,6 +68,42 @@ resolve_rnb_sample_ann <- function(rnb_set, sample_ann, sample_ids) {
 #' @param sample_ann Optional \code{data.frame} of sample annotation.
 #' @return a \code{methylTFRdeviations} object with bias-corrected
 #' deviations and Z-scores.
+#' @examples
+#' # A minimal end-to-end run on the BATF example data bundled with the
+#' # package. The bundled calls are wrapped in an RnBiseqSet so that the
+#' # example exercises the same code path as a preprocessed RnBeads object.
+#' 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"))
+#'
+#' if (requireNamespace("RnBeads", quietly = TRUE) &&
+#'     requireNamespace("RnBeads.hg38", quietly = TRUE)) {
+#'     sites <- data.frame(
+#'         chr = as.character(GenomicRanges::seqnames(msites)),
+#'         start = GenomicRanges::start(msites),
+#'         strand = "*",
+#'         stringsAsFactors = FALSE
+#'     )
+#'     rnb_set <- RnBeads::RnBiseqSet(
+#'         pheno = data.frame(sampleName = "sample_1"),
+#'         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)
+#' }
+#' @seealso \code{\link{run_methyltfr}} for running methylTFR from
+#' per-sample BED files.
 #' @importFrom logger log_info
 #' @importFrom methods is
 #' @export
@@ -151,11 +187,57 @@ rnb_annotation_target <- function(rnb_set) {
 }
 
 
+#' @title rnb_probe_ids
+#' @description Read the probe identifiers of an RnBeads object.
+#' @details The \code{sites} slot of an \code{RnBSet} is an index matrix
+#' whose row names carry the probe identifiers.
+#' @param rnb_set An \code{RnBSet} object.
+#' @return A character vector of probe identifiers, or \code{NULL}.
+#' @keywords internal
+rnb_probe_ids <- function(rnb_set) {
+    ids <- tryCatch(rownames(rnb_set@sites), error = function(e) NULL)
+    if (length(ids) == 0) {
+        ids <- tryCatch(
+            rownames(RnBeads::meth(rnb_set, type = "sites")),
+            error = function(e) NULL
+        )
+    }
+    if (length(ids) == 0) NULL else as.character(ids)
+}
+
+
+#' @title rnb_annotation_by_id
+#' @description Subset a genome-wide annotation track to the probes of an
+#' RnBeads object, matching on probe identifier.
+#' @details The \code{sites} slot holds a three-column index matrix, not a
+#' row index into the annotation, so the track is matched by identifier
+#' rather than by position.
+#' @param rnb_set An \code{RnBSet} object.
+#' @param target Character scalar naming the annotation target.
+#' @param assembly Character scalar naming the genome assembly.
+#' @return A \code{data.frame} of annotation rows, or \code{NULL} when the
+#' identifiers cannot be matched.
+#' @keywords internal
+rnb_annotation_by_id <- function(rnb_set, target, assembly) {
+    ids <- rnb_probe_ids(rnb_set)
+    if (is.null(ids)) {
+        return(NULL)
+    }
+    track <- RnBeads::rnb.get.annotation(target, assembly)
+    full <- RnBeads::rnb.annotation2data.frame(track)
+    idx <- match(ids, rownames(full))
+    if (anyNA(idx)) {
+        return(NULL)
+    }
+    full[idx, , drop = FALSE]
+}
+
+
 #' @title rnb_annotation_table
 #' @description Look up the site or probe annotation of an RnBeads object.
 #' @details The annotation stored in the object is used when available,
 #' otherwise the genome-wide track registered for \code{target} is
-#' subset to the sites of the object.
+#' matched to the probes of the object by identifier.
 #' @param rnb_set An \code{RnBSet} object.
 #' @param target Character scalar naming the annotation target.
 #' @param assembly Character scalar naming the genome assembly.
@@ -169,10 +251,7 @@ rnb_annotation_table <- function(rnb_set, target, assembly) {
         ann <- pull(RnBeads::annotation(rnb_set, type = "sites"))
     }
     if (is.null(ann) && target != "sites") {
-        ann <- pull({
-            track <- RnBeads::rnb.get.annotation(target, assembly)
-            RnBeads::rnb.annotation2data.frame(track)[rnb_set@sites, ]
-        })
+        ann <- pull(rnb_annotation_by_id(rnb_set, target, assembly))
     }
     if (is.null(ann) || nrow(ann) == 0) {
         stop(
diff --git a/man/rnb_annotation_by_id.Rd b/man/rnb_annotation_by_id.Rd
new file mode 100644
index 0000000..87880af
--- /dev/null
+++ b/man/rnb_annotation_by_id.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/rnbeads_interface.R
+\name{rnb_annotation_by_id}
+\alias{rnb_annotation_by_id}
+\title{rnb_annotation_by_id}
+\usage{
+rnb_annotation_by_id(rnb_set, target, assembly)
+}
+\arguments{
+\item{rnb_set}{An \code{RnBSet} object.}
+
+\item{target}{Character scalar naming the annotation target.}
+
+\item{assembly}{Character scalar naming the genome assembly.}
+}
+\value{
+A \code{data.frame} of annotation rows, or \code{NULL} when the
+identifiers cannot be matched.
+}
+\description{
+Subset a genome-wide annotation track to the probes of an
+RnBeads object, matching on probe identifier.
+}
+\details{
+The \code{sites} slot holds a three-column index matrix, not a
+row index into the annotation, so the track is matched by identifier
+rather than by position.
+}
+\keyword{internal}
diff --git a/man/rnb_annotation_table.Rd b/man/rnb_annotation_table.Rd
index ccb9b89..243429c 100644
--- a/man/rnb_annotation_table.Rd
+++ b/man/rnb_annotation_table.Rd
@@ -22,6 +22,6 @@ Look up the site or probe annotation of an RnBeads object.
 \details{
 The annotation stored in the object is used when available,
 otherwise the genome-wide track registered for \code{target} is
-subset to the sites of the object.
+matched to the probes of the object by identifier.
 }
 \keyword{internal}
diff --git a/man/rnb_probe_ids.Rd b/man/rnb_probe_ids.Rd
new file mode 100644
index 0000000..583fa62
--- /dev/null
+++ b/man/rnb_probe_ids.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/rnbeads_interface.R
+\name{rnb_probe_ids}
+\alias{rnb_probe_ids}
+\title{rnb_probe_ids}
+\usage{
+rnb_probe_ids(rnb_set)
+}
+\arguments{
+\item{rnb_set}{An \code{RnBSet} object.}
+}
+\value{
+A character vector of probe identifiers, or \code{NULL}.
+}
+\description{
+Read the probe identifiers of an RnBeads object.
+}
+\details{
+The \code{sites} slot of an \code{RnBSet} is an index matrix
+whose row names carry the probe identifiers.
+}
+\keyword{internal}

From 21a96cc70e1e4ddd5259b0af51c9b4ac5b549bba Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=E2=80=9Cigunduz=E2=80=9D?= <irembgunduz@gmail.com>
Date: Sun, 30 Aug 2026 15:31:12 +0200
Subject: [PATCH 3/4] update

---
 man/run_methylTFR_RnBeads.Rd | 39 ++++++++++++++++++++++++++++++++++++
 1 file changed, 39 insertions(+)

diff --git a/man/run_methylTFR_RnBeads.Rd b/man/run_methylTFR_RnBeads.Rd
index a5a8b13..69e1f51 100644
--- a/man/run_methylTFR_RnBeads.Rd
+++ b/man/run_methylTFR_RnBeads.Rd
@@ -56,3 +56,42 @@ sets (\code{RnBiseqSet}) are filtered by coverage, array sets
 (\code{RnBeadSet}) by detection p-value. The number of sites retained
 per sample is reported through \pkg{logger}.
 }
+\examples{
+# A minimal end-to-end run on the BATF example data bundled with the
+# package. The bundled calls are wrapped in an RnBiseqSet so that the
+# example exercises the same code path as a preprocessed RnBeads object.
+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"))
+
+if (requireNamespace("RnBeads", quietly = TRUE) &&
+    requireNamespace("RnBeads.hg38", quietly = TRUE)) {
+    sites <- data.frame(
+        chr = as.character(GenomicRanges::seqnames(msites)),
+        start = GenomicRanges::start(msites),
+        strand = "*",
+        stringsAsFactors = FALSE
+    )
+    rnb_set <- RnBeads::RnBiseqSet(
+        pheno = data.frame(sampleName = "sample_1"),
+        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)
+}
+}
+\seealso{
+\code{\link{run_methyltfr}} for running methylTFR from
+per-sample BED files.
+}

From 818940a4c63fab232a300d0ae96a5588152d02c0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=E2=80=9Cigunduz=E2=80=9D?= <irembgunduz@gmail.com>
Date: Sun, 30 Aug 2026 16:05:14 +0200
Subject: [PATCH 4/4] version bump for rnbeads wrapper update

---
 DESCRIPTION | 2 +-
 NEWS.md     | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 0e491d1..872da72 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: methylTFR
 Title: Quantification of DNA Methylation Signatures in TFBS
-Version: 0.99.6
+Version: 0.99.7
 Date: 2026-08-19
 Authors@R: c(
     person("Irem B.", "Gündüz", , "irembgunduz@gmail.com", role = c("aut", "cre"),
diff --git a/NEWS.md b/NEWS.md
index b6fa53f..f9b6f94 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,4 +1,4 @@
-# methylTFR 0.99.6
+# methylTFR 0.99.7
 
 NEW FEATURES