@@ -69,10 +69,12 @@ setValidity2("BSseq", function(object) {

 # Exported functions -----------------------------------------------------------

+# TODO: Document that if duplicate loci are detected then the output `M` and

+#       `Cov` will be numeric even if input `M` and `Cov` were integer.

 # TODO: BSseq() is arguably a bad constructor. It doesn't return a valid BSseq

 #       object when called without any arguments. It also does some pretty

-#       complicated parsing of the inputs. Still, I think we're stuck with it

-#       because it's been around for a long time.

+#       complicated parsing of the inputs. But we're stuck with it because it's

+#       been around for a long time.

 BSseq <- function(M = NULL, Cov = NULL, coef = NULL, se.coef = NULL,

                   trans = NULL, parameters = NULL, pData = NULL,

                   gr = NULL, pos = NULL, chr = NULL, sampleNames = NULL,

@@ -171,27 +173,15 @@ BSseq <- function(M = NULL, Cov = NULL, coef = NULL, se.coef = NULL,

         if (!is.null(coef) || !is.null(se.coef)) {

             stop("Cannot collapse when 'coef' or 'se.coef' are present.")

         }

-        # NOTE: reduce() sorts the output

-        # TODO: Clarify above comment.

         unique_gr <- unique(gr)

-        # TODO: type = "equal"?

-        ol <- findOverlaps(unique_gr, gr)

-        gr <- unique(gr)

+        ol <- findOverlaps(unique_gr, gr, type = "equal")

+        gr <- unique_gr

+        # TODO: Can we avoid making `idx`?

         idx <- as.list(ol)

         M <- .collapseMatrixLike(M, idx, MARGIN = 2L)

         Cov <- .collapseMatrixLike(Cov, idx, MARGIN = 2L)

     }

-    # Sort loci

-    if (is.unsorted(gr)) {

-        o <- order(gr)

-        gr <- gr[o]

-        M <- M[o, , drop = FALSE]

-        Cov <- Cov[o, , drop = FALSE]

-        coef <- coef[o, , drop = FALSE]

-        se.coef <- se.coef[o, , drop = FALSE]

-    }

-

     # Construct BSseq object

     assays <- SimpleList(M = M, Cov = Cov, coef = coef, se.coef = se.coef)

     assays <- assays[!S4Vectors:::sapply_isNULL(assays)]

@@ -13,16 +13,42 @@ setGeneric(

 # Internal methods -------------------------------------------------------------

+# TODO: Re-write in C/C++ to avoid allocating `ans` multiple times

 setMethod(".collapseMatrixLike", "matrix", function(x, idx, MARGIN) {

     if (MARGIN == 1L) {

-        # TODO: Need colnames?

-        return(do.call(cbind, lapply(idx, function(j) rowSums2(x, cols = j))))

+        ans <- matrix(

+            # NOTE: rowSums2() always returns numeric

+            data = NA_real_,

+            nrow = nrow(x),

+            ncol = length(idx),

+            dimnames = list(rownames(x), NULL))

+        l <- lengths(idx)

+        stopifnot(min(l) > 0)

+        l1 <- which(l == 1)

+        ans[, l1] <- x[, unlist(idx[l1])]

+        lmore <- which(l != 1)

+        ans[, lmore] <- do.call(cbind, lapply(idx[lmore], function(j) {

+            rowSums2(x, cols = j)

+        }))

     } else if (MARGIN == 2L) {

-        # TODO: Need rownames?

-        do.call(rbind, lapply(idx, function(i) colSums2(x, rows = i)))

+        ans <- matrix(

+            # NOTE: colSums2() always returns numeric

+            data = NA_real_,

+            nrow = length(idx),

+            ncol = ncol(x),

+            dimnames = list(NULL, colnames(x)))

+        l <- lengths(idx)

+        stopifnot(min(l) > 0)

+        l1 <- which(l == 1)

+        ans[l1, ] <- x[unlist(idx[l1]), ]

+        lmore <- which(l != 1)

+        ans[lmore, ] <- do.call(rbind, lapply(idx[lmore], function(i) {

+            colSums2(x, rows = i)

+        }))

     } else {

         stop("'MARGIN' must be 1 or 2")

     }

+    ans

 })

 setMethod(

@@ -3,6 +3,13 @@

 \encoding{UTF-8}

 \section{Version 1.17.x}{

   \itemize{

+    \item{\code{BSseq()} will no longer reorder inputs. Previously, the

+    returned \emph{BSseq} object was ordered by ordering the loci, although

+    this behaviour was not documented. \code{BSseq()} may still filter out loci

+    if \code{rmZeroCov = FALSE} or collapse loci if

+    \code{strandCollapse = FALSE} or duplicate loci are detected, but the

+    relative order of loci in the output will match that of the input.

+    }

     \item{Fix bug with \code{maxGap} argument of \code{BSmooth()}. The bug

     meant that the 'maximum gap between two methylation loci' was incorrectly

     set to \code{2 * maxGap + 1} instead of \code{maxGap}. This likely did not

@@ -47,7 +47,7 @@ BSseq(M = NULL, Cov = NULL, coef = NULL, se.coef = NULL,

   In case one or more methylation loci appears multiple times, the

   \code{M} and \code{Cov} matrices are summed over rows linked to the

-  same methylation loci.  See the example below.

+  same methylation loci. See the example below.

   Users should never have to specify \code{coef}, \code{se.coef},

   \code{trans}, and \code{parameters}, this is for internal use (they

@@ -29,14 +29,18 @@ test_that("BSseq", {

     BStest2 <- BSseq(pos = 3:1, chr = rep("chr1", 3), M = M[3:1, ],

                      Cov = M[3:1, ] + 2)

-    expect_equal(BStest, BStest2)

+    # NOTE: BSseq() does not reorder inputs

+    expect_true(!isTRUE(all.equal(BStest, BStest2)))

+    expect_equal(BStest, sort(BStest2))

     M2 <- rbind(M[3:1, ], c(1, 1, 1))

     Cov2 <- rbind(M[3:1, ] + 2, c(1, 1, 1))

     M2[3, ] <- M2[3, ] - 1

     Cov2[3, ] <- Cov2[3, ] - 1

     BStest3 <- suppressWarnings(

         BSseq(pos = c(3:1,1), chr = rep("chr1", 4), M = M2, Cov = Cov2))

-    expect_equal(BStest, BStest3)

+    # NOTE: BSseq() does not reorder inputs

+    expect_true(!isTRUE(all.equal(BStest, BStest3)))

+    expect_equal(BStest, sort(BStest3))

 })

 test_that("BSseq with HDF5Backend", {