@@ -1,7 +1,7 @@

 Package: NanoMethViz

 Type: Package

 Title: Visualise methlation data from Oxford Nanopore sequencing

-Version: 1.99.3

+Version: 1.99.4

 Authors@R: c(

     person("Shian", "Su", email = "[email protected]", role = c("cre", "aut")))

 Description: NanoMethViz is a toolkit for visualising methylation data from 

@@ -36,7 +36,7 @@ Imports:

     GenomicRanges,

     ggthemes,

     glue,

-    limma,

+    limma (>= 3.44.0),

     patchwork,

     purrr,

     rlang,

@@ -3,6 +3,9 @@

 #' @param x the NanoMethResult object.

 #' @param regions a table of regions containing at least columns chr, strand,

 #'   start and end. Any additiona columns can be used for grouping.

+#' @param binary_threshold the modification probability such that calls with

+#'   modification probability above the threshold are considered methylated, and

+#'   those with probability equal or below are considered unmethylated.

 #' @param group_col the column to group aggregated trends by. This column can

 #'   be in from the regions table or samples(x).

 #' @param flank the number of flanking bases to add to each side of each region.

@@ -24,6 +27,7 @@

 plot_agg_regions <- function(

     x,

     regions,

+    binary_threshold = 0.5,

     group_col = NULL,

     flank = 2000,

     stranded = TRUE,

@@ -41,7 +45,7 @@ plot_agg_regions <- function(

     # query methylation data

     methy_data <- .get_methy_data(x, regions, flank)

     methy_data <- .scale_methy_data(methy_data, stranded, flank)

-    methy_data <- .pos_avg(methy_data)

+    methy_data <- .pos_avg(methy_data, binary_threshold = binary_threshold)

     # unnest and annotate after position averaging to reduce data burden

     methy_data <- .unnest_anno(methy_data, NanoMethViz::samples(x))

     methy_data <- .bin_avg(methy_data, group_col = group_col)

@@ -132,12 +136,12 @@ plot_agg_regions <- function(

     methy_data

 }

-.pos_avg <- function(x) {

+.pos_avg <- function(x, binary_threshold = 0.5) {

     average_by_pos <- function(x) {

           x %>%

               dplyr::group_by(.data$sample, .data$chr, .data$strand, .data$rel_pos) %>%

               dplyr::summarise(

-                  methy_prop = mean(.data$statistic > 0),

+                  methy_prop = mean(e1071::sigmoid(.data$statistic) > binary_threshold),

                   .groups = "drop")

       }

@@ -22,6 +22,9 @@ setGeneric("plot_gene", function(x, gene, ...) {

 #'   values for left and right window size. Values indicate proportion of gene

 #'   length.

 #' @param anno_regions the data.frame of regions to annotate.

+#' @param binary_threshold the modification probability such that calls with

+#'   modification probability above the threshold are set to 1 and probabilities

+#'   equalt to or below the threshold are set to 0.

 #' @param spaghetti whether or not individual reads should be shown.

 #' @param span the span for loess smoothing.

 #' @param gene_anno whether or not gene annotation tracks are plotted.

@@ -40,6 +43,7 @@ setMethod("plot_gene", signature(x = "NanoMethResult", gene = "character"),

         gene,

         window_prop = 0.3,

         anno_regions = NULL,

+        binary_threshold = NULL,

         spaghetti = FALSE,

         span = NULL,

         gene_anno = TRUE

@@ -49,6 +53,7 @@ setMethod("plot_gene", signature(x = "NanoMethResult", gene = "character"),

             gene,

             window_prop = window_prop,

             anno_regions = anno_regions,

+            binary_threshold = binary_threshold,

             spaghetti = spaghetti,

             span = span,

             gene_anno = gene_anno

@@ -61,6 +66,7 @@ setMethod("plot_gene", signature(x = "NanoMethResult", gene = "character"),

     gene,

     window_prop,

     anno_regions,

+    binary_threshold,

     spaghetti,

     span,

     gene_anno

@@ -91,6 +97,7 @@ setMethod("plot_gene", signature(x = "NanoMethResult", gene = "character"),

             window_prop = window_prop,

             sample_anno = sample_anno,

             anno_regions = anno_regions,

+            binary_threshold = binary_threshold,

             spaghetti = spaghetti,

             span = span

         )

@@ -14,6 +14,7 @@ plot_grange <- function(

     x,

     grange,

     anno_regions = NULL,

+    binary_threshold = NULL,

     spaghetti = FALSE,

     span = NULL,

     window_prop = 0.3

@@ -34,6 +35,7 @@ plot_grange <- function(

         start = start,

         end = end,

         anno_regions = anno_regions,

+        binary_threshold = binary_threshold,

         spaghetti = spaghetti,

         span = span,

         window_prop = window_prop

@@ -38,16 +38,25 @@ plot_mds <- function(x, top = 500, plot_dims = c(1, 2), labels = colnames(x), gr

     var_exp1 <- round(100 * mds_res$var.explained[plot_dims[1]])

     var_exp2 <- round(100 * mds_res$var.explained[plot_dims[2]])

-    plot_data <- data.frame(

-        dim1 = mds_res$eigen.vectors[, plot_dims[1]],

-        dim2 = mds_res$eigen.vectors[, plot_dims[2]]

-    )

+    if ("eigen.vectors" %in% names(mds_res)) {

+        plot_data <- data.frame(

+            dim1 = mds_res$eigen.vectors[, plot_dims[1]],

+            dim2 = mds_res$eigen.vectors[, plot_dims[2]]

+        )

+        xlabel <- glue::glue("Leading logFC Dim {plot_dims[1]} ({var_exp1}%)")

+        ylabel <- glue::glue("Leading logFC Dim {plot_dims[2]} ({var_exp2}%)")

+    } else {

+        plot_data <- data.frame(

+            dim1 = mds_res$cmdscale.out[, plot_dims[1]],

+            dim2 = mds_res$cmdscale.out[, plot_dims[2]]

+        )

+        xlabel <- glue::glue("Leading logFC Dim {plot_dims[1]}")

+        ylabel <- glue::glue("Leading logFC Dim {plot_dims[2]}")

+    }

     plot_data$labels <- labels

     plot_data$groups <- groups

-    xlabel <- glue::glue("Leading logFC Dim {plot_dims[1]} ({var_exp1}%)")

-    ylabel <- glue::glue("Leading logFC Dim {plot_dims[2]} ({var_exp2}%)")

     if (!is.null(groups)) {

         p <- ggplot2::ggplot(plot_data, ggplot2::aes_string(x = "dim1", y = "dim2", col = "groups"))

@@ -8,6 +8,7 @@ plot_methylation_internal <- function(

     title,

     palette_col = ggplot2::scale_colour_brewer(palette = "Set1"),

     anno_regions = NULL,

+    binary_threshold = NULL,

     spaghetti = FALSE,

     span = NULL

 ) {

@@ -26,11 +27,21 @@ plot_methylation_internal <- function(

     }

     # extract group information and convert probabilities

-    plot_data <- methy_data %>%

-        dplyr::inner_join(sample_anno, by = "sample") %>%

-        dplyr::mutate(

-            mod_prob = e1071::sigmoid(.data$statistic)

-        )

+    if (is.null(binary_threshold)) {

+        plot_data <- methy_data %>%

+            dplyr::inner_join(sample_anno, by = "sample") %>%

+            dplyr::mutate(

+                mod_prob = e1071::sigmoid(.data$statistic)

+            )

+    } else {

+        plot_data <- methy_data %>%

+            dplyr::inner_join(sample_anno, by = "sample") %>%

+            dplyr::mutate(

+                mod_prob = as.numeric(

+                    e1071::sigmoid(.data$statistic) > binary_threshold

+                )

+            )

+    }

     # set up plot

     p <- ggplot(plot_data, aes(x = .data$pos, col = .data$group))

@@ -99,7 +110,8 @@ plot_feature <- function(

     sample_anno,

     anno_regions = NULL,

     window_prop = c(0.3, 0.3),

-    spaghetti = TRUE,

+    binary_threshold = NULL,

+    spaghetti = FALSE,

     span = NULL

 ) {

     chr <- feature$chr

@@ -135,6 +147,7 @@ plot_feature <- function(

         xlim = xlim,

         title = title,

         anno_regions = anno_regions,

+        binary_threshold = binary_threshold,

         spaghetti = spaghetti,

         sample_anno = sample_anno,

         span = span

@@ -21,6 +21,9 @@ setGeneric("plot_region", function(x, chr, start, end, ...) {

 #' @rdname plot_region

 #'

 #' @param anno_regions the data.frame of regions to be annotated.

+#' @param binary_threshold the modification probability such that calls with

+#'   modification probability above the threshold are set to 1 and probabilities

+#'   equalt to or below the threshold are set to 0.

 #' @param spaghetti whether or not individual reads should be shown.

 #' @param span the span for loess smoothing.

 #' @param window_prop the size of flanking region to plot. Can be a vector of two

@@ -47,6 +50,7 @@ setMethod("plot_region",

         start,

         end,

         anno_regions = NULL,

+        binary_threshold = NULL,

         spaghetti = FALSE,

         span = NULL,

         window_prop = 0

@@ -57,6 +61,7 @@ setMethod("plot_region",

             start = start,

             end = end,

             anno_regions = anno_regions,

+            binary_threshold = binary_threshold,

             spaghetti = spaghetti,

             span = span,

             window_prop = window_prop

@@ -80,6 +85,7 @@ setMethod("plot_region",

         start,

         end,

         anno_regions = NULL,

+        binary_threshold = NULL,

         spaghetti = FALSE,

         span = NULL,

         window_prop = 0

@@ -91,6 +97,7 @@ setMethod("plot_region",

             start = start,

             end = end,

             anno_regions = anno_regions,

+            binary_threshold = binary_threshold,

             spaghetti = spaghetti,

             span = span,

             window_prop = window_prop

@@ -105,6 +112,7 @@ setMethod("plot_region",

     start,

     end,

     anno_regions = NULL,

+    binary_threshold = NULL,

     spaghetti = FALSE,

     span = NULL,

     window_prop = 0

@@ -150,6 +158,7 @@ setMethod("plot_region",

         xlim = xlim,

         title = title,

         anno_regions = anno_regions,

+        binary_threshold = binary_threshold,

         spaghetti = spaghetti,

         sample_anno = sample_anno,

         span = span

@@ -7,6 +7,7 @@

 plot_agg_regions(

   x,

   regions,

+  binary_threshold = 0.5,

   group_col = NULL,

   flank = 2000,

   stranded = TRUE,

@@ -20,6 +21,10 @@ plot_agg_regions(

 \item{regions}{a table of regions containing at least columns chr, strand,

 start and end. Any additiona columns can be used for grouping.}

+\item{binary_threshold}{the modification probability such that calls with

+modification probability above the threshold are considered methylated, and

+those with probability equal or below are considered unmethylated.}

+

 \item{group_col}{the column to group aggregated trends by. This column can

 be in from the regions table or samples(x).}

@@ -12,6 +12,7 @@ plot_gene(x, gene, ...)

   gene,

   window_prop = 0.3,

   anno_regions = NULL,

+  binary_threshold = NULL,

   spaghetti = FALSE,

   span = NULL,

   gene_anno = TRUE

@@ -30,6 +31,10 @@ length.}

 \item{anno_regions}{the data.frame of regions to annotate.}

+\item{binary_threshold}{the modification probability such that calls with

+modification probability above the threshold are set to 1 and probabilities

+equalt to or below the threshold are set to 0.}

+

 \item{spaghetti}{whether or not individual reads should be shown.}

 \item{span}{the span for loess smoothing.}

@@ -8,6 +8,7 @@ plot_grange(

   x,

   grange,

   anno_regions = NULL,

+  binary_threshold = NULL,

   spaghetti = FALSE,

   span = NULL,

   window_prop = 0.3

@@ -20,6 +21,10 @@ plot_grange(

 \item{anno_regions}{the data.frame of regions to be annotated.}

+\item{binary_threshold}{the modification probability such that calls with

+modification probability above the threshold are set to 1 and probabilities

+equalt to or below the threshold are set to 0.}

+

 \item{spaghetti}{whether or not individual reads should be shown.}

 \item{span}{the span for loess smoothing.}

@@ -14,6 +14,7 @@ plot_region(x, chr, start, end, ...)

   start,

   end,

   anno_regions = NULL,

+  binary_threshold = NULL,

   spaghetti = FALSE,

   span = NULL,

   window_prop = 0

@@ -25,6 +26,7 @@ plot_region(x, chr, start, end, ...)

   start,

   end,

   anno_regions = NULL,

+  binary_threshold = NULL,

   spaghetti = FALSE,

   span = NULL,

   window_prop = 0

@@ -43,6 +45,10 @@ plot_region(x, chr, start, end, ...)

 \item{anno_regions}{the data.frame of regions to be annotated.}

+\item{binary_threshold}{the modification probability such that calls with

+modification probability above the threshold are set to 1 and probabilities

+equalt to or below the threshold are set to 0.}

+

 \item{spaghetti}{whether or not individual reads should be shown.}

 \item{span}{the span for loess smoothing.}

@@ -5,6 +5,11 @@

 using namespace Rcpp;

+#ifdef RCPP_USE_GLOBAL_ROSTREAM

+Rcpp::Rostream<true>&  Rcpp::Rcout = Rcpp::Rcpp_cout_get();

+Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();

+#endif

+

 // convert_methy_to_dss_cpp

 std::vector<std::string> convert_methy_to_dss_cpp(std::string input, std::string output_dir);

 RcppExport SEXP _NanoMethViz_convert_methy_to_dss_cpp(SEXP inputSEXP, SEXP output_dirSEXP) {