@@ -3,7 +3,9 @@ CHANGES IN VERSION 1.5.1

 SIGNIFICANT USER-LEVEL CHANGES

-    o Column 'peakScores' is now bounded between 0 and 1000 and approximately uniformly distributed.

+    o New column 'maxPostInPeak' containing the maximum posterior within each peak.

+

+    o Score in exported BED files is calculated as -10*log10(maxPostInPeak).

     o 'changeFDR()' was renamed to 'adjustSensitivity()'.

@@ -75,7 +75,11 @@ adjustSensitivity.multivariate <- function(model, sensitivity, invert=FALSE) {

     if (is(model, class.multivariate.hmm)) {

         ## Calculate states

         ptm <- startTimedMessage("Calculating states from maximum-posterior in each peak ...")

-        p <- getMaxPostInPeaks(model$bins$state, model$bins$posteriors)

+        if (is.null(model$bins$maxPostInPeak)) {

+            p <- getMaxPostInPeaks(model$bins$state, model$bins$posteriors)

+        } else {

+            p <- model$bins$maxPostInPeak

+        }

         p.thresholded <- matrix(FALSE, ncol=ncol(p), nrow=nrow(p))

         for (icol in 1:ncol(p)) {

             if (!invert) {

@@ -92,25 +96,23 @@ adjustSensitivity.multivariate <- function(model, sensitivity, invert=FALSE) {

             model$bins$combination <- factor(mapping[as.character(model$bins$state)], levels=mapping[as.character(levels(model$bins$state))])

         }

         stopTimedMessage(ptm)

-        ## Redo peakScores

-        ptm <- startTimedMessage("Re-estimating peakScores ...")

-        p <- getMaxPostInPeaks(model$bins$state, model$bins$posteriors)

-        model$bins$peakScores <- calculatePeakScores(p)

-        rm(p)

+        ## Redo maxPostInPeak

+        ptm <- startTimedMessage("Re-estimating maximum posterior in peaks ...")

+        model$bins$maxPostInPeak <- getMaxPostInPeaks(model$bins$state, model$bins$posteriors)

         stopTimedMessage(ptm)

         ## Redo segmentation

         model$segments <- multivariateSegmentation(model$bins, column2collapseBy='state')

         ## Redo peaks

         ptm <- startTimedMessage("Recalculating peaks ...")

         model$peaks <- list()

-        for (i1 in 1:ncol(model$segments$peakScores)) {

-            mask <- model$segments$peakScores[,i1] > 0

+        for (i1 in 1:ncol(model$segments$maxPostInPeak)) {

+            mask <- model$segments$maxPostInPeak[,i1] > 0

             peaks <- model$segments[mask]

             mcols(peaks) <- NULL

-            peaks$peakScores <- model$segments$peakScores[mask,i1]

+            peaks$maxPostInPeak <- model$segments$maxPostInPeak[mask,i1]

             model$peaks[[i1]] <- peaks

         }

-        names(model$peaks) <- colnames(model$segments$peakScores)

+        names(model$peaks) <- colnames(model$segments$maxPostInPeak)

         stopTimedMessage(ptm)

     } else if (is(model, class.combined.multivariate.hmm)) {

@@ -118,7 +120,11 @@ adjustSensitivity.multivariate <- function(model, sensitivity, invert=FALSE) {

         mapping.df <- stateBrewer(model$info[,setdiff(names(model$info), 'ID')], mode='full')

         mapping <- mapping.df$combination

         names(mapping) <- mapping.df$state

-        p <- getMaxPostInPeaks(model$bins$state, model$bins$posteriors)

+        if (is.null(model$bins$maxPostInPeak)) {

+            p <- getMaxPostInPeaks(model$bins$state, model$bins$posteriors)

+        } else {

+            p <- model$bins$maxPostInPeak

+        }

         p.thresholded <- matrix(FALSE, ncol=ncol(p), nrow=nrow(p))

         for (icol in 1:ncol(p)) {

             if (!invert) {

@@ -137,11 +143,9 @@ adjustSensitivity.multivariate <- function(model, sensitivity, invert=FALSE) {

         multiHMM$bins$combination <- mapping[as.character(states)]

         multiHMM$bins$state <- factor(states)

         multiHMM$bins$posteriors <- model$bins$posteriors

-        ## Redo peakScores

-        ptm <- startTimedMessage("Re-estimating peakScores ...")

-        p <- getMaxPostInPeaks(multiHMM$bins$state, multiHMM$bins$posteriors)

-        multiHMM$bins$peakScores <- calculatePeakScores(p)

-        rm(p)

+        ## Redo maxPostInPeak

+        ptm <- startTimedMessage("Re-estimating maximum posterior in peaks ...")

+        multiHMM$bins$maxPostInPeak <- getMaxPostInPeaks(multiHMM$bins$state, multiHMM$bins$posteriors)

         stopTimedMessage(ptm)

         multiHMM$mapping <- mapping

         multiHMM$peaks <- model$peaks

@@ -149,14 +153,14 @@ adjustSensitivity.multivariate <- function(model, sensitivity, invert=FALSE) {

         ## Redo peaks

         ptm <- startTimedMessage("Recalculating peaks ...")

         model$peaks <- list()

-        for (i1 in 1:ncol(model$segments$peakScores)) {

-            mask <- model$segments$peakScores[,i1] > 0

+        for (i1 in 1:ncol(model$segments$maxPostInPeak)) {

+            mask <- model$segments$maxPostInPeak[,i1] > 0

             peaks <- model$segments[mask]

             mcols(peaks) <- NULL

-            peaks$peakScores <- model$segments$peakScores[mask,i1]

+            peaks$maxPostInPeak <- model$segments$maxPostInPeak[mask,i1]

             model$peaks[[i1]] <- peaks

         }

-        names(model$peaks) <- colnames(model$segments$peakScores)

+        names(model$peaks) <- colnames(model$segments$maxPostInPeak)

         stopTimedMessage(ptm)

     } else {

         stop("Supply either a uniHMM, multiHMM or combinedMultiHMM object.")

@@ -187,7 +191,9 @@ adjustSensitivity.univariate <- function(model, sensitivity, invert=FALSE) {

     if (is.null(model$bins$posterior.modified)) stop("Cannot recalculate states because column 'posterior.modified' is missing.")

     ## Calculate states

     ptm <- startTimedMessage("Calculating states from maximum-posterior in each peak ...")

-    model$bins$maxPostInPeak <- getMaxPostInPeaks.univariate(model$bins$state, model$bins$posterior.modified)

+    if (is.null(model$bins$maxPostInPeak)) {

+        model$bins$maxPostInPeak <- getMaxPostInPeaks.univariate(model$bins$state, model$bins$posterior.modified)

+    }

     states <- model$bins$state

     states[model$bins$state == 'modified'] <- 'unmodified'

     if (!invert) {

@@ -198,17 +204,15 @@ adjustSensitivity.univariate <- function(model, sensitivity, invert=FALSE) {

     model$bins$state <- states

     model$bins$maxPostInPeak <- NULL

     stopTimedMessage(ptm)

-    ## Redo peakScores

-    ptm <- startTimedMessage("Re-estimating peakScores ...")

+    ## Redo maxPostInPeak

+    ptm <- startTimedMessage("Re-estimating maximum posterior in peaks ...")

     model$bins$maxPostInPeak <- getMaxPostInPeaks.univariate(model$bins$state, model$bins$posterior.modified)

-    model$bins$peakScores <- calculatePeakScores.univariate(model$bins$maxPostInPeak)

-    model$bins$maxPostInPeak <- NULL

     stopTimedMessage(ptm)

     ## Redo segmentation

     ptm <- startTimedMessage("Making segmentation ...")

     gr <- model$bins

     df <- as.data.frame(gr)

-    red.df <- suppressMessages(collapseBins(df, column2collapseBy='state', columns2average=c('score'), columns2drop=c('width',grep('posteriors', names(df), value=TRUE), 'counts', 'counts.rpkm', 'posterior.modified')))

+    red.df <- suppressMessages(collapseBins(df, column2collapseBy='state', columns2drop=c('width',grep('posteriors', names(df), value=TRUE), 'counts', 'counts.rpkm', 'posterior.modified')))

     model$segments <- methods::as(red.df, 'GRanges')

     seqlengths(model$segments) <- seqlengths(model$bins)[seqlevels(model$segments)]

     ## Redo peaks

@@ -317,25 +321,23 @@ changePostCutoff.multivariate <- function(model, post.cutoff) {

             model$bins$combination <- factor(mapping[as.character(model$bins$state)], levels=mapping[as.character(levels(model$bins$state))])

         }

         stopTimedMessage(ptm)

-        ## Redo peakScores

-        ptm <- startTimedMessage("Re-estimating peakScores ...")

-        p <- getMaxPostInPeaks(model$bins$state, model$bins$posteriors)

-        model$bins$peakScores <- calculatePeakScores(p)

-        rm(p)

+        ## Redo maxPostInPeak

+        ptm <- startTimedMessage("Re-estimating maximum posterior in peaks ...")

+        model$bins$maxPostInPeak <- getMaxPostInPeaks(model$bins$state, model$bins$posteriors)

         stopTimedMessage(ptm)

         ## Redo segmentation

         model$segments <- multivariateSegmentation(model$bins, column2collapseBy='state')

         ## Redo peaks

         ptm <- startTimedMessage("Recalculating peaks ...")

         model$peaks <- list()

-        for (i1 in 1:ncol(model$segments$peakScores)) {

-            mask <- model$segments$peakScores[,i1] > 0

+        for (i1 in 1:ncol(model$segments$maxPostInPeak)) {

+            mask <- model$segments$maxPostInPeak[,i1] > 0

             peaks <- model$segments[mask]

             mcols(peaks) <- NULL

-            peaks$peakScores <- model$segments$peakScores[mask,i1]

+            peaks$maxPostInPeak <- model$segments$maxPostInPeak[mask,i1]

             model$peaks[[i1]] <- peaks

         }

-        names(model$peaks) <- colnames(model$segments$peakScores)

+        names(model$peaks) <- colnames(model$segments$maxPostInPeak)

         stopTimedMessage(ptm)

     } else if (is(model, class.combined.multivariate.hmm)) {

         ptm <- startTimedMessage("Calculating states from posteriors ...")

@@ -357,11 +359,9 @@ changePostCutoff.multivariate <- function(model, post.cutoff) {

         multiHMM$bins$combination <- mapping[as.character(states)]

         multiHMM$bins$state <- factor(states)

         multiHMM$bins$posteriors <- post

-        ## Redo peakScores

-        ptm <- startTimedMessage("Re-estimating peakScores ...")

-        p <- getMaxPostInPeaks(multiHMM$bins$state, multiHMM$bins$posteriors)

-        multiHMM$bins$peakScores <- calculatePeakScores(p)

-        rm(p)

+        ## Redo maxPostInPeak

+        ptm <- startTimedMessage("Re-estimating maximum posterior in peaks ...")

+        multiHMM$bins$maxPostInPeak <- getMaxPostInPeaks(multiHMM$bins$state, multiHMM$bins$posteriors)

         stopTimedMessage(ptm)

         multiHMM$mapping <- mapping

         multiHMM$peaks <- model$peaks

@@ -369,14 +369,14 @@ changePostCutoff.multivariate <- function(model, post.cutoff) {

         ## Redo peaks

         ptm <- startTimedMessage("Recalculating peaks ...")

         model$peaks <- list()

-        for (i1 in 1:ncol(model$segments$peakScores)) {

-            mask <- model$segments$peakScores[,i1] > 0

+        for (i1 in 1:ncol(model$segments$maxPostInPeak)) {

+            mask <- model$segments$maxPostInPeak[,i1] > 0

             peaks <- model$segments[mask]

             mcols(peaks) <- NULL

-            peaks$peakScores <- model$segments$peakScores[mask,i1]

+            peaks$maxPostInPeak <- model$segments$maxPostInPeak[mask,i1]

             model$peaks[[i1]] <- peaks

         }

-        names(model$peaks) <- colnames(model$segments$peakScores)

+        names(model$peaks) <- colnames(model$segments$maxPostInPeak)

         stopTimedMessage(ptm)

     } else {

         stop("Supply either a uniHMM, multiHMM or combinedMultiHMM object.")

@@ -416,11 +416,9 @@ changePostCutoff.univariate <- function(model, post.cutoff) {

     states <- state.labels[states]

     model$bins$state <- states

     stopTimedMessage(ptm)

-    ## Redo peakScores

-    ptm <- startTimedMessage("Re-estimating peakScores ...")

+    ## Redo maxPostInPeak

+    ptm <- startTimedMessage("Re-estimating maximum posterior in peaks ...")

     model$bins$maxPostInPeak <- getMaxPostInPeaks.univariate(model$bins$state, model$bins$posterior.modified)

-    model$bins$peakScores <- calculatePeakScores.univariate(model$bins$maxPostInPeak)

-    model$bins$maxPostInPeak <- NULL

     stopTimedMessage(ptm)

     ## Redo segmentation

     ptm <- startTimedMessage("Making segmentation ...")

@@ -498,11 +498,9 @@ runMultivariate <- function(binned.data, stepbins, info, comb.states, use.states

     stopTimedMessage(ptm)

     if (get.posteriors) {

-        ptm <- startTimedMessage("Calculating peak scores ...")

+        ptm <- startTimedMessage("Getting maximum posterior in peaks ...")

         result$bins$maxPostInPeak <- getMaxPostInPeaks(result$bins$state, result$bins$posteriors)

-        result$bins$peakScores <- calculatePeakScores(result$bins$maxPostInPeak)

-        result$bins$maxPostInPeak <- NULL

-        result$bins$differential.score <- differentialScoreSum(result$bins$peakScores, result$info)

+        result$bins$differential.score <- differentialScoreSum(result$bins$maxPostInPeak, result$info)

         stopTimedMessage(ptm)

     }

     if (!keep.posteriors) {

@@ -527,14 +525,14 @@ runMultivariate <- function(binned.data, stepbins, info, comb.states, use.states

     ## Peaks ##

     ptm <- startTimedMessage("Obtaining peaks ...")

     result$peaks <- list()

-    for (i1 in 1:ncol(result$segments$peakScores)) {

-        mask <- result$segments$peakScores[,i1] > 0

+    for (i1 in 1:ncol(result$segments$maxPostInPeak)) {

+        mask <- result$segments$maxPostInPeak[,i1] > 0

         peaks <- result$segments[mask]

         mcols(peaks) <- NULL

-        peaks$peakScores <- result$segments$peakScores[mask,i1]

+        peaks$maxPostInPeak <- result$segments$maxPostInPeak[mask,i1]

         result$peaks[[i1]] <- peaks

     }

-    names(result$peaks) <- colnames(result$segments$peakScores)

+    names(result$peaks) <- colnames(result$segments$maxPostInPeak)

     stopTimedMessage(ptm)

     return(result)

@@ -560,8 +560,6 @@ callPeaksUnivariateAllChr <- function(binned.data, input.data=NULL, eps=0.01, in

     ## Peak score as maximum posterior in that peak ##

     result$bins$maxPostInPeak <- getMaxPostInPeaks.univariate(result$bins$state, result$bins$posterior.modified)

-    result$bins$peakScores <- calculatePeakScores.univariate(result$bins$maxPostInPeak)

-    result$bins$maxPostInPeak <- NULL

     ## Segmentation ##

     ptm <- startTimedMessage("Making segmentation ...")

@@ -92,8 +92,8 @@ combineMultivariates <- function(hmms, mode) {

         counts[[1]] <- hmm$bins$counts.rpkm

         posteriors <- list()

         posteriors[[1]] <- hmm$bins$posteriors

-        peakScores <- list()

-        peakScores[[1]] <- hmm$bins$peakScores

+        maxPostInPeak <- list()

+        maxPostInPeak[[1]] <- hmm$bins$maxPostInPeak

         peaks <- list()

         peaks[[1]] <- hmm$peaks

         binstates <- list()

@@ -108,7 +108,7 @@ combineMultivariates <- function(hmms, mode) {

                 combs[[i1]] <- hmm$bins$combination

                 counts[[i1]] <- hmm$bins$counts.rpkm

                 posteriors[[i1]] <- hmm$bins$posteriors

-                peakScores[[i1]] <- hmm$bins$peakScores

+                maxPostInPeak[[i1]] <- hmm$bins$maxPostInPeak

                 peaks[[i1]] <- hmm$peaks

                 binstates[[i1]] <- dec2bin(hmm$bins$state, colnames=hmm$info$ID)

                 stopTimedMessage(ptm)

@@ -116,7 +116,7 @@ combineMultivariates <- function(hmms, mode) {

         }

         ptm <- startTimedMessage("Concatenating conditions ...")

         posteriors <- do.call(cbind, posteriors)

-        peakScores <- do.call(cbind, peakScores)

+        maxPostInPeak <- do.call(cbind, maxPostInPeak)

         infos <- do.call(rbind, infos)

         conditions <- unique(infos$condition)

         states <- factor(bin2dec(do.call(cbind, binstates)))

@@ -131,7 +131,7 @@ combineMultivariates <- function(hmms, mode) {

         infos <- list()

         counts <- list()

         posteriors <- list()

-        peakScores <- list()

+        maxPostInPeak <- list()

         peaks <- list()

         binstates <- list()

         for (i1 in 1:length(hmms)) {

@@ -140,7 +140,7 @@ combineMultivariates <- function(hmms, mode) {

             infos[[i1]] <- hmm$info

             counts[[i1]] <- hmm$bins$counts.rpkm

             posteriors[[i1]] <- hmm$bins$posteriors

-            peakScores[[i1]] <- hmm$bins$peakScores

+            maxPostInPeak[[i1]] <- hmm$bins$maxPostInPeak

             peaks[[i1]] <- hmm$peaks

             binstates[[i1]] <- dec2bin(hmm$bins$state, colnames=hmm$info$ID)

             stopTimedMessage(ptm)

@@ -158,8 +158,8 @@ combineMultivariates <- function(hmms, mode) {

         counts <- counts[,infos$ID]

         posteriors <- do.call(cbind, posteriors)

         posteriors <- posteriors[,infos$ID]

-        peakScores <- do.call(cbind, peakScores)

-        peakScores <- peakScores[,infos$ID]

+        maxPostInPeak <- do.call(cbind, maxPostInPeak)

+        maxPostInPeak <- maxPostInPeak[,infos$ID]

         binstates <- do.call(cbind, binstates)

         binstates <- binstates[,infos$ID]

         states <- factor(bin2dec(binstates))

@@ -199,7 +199,7 @@ combineMultivariates <- function(hmms, mode) {

         conditions <- unique(infos$condition)

         counts <- hmm$bins$counts.rpkm

         posteriors <- hmm$bins$posteriors

-        peakScores <- hmm$bins$peakScores

+        maxPostInPeak <- hmm$bins$maxPostInPeak

         peaks <- hmm$peaks

         states <- hmm$bins$state

         combs <- list()

@@ -227,8 +227,8 @@ combineMultivariates <- function(hmms, mode) {

         counts[[1]] <- hmm$bins$counts.rpkm

         posteriors <- list()

         posteriors[[1]] <- hmm$bins$posteriors

-        peakScores <- list()

-        peakScores[[1]] <- hmm$bins$peakScores

+        maxPostInPeak <- list()

+        maxPostInPeak[[1]] <- hmm$bins$maxPostInPeak

         peaks <- list()

         peaks[[1]] <- hmm$peaks

         binstates <- list()

@@ -242,7 +242,7 @@ combineMultivariates <- function(hmms, mode) {

                 infos[[i1]] <- hmm$info

                 counts[[i1]] <- hmm$bins$counts.rpkm

                 posteriors[[i1]] <- hmm$bins$posteriors

-                peakScores[[i1]] <- hmm$bins$peakScores

+                maxPostInPeak[[i1]] <- hmm$bins$maxPostInPeak

                 peaks[[i1]] <- hmm$peaks

                 binstates[[i1]] <- dec2bin(hmm$bins$state, colnames=hmm$info$ID)

                 stopTimedMessage(ptm)

@@ -251,7 +251,7 @@ combineMultivariates <- function(hmms, mode) {

         ptm <- startTimedMessage("Concatenating mark-conditions ...")

         counts <- do.call(cbind, counts)

         posteriors <- do.call(cbind, posteriors)

-        peakScores <- do.call(cbind, peakScores)

+        maxPostInPeak <- do.call(cbind, maxPostInPeak)

         binstates <- do.call(cbind, binstates)

         infos <- do.call(rbind, infos)

         conditions <- unique(infos$condition)

@@ -297,12 +297,12 @@ combineMultivariates <- function(hmms, mode) {

     ptm <- startTimedMessage("Transferring counts and posteriors ...")

     bins$counts.rpkm <- counts

     bins$posteriors <- posteriors

-    bins$peakScores <- peakScores

+    bins$maxPostInPeak <- maxPostInPeak

     stopTimedMessage(ptm)

     ## Add differential score ##

     ptm <- startTimedMessage("Adding differential score ...")

-    bins$differential.score <- differentialScoreSum(bins$peakScores, infos)

+    bins$differential.score <- differentialScoreSum(bins$maxPostInPeak, infos)

     stopTimedMessage(ptm)

     ### Redo the segmentation for all conditions combined

@@ -193,8 +193,11 @@ exportUnivariatePeaks <- function(hmm.list, filename, header=TRUE, separate.file

             }

             cat(paste0("track name=\"",trackname.string,"\" description=\"",trackname.string,"\" visibility=1 itemRgb=On priority=",priority,"\n"), file=filename.gz, append=TRUE)

         }

-        if (is.null(peaks$peakScores)) {

+        if (is.null(peaks$maxPostInPeak)) {

             peaks$peakScores <- 0

+        } else {

+            peaks$peakScores <- suppressWarnings( -10*log10(1-peaks$maxPostInPeak) )

+            peaks$peakScores[is.nan(peaks$peakScores) | peaks$peakScores > 1000] <- 1000

         }

         df <- as.data.frame(peaks)

         df$peakNumber <- paste0('peak_', 1:nrow(df))

@@ -434,12 +437,19 @@ exportMultivariatePeaks <- function(hmm, filename, trackname=NULL, header=TRUE,

         cat("", file=filename.gz)

     }

-    ## Export peaks on a per bin basis to obtain proper posterior scores

+    ## Export peaks

     for (imod in 1:length(hmm$info$ID)) {

         ID <- hmm$info$ID[imod]

-        ## Select only segments with peaks

+        ## Add peak score

         peaks <- insertchr(hmm$peaks[[ID]])

+        if (is.null(peaks$maxPostInPeak)) {

+            peaks$peakScores <- 0

+        } else {

+            peaks$peakScores <- suppressWarnings( -10*log10(1-peaks$maxPostInPeak) )

+            peaks$peakScores[is.nan(peaks$peakScores) | peaks$peakScores > 1000] <- 1000

+        }

+        

         peaks.df <- as.data.frame(peaks)

         peaks.df$peakNumber <- paste0('peak_', 1:nrow(peaks.df))

@@ -606,12 +616,20 @@ exportCombinedMultivariatePeaks <- function(hmm, filename, trackname=NULL, heade

         cat("", file=filename.gz)

     }

+    ## Export peaks

     for (imod in 1:length(hmm$info$ID)) {

         ID <- hmm$info$ID[imod]

         cond <- hmm$info$condition[imod]

-        ## Select only segments with peaks

+        ## Add peak score

         peaks <- insertchr(hmm$peaks[[ID]])

+        if (is.null(peaks$maxPostInPeak)) {

+            peaks$peakScores <- 0

+        } else {

+            peaks$peakScores <- suppressWarnings( -10*log10(1-peaks$maxPostInPeak) )

+            peaks$peakScores[is.nan(peaks$peakScores) | peaks$peakScores > 1000] <- 1000

+        }

+        

         peaks.df <- as.data.frame(peaks)

         peaks.df$peakNumber <- paste0('peak_', 1:nrow(peaks.df))

@@ -2,7 +2,7 @@ getMaxPostInPeaks <- function(states, posteriors) {

     binstates <- dec2bin(states, colnames = colnames(posteriors))

     rownames(binstates) <- NULL

-    peakScores <- matrix(NA, nrow=nrow(binstates), ncol=ncol(binstates), dimnames=list(NULL, colnames(binstates)))

+    maxPostInPeak <- matrix(NA, nrow=nrow(binstates), ncol=ncol(binstates), dimnames=list(NULL, colnames(binstates)))

     for (icol in 1:ncol(binstates)) {

         r.bin <- rle(binstates[,icol])

         r <- r.bin

@@ -16,9 +16,9 @@ getMaxPostInPeaks <- function(states, posteriors) {

         r <- r.bin

         r$values[r$values == TRUE] <- df$x

         r$values[r$values == FALSE] <- 0

-        peakScores[,icol] <- inverse.rle(r)

+        maxPostInPeak[,icol] <- inverse.rle(r)

     }

-    return(peakScores)

+    return(maxPostInPeak)

 }

@@ -76,21 +76,21 @@ removeCondition <- function(model, conditions) {

         }

         model$bins$posteriors <- posteriors

     }

-    if (!is.null(model$bins$peakScores)) {

-        peakScores <- model$bins$peakScores

+    if (!is.null(model$bins$maxPostInPeak)) {

+        maxPostInPeak <- model$bins$maxPostInPeak

         removeconds <- paste0(paste0('-', conditions, '-'), collapse='|')

-        keepconds <- grep(removeconds, colnames(peakScores), invert=TRUE, value=TRUE)

-        peakScores <- peakScores[,keepconds]

-        if (class(peakScores) != "matrix") {

-            peakScores <- matrix(peakScores, ncol=1, dimnames=list(NULL, keepconds))

+        keepconds <- grep(removeconds, colnames(maxPostInPeak), invert=TRUE, value=TRUE)

+        maxPostInPeak <- maxPostInPeak[,keepconds]

+        if (class(maxPostInPeak) != "matrix") {

+            maxPostInPeak <- matrix(maxPostInPeak, ncol=1, dimnames=list(NULL, keepconds))

         }

-        model$bins$peakScores <- peakScores

+        model$bins$maxPostInPeak <- maxPostInPeak

     }

     stopTimedMessage(ptm)

     # Redo differential score

     if (!is.null(model$bins$differential.score)) {

         ptm <- startTimedMessage("Differential score ...")

-        model$bins$differential.score <- differentialScoreSum(model$bins$peakScores, model$info)

+        model$bins$differential.score <- differentialScoreSum(model$bins$maxPostInPeak, model$info)

         stopTimedMessage(ptm)

     }

     # Redo transition frequencies and groups

@@ -12,16 +12,16 @@ multivariateSegmentation <- function(bins, column2collapseBy='state') {

     df <- as.data.frame(bins)

     ind.readcols <- grep('^counts', names(df))

     ind.postcols <- grep('^posteriors', names(df))

-    ind.peakcols <- grep('^peakScores', names(df))

+    ind.peakcols <- grep('^maxPostInPeak', names(df))

     ind.widthcol <- grep('width', names(df))

     ind.scorecol <- grep('differential.score', names(df))

     red.df <- suppressMessages(collapseBins(df, column2collapseBy=column2collapseBy, columns2average=c(ind.scorecol), columns2drop=c(ind.readcols, ind.widthcol, ind.postcols)))

     names(red.df) <- sub('^mean.','', names(red.df))

-    red.peakscores <- as.matrix(red.df[grep('^peakScores', names(red.df))])

-    red.df[grep('^peakScores', names(red.df))] <- NULL

+    red.maxPostInPeak <- as.matrix(red.df[grep('^maxPostInPeak', names(red.df))])

+    red.df[grep('^maxPostInPeak', names(red.df))] <- NULL

     segments <- methods::as(red.df, 'GRanges')

-    segments$peakScores <- red.peakscores

-    colnames(segments$peakScores) <- colnames(bins$peakScores)

+    segments$maxPostInPeak <- red.maxPostInPeak

+    colnames(segments$maxPostInPeak) <- colnames(bins$maxPostInPeak)

     ## Reorder properly to match the order in bins

     diffscore.temp <- segments$differential.score

     segments$differential.score <- NULL

Binary files a/vignettes/chromstaR.pdf and b/vignettes/chromstaR.pdf differ