From e89a6fd854bb73c4e588333f629f9ecabd5a296f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9Cigunduz=E2=80=9D?= Date: Thu, 27 Aug 2026 14:35:44 +0200 Subject: [PATCH 01/11] version bump --- DESCRIPTION | 2 +- NEWS.md | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 3cdb343..db4344d 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.5 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 8091345..953597f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,4 @@ -# methylTFR 0.99.4 +# methylTFR 0.99.5 NEW FEATURES From 685dde3167748f8016f9854ee675e75d29571bac Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9Cigunduz=E2=80=9D?= Date: Thu, 27 Aug 2026 15:00:28 +0200 Subject: [PATCH 02/11] styling --- R/compute_deviations.R | 13 +++++++++---- R/memory_helpers.R | 6 +++++- R/variability.R | 5 ++++- vignettes/memTcells.Rmd | 6 ++++-- vignettes/methylTFR.Rmd | 9 ++++++--- 5 files changed, 28 insertions(+), 11 deletions(-) diff --git a/R/compute_deviations.R b/R/compute_deviations.R index a659636..2458a6b 100644 --- a/R/compute_deviations.R +++ b/R/compute_deviations.R @@ -55,8 +55,10 @@ #' ) #' @export computeDeviation <- function( - motif, msites, tf_bindsites, gcfreqs, - enhancer = NULL, ignoreStrand = TRUE, + motif, msites, + tf_bindsites, gcfreqs, + enhancer = NULL, + ignoreStrand = TRUE, binMsites ) { if (!is.logical(ignoreStrand)) { @@ -73,13 +75,16 @@ computeDeviation <- 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") + stop("Please provide a valid tf binding sites + as GRangesList") } if (!is.null(enhancer) && !is(enhancer, "GRanges")) { stop("Please provide a valid enhancer regions") } 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/memory_helpers.R b/R/memory_helpers.R index 007b43d..b9473bb 100644 --- a/R/memory_helpers.R +++ b/R/memory_helpers.R @@ -19,7 +19,11 @@ create_sink <- function( 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/variability.R b/R/variability.R index dc8a531..43b84c8 100644 --- a/R/variability.R +++ b/R/variability.R @@ -120,7 +120,10 @@ calibrateDeviations <- function(devs, method = c("robust", "gaussian")) { #' @export computeZScoreVariability <- function( object, - method = c("robust", "gaussian"), + method = c( + "robust", + "gaussian" + ), bootstrap = FALSE, niterations = 1000L, conf_level = 0.95, diff --git a/vignettes/memTcells.Rmd b/vignettes/memTcells.Rmd index a6619a5..1e13f2d 100644 --- a/vignettes/memTcells.Rmd +++ b/vignettes/memTcells.Rmd @@ -1,11 +1,13 @@ --- -title: "Case study: TF activity in memory vs. naive T cells" +title: "Case study: TF activity +in memory vs. naive T cells" author: "Irem B. Gündüz, Sarath Kumar Murugan, Fabian Muller" date: "`r Sys.Date()`" package: methylTFR output: BiocStyle::html_document vignette: > - %\VignetteIndexEntry{Case study: TF activity in memory vs. naive T cells} + %\VignetteIndexEntry{Case study: + TF activity in memory vs. naive T cells} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- diff --git a/vignettes/methylTFR.Rmd b/vignettes/methylTFR.Rmd index 0b37c6c..1939a8e 100644 --- a/vignettes/methylTFR.Rmd +++ b/vignettes/methylTFR.Rmd @@ -1,11 +1,13 @@ --- -title: "methylTFR: DNA Methylation Signatures in Transcription Factor Binding Sites" +title: "methylTFR: DNA Methylation Signatures +in TFBS" author: "Irem B. Gündüz, Sarath Kumar Murugan, Fabian Muller" date: "`r Sys.Date()`" package: methylTFR output: BiocStyle::html_document vignette: > - %\VignetteIndexEntry{methylTFR: DNA Methylation Signatures in TFBS} + %\VignetteIndexEntry{methylTFR: DNA Methylation + Signatures in TFBS} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- @@ -104,7 +106,8 @@ the multi-sample section below. To build annotations for a different assembly, motif set, or region restriction, use -[`methylTFRAnnotationBuilder`](https://github.com/EpigenomeInformatics/methylTFRAnnotationBuilder). +[`methylTFRAnnotationBuilder`]( + https://github.com/EpigenomeInformatics/methylTFRAnnotationBuilder). The examples in this section instead use a small `BATF` subset bundled with `methylTFR`, so they run without any annotation package. From 60b6ede8774e918a504dce848e8fdb3a78a72459 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9Cigunduz=E2=80=9D?= Date: Thu, 27 Aug 2026 15:28:14 +0200 Subject: [PATCH 03/11] style --- R/compute_deviations.R | 2 +- R/expected_deviations.R | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/R/compute_deviations.R b/R/compute_deviations.R index 2458a6b..0f8b322 100644 --- a/R/compute_deviations.R +++ b/R/compute_deviations.R @@ -76,7 +76,7 @@ computeDeviation <- function( !any(c(!is(tf_bindsites, "GRangesList") || !is.list(tf_bindsites)))) { stop("Please provide a valid tf binding sites - as GRangesList") + as GRangesList") } if (!is.null(enhancer) && !is(enhancer, "GRanges")) { stop("Please provide a valid enhancer regions") diff --git a/R/expected_deviations.R b/R/expected_deviations.R index 48731b5..1855ab1 100644 --- a/R/expected_deviations.R +++ b/R/expected_deviations.R @@ -37,7 +37,7 @@ addGCBintoMethylome <- function( ) { if (!is.logical(ignoreStrand)) { warning("Found invalid strand option, - using the default") + using the default") ignoreStrand <- TRUE } if (is.null(msites) || !is(msites, "GRanges")) { @@ -86,7 +86,7 @@ computeExpectations <- function(binMsites, gcfreq) { } if (!is.matrix(gcfreq)) { stop("Please provide a valid - GC bin frequency table as a matrix") + GC bin frequency table as a matrix") } exp.data <- t(gcfreq) %*% binMsites[, 2] mpos <- round(seq(-floor(length(exp.data) / 2), From 891961658924e2780d9a82538dba9fa144e1bbdc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9Cigunduz=E2=80=9D?= Date: Fri, 28 Aug 2026 13:27:34 +0200 Subject: [PATCH 04/11] Revert "style" This reverts commit 60b6ede8774e918a504dce848e8fdb3a78a72459. --- R/compute_deviations.R | 2 +- R/expected_deviations.R | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/R/compute_deviations.R b/R/compute_deviations.R index 0f8b322..2458a6b 100644 --- a/R/compute_deviations.R +++ b/R/compute_deviations.R @@ -76,7 +76,7 @@ computeDeviation <- function( !any(c(!is(tf_bindsites, "GRangesList") || !is.list(tf_bindsites)))) { stop("Please provide a valid tf binding sites - as GRangesList") + as GRangesList") } if (!is.null(enhancer) && !is(enhancer, "GRanges")) { stop("Please provide a valid enhancer regions") diff --git a/R/expected_deviations.R b/R/expected_deviations.R index 1855ab1..48731b5 100644 --- a/R/expected_deviations.R +++ b/R/expected_deviations.R @@ -37,7 +37,7 @@ addGCBintoMethylome <- function( ) { if (!is.logical(ignoreStrand)) { warning("Found invalid strand option, - using the default") + using the default") ignoreStrand <- TRUE } if (is.null(msites) || !is(msites, "GRanges")) { @@ -86,7 +86,7 @@ computeExpectations <- function(binMsites, gcfreq) { } if (!is.matrix(gcfreq)) { stop("Please provide a valid - GC bin frequency table as a matrix") + GC bin frequency table as a matrix") } exp.data <- t(gcfreq) %*% binMsites[, 2] mpos <- round(seq(-floor(length(exp.data) / 2), From 028e789e81a828c172223a0b4469a48184b5a45e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9Cigunduz=E2=80=9D?= Date: Fri, 28 Aug 2026 14:16:37 +0200 Subject: [PATCH 05/11] attempt to fix bioc notes --- .gitignore | 1 + R/compute_deviations.R | 19 ++++----- R/differential_analysis.R | 12 +++--- R/expected_deviations.R | 24 ++++++----- R/memory_helpers.R | 4 +- R/methyltfr_core.R | 28 ++++++------- R/plot_helpers.R | 10 ++--- R/plots.R | 16 ++++---- R/rnbeads_interface.R | 14 +++---- R/run_methyltfr.R | 55 +++++++++++++++++++++++--- R/variability.R | 18 ++++----- man/run_methyltfr.Rd | 44 +++++++++++++++++++++ tests/testthat/test_computeDeviation.R | 4 +- 13 files changed, 168 insertions(+), 81 deletions(-) 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/R/compute_deviations.R b/R/compute_deviations.R index 2458a6b..47e9760 100644 --- a/R/compute_deviations.R +++ b/R/compute_deviations.R @@ -55,11 +55,11 @@ #' ) #' @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") @@ -69,14 +69,15 @@ computeDeviation <- function( 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") + 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") + stop("Please provide a valid tf binding sites as GRangesList") } if (!is.null(enhancer) && !is(enhancer, "GRanges")) { stop("Please provide a valid enhancer regions") diff --git a/R/differential_analysis.R b/R/differential_analysis.R index 92ebab9..ba39ef7 100644 --- a/R/differential_analysis.R +++ b/R/differential_analysis.R @@ -53,12 +53,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"))) { 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 b9473bb..602f680 100644 --- a/R/memory_helpers.R +++ b/R/memory_helpers.R @@ -12,8 +12,8 @@ #' @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)) { diff --git a/R/methyltfr_core.R b/R/methyltfr_core.R index a551414..02b6827 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") @@ -147,16 +147,16 @@ 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.") 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..3941e9d 100644 --- a/R/plots.R +++ b/R/plots.R @@ -58,9 +58,9 @@ #' #' @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( @@ -211,11 +211,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")) { diff --git a/R/rnbeads_interface.R b/R/rnbeads_interface.R index ee8df49..59051bc 100644 --- a/R/rnbeads_interface.R +++ b/R/rnbeads_interface.R @@ -81,11 +81,11 @@ #' } #' @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( @@ -268,8 +268,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..fa16671 100644 --- a/R/run_methyltfr.R +++ b/R/run_methyltfr.R @@ -47,14 +47,57 @@ #' 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", diff --git a/R/variability.R b/R/variability.R index 43b84c8..69e6543 100644 --- a/R/variability.R +++ b/R/variability.R @@ -119,15 +119,15 @@ 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")) { diff --git a/man/run_methyltfr.Rd b/man/run_methyltfr.Rd index 8ef2393..d1b82e2 100644 --- a/man/run_methyltfr.Rd +++ b/man/run_methyltfr.Rd @@ -76,6 +76,50 @@ calculate the deviation in transcription factor footprint base for all given motifs per raw samples } +\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) +} \seealso{ \code{\link{run_methylTFR_RnBeads}} for running methylTFR directly on a preprocessed RnBeads object. diff --git a/tests/testthat/test_computeDeviation.R b/tests/testthat/test_computeDeviation.R index d76976b..e749273 100644 --- a/tests/testthat/test_computeDeviation.R +++ b/tests/testthat/test_computeDeviation.R @@ -26,7 +26,7 @@ test_that("computeDeviation", { ) expect_error( computeDeviation(motif = "BATF", NULL, tf_bindsites, gcfreqs, enhancer = NULL, ignoreStrand = TRUE, bin_meth), - "Please provide a valid methylation\n sites with read_methylome function" + "Please provide a valid methylation sites with read_methylome function" ) expect_error( computeDeviation(motif = "BATF", msites, NULL, gcfreqs, enhancer = NULL, ignoreStrand = TRUE, bin_meth), @@ -34,7 +34,7 @@ test_that("computeDeviation", { ) expect_error( computeDeviation("BATF", msites, tf_bindsites, NULL, enhancer = NULL, ignoreStrand = TRUE, bin_meth), - "Please provide a valid\n GC bin frequency table as a matrix" + "Please provide a valid GC bin frequency table as a matrix" ) # Check that the function returns a warning when given incorrect options From 94637f0857dcdf2061bd57d00c45bbb28e8f792d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9Cigunduz=E2=80=9D?= Date: Fri, 28 Aug 2026 14:32:28 +0200 Subject: [PATCH 06/11] update pkgdown --- R/methyltfr_core.R | 263 +++++++++++------- docs/404.html | 2 +- docs/LICENSE-text.html | 2 +- docs/articles/index.html | 4 +- docs/articles/memTcells.html | 82 +++--- docs/articles/methylTFR.html | 142 +++------- docs/authors.html | 6 +- docs/index.html | 2 +- docs/news/index.html | 4 +- docs/pkgdown.yml | 2 +- docs/reference/addGCBintoMethylome.html | 2 +- docs/reference/anova_helper.html | 2 +- docs/reference/calibrateDeviations.html | 2 +- .../cbind-methylTFRdeviations-method.html | 2 +- docs/reference/check_annotation_inputs.html | 2 +- docs/reference/check_run_options.html | 2 +- docs/reference/cleanMem.html | 2 +- docs/reference/computeColZScore.html | 2 +- docs/reference/computeDeviation.html | 2 +- docs/reference/computeExpectations.html | 2 +- docs/reference/computeExpectedFootprint.html | 2 +- docs/reference/computeFootprint.html | 2 +- docs/reference/computeRowZScore.html | 2 +- docs/reference/computeZScoreVariability.html | 4 +- docs/reference/create_sink.html | 2 +- docs/reference/dev_helper.html | 2 +- ...ionZScores-methylTFRdeviations-method.html | 2 +- docs/reference/deviationZScores.html | 2 +- ...deviations-methylTFRdeviations-method.html | 2 +- docs/reference/deviations.html | 2 +- .../differential_deviation_test.html | 2 +- docs/reference/granges_helper.html | 2 +- docs/reference/index.html | 2 +- docs/reference/kw_helper.html | 2 +- docs/reference/methylTFRdeviations-class.html | 2 +- docs/reference/methyltfr_core.html | 2 +- docs/reference/plotExpectedFootprint.html | 2 +- docs/reference/plotMotifFootprint.html | 2 +- .../rbind-methylTFRdeviations-method.html | 2 +- docs/reference/read_methylome.html | 2 +- docs/reference/read_sample_annotation.html | 2 +- docs/reference/rnb_column.html | 2 +- docs/reference/rnb_has_coverage.html | 2 +- docs/reference/rnb_sample_ids.html | 2 +- docs/reference/rnb_sample_msites.html | 2 +- docs/reference/rnb_sites_to_granges.html | 2 +- docs/reference/run_methylTFR_RnBeads.html | 2 +- docs/reference/run_methyltfr.html | 56 +++- docs/reference/set_grid.html | 2 +- docs/reference/t_helper.html | 2 +- docs/reference/wilcoxon_helper.html | 2 +- docs/reference/write_block_to_sink.html | 2 +- docs/search.json | 2 +- vignettes/methylTFR.Rmd | 2 +- 54 files changed, 345 insertions(+), 308 deletions(-) diff --git a/R/methyltfr_core.R b/R/methyltfr_core.R index 02b6827..1322e97 100644 --- a/R/methyltfr_core.R +++ b/R/methyltfr_core.R @@ -107,6 +107,157 @@ 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 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 threads Thread count for parallel processing. +#' @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 parallel mclapply +#' @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, threads, 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 <- 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, 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 @@ -147,48 +298,11 @@ 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) # Split the motifs into chunks numChunks <- ceiling(length(motifs) / chunkSize) @@ -213,67 +327,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, + threads = threads, 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/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