Merge branch 'clustering_refactoring2' into devel

Author: Charles Plessy

DESCRIPTION

@@ -15,8 +15,6 @@ Imports:
     BSgenome,
     CAGEfightR,
     data.table,
-    DelayedArray,
-    DelayedMatrixStats,
     formula.tools,
     GenomeInfoDb,
     GenomicAlignments,
@@ -81,8 +79,8 @@ Collate:
     'AggregationMethods.R'
     'CAGEfightR.R'
     'CAGEr-package.R'
+    'Paraclu.R'
     'CorrelationMethods.R'
-    'CumulativeDistributionFunctions.R'
     'GetMethods.R'
     'CumulativeDistributionMethods.R'
     'ExportMethods.R'

NAMESPACE

@@ -19,7 +19,6 @@ export(CTSScumulativesCC)
 export(CTSScumulativesTagClusters)
 export(CTSSnormalizedTpmDF)
 export(CTSSnormalizedTpmGR)
-export(CTSStagCountDA)
 export(CTSStagCountDF)
 export(CTSStagCountGR)
 export(CTSStagCountSE)
@@ -37,6 +36,7 @@ export(consensusClustersGR)
 export(consensusClustersSE)
 export(consensusClustersTpm)
 export(cumulativeCTSSdistribution)
+export(distclu)
 export(exportToTrack)
 export(expressionClasses)
 export(findStrandInvaders)
@@ -54,6 +54,7 @@ export(librarySizes)
 export(mergeCAGEsets)
 export(mergeSamples)
 export(normalizeTagCount)
+export(paraclu)
 export(plotAnnot)
 export(plotCorrelation)
 export(plotCorrelation2)
@@ -79,8 +80,6 @@ exportMethods(flagByUpstreamSequences)
 exportMethods(quickEnhancers)
 exportMethods(removeStrandInvaders)
 import(BiocGenerics)
-import(DelayedArray)
-import(DelayedMatrixStats)
 import(MultiAssayExperiment)
 import(SummarizedExperiment)
 import(methods)
@@ -112,35 +111,34 @@ importFrom(GenomicRanges,findOverlaps)
 importFrom(GenomicRanges,from_GPos_to_GRanges)
 importFrom(GenomicRanges,granges)
 importFrom(GenomicRanges,promoters)
-importFrom(IRanges,IPos)
 importFrom(IRanges,IRanges)
 importFrom(IRanges,IRangesList)
 importFrom(IRanges,RleList)
-importFrom(IRanges,Views)
 importFrom(IRanges,extractList)
 importFrom(IRanges,reduce)
-importFrom(IRanges,viewApply)
 importFrom(KernSmooth,bkde2D)
 importFrom(Rsamtools,"bamFlag<-")
 importFrom(Rsamtools,ScanBamParam)
 importFrom(Rsamtools,scanBamFlag)
 importFrom(S4Vectors,"%in%")
 importFrom(S4Vectors,DataFrame)
 importFrom(S4Vectors,List)
+importFrom(S4Vectors,Pairs)
 importFrom(S4Vectors,Rle)
 importFrom(S4Vectors,SimpleList)
 importFrom(S4Vectors,decode)
 importFrom(S4Vectors,endoapply)
+importFrom(S4Vectors,first)
 importFrom(S4Vectors,mcols)
 importFrom(S4Vectors,queryHits)
 importFrom(S4Vectors,runLength)
 importFrom(S4Vectors,runValue)
+importFrom(S4Vectors,second)
 importFrom(S4Vectors,subjectHits)
 importFrom(S4Vectors,unstrsplit)
 importFrom(SummarizedExperiment,SummarizedExperiment)
 importFrom(SummarizedExperiment,assay)
 importFrom(SummarizedExperiment,assays)
-importFrom(SummarizedExperiment,rowRanges)
 importFrom(VGAM,zeta)
 importFrom(data.table,data.table)
 importFrom(data.table,fread)

NEWS.md

@@ -10,6 +10,9 @@ NEW FEATURES
     `legend`, `xlab`, `ylab`, `xlim` and `ylim` arguments were removed as this
     can be controlled via `ggplot2` functions.
 -   New `TSSlogo` function wrapping the `ggseqlogo` package.
+-   New `distclu` and `paraclu` functions that can run directly on `CTSS`
+    objects.  You can use them to test parameters before running the whole
+    `CAGEexp` object through `clusterCTSS`.
 
 # Changes in version 2.11.1
 
@@ -22,6 +25,16 @@ BUG FIXES
 
 # Changes in version 2.8.0 (25/10/2023)
 
+BACKWARDS-INCOMPATIBLE CHANGES
+
+-   The `CTSStagCountDA` function is removed.
+-   The _dominant peak_ in `TagClusters` objects is now a `GRanges` object like
+    in `ConsensusClusters`.
+-   The `custom` method for tag clustering is removed.  It was obsoleted by
+    the newer `CustomConsensusClusters` function.
+-   The `exportToTrack` function now exports scores of _tag clusters_ and
+    _consensus clusters_ instead of setting them to zero.
+
 BUG FIXES
 
 -   Correct quantile positions, which were shifted by one base.  This bug may
@@ -41,7 +54,7 @@ BUG FIXES
 NEW FEATURES
 
 -   Allow URLs to files in `getCTSS()` (Fixes #50).
--   Accelerated the computation of quantile position by ~20 times.
+-   Accelerated the computation of cumulative sums and quantile position ~10×.
 -   New `resetCAGEexp()` function.
 -   New `flagByUpstreamSequences()` function.
 -   The `annotateCTSS` and `annotateConsensusClusters` function gain a

R/AggregationMethods.R

@@ -144,35 +144,14 @@ setMethod( "aggregateTagClusters", "CAGEr"
 
   # CTSS with score that is sum of all samples
   ctss <- CTSScoordinatesGR(CAGEexp_obj)
-  score(ctss) <- rowSums(CTSSnormalizedTpmDF(CAGEexp_obj) |> DelayedArray::DelayedArray() )
+  score(ctss) <- rowSums.RleDataFrame(CTSSnormalizedTpmDF(CAGEexp_obj))
 
   # Stop if some TCs do not overlap with any CTSS, because the rest of the code
   # is not robust against that.
   if (any(countOverlaps(clusters.gr, ctss) == 0))
     stop("Some TCs do not overlap any CTSS!")
 
-  # See `benchmarks/dominant_ctss.md`.
-  o <- findOverlaps(clusters.gr, ctss)
-  
-  rl <- rle(queryHits(o))$length
-  cluster_start_idx <- cumsum(c(1, head(rl, -1))) # Where each run starts
-  grouped_scores <- extractList(score(ctss), o)
-  # grouped_pos    <- extractList(pos(ctss), o)
-  
-  find.dominant.idx <- function (x) {
-    # which.max is breaking ties by taking the last, but this will give slightly
-    # different biases on plus an minus strands.
-    w <- which(x == max(x))
-    w[ceiling(length(w)/2)]
-  }
-  local_max_idx <- sapply(grouped_scores, find.dominant.idx) -1  # Start at zero
-  global_max_ids <- cluster_start_idx + local_max_idx
-  # start(clusters.gr) <- min(grouped_pos)
-  # end  (clusters.gr) <- max(grouped_pos)
-  score(clusters.gr) <- sum(grouped_scores)
-  clusters.gr$dominant_ctss     <- granges(ctss)[subjectHits(o)][global_max_ids]
-  clusters.gr$tpm.dominant_ctss <- score(ctss)  [subjectHits(o)][global_max_ids]
-  clusters.gr
+  clusters.gr <- .ctss_summary_for_clusters(ctss, clusters.gr, removeSingletons = FALSE)
 
   names(clusters.gr) <- as.character(clusters.gr)
   .ConsensusClusters(clusters.gr)
@@ -190,7 +169,7 @@ setMethod( ".CCtoSE"
       rowRanges(se)$cluster <- ranges2names(rowRanges(se), consensus.clusters)
 
     if (tpmThreshold > 0)
-      se <- se[rowSums(DelayedArray(assays(se)[["normalizedTpmMatrix"]])) > tpmThreshold,]
+      se <- se[rowSums.RleDataFrame(assays(se)[["normalizedTpmMatrix"]]) > tpmThreshold,]
 
     .rowsumAsMatrix <- function(DF, names) {
       # First, remove CTSS that do not match clusters

R/CAGEr.R

@@ -186,10 +186,15 @@ setMethod(".filterCtss", "RangedSummarizedExperiment", function (object, thresho
   assay <- ifelse(thresholdIsTpm, "normalizedTpmMatrix", "counts")
   if(assay == "normalizedTpmMatrix" & is.null(assays(object)[[assay]]))
     stop("Normalise the CAGEr object first with ", sQuote("normalizeTagCount()"), ".")
-  .filterCtss(DelayedArray(assays(object)[[assay]]), threshold, nrPassThreshold, thresholdIsTpm)
+  .filterCtss(assays(object)[[assay]], threshold, nrPassThreshold, thresholdIsTpm)
 })
 
-setMethod(".filterCtss", "DelayedArray", function (object, threshold, nrPassThreshold, thresholdIsTpm) {
-  nr.pass.threshold <- rowSums(object >= threshold)
-  Rle(nr.pass.threshold >= min(nrPassThreshold, ncol(object)))
+setMethod(".filterCtss", "DataFrame", function (object, threshold, nrPassThreshold, thresholdIsTpm) {
+  nr.pass.threshold <- rowSums.RleDataFrame(lapply(object, \(x) x > threshold) |> DataFrame())
+  nr.pass.threshold >= min(nrPassThreshold, ncol(object))
 })
+
+setMethod(".filterCtss", "matrix", function (object, threshold, nrPassThreshold, thresholdIsTpm) {
+  nr.pass.threshold <- rowSums(object > threshold)
+  nr.pass.threshold >= min(nrPassThreshold, ncol(object))
+})

R/CTSS.R

@@ -105,30 +105,6 @@ setMethod("coerce", c("GRanges", "CTSS"),
   as(gp, "CTSS")
 })
 
-
-#' @description The `CTSS.chr` class represents a `CTSS` object that is 
-#' guaranteed to be only on a single chromosome.  It is used internally by
-#' _CAGEr_ for type-safe polymorphic dispatch.
-#' 
-#' @rdname CTSS-class
-#' 
-#' @import methods
-#' @importFrom GenomicRanges GPos
-#' 
-#' @examples
-#' 
-#' # (internal use) Transform CTSS to CTSS.chr object
-#' ctss.chr <- as(CTSScoordinatesGR(exampleCAGEexp), "CTSS.chr")
-
-setClass( "CTSS.chr"
-        , contains = "CTSS"
-        , validity =
-  function(object) 
-    if (length(seqlevelsInUse(object)) > 1)
-      return("Mutiple sequnames found: CTSS.chr objects should be only on a single chromosome.")
-)
-
-
 #' ConsensusClusters
 #' 
 #' The \code{ConsensusClusters} class represents consensus clusters.