@@ -1,8 +1,8 @@

 Package: chromstaR

 Type: Package

 Title: Combinatorial and Differential Chromatin State Analysis for ChIP-Seq Data

-Version: 1.1.0

-Date: 2014-05-03

+Version: 1.1.1

+Date: 2016-10-19

 Author: Aaron Taudt, Maria Colome Tatche, Matthias Heinig, Minh Anh Nguyen

 Maintainer: Aaron Taudt <aaron.taudt@gmail.com>

 Description: This package implements functions for combinatorial and differential analysis of ChIP-seq data. It includes uni- and multivariate peak-calling, export to genome browser viewable files, and functions for enrichment analyses.

@@ -6,17 +6,17 @@ export(binReads)

 export(callPeaksMultivariate)

 export(callPeaksReplicates)

 export(callPeaksUnivariate)

-export(changePostCutoff)

+export(changeFDR)

 export(collapseBins)

 export(combineMultivariates)

 export(dec2bin)

-export(exportBinnedData)

-export(exportCombinedMultivariate)

+export(exportCombinations)

+export(exportCounts)

 export(exportGRangesAsBedFile)

-export(exportMultivariate)

-export(exportUnivariates)

+export(exportPeaks)

 export(fixedWidthBins)

 export(genomicFrequencies)

+export(getCombinations)

 export(getDistinctColors)

 export(getStateColors)

 export(heatmapCombinations)

@@ -24,9 +24,9 @@ export(heatmapCountCorrelation)

 export(heatmapTransitionProbs)

 export(loadHmmsFromFiles)

 export(plotEnrichCountHeatmap)

+export(plotEnrichment)

 export(plotExpression)

 export(plotFoldEnrichHeatmap)

-export(plotFoldEnrichment)

 export(plotHistogram)

 export(readBamFileAsGRanges)

 export(readBedFileAsGRanges)

new file mode 100644

@@ -0,0 +1,26 @@

+CHANGES IN VERSION 1.1.1

+------------------------

+

+NEW FEATURES

+

+    o Selection of peaks can be done with 'changeFDR'.

+

+    o Peak calls are available in each chromstaR-object as list entry '$peaks'.

+

+SIGNIFICANT USER-LEVEL CHANGES

+

+    o 'plotFoldEnrichment' renamed to 'plotEnrichment'.

+

+    o 'exportBinnedData', 'exportUnivariates', 'exportMultivariates', 'exportCombinedMultivariates' replaced by 'exportCounts', 'exportPeaks', 'exportCombinations'.

+

+BUG FIXES

+

+    o Proper computation of fold enrichments, with < 1 indicating depletion and > 1 indicating enrichment.

+

+DEPRECATED AND DEFUNCT

+

+    o 'changePostCutoff'.

+

+    o 'plotFoldEnrichment'.

+

+    o 'exportBinnedData', 'exportUnivariates', 'exportMultivariates', 'exportCombinedMultivariates'.

@@ -11,8 +11,7 @@

 #' @param per.chrom If \code{per.chrom=TRUE} chromosomes will be treated separately. This tremendously speeds up the calculation but results might be noisier as compared to \code{per.chrom=FALSE}, where all chromosomes are concatenated for the HMM.

 #' @param chromosomes A vector specifying the chromosomes to use from the models in \code{hmms}. The default (\code{NULL}) uses all available chromosomes.

 #' @param eps Convergence threshold for the Baum-Welch algorithm.

-#' @param post.cutoff False discovery rate. The default \code{NULL} means that the state with maximum posterior probability will be chosen, irrespective of its absolute probability.

-#' @param keep.posteriors If set to \code{TRUE}, posteriors will be available in the output. This is useful to change the post.cutoff later, but increases the necessary disk space to store the result immense.

+#' @param keep.posteriors If set to \code{TRUE}, posteriors will be available in the output. This can be useful to change the posterior cutoff later, but increases the necessary disk space to store the result immensely.

 #' @param num.threads Number of threads to use. Setting this to >1 may give increased performance.

 #' @param max.time The maximum running time in seconds for the Baum-Welch algorithm. If this time is reached, the Baum-Welch will terminate after the current iteration finishes. The default \code{NULL} is no limit.

 #' @param max.iter The maximum number of iterations for the Baum-Welch algorithm. The default \code{NULL} is no limit.

@@ -50,7 +49,7 @@

 #'heatmapTransitionProbs(multimodel)

 #'heatmapCountCorrelation(multimodel)

 #'

-callPeaksMultivariate <- function(hmms, use.states, max.states=NULL, per.chrom=TRUE, chromosomes=NULL, eps=0.01, post.cutoff=NULL, keep.posteriors=TRUE, num.threads=1, max.time=NULL, max.iter=NULL, keep.densities=FALSE, verbosity=1) {

+callPeaksMultivariate <- function(hmms, use.states, max.states=NULL, per.chrom=TRUE, chromosomes=NULL, eps=0.01, keep.posteriors=FALSE, num.threads=1, max.time=NULL, max.iter=NULL, keep.densities=FALSE, verbosity=1) {

     ## Intercept user input

     if (!is.null(use.states)) {

@@ -71,9 +70,6 @@ callPeaksMultivariate <- function(hmms, use.states, max.states=NULL, per.chrom=T

     if (check.logical(keep.posteriors)!=0) stop("argument 'keep.posteriors' expects a logical (TRUE or FALSE)")

     if (check.logical(keep.densities)!=0) stop("argument 'keep.densities' expects a logical (TRUE or FALSE)")

     if (check.integer(verbosity)!=0) stop("argument 'verbosity' expects an integer")

-    if (!is.null(post.cutoff)) {

-        if (post.cutoff>1 | post.cutoff<0) stop("argument 'post.cutoff' has to be between 0 and 1 if specified")

-    }

     if (length(hmms)==0) {

         stop("argument 'hmms' is of length=0. Cannot call multivariate peaks with no models.")

     }

@@ -90,9 +86,6 @@ callPeaksMultivariate <- function(hmms, use.states, max.states=NULL, per.chrom=T

             result$bins$state <- factor(c(0,0,1))[hmm$bins$state]

             result$bins$posteriors <- matrix(hmm$bins$posterior.modified, ncol=1, dimnames=list(NULL, result$info$ID))

             result$bins$differential.score <- 0

-            if (!is.null(post.cutoff)) {

-                result$bins$state <- factor(as.integer(result$bins$posteriors >= post.cutoff))

-            }

         ## Add combinations

             mapping <- NULL

             if (!is.null(use.states)) {

@@ -125,8 +118,6 @@ callPeaksMultivariate <- function(hmms, use.states, max.states=NULL, per.chrom=T

             # Distributions

             result$distributions <- list(hmm$distributions)

             names(result$distributions) <- result$info$ID

-            # post.cutoff

-            result$post.cutoff <- post.cutoff

         ## Convergence info

             convergenceInfo <- list(eps=NA, loglik=NA, loglik.delta=NA, num.iterations=NA, time.sec=NA)

             result$convergenceInfo <- convergenceInfo

@@ -167,7 +158,7 @@ callPeaksMultivariate <- function(hmms, use.states, max.states=NULL, per.chrom=T

             ptm <- startTimedMessage("Running multivariate ...")

             models <- foreach (chrom = chromosomes, .packages='chromstaR') %dopar% {

                 bins <- p$bins[seqnames(p$bins)==chrom]

-                model <- runMultivariate(bins=bins, info=p$info, comb.states=p$comb.states, use.states=p$use.states, distributions=p$distributions, weights=p$weights, correlationMatrix=p$correlationMatrix, correlationMatrixInverse=p$correlationMatrixInverse, determinant=p$determinant, max.iter=max.iter, max.time=max.time, eps=eps, num.threads=1, post.cutoff=post.cutoff, keep.posteriors=keep.posteriors, keep.densities=keep.densities, verbosity=verbosity)

+                model <- runMultivariate(bins=bins, info=p$info, comb.states=p$comb.states, use.states=p$use.states, distributions=p$distributions, weights=p$weights, correlationMatrix=p$correlationMatrix, correlationMatrixInverse=p$correlationMatrixInverse, determinant=p$determinant, max.iter=max.iter, max.time=max.time, eps=eps, num.threads=1, keep.posteriors=keep.posteriors, keep.densities=keep.densities, verbosity=verbosity)

                 model

             }

             stopTimedMessage(ptm)

@@ -176,7 +167,7 @@ callPeaksMultivariate <- function(hmms, use.states, max.states=NULL, per.chrom=T

             for (chrom in chromosomes) {

                 ptm <- startTimedMessage("Chromosome = ", chrom, "\n")

                 bins <- p$bins[seqnames(p$bins)==chrom]

-                model <- suppressMessages( runMultivariate(bins=bins, info=p$info, comb.states=p$comb.states, use.states=p$use.states, distributions=p$distributions, weights=p$weights, correlationMatrix=p$correlationMatrix, correlationMatrixInverse=p$correlationMatrixInverse, determinant=p$determinant, max.iter=max.iter, max.time=max.time, eps=eps, num.threads=1, post.cutoff=post.cutoff, keep.posteriors=keep.posteriors, keep.densities=keep.densities, verbosity=verbosity) )

+                model <- suppressMessages( runMultivariate(bins=bins, info=p$info, comb.states=p$comb.states, use.states=p$use.states, distributions=p$distributions, weights=p$weights, correlationMatrix=p$correlationMatrix, correlationMatrixInverse=p$correlationMatrixInverse, determinant=p$determinant, max.iter=max.iter, max.time=max.time, eps=eps, num.threads=1, keep.posteriors=keep.posteriors, keep.densities=keep.densities, verbosity=verbosity) )

                 message("Time spent for chromosome = ", chrom, ":", appendLF=FALSE)

                 stopTimedMessage(ptm)

                 models[[chrom]] <- model

@@ -191,7 +182,7 @@ callPeaksMultivariate <- function(hmms, use.states, max.states=NULL, per.chrom=T

     ## Run multivariate for all chromosomes

     } else {

-        model <- runMultivariate(bins=p$bins, info=p$info, comb.states=p$comb.states, use.states=p$use.states, distributions=p$distributions, weights=p$weights, correlationMatrix=p$correlationMatrix, correlationMatrixInverse=p$correlationMatrixInverse, determinant=p$determinant, max.iter=max.iter, max.time=max.time, eps=eps, num.threads=num.threads, post.cutoff=post.cutoff, keep.posteriors=keep.posteriors, keep.densities=keep.densities, verbosity=verbosity)

+        model <- runMultivariate(bins=p$bins, info=p$info, comb.states=p$comb.states, use.states=p$use.states, distributions=p$distributions, weights=p$weights, correlationMatrix=p$correlationMatrix, correlationMatrixInverse=p$correlationMatrixInverse, determinant=p$determinant, max.iter=max.iter, max.time=max.time, eps=eps, num.threads=num.threads, keep.posteriors=keep.posteriors, keep.densities=keep.densities, verbosity=verbosity)

     }

@@ -200,7 +191,7 @@ callPeaksMultivariate <- function(hmms, use.states, max.states=NULL, per.chrom=T

 }

-runMultivariate <- function(bins, info, comb.states, use.states, distributions, weights, correlationMatrix, correlationMatrixInverse, determinant, max.iter, max.time, eps, num.threads, post.cutoff, keep.posteriors, keep.densities, transitionProbs.initial=NULL, startProbs.initial=NULL, verbosity=1) {

+runMultivariate <- function(bins, info, comb.states, use.states, distributions, weights, correlationMatrix, correlationMatrixInverse, determinant, max.iter, max.time, eps, num.threads, keep.posteriors, keep.densities, transitionProbs.initial=NULL, startProbs.initial=NULL, verbosity=1) {

     ptm <- startTimedMessage("Starting multivariate HMM with ", length(comb.states), " combinatorial states")

     message("")

@@ -273,6 +264,12 @@ runMultivariate <- function(bins, info, comb.states, use.states, distributions,

         result$info <- info

     ## Bin coordinates, posteriors and states

         result$bins <- bins

+        if (!is.null(use.states$state)) {

+            state.levels <- levels(use.states$state)

+        } else {

+            state.levels <- hmm$comb.states

+        }

+        result$bins$state <- factor(hmm$states, levels=state.levels)

         if (get.posteriors) {

             ptm <- startTimedMessage("Transforming posteriors to `per sample` representation ...")

             hmm$posteriors <- matrix(hmm$posteriors, ncol=hmm$max.states)

@@ -281,20 +278,11 @@ runMultivariate <- function(bins, info, comb.states, use.states, distributions,

             post.per.track <- hmm$posteriors %*% binstates

             colnames(post.per.track) <- result$info$ID

             result$bins$posteriors <- post.per.track

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

+            result$bins$peakScores <- getPeakScores(result$bins)

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

             stopTimedMessage(ptm)

         }

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

-        if (!is.null(use.states$state)) {

-            state.levels <- levels(use.states$state)

-        } else {

-            state.levels <- hmm$comb.states

-        }

-        if (!is.null(post.cutoff)) {

-            result$bins$state <- factor(bin2dec(result$bins$posteriors >= post.cutoff), levels=state.levels)

-        } else {

-            result$bins$state <- factor(hmm$states, levels=state.levels)

-        }

         stopTimedMessage(ptm)

         if (keep.densities) {

             result$bins$densities <- matrix(hmm$densities, ncol=hmm$max.states)

@@ -314,6 +302,16 @@ runMultivariate <- function(bins, info, comb.states, use.states, distributions,

         if (!keep.posteriors) {

             result$bins$posteriors <- NULL

         }

+    ## Peaks

+        result$peaks <- list()

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

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

+            peaks <- result$segments[mask]

+            mcols(peaks) <- NULL

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

+            result$peaks[[i1]] <- peaks

+        }

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

     ## Parameters

         result$mapping <- mapping

         combinations <- mapping[as.character(comb.states)]

@@ -337,8 +335,6 @@ runMultivariate <- function(bins, info, comb.states, use.states, distributions,

         # Distributions

         result$distributions <- distributions

         names(result$distributions) <- result$info$ID

-        # post.cutoff

-        result$post.cutoff <- post.cutoff

     ## Convergence info

         convergenceInfo <- list(eps=eps, loglik=hmm$loglik, loglik.delta=hmm$loglik.delta, num.iterations=hmm$num.iterations, time.sec=hmm$time.sec)

         result$convergenceInfo <- convergenceInfo

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

                                                         ranges=ranges(binned.data),

                                                         counts=hmm$counts,

                                                         state=states) 

-        result$bins$score <- apply(hmm$posteriors, 1, max)

         result$bins$posteriors <- hmm$posteriors

         if (!control) {

             result$bins$posterior.modified <- hmm$posteriors[,'modified']

@@ -430,17 +429,22 @@ callPeaksUnivariateAllChr <- function(binned.data, input.data=NULL, eps=0.01, in

         }

         seqlengths(result$bins) <- seqlengths(binned.data)

         stopTimedMessage(ptm)

+    ## Peak score as maximum posterior in that peak

+        result$bins$peakScores <- getPeakScore.univariate(result$bins)

     ## Segmentation

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

         df <- as.data.frame(result$bins)

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

-        red.gr <- GRanges(seqnames=red.df[,1], ranges=IRanges(start=red.df[,2], end=red.df[,3]), strand=red.df[,4], state=red.df[,'state'], score=red.df[,'max.score'])

-        result$segments <- red.gr

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

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

         seqlengths(result$segments) <- seqlengths(binned.data)[seqlevels(result$segments)]

         if (!keep.posteriors) {

             result$bins$posteriors <- NULL

         }

         stopTimedMessage(ptm)

+    ## Peaks

+        result$peaks <- result$segments[result$segments$state == 'modified']

+        result$peaks$state <- NULL

+        result$segments <- NULL

     ## Parameters

         # Weights

         result$weights <- hmm$weights

new file mode 100644

@@ -0,0 +1,199 @@

+#' Change the false discovery rate of a Hidden Markov Model

+#'

+#' Adjusts the peak calls of a \code{\link{uniHMM}} or \code{\link{multiHMM}} object with the given false discovery rate.

+#'

+#' Each peak has a peak score associated (between 0 and 1). The false discovery rate is a simple cutoff on 1-peakScore, e.g. for an \code{fdr = 0.01} only peaks with a peak score \code{>= 0.99} (1-fdr) will be selected.

+#' 

+#' @author Aaron Taudt

+#' @param model A \code{\link{uniHMM}} or \code{\link{multiHMM}} object with posteriors.

+#' @param fdr A vector of values between 0 and 1 for each column in \code{model$bins$peakScores}. If only one value is given, it will be reused for all columns. Values close to 0 will yield more stringent peak calls with lower false positive but higher false negative rate.

+#' @param invert Select peaks below (\code{FDR}) or above (\code{TRUE}) the given \code{fdr}. This is useful to select low confidence peaks.

+#' @return The input object is returned with adjusted peak calls.

+#' @export

+#' @examples

+#'## Get an example BAM file

+#'file <- system.file("extdata", "euratrans",

+#'                       "lv-H3K27me3-BN-male-bio2-tech1.bam",

+#'                        package="chromstaRData")

+#'## Bin the file into bin size 1000bp

+#'data(rn4_chrominfo)

+#'binned <- binReads(file, assembly=rn4_chrominfo, binsizes=1000,

+#'                   chromosomes='chr12')

+#'## Fit the univariate Hidden Markov Model

+#'hmm <- callPeaksUnivariate(binned, max.time=60, eps=1)

+#'## Compare fits with different fdrs

+#'plotHistogram(hmm) + ylim(0,0.25) + ylim(0,0.3)

+#'plotHistogram(changeFDR(hmm, fdr=0.01)) + ylim(0,0.3)

+#'plotHistogram(changeFDR(hmm, fdr=1e-12)) + ylim(0,0.3)

+#'

+changeFDR <- function(model, fdr=0.01, invert=FALSE) {

+

+    if (!is.numeric(fdr)) {

+        warning("'fdr' is not numeric. Nothing done.")

+        return(model)

+    }

+    if (is(model, class.univariate.hmm)) {

+        model.new <- changeFDR.univariate(model=model, fdr=fdr, invert=invert)

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

+        model.new <- changeFDR.multivariate(model=model, fdr=fdr, invert=invert)

+    }

+    return(model.new)

+    

+}

+

+

+changeFDR.multivariate <- function(model, fdr, invert=FALSE) {

+  

+    ## Make fdr vector

+    if (is.null(model$bins$peakScores)) stop("Cannot recalculate states because peakScores are missing.")

+    numcol <- ncol(model$bins$peakScores)

+    if (length(fdr) == 1) {

+        fdr <- rep(fdr, numcol)

+    }

+    if (length(fdr) != numcol) {

+        stop("Need ", numcol, " values in 'fdr' but only ", length(fdr), " are provided.")

+    }

+    

+    for (fdr.i in fdr) {

+        if (fdr.i < 0 | fdr.i > 1) {

+            stop("Values for 'fdr' need to be inside the interval [0,1].")

+        }

+    }

+    names(fdr) <- colnames(model$bins$peakScores)

+

+    # Check if replicates have the same fdr value

+    reps <- sub('-rep.*', '', model$info$ID)

+    if (any(sapply(split(fdr, reps), function(x) { Reduce('&', x==x[1]) }) == FALSE)) {

+        stop("Replicates must have the same fdr value.")

+    }

+    

+    ## fdr threshold

+    threshold <- 1-fdr

+

+    ### Multivariate HMM ###

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

+        ## Calculate states

+        ptm <- startTimedMessage("Calculating states from peakScores ...")

+        p <- model$bins$peakScores

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

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

+            if (!invert) {

+                p.thresholded[,icol] <- p[,icol] >= threshold[icol]

+            } else {

+                p.thresholded[,icol] <- p[,icol] < threshold[icol] & p[,icol] > 0

+            }

+        }

+        states <- factor(bin2dec(p.thresholded), levels=levels(model$bins$state))

+        model$bins$state <- states

+        ## Combinations

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

+            mapping <- model$mapping

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

+        }

+        stopTimedMessage(ptm)

+        ## Redo segmentation

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

+        ## Redo peaks

+        model$peaks <- list()

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

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

+            peaks <- model$segments[mask]

+            mcols(peaks) <- NULL

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

+            model$peaks[[i1]] <- peaks

+        }

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

+        

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

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

+        mapping <- mapping.df$combination

+        names(mapping) <- mapping.df$state

+        p <- model$bins$peakScores

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

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

+            if (!invert) {

+                p.thresholded[,icol] <- p[,icol] >= threshold[icol]

+            } else {

+                p.thresholded[,icol] <- p[,icol] < threshold[icol] & p[,icol] > 0

+            }

+        }

+        states <- factor(bin2dec(p.thresholded), levels=mapping.df$state)

+        ## Make fake multiHMM

+        multiHMM <- list()

+        class(multiHMM) <- class.multivariate.hmm

+        multiHMM$info <- model$info

+        multiHMM$bins <- model$bins

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

+        multiHMM$bins$state <- factor(states)

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

+        multiHMM$bins$peakScores <- p

+        multiHMM$mapping <- mapping

+        multiHMM$peaks <- model$peaks

+        model <- combineMultivariates(list(multiHMM), mode='full')

+        ## Apply FDR on peaks ##

+        for (i1 in 1:length(model$peaks)) {

+            peaks <- model$peaks[[i1]]

+            if (!invert) {

+                model$peaks[[i1]] <- peaks[peaks$peakScores >= threshold[names(model$peaks)[i1]]]

+            } else {

+                model$peaks[[i1]] <- peaks[peaks$peakScores < threshold[names(model$peaks)[i1]]]

+            }

+        }

+    } else {

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

+    }

+    ## Return model

+    model$fdr <- fdr

+    return(model)

+    

+

+    

+}

+

+

+changeFDR.univariate <- function(model, fdr, invert=FALSE) {

+  

+    fdr <- suppressWarnings( as.numeric(fdr) )

+    if (!is.numeric(fdr) | is.na(fdr)) {

+        warning("Not changing false discovery rate because given 'fdr' is not numeric.")

+        return(model)

+    } else if (fdr < 0 | fdr > 1) {

+        warning("fdr has to be inside the interval [0,1]. Nothing done.")

+        return(model)

+    }

+

+    ## fdr threshold

+    threshold <- 1-fdr

+

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

+    ## Calculate states

+    ptm <- startTimedMessage("Calculating states from peakScores ...")

+    states <- model$bins$state

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

+    if (!invert) {

+        states[ model$bins$peakScores >= threshold ] <- 'modified'

+    } else {

+        states[ model$bins$peakScores < threshold & model$bins$peakScores > 0 ] <- 'modified'

+    }

+    model$bins$state <- states

+    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')))

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

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

+    ## Redo peaks

+    model$peaks <- model$segments[model$segments$state == 'modified']

+    model$peaks$state <- NULL

+    model$segments <- NULL

+    ## Redo weights

+    model$weights <- table(model$bins$state) / length(model$bins)

+    stopTimedMessage(ptm)

+

+    ## Return model

+    model$fdr <- fdr

+    return(model)

+  

+}

@@ -1,4 +1,4 @@

-#' Change the false discovery rate of a Hidden Markov Model

+#' Change the posterior cutoff of a Hidden Markov Model

 #'

 #' Adjusts the peak calls of a \code{\link{uniHMM}} or \code{\link{multiHMM}} object with the given posterior cutoff.

 #'

@@ -8,25 +8,24 @@

 #' @param model A \code{\link{uniHMM}} or \code{\link{multiHMM}} object with posteriors.

 #' @param post.cutoff A vector of posterior cutoff values between 0 and 1 the same length as \code{ncol(model$bins$posteriors)}. If only one value is given, it will be reused for all columns. Values close to 1 will yield more stringent peak calls with lower false positive but higher false negative rate.

 #' @return The input object is returned with adjusted peak calls.

-#' @export

-#' @examples

-#'## Get an example BAM file

-#'file <- system.file("extdata", "euratrans",

-#'                       "lv-H3K27me3-BN-male-bio2-tech1.bam",

-#'                        package="chromstaRData")

-#'## Bin the file into bin size 1000bp

-#'data(rn4_chrominfo)

-#'binned <- binReads(file, assembly=rn4_chrominfo, binsizes=1000,

-#'                   chromosomes='chr12')

-#'## Fit the univariate Hidden Markov Model

-#'# !Keep posteriors to change the post.cutoff later!

-#'hmm <- callPeaksUnivariate(binned, max.time=60, eps=1,

-#'                           keep.posteriors=TRUE)

-#'## Compare fits with different post.cutoffs

-#'plotHistogram(changePostCutoff(hmm, post.cutoff=0.01)) + ylim(0,0.25)

-#'plotHistogram(hmm) + ylim(0,0.25)

-#'plotHistogram(changePostCutoff(hmm, post.cutoff=0.99)) + ylim(0,0.25)

-#'

+# #' @examples

+# #'## Get an example BAM file

+# #'file <- system.file("extdata", "euratrans",

+# #'                       "lv-H3K27me3-BN-male-bio2-tech1.bam",

+# #'                        package="chromstaRData")

+# #'## Bin the file into bin size 1000bp

+# #'data(rn4_chrominfo)

+# #'binned <- binReads(file, assembly=rn4_chrominfo, binsizes=1000,

+# #'                   chromosomes='chr12')

+# #'## Fit the univariate Hidden Markov Model

+# #'# !Keep posteriors to change the post.cutoff later!

+# #'hmm <- callPeaksUnivariate(binned, max.time=60, eps=1,

+# #'                           keep.posteriors=TRUE)

+# #'## Compare fits with different post.cutoffs

+# #'plotHistogram(changePostCutoff(hmm, post.cutoff=0.01)) + ylim(0,0.25)

+# #'plotHistogram(hmm) + ylim(0,0.25)

+# #'plotHistogram(changePostCutoff(hmm, post.cutoff=0.99)) + ylim(0,0.25)

+# #'

 changePostCutoff <- function(model, post.cutoff=0.5) {

     if (!is.numeric(post.cutoff)) {

@@ -80,7 +79,7 @@ changePostCutoff.multivariate <- function(model, post.cutoff) {

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

             post.thresholded[,icol] <- post[,icol] >= threshold[icol]

         }

-        states <- factor(bin2dec(model$bins$posteriors >= threshold), levels=levels(model$bins$state))

+        states <- factor(bin2dec(post.thresholded), levels=levels(model$bins$state))

         model$bins$state <- states

         ## Combinations

         if (!is.null(model$bins$combination)) {

@@ -90,6 +90,10 @@ combineMultivariates <- function(hmms, mode) {

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

         posteriors <- list()

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

+        peakScores <- list()

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

+        peaks <- list()

+        peaks[[1]] <- hmm$peaks

         binstates <- list()

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

         stopTimedMessage(ptm)

@@ -102,6 +106,8 @@ combineMultivariates <- function(hmms, mode) {

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

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

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

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

+                peaks[[i1]] <- hmm$peaks

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

                 stopTimedMessage(ptm)

             }

@@ -109,12 +115,13 @@ combineMultivariates <- function(hmms, mode) {

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

         counts <- do.call(cbind, counts)

         posteriors <- do.call(cbind, posteriors)

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

         infos <- do.call(rbind, infos)

         conditions <- unique(infos$condition)

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

         names(states) <- NULL

         names(combs) <- conditions

-        combs.df <- as(combs,'DataFrame')

+        combs.df <- methods::as(combs,'DataFrame')

         stopTimedMessage(ptm)

     } else if (mode == 'differential') {

@@ -122,6 +129,8 @@ combineMultivariates <- function(hmms, mode) {

         infos <- list()

         counts <- list()

         posteriors <- list()

+        peakScores <- list()

+        peaks <- list()

         binstates <- list()

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

             ptm <- startTimedMessage("Processing HMM ",i1," ...")

@@ -129,12 +138,15 @@ combineMultivariates <- function(hmms, mode) {

             infos[[i1]] <- hmm$info

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

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

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

+            peaks[[i1]] <- hmm$peaks

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

             stopTimedMessage(ptm)

         }

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

         counts <- do.call(cbind, counts)

         posteriors <- do.call(cbind, posteriors)

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

         infos <- do.call(rbind, infos)

         conditions <- unique(infos$condition)

         binstates <- do.call(cbind, binstates)

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

             combs[[condition]] <- mapping[as.character(states.cond)]

         }

         combs.df <- as.data.frame(combs) # get factors instead of characters

-        combs.df <- as(combs.df, 'DataFrame')

+        combs.df <- methods::as(combs.df, 'DataFrame')

         stopTimedMessage(ptm)

     } else if (mode == 'full') {

@@ -175,6 +187,8 @@ combineMultivariates <- function(hmms, mode) {

         conditions <- unique(infos$condition)

         counts <- hmm$bins$counts

         posteriors <- hmm$bins$posteriors

+        peakScores <- hmm$bins$peakScores

+        peaks <- hmm$peaks

         states <- hmm$bins$state

         combs <- list()

         for (condition in conditions) {

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

             stopTimedMessage(ptm)

         }

         combs.df <- as.data.frame(combs) # get factors instead of characters

-        combs.df <- as(combs.df, 'DataFrame')

+        combs.df <- methods::as(combs.df, 'DataFrame')

     } else if (mode == 'replicate') {

         ## Load first HMM for coordinates

@@ -201,6 +215,10 @@ combineMultivariates <- function(hmms, mode) {

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

         posteriors <- list()

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

+        peakScores <- list()

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

+        peaks <- list()

+        peaks[[1]] <- hmm$peaks

         binstates <- list()

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

         stopTimedMessage(ptm)

@@ -212,6 +230,8 @@ combineMultivariates <- function(hmms, mode) {

                 infos[[i1]] <- hmm$info

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

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

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

+                peaks[[i1]] <- hmm$peaks

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

                 stopTimedMessage(ptm)

             }

@@ -219,6 +239,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)

         binstates <- do.call(cbind, binstates)

         infos <- do.call(rbind, infos)

         conditions <- unique(infos$condition)

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

             combs[[condition]] <- mapping[as.character(states)]

         }

         names(combs) <- conditions

-        combs.df <- as(combs,'DataFrame')

+        combs.df <- methods::as(combs,'DataFrame')

         stopTimedMessage(ptm)

     } else {

         stop("Unknown mode '", mode, "'.")

@@ -261,9 +282,10 @@ combineMultivariates <- function(hmms, mode) {

     ## Transferring counts and posteriors

     bins$counts <- counts

     bins$posteriors <- posteriors

+    bins$peakScores <- peakScores

     ## Add differential score ##

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

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

     ### Redo the segmentation for all conditions combined

     ptm <- startTimedMessage("Redoing segmentation for all conditions combined ...")

@@ -273,20 +295,23 @@ combineMultivariates <- function(hmms, mode) {

     ### Redo the segmentation for each condition separately

     ptm <- startTimedMessage("Redoing segmentation for each condition separately ...")

-    segments.separate <- list()

+    segments.per.condition <- list()

     for (cond in names(combs)) {

         bins.cond <- bins

         mcols(bins.cond) <- mcols(bins)[paste0('combination.',cond)]

         df <- as.data.frame(bins.cond)

         names(df)[6] <- cond

         segments.cond <- suppressMessages( collapseBins(df, column2collapseBy=cond, columns2drop=c('width', grep('posteriors', names(df), value=TRUE))) )

-        segments.cond <- as(segments.cond, 'GRanges')

+        segments.cond <- methods::as(segments.cond, 'GRanges')

         names(mcols(segments.cond)) <- 'combination'

         seqlengths(segments.cond) <- seqlengths(bins)[seqlevels(segments.cond)]

-        segments.separate[[cond]] <- segments.cond

+        segments.per.condition[[cond]] <- segments.cond

     }

     stopTimedMessage(ptm)

+    ### Flatten peak list ###

+    peaks <- unlist(peaks)

+    

     ### Make return object

     hmm <- list()

     class(hmm) <- class.combined.multivariate.hmm

@@ -294,7 +319,8 @@ combineMultivariates <- function(hmms, mode) {

     rownames(hmm$info) <- NULL

     hmm$bins <- bins

     hmm$segments <- segments

-    hmm$segments.separate <- segments.separate

+    hmm$segments.per.condition <- segments.per.condition

+    hmm$peaks <- peaks

     hmm$frequencies <- freqs$table

     return(hmm)

@@ -41,6 +41,19 @@ NULL

 NULL

+#' Gene coordinates for rn4

+#'

+#' A data.frame containing gene coordinates and biotypes of the rn4 assembly.

+#'

+#' @docType data

+#' @name genes_rn4

+#' @format A data.frame.

+#' @examples

+#'data(genes_rn4)

+#'head(genes_rn4)

+NULL

+

+

 #' chromstaR objects

 #'

 #' @description

@@ -136,7 +149,7 @@ NULL

 #' A \code{list()} with the following entries:

 #' \item{bins}{A \code{\link[GenomicRanges]{GRanges}} object containing genomic bin coordinates and human readable combinations for the combined \code{\link{multiHMM}} objects.}

 #' \item{segments}{Same as \code{bins}, but consecutive bins with the same state are collapsed into segments.}

-#' \item{segments.separate}{A \code{list} with segments for each condition separately.}

+#' \item{segments.per.condition}{A \code{list} with segments for each condition separately.}

 #' @seealso \code{\link{combineMultivariates}}, \code{\link{uniHMM}}, \code{\link{multiHMM}}

 #' @name combinedMultiHMM

 #' @aliases combinedHMM

@@ -18,7 +18,7 @@

 #' \describe{

 #'   \item{\code{\link{heatmapCountCorrelation}}}{Heatmap of read count correlations.}

 #'   \item{\code{\link{plotEnrichCountHeatmap}}}{Heatmap of read counts around annotation.}

-#'   \item{\code{\link{plotFoldEnrichment}}}{Enrichment of combinatorial states around annotation.}

+#'   \item{\code{\link{plotEnrichment}}}{Enrichment of combinatorial states around annotation.}

 #'   \item{\code{\link{plotFoldEnrichHeatmap}}}{Enrichment of combinatorial states at multiple annotations.}

 #'   \item{\code{\link{plotExpression}}}{Boxplot of expression values that overlap combinatorial states.}

 #' }

@@ -32,17 +32,17 @@

 #'### Make the enrichment plots ###

 #'# We expect promoter [H3K4me3] and bivalent-promoter signatures [H3K4me3+H3K27me3]

 #'# to be enriched at transcription start sites.

-#'    plotFoldEnrichment(hmm = model, annotation = genes, bp.around.annotation = 15000) +

+#'    plotEnrichment(hmm = model, annotation = genes, bp.around.annotation = 15000) +

 #'    ggtitle('Fold enrichment around genes') +

 #'    xlab('distance from gene body')

 #'  

 #'# Plot enrichment only at TSS. We make use of the fact that TSS is the start of a gene.

-#'    plotFoldEnrichment(model, genes, region = 'start') +

+#'    plotEnrichment(model, genes, region = 'start') +

 #'    ggtitle('Fold enrichment around TSS') +

 #'    xlab('distance from TSS in [bp]')

 #'# Note: If you want to facet the plot because you have many combinatorial states you

 #'# can do that with

-#'    plotFoldEnrichment(model, genes, region = 'start') +

+#'    plotEnrichment(model, genes, region = 'start') +

 #'    facet_wrap(~ combination)

 #'  

 #'# Another form of visualization that shows every TSS in a heatmap

@@ -65,6 +65,7 @@ NULL

 #' @param annotations A \code{list()} with \code{\link{GRanges}} objects containing coordinates of multiple annotations The names of the list entries will be used to name the return values.

 #' @param plot A logical indicating whether the plot or an array with the fold enrichment values is returned.

 #' @importFrom S4Vectors subjectHits queryHits

+#' @importFrom IRanges subsetByOverlaps

 #' @importFrom reshape2 melt

 #' @export

 plotFoldEnrichHeatmap <- function(hmm, annotations, what="combinations", combinations=NULL, marks=NULL, plot=TRUE) {

@@ -89,7 +90,8 @@ plotFoldEnrichHeatmap <- function(hmm, annotations, what="combinations", combina

         mark.levels <- marks

     }

     genome <- sum(as.numeric(width(bins)))

-    feature.lengths <- lapply(annotations, function(x) { sum(as.numeric(width(x))) })

+    annotationsAtBins <- lapply(annotations, function(x) { IRanges::subsetByOverlaps(x, bins) })

+    feature.lengths <- lapply(annotationsAtBins, function(x) { sum(as.numeric(width(x))) })

     if (what == 'peaks') {

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

@@ -98,16 +100,17 @@ plotFoldEnrichHeatmap <- function(hmm, annotations, what="combinations", combina

     ## Fold enrichment

     ggplts <- list()

     folds <- list()

+    maxfolds <- list()

     for (condition in conditions) {

         if (what == 'combinations') {

             bins$combination <- mcols(bins)[,paste0('combination.', condition)]

-            fold <- array(NA, dim=c(length(annotations), length(comb.levels)), dimnames=list(annotation=names(annotations), combination=comb.levels))

+            fold <- array(NA, dim=c(length(annotationsAtBins), length(comb.levels)), dimnames=list(annotation=names(annotationsAtBins), combination=comb.levels))

             for (icomb in 1:length(comb.levels)) {

                 mask <- bins$combination == comb.levels[icomb]

                 bins.mask <- bins[mask]

                 combstate.length <- sum(as.numeric(width(bins.mask)))

-                for (ifeat in 1:length(annotations)) {

-                    feature <- annotations[[ifeat]]

+                for (ifeat in 1:length(annotationsAtBins)) {

+                    feature <- annotationsAtBins[[ifeat]]

                     ind <- findOverlaps(bins.mask, feature)

                     binsinfeature <- bins.mask[unique(S4Vectors::queryHits(ind))]

@@ -121,7 +124,7 @@ plotFoldEnrichHeatmap <- function(hmm, annotations, what="combinations", combina

             }

         } else if (what == 'peaks') {

-            fold <- array(NA, dim=c(length(annotations), length(mark.levels)), dimnames=list(annotation=names(annotations), mark=mark.levels))

+            fold <- array(NA, dim=c(length(annotationsAtBins), length(mark.levels)), dimnames=list(annotation=names(annotationsAtBins), mark=mark.levels))

             for (imark in 1:length(mark.levels)) {

                 mark <- mark.levels[imark]

                 colmask <- hmm$info$mark == mark

@@ -133,8 +136,8 @@ plotFoldEnrichHeatmap <- function(hmm, annotations, what="combinations", combina

                 mask <- binstates.cond

                 bins.mask <- bins[mask]

                 combstate.length <- sum(as.numeric(width(bins.mask)))

-                for (ifeat in 1:length(annotations)) {

-                    feature <- annotations[[ifeat]]

+                for (ifeat in 1:length(annotationsAtBins)) {

+                    feature <- annotationsAtBins[[ifeat]]

                     ind <- findOverlaps(bins.mask, feature)

                     binsinfeature <- bins.mask[unique(S4Vectors::queryHits(ind))]

@@ -151,16 +154,22 @@ plotFoldEnrichHeatmap <- function(hmm, annotations, what="combinations", combina

         if (plot) {

             df <- reshape2::melt(fold, value.name='foldEnrichment')

+            maxfolds[[condition]] <- max(df$foldEnrichment, na.rm=TRUE)

             if (what == 'combinations') {

-                ggplt <- ggplot(df) + geom_tile(aes_string(x='combination', y='annotation', fill='foldEnrichment')) + theme(axis.text.x = element_text(angle=90, hjust=1)) + scale_fill_gradient(low='white', high='blue')

+                ggplt <- ggplot(df) + geom_tile(aes_string(x='combination', y='annotation', fill='foldEnrichment'))

             } else if (what == 'peaks') {

-                ggplt <- ggplot(df) + geom_tile(aes_string(x='mark', y='annotation', fill='foldEnrichment')) + theme(axis.text.x = element_text(angle=90, hjust=1)) + scale_fill_gradient(low='white', high='blue')

+                ggplt <- ggplot(df) + geom_tile(aes_string(x='mark', y='annotation', fill='foldEnrichment'))

             }

+            ggplt <- ggplt + theme(axis.text.x = element_text(angle=90, hjust=1, vjust=0.5))

             ggplts[[condition]] <- ggplt

         } else {

             folds[[condition]] <- fold

         }

     }

+    if (plot) {

+        maxfolds <- unlist(maxfolds)

+        ggplts <- lapply(ggplts, function(ggplt) { ggplt + scale_fill_gradientn(colors=grDevices::colorRampPalette(c("blue","white","red"))(20), values=c(seq(0,1,length.out=10), seq(1,max(maxfolds, na.rm=TRUE),length.out=10)), rescaler=function(x,...) {x}, oob=identity, limits=c(0,max(maxfolds, na.rm=TRUE))) })

+    }

     if (class(hmm) == class.multivariate.hmm) {

         if (plot) {

@@ -183,6 +192,7 @@ plotFoldEnrichHeatmap <- function(hmm, annotations, what="combinations", combina

 #' @inheritParams enrichmentAtAnnotation

 #' @param max.rows An integer specifying the number of randomly subsampled rows that are plotted from the \code{annotation} object. This is necessary to avoid crashing for heatmaps with too many rows.

 #' @importFrom reshape2 melt

+#' @importFrom IRanges subsetByOverlaps

 #' @export

 plotEnrichCountHeatmap <- function(hmm, annotation, bp.around.annotation=10000, max.rows=1000) {

@@ -190,24 +200,27 @@ plotEnrichCountHeatmap <- function(hmm, annotation, bp.around.annotation=10000,

     ## Variables

     bins <- hmm$bins

     if (class(hmm) == class.combined.multivariate.hmm) {

+        conditions <- sub('combination.', '', grep('combination', names(mcols(bins)), value=TRUE))

+        comb.levels <- levels(mcols(bins)[,paste0('combination.', conditions[1])])

+        ## Create new column combination with all conditions combined

+        combinations <- list()

+        for (condition in conditions) {

+            combinations[[condition]] <- paste0(condition, ":", mcols(bins)[,paste0('combination.', condition)])

+        }

+        combinations$sep <- ', '

+        bins$combination <- factor(do.call(paste, combinations))

     } else if (class(hmm) == class.multivariate.hmm) {

-        # Rename 'combination' to 'combination.' for coherence with combinedMultiHMM

-        names(mcols(bins))[grep('combination', names(mcols(bins)))] <- 'combination.'

+        comb.levels <- levels(bins$combination)

     }

-    conditions <- sub('combination.', '', grep('combination', names(mcols(bins)), value=TRUE))

-    comb.levels <- levels(mcols(bins)[,paste0('combination.', conditions[1])])

     binsize <- width(bins)[1]

     around <- round(bp.around.annotation/binsize)

-    ## Create new column combination with all conditions combined

-    combinations <- list()

-    for (condition in conditions) {

-        combinations[[condition]] <- paste0(condition, ":", mcols(bins)[,paste0('combination.', condition)])

-    }

-    combinations$sep <- ', '

-    bins$combination <- factor(do.call(paste, combinations))

+    

+    ## Get RPKM values

+    bins$counts <- sweep(bins$counts, MARGIN = 2, STATS = colSums(bins$counts), FUN = '/')

+    bins$counts <- bins$counts * 1e6 * 1000/mean(width(bins))