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

 Package: chromstaR

 Type: Package

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

-Version: 1.1.1

+Version: 1.1.2

 Date: 2016-10-19

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

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

@@ -25,7 +25,8 @@ Imports:

     reshape2,

     Rsamtools,

     GenomicAlignments,

-    bamsignals

+    bamsignals,

+    mvtnorm

 Suggests:

     knitr,

     BiocStyle,

@@ -36,5 +37,5 @@ License: Artistic-2.0

 LazyLoad: yes

 VignetteBuilder: knitr

 Packaged: 2014-05-03 00:00:00 CET+1; Taudt

-biocViews: Software, DifferentialPeakCalling, HiddenMarkovModel, ChIPSeq, MultipleComparison

+biocViews: Software, DifferentialPeakCalling, HiddenMarkovModel, ChIPSeq, HistoneModification, MultipleComparison

 RoxygenNote: 5.0.1

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

 # Generated by roxygen2: do not edit by hand

+S3method(print,combinedMultiHMM)

+S3method(print,multiHMM)

+S3method(print,uniHMM)

 export(Chromstar)

 export(bin2dec)

 export(binReads)

@@ -7,6 +10,7 @@ export(callPeaksMultivariate)

 export(callPeaksReplicates)

 export(callPeaksUnivariate)

 export(changeFDR)

+export(changePostCutoff)

 export(collapseBins)

 export(combineMultivariates)

 export(dec2bin)

@@ -31,6 +35,7 @@ export(plotHistogram)

 export(readBamFileAsGRanges)

 export(readBedFileAsGRanges)

 export(readCustomBedFile)

+export(removeCondition)

 export(stateBrewer)

 export(transitionFrequencies)

 export(unis2pseudomulti)

@@ -57,9 +62,11 @@ importFrom(S4Vectors,runmean)

 importFrom(S4Vectors,subjectHits)

 importFrom(bamsignals,bamCount)

 importFrom(grDevices,col2rgb)

+importFrom(grDevices,rgb2hsv)

 importFrom(graphics,hist)

 importFrom(methods,as)

 importFrom(methods,is)

+importFrom(mvtnorm,rmvnorm)

 importFrom(reshape2,melt)

 importFrom(stats,aggregate)

 importFrom(stats,cutree)

@@ -1,3 +1,15 @@

+CHANGES IN VERSION 1.1.2

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

+

+NEW FEATURES

+

+    o Proper print() methods for all objects.

+

+BUG FIXES

+

+    o Fixed a bug where chromosomes with a single bin where making problems.

+

+

 CHANGES IN VERSION 1.1.1

 ------------------------

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

-#' Wrapper function for the \code{\link{chromstaR}} package

+#' Wrapper function for the \pkg{\link{chromstaR}} package

 #' 

 #' This function performs \code{\link[chromstaR:binReads]{binning}}, \code{\link[chromstaR:callPeaksUnivariate]{univariate peak calling}} and \code{\link[chromstaR:callPeaksMultivariate]{multivariate peak calling}} from a list of input files.

 #' 

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

 #' Convert aligned reads from .bam or .bed(.gz) files into read counts in equidistant windows (bins). This function uses \code{\link[GenomicRanges]{countOverlaps}} to calculate the read counts, or alternatively \code{\link[bamsignals]{bamProfile}} if option \code{use.bamsignals} is set (only effective for .bam files).

 #'

 #' @aliases binning

-#' @param file A file with aligned reads. Alternatively a \code{\link{GRanges}} with aligned reads if format is set to 'GRanges'.

+#' @param file A file with aligned reads. Alternatively a \code{\link{GRanges}} with aligned reads.

 #' @param experiment.table An \code{\link{experiment.table}} containing the supplied \code{file}. This is necessary to uniquely identify the file in later steps of the workflow. Set to \code{NULL} if you don't have it (not recommended).

 #' @inheritParams readBamFileAsGRanges

 #' @inheritParams readBedFileAsGRanges

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

 #' @param variable.width.reference A BAM file that is used as reference to produce variable width bins. See \code{\link{variableWidthBins}} for details.

 #' @param chromosomes If only a subset of the chromosomes should be binned, specify them here.

 #' @param use.bamsignals If \code{TRUE} the \pkg{\link[bamsignals]{bamsignals}} package is used for parsing of BAM files. This gives tremendous speed advantage for only one binsize but linearly increases for multiple binsizes, while \code{use.bamsignals=FALSE} has a binsize dependent runtime and might be faster if many binsizes are calculated.

+#' @param format One of \code{c('bed','bam','GRanges')} or \code{NULL} if the format is to be determined automatically.

 #' @return If only one bin size was specified for option \code{binsizes}, the function returns a single \code{\link{GRanges}} object with meta data column 'counts' that contains the read count. If multiple \code{binsizes} were specified , the function returns a named \code{list()} of \link{GRanges} objects.

 #' @importFrom Rsamtools BamFile indexBam

 #' @importFrom bamsignals bamCount

@@ -35,14 +36,21 @@

 #'                   chromosomes='chr12')

 #'print(binned)

 #'

-binReads <- function(file, experiment.table=NULL, assembly, bamindex=file, chromosomes=NULL, pairedEndReads=FALSE, min.mapq=10, remove.duplicate.reads=TRUE, max.fragment.width=1000, blacklist=NULL, binsizes=1000, reads.per.bin=NULL, bins=NULL, variable.width.reference=NULL, use.bamsignals=TRUE) {

+binReads <- function(file, experiment.table=NULL, assembly, bamindex=file, chromosomes=NULL, pairedEndReads=FALSE, min.mapq=10, remove.duplicate.reads=TRUE, max.fragment.width=1000, blacklist=NULL, binsizes=1000, reads.per.bin=NULL, bins=NULL, variable.width.reference=NULL, use.bamsignals=TRUE, format=NULL) {

     ## Determine format

-    if (is.character(file)) {

-        file.clean <- sub('\\.gz$','', file)

-        format <- rev(strsplit(file.clean, '\\.')[[1]])[1]

-    } else if (class(file)=='GRanges') {

-        format <- 'GRanges'

+    if (is.null(format)) {

+        if (is.character(file)) {

+            file.clean <- sub('\\.gz$','', file)

+            format <- rev(strsplit(file.clean, '\\.')[[1]])[1]

+        } else if (class(file)=='GRanges') {

+            format <- 'GRanges'

+        } else {

+            stop("Could not determine format automatically. Please specify it via the 'format' parameter.")

+        }

+    }

+    if (! format %in% c('bed','bam','GRanges')) {

+        stop("Unknown format. Needs to be one of 'bed', 'bam' or 'GRanges'.")

     }

     if (format=='bed') {

@@ -111,6 +119,8 @@ binReads <- function(file, experiment.table=NULL, assembly, bamindex=file, chrom

         ## GRanges (1-based)

         data <- file

         chrom.lengths <- seqlengths(data)

+    } else {

+        stop("Unknown file format.")

     }

     ## Select chromosomes to bin

@@ -78,6 +78,11 @@ callPeaksMultivariate <- function(hmms, use.states, max.states=NULL, per.chrom=T

         ### Make return object ###

             result <- list()

             result$info <- hmm$info

+            if (is.null(result$info)) {

+                n <- 1

+                result$info <- data.frame(file=rep(NA, n), mark=1:n, condition=1:n, replicate=1, pairedEndReads=rep(NA, n), controlFiles=rep(NA, n))

+                result$info$ID <- paste0(result$info$mark, '-', result$info$condition, '-rep', result$info$replicate)

+            }

         ## Bin coordinates, posteriors and states

             result$bins <- hmm$bins

             result$bins$score <- NULL

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

     p <- prepareMultivariate(hmms, use.states=use.states, max.states=max.states, chromosomes=chromosomes)

     if (is.null(chromosomes)) {

-        chromosomes <- seqlevels(p$bins)

+        chromosomes <- intersect(seqlevels(p$bins), unique(seqnames(p$bins)))

     }

     ## Run multivariate per chromosome

@@ -167,7 +172,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, 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)

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

                 stopTimedMessage(ptm)

                 models[[chrom]] <- model

@@ -396,6 +401,11 @@ prepareMultivariate = function(hmms, use.states=NULL, max.states=NULL, chromosom

     }

     info <- do.call(rbind, info)

     rownames(info) <- NULL

+    if (is.null(info)) {

+        n <- nummod

+        info <- data.frame(file=rep(NA, n), mark=1:n, condition=1:n, replicate=rep(1, n), pairedEndReads=rep(NA, n), controlFiles=rep(NA, n))

+        info$ID <- paste0(info$mark, '-', info$condition, '-rep', info$replicate)

+    }

     bins$counts <- counts

     colnames(bins$counts) <- info$ID

     maxcounts <- max(bins$counts)

@@ -507,7 +517,8 @@ prepareMultivariate = function(hmms, use.states=NULL, max.states=NULL, chromosom

         }

         temp <- tryCatch({

             if (nrow(z.temp) > 100) {

-                correlationMatrix[,,istate] <- cor(z.temp)

+                correlationMatrix[,,istate] <- suppressWarnings( cor(z.temp) )

+                correlationMatrix[,,istate][is.na(correlationMatrix[,,istate])] <- 0

                 determinant[istate] <- det( correlationMatrix[,,istate] )

                 correlationMatrixInverse[,,istate] <- chol2inv(chol(correlationMatrix[,,istate]))

             } else {

@@ -516,12 +527,10 @@ prepareMultivariate = function(hmms, use.states=NULL, max.states=NULL, chromosom

                 correlationMatrixInverse[,,istate] <- diag(nummod) # solve(diag(x)) == diag(x)

             }

             0

-        }, warning = function(war) {

-            1

         }, error = function(err) {

-            1

+            2

         })

-        if (temp!=0) {

+        if (temp==2) {

             correlationMatrix[,,istate] <- diag(nummod)

             determinant[istate] <- 1

             correlationMatrixInverse[,,istate] <- diag(nummod)

@@ -572,3 +581,13 @@ prepareMultivariate = function(hmms, use.states=NULL, max.states=NULL, chromosom

     return(out)

 }

+

+### Get real transition probabilities ###

+transProbs <- function(model) {

+    bins <- model$bins

+    df <- data.frame(from = bins$combination[-length(bins)], to = bins$combination[-1])

+    t <- table(df)

+    t <- t[rownames(model$transitionProbs), colnames(model$transitionProbs)]

+    t <- sweep(t, MARGIN = 1, STATS = rowSums(t), FUN = '/')

+    return(t)

+}

@@ -64,8 +64,8 @@ callPeaksReplicates <- function(hmm.list, max.states=32, force.equal=FALSE, eps=

         ## Univariate replicateInfo

         ids <- sapply(hmms, function(x) { x$info$ID })

-        weight.univariate <- unlist(lapply(hmms, function(x) { x$weights['modified'] }))

-        total.count <- unlist(lapply(hmms, function(x) { sum(x$bins$counts) }))

+        weight.univariate <- sapply(hmms, function(x) { x$weights['modified'] })

+        total.count <- sapply(hmms, function(x) { sum(x$bins$counts) })

         info.df <- data.frame(total.count=total.count, weight.univariate=weight.univariate)

         if (!is.null(unlist(ids))) {

             rownames(info.df) <- ids

@@ -57,7 +57,7 @@ callPeaksUnivariate <- function(binned.data, input.data=NULL, prefit.on.chr=NULL

     }

     if (!is.null(input.data)) {

         if (class(input.data) == 'character') { 

-            message("Loading input file ",input.data)

+            message("Loading input file(s) ", paste0(input.data, collapse=', '))

             input.datas <- loadHmmsFromFiles(input.data)

             input.data <- input.datas[[1]]

             if (length(input.datas) > 1) {

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

 #' @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.

-# #' @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)

-# #'

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

+#'

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

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

@@ -108,10 +109,27 @@ changePostCutoff.multivariate <- function(model, post.cutoff) {

         multiHMM$bins$state <- factor(states)

         multiHMM$bins$posteriors <- post

         multiHMM$mapping <- mapping

+        multiHMM$peaks <- model$peaks

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

     } else {

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

     }

+    

+    ## Redo peaks

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

+    binstates <- dec2bin(model$segments$state, colnames=colnames(post))

+    segments <- model$segments

+    mcols(segments) <- NULL

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

+        segments$peakScores <- model$segments$peakScores[,icol]

+        segments$state <- binstates[,icol]

+        red.df <- suppressMessages( collapseBins(as.data.frame(segments), column2collapseBy = 'state') )

+        red.df <- red.df[red.df$state == 1,]

+        red.df$state <- NULL

+        model$peaks[[icol]] <- methods::as(red.df, 'GRanges')

+    }

+    stopTimedMessage(ptm)

+    

     ## Return model

     model$post.cutoff <- threshold

     return(model)

@@ -150,10 +168,11 @@ changePostCutoff.univariate <- function(model, post.cutoff) {

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

-    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[,'mean.score'])

-    model$segments <- red.gr

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

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

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

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

+    model$peaks$state <- NULL

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

     stopTimedMessage(ptm)

     ## Return model

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

             stopTimedMessage(ptm)

         }

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

+        # Slightly more complicated selection procedure for conditions in case one mark is missing

+        conds.help <- lapply(infos, function(x) { unique(x$condition) })

+        conditions <- conds.help[[which.max(sapply(conds.help, length))]]

+        infos <- do.call(rbind, infos)

+        infos$condition <- factor(infos$condition, levels=conditions)

+        infos <- infos[order(infos$condition, infos$mark, infos$replicate),]

+        infos$condition <- as.character(infos$condition)

+        # Reorder everything according to conditions

         counts <- do.call(cbind, counts)

+        counts <- counts[,infos$ID]

         posteriors <- do.call(cbind, posteriors)

+        posteriors <- posteriors[,infos$ID]

         peakScores <- do.call(cbind, peakScores)

-        infos <- do.call(rbind, infos)

-        conditions <- unique(infos$condition)

+        peakScores <- peakScores[,infos$ID]

         binstates <- do.call(cbind, binstates)

+        binstates <- binstates[,infos$ID]

         states <- factor(bin2dec(binstates))

         names(states) <- NULL

         stopTimedMessage(ptm)

@@ -3,8 +3,9 @@

 #' Plotting functions for enrichment analysis of \code{\link{multiHMM}} or \code{\link{combinedMultiHMM}} objects with any annotation of interest, specified as a \code{\link[GenomicRanges]{GRanges}} object.

 #' 

 #' @name enrichment_analysis

-#' @param combinations A vector with combinations for which the fold enrichment will be calculated. If \code{NULL} all combinations will be considered.

+#' @param combinations A vector with combinations for which the enrichment will be calculated. If \code{NULL} all combinations will be considered.

 #' @param marks A vector with marks for which the enrichment is plotted. If \code{NULL} all marks will be considered.

+#' @param logscale Set to \code{TRUE} to plot fold enrichment on log-scale. Ignored if \code{statistic = 'fraction'}.

 #' @return A \code{\link[ggplot2:ggplot]{ggplot}} object containing the plot or a list() with \code{\link[ggplot2:ggplot]{ggplot}} objects if several plots are returned. For \code{plotFoldEnrichHeatmap} a named array with fold enrichments if \code{plot=FALSE}.

 #' @author Aaron Taudt

 #' @seealso \code{\link{plotting}}

@@ -54,6 +55,7 @@

 #'# Fold enrichment with different biotypes, showing that protein coding genes are

 #'# enriched with (bivalent) promoter combinations [H3K4me3] and [H3K4me3+H3K27me3],

 #'# while rRNA is enriched with the empty [] and repressive combinations [H3K27me3].

+#'    tss <- resize(genes, width = 3, fix = 'start')

 #'    biotypes <- split(tss, tss$biotype)

 #'    plotFoldEnrichHeatmap(model, annotations=biotypes) + coord_flip()

 #'

@@ -64,11 +66,12 @@ NULL

 #' @param hmm A \code{\link{combinedMultiHMM}} or \code{\link{multiHMM}} object or a file that contains such an object.

 #' @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.

+#' @param what One of \code{c('combinations','peaks','counts','transitions')} specifying on which feature the statistic is calculated.

 #' @importFrom S4Vectors subjectHits queryHits

 #' @importFrom IRanges subsetByOverlaps

 #' @importFrom reshape2 melt

 #' @export

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

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

     hmm <- loadHmmsFromFiles(hmm, check.class=c(class.multivariate.hmm, class.combined.multivariate.hmm))[[1]]

     ## Variables

@@ -76,9 +79,17 @@ plotFoldEnrichHeatmap <- function(hmm, annotations, what="combinations", combina

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

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

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

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

+        names(mcols(bins))[grep('combination', names(mcols(bins)))] <- paste0('combination.', unique(hmm$info$condition))

     }

     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

+    combs <- list()

+    for (condition in conditions) {

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

+    }

+    combs$sep <- ', '

+    bins$transitions <- factor(do.call(paste, combs))

     if (is.null(combinations)) {

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

     } else {

@@ -93,85 +104,137 @@ plotFoldEnrichHeatmap <- function(hmm, annotations, what="combinations", combina

     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)

-    }

-    

-    ## 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(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(annotationsAtBins)) {

-                    feature <- annotationsAtBins[[ifeat]]

-                    ind <- findOverlaps(bins.mask, feature)

+    minfolds <- list()

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

-                    sum.binsinfeature <- sum(as.numeric(width(binsinfeature)))

-    

-                    featuresinbins <- feature[unique(S4Vectors::subjectHits(ind))]

-                    sum.featuresinbins <- sum(as.numeric(width(featuresinbins)))

-    

-                    fold[ifeat,icomb] <- sum.binsinfeature / combstate.length / feature.lengths[[ifeat]] * genome

-                }

-            }

-    

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

-            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

-                colmask <- colmask & (!duplicated(paste0(hmm$info$mark, hmm$info$condition))) # remove replicates

-                if (condition != "") {

-                    colmask <- colmask & (hmm$info$condition == condition)

-                }

-                binstates.cond <- binstates[,colmask]

-                mask <- binstates.cond

-                bins.mask <- bins[mask]

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

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

-                    feature <- annotationsAtBins[[ifeat]]

-                    ind <- findOverlaps(bins.mask, feature)

-    

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

-                    sum.binsinfeature <- sum(as.numeric(width(binsinfeature)))

-    

-                    featuresinbins <- feature[unique(S4Vectors::subjectHits(ind))]

-                    sum.featuresinbins <- sum(as.numeric(width(featuresinbins)))

-    

-                    fold[ifeat,imark] <- sum.binsinfeature / combstate.length / feature.lengths[[ifeat]] * genome

-                }

+    if (what == 'peaks') {

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

+    }

+    if (what == 'transitions') {

+        bins$combination <- bins$transitions

+        comb.levels <- names(sort(table(bins$combination), decreasing = TRUE))

+        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(annotationsAtBins)) {

+                feature <- annotationsAtBins[[ifeat]]

+                ind <- findOverlaps(bins.mask, feature)

+

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

+                sum.binsinfeature <- sum(as.numeric(width(binsinfeature)))

+

+                featuresinbins <- feature[unique(S4Vectors::subjectHits(ind))]

+                sum.featuresinbins <- sum(as.numeric(width(featuresinbins)))

+

+                fold[ifeat,icomb] <- sum.binsinfeature / combstate.length / feature.lengths[[ifeat]] * genome

             }

-          

         }

-      

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

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

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

+            if (logscale) {

+                df$foldEnrichment <- log(df$foldEnrichment)

             }

+            foldsnoinf <- setdiff(df$foldEnrichment, c(Inf, -Inf))

+            maxfolds[[1]] <- max(foldsnoinf, na.rm=TRUE)

+            minfolds[[1]] <- min(foldsnoinf, na.rm=TRUE)

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

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

-            ggplts[[condition]] <- ggplt

+            ggplts[[1]] <- ggplt

         } else {

-            folds[[condition]] <- fold

+            folds[[1]] <- fold

+        }

+    }

+    

+    ## Fold enrichment

+    if (what %in% c('peaks', 'combinations')) {

+        for (condition in conditions) {

+            if (what == 'combinations') {

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

+                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(annotationsAtBins)) {

+                        feature <- annotationsAtBins[[ifeat]]

+                        ind <- findOverlaps(bins.mask, feature)

+        

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

+                        sum.binsinfeature <- sum(as.numeric(width(binsinfeature)))

+        

+                        featuresinbins <- feature[unique(S4Vectors::subjectHits(ind))]

+                        sum.featuresinbins <- sum(as.numeric(width(featuresinbins)))

+        

+                        fold[ifeat,icomb] <- sum.binsinfeature / combstate.length / feature.lengths[[ifeat]] * genome

+                    }

+                }

+        

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

+                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

+                    colmask <- colmask & (!duplicated(paste0(hmm$info$mark, hmm$info$condition))) # remove replicates

+                    if (condition != "") {

+                        colmask <- colmask & (hmm$info$condition == condition)

+                    }

+                    binstates.cond <- binstates[,colmask]

+                    mask <- binstates.cond

+                    bins.mask <- bins[mask]

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

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

+                        feature <- annotationsAtBins[[ifeat]]

+                        ind <- findOverlaps(bins.mask, feature)

+        

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

+                        sum.binsinfeature <- sum(as.numeric(width(binsinfeature)))

+        

+                        featuresinbins <- feature[unique(S4Vectors::subjectHits(ind))]

+                        sum.featuresinbins <- sum(as.numeric(width(featuresinbins)))

+        

+                        fold[ifeat,imark] <- sum.binsinfeature / combstate.length / feature.lengths[[ifeat]] * genome

+                    }

+                }

+              

+            }

+          

+            if (plot) {

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

+                if (logscale) {

+                    df$foldEnrichment <- log(df$foldEnrichment)

+                }

+                foldsnoinf <- setdiff(df$foldEnrichment, c(Inf, -Inf))

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

+                minfolds[[condition]] <- min(foldsnoinf, na.rm=TRUE)

+                if (what == 'combinations') {

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

+                }

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

+        maxfold <- max(unlist(maxfolds), na.rm=TRUE)

+        minfold <- min(unlist(minfolds), na.rm=TRUE)

+        limits <- max(abs(maxfold), abs(minfold))

+        if (logscale) {

+            ggplts <- lapply(ggplts, function(ggplt) { ggplt + scale_fill_gradientn(name='log(observed/expected)', colors=grDevices::colorRampPalette(c("blue","white","red"))(20), values=c(seq(-limits,0,length.out=10), seq(0,limits,length.out=10)), rescaler=function(x,...) {x}, oob=identity, limits=c(-limits, limits)) })

+            ggplts <- lapply(ggplts, function(ggplt) { ggplt$data$foldEnrichment[ggplt$data$foldEnrichment == -Inf] <- -limits; ggplt })

+        } else {

+            ggplts <- lapply(ggplts, function(ggplt) { ggplt + scale_fill_gradientn(name='observed/expected', colors=grDevices::colorRampPalette(c("blue","white","red"))(20), values=c(seq(0,1,length.out=10), seq(1,maxfold,length.out=10)), rescaler=function(x,...) {x}, oob=identity, limits=c(0,maxfold)) })

+        }

     }

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

+    if (class(hmm) == class.multivariate.hmm | what == 'transitions') {

         if (plot) {

             return(ggplts[[1]])

         } else {

@@ -203,12 +266,12 @@ plotEnrichCountHeatmap <- function(hmm, annotation, bp.around.annotation=10000,

         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()

+        combs <- list()

         for (condition in conditions) {

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

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

         }

-        combinations$sep <- ', '

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

+        combs$sep <- ', '

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

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

         comb.levels <- levels(bins$combination)

     }

@@ -254,18 +317,19 @@ plotEnrichCountHeatmap <- function(hmm, annotation, bp.around.annotation=10000,

     rownames(ext.index) <- 1:nrow(ext.index)

     stopTimedMessage(ptm)

-    ## Go through combinations and then tracks to get the read counts

+    ## Go through combinations and then IDs to get the read counts

     ptm <- startTimedMessage("Getting read counts")

     counts <- list()

-    for (combination in levels(bins$combination)) {

+    combinations <- names(sort(table(bins$combination[index]), decreasing = TRUE))

+    for (combination in combinations) {

         counts[[combination]] <- list()

         index.combination <- which(bins$combination[index]==combination)

         ext.index.combination <- ext.index[index.combination,]

         if (is.null(dim(ext.index.combination))) {

             ext.index.combination <- array(ext.index.combination, dim=c(1,dim(ext.index)[[2]]), dimnames=list(anno=rownames(ext.index)[index.combination], position=dimnames(ext.index)[[2]]))

         }

-        for (ntrack in colnames(bins$counts)) {

-            counts[[combination]][[ntrack]] <- array(bins$counts[ext.index.combination,ntrack], dim=dim(ext.index.combination), dimnames=dimnames(ext.index.combination))

+        for (nID in colnames(bins$counts)) {

+            counts[[combination]][[nID]] <- array(bins$counts[ext.index.combination,nID], dim=dim(ext.index.combination), dimnames=dimnames(ext.index.combination))

         }

     }

     stopTimedMessage(ptm)

@@ -287,16 +351,16 @@ plotEnrichCountHeatmap <- function(hmm, annotation, bp.around.annotation=10000,

     comb2keep <- names(num.comb)[num.comb/sum(num.comb) > 0.005]

     counts <- counts[comb2keep]

     df <- reshape2::melt(counts)

-    names(df) <- c('id','position','RPKM','track','combination')

+    names(df) <- c('id','position','RPKM','ID','combination')

     df$id <- factor(df$id, levels=rev(unique(df$id)))

     df$combination <- factor(df$combination, levels=unique(df$combination))

-    df$track <- factor(df$track, levels=colnames(bins$counts))

+    df$ID <- factor(df$ID, levels=hmm$info$ID)

     ## Plot as heatmap

     ggplt <- ggplot(df) + geom_tile(aes_string(x='position', y='id', color='combination'))

     ggplt <- ggplt + scale_color_manual(values = getDistinctColors(length(unique(df$combination))))

     ggplt <- ggplt + geom_tile(aes_string(x='position', y='id', fill='RPKM'), alpha=0.6)

-    ggplt <- ggplt + facet_wrap( ~ track, nrow=1) + custom_theme

+    ggplt <- ggplt + facet_wrap( ~ ID, nrow=1) + custom_theme

     ggplt <- ggplt + xlab('distance from annotation in [bp]') + ylab('')

     ggplt <- ggplt + scale_fill_continuous(trans='log1p', low='white', high='black')

     # Insert horizontal lines

@@ -315,7 +379,7 @@ plotEnrichCountHeatmap <- function(hmm, annotation, bp.around.annotation=10000,

 #' @inheritParams enrichmentAtAnnotation

 #' @importFrom reshape2 melt

 #' @export

-plotEnrichment <- function(hmm, annotation, bp.around.annotation=10000, region=c("start","inside","end"), num.intervals=20, what='combinations', combinations=NULL, marks=NULL, statistic='fold') {

+plotEnrichment <- function(hmm, annotation, bp.around.annotation=10000, region=c("start","inside","end"), num.intervals=20, what='combinations', combinations=NULL, marks=NULL, statistic='fold', logscale=TRUE) {

     ## Check user input

     if ((!what %in% c('combinations','peaks','counts')) | length(what) > 1) {

@@ -326,11 +390,17 @@ plotEnrichment <- function(hmm, annotation, bp.around.annotation=10000, region=c

     }

     ## Variables

-    hmm <- loadHmmsFromFiles(hmm, check.class=c(class.multivariate.hmm, class.combined.multivariate.hmm))[[1]]

-    hmm$bins$counts <- rpkm.matrix(hmm$bins$counts, binsize=mean(width(hmm$bins)))

+    hmm <- loadHmmsFromFiles(hmm, check.class=c(class.univariate.hmm, class.multivariate.hmm, class.combined.multivariate.hmm))[[1]]

     bins <- hmm$bins

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

+    if (class(hmm) == class.univariate.hmm) {

+        bins$counts <- rpkm.vector(hmm$bins$counts, binsize=mean(width(hmm$bins)))

+        mcols(bins)['combination.'] <- bins$state

+        bins$state <- c('zero-inflation' = 0, 'unmodified' = 0, 'modified' = 1)[bins$state]

+        hmm$info <- data.frame(file=NA, mark=1, condition=1, replicate=1, pairedEndReads=NA, controlFiles=NA, ID='1-1-rep1')

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

+        bins$counts <- rpkm.matrix(hmm$bins$counts, binsize=mean(width(hmm$bins)))

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

+        bins$counts <- rpkm.matrix(hmm$bins$counts, binsize=mean(width(hmm$bins)))

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

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

     }

@@ -348,10 +418,12 @@ plotEnrichment <- function(hmm, annotation, bp.around.annotation=10000, region=c

     if (what %in% c('peaks','counts')) {

         ### Get fold enrichment

-        enrich <- enrichmentAtAnnotation(hmm$bins, hmm$info, annotation, bp.around.annotation=bp.around.annotation, region=region, what=what, num.intervals=num.intervals, statistic=statistic)

+        enrich <- enrichmentAtAnnotation(bins, hmm$info, annotation, bp.around.annotation=bp.around.annotation, region=region, what=what, num.intervals=num.intervals, statistic=statistic)

     }

     ggplts <- list()

     maxfolds <- list()

+    minfolds <- list()

+    maxcols <- list()

     for (condition in conditions) {

         if (what == 'combinations') {

             ### Get fold enrichment

@@ -382,37 +454,59 @@ plotEnrichment <- function(hmm, annotation, bp.around.annotation=10000, region=c

         df$position[df$L1 == 'inside'] <- df$position[df$L1 == 'inside'] * bp.around.annotation

         if (what == 'combinations') {

             df <- df[df$combination %in% comb.levels,]

+            df$combination <- factor(df$combination, levels=comb.levels)

         } else if (what %in% c('peaks','counts')) {

-            df$mark <- sub("-.*", "", df$track)

+            df$mark <- sub("-.*", "", df$ID)

             df <- df[df$mark %in% mark.levels, ]

-            df$condition <- sapply(strsplit(as.character(df$track), '-'), '[[', 2)

+            df$condition <- sapply(strsplit(as.character(df$ID), '-'), '[[', 2)

             if (condition != "") {

                 df <- df[df$condition == condition, ]

             }

         }

+        if (logscale & !(statistic=='fraction')) {

+            df$value <- log(df$value)

+        }

         ### Plot

         if (what == 'combinations') {

             ggplt <- ggplot(df) + geom_line(aes_string(x='position', y='value', col='combination'), size=2)

             if (statistic == 'fold') {

-                ggplt <- ggplt + ylab('fold enrichment')

-                ggplt <- ggplt + geom_hline(yintercept=1, lty=2)

+                if (logscale) {

+                    ggplt <- ggplt + ylab('log(observed/expected)')

+                    ggplt <- ggplt + geom_hline(yintercept=0, lty=2)

+                } else {

+                    ggplt <- ggplt + ylab('observed/expected')

+                    ggplt <- ggplt + geom_hline(yintercept=1, lty=2)

+                }

             } else if (statistic == 'fraction') {

                 ggplt <- ggplt + ylab('fraction')

             }

-            ggplt <- ggplt + scale_color_manual(values = getDistinctColors(length(unique(df$combination))))

+            # ggplt <- ggplt + scale_color_manual(values = getDistinctColors(length(unique(df$combination))))

+            maxcols[[condition]] <- length(unique(df$combination))

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

             ggplt <- ggplot(df) + geom_line(aes_string(x='position', y='value', col='mark'), size=2)

             if (statistic == 'fold') {

-                ggplt <- ggplt + ylab('fold enrichment')

-                ggplt <- ggplt + geom_hline(yintercept=1, lty=2)

+                if (logscale) {

+                    ggplt <- ggplt + ylab('log(observed/expected)')

+                    ggplt <- ggplt + geom_hline(yintercept=0, lty=2)

+                } else {

+                    ggplt <- ggplt + ylab('observed/expected')

+                    ggplt <- ggplt + geom_hline(yintercept=1, lty=2)

+                }

             } else if (statistic == 'fraction') {

                 ggplt <- ggplt + ylab('fraction')

             }

-            ggplt <- ggplt + scale_color_manual(values = getDistinctColors(length(unique(df$mark))))

+            # ggplt <- ggplt + scale_color_manual(values = getDistinctColors(length(unique(df$mark))))

+            maxcols[[condition]] <- length(unique(df$mark))

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

-            ggplt <- ggplot(df) + geom_line(aes_string(x='position', y='value', col='track'), size=2) + ylab('RPKM')

-            ggplt <- ggplt + scale_color_manual(values = getDistinctColors(length(unique(df$track))))

+            ggplt <- ggplot(df) + geom_line(aes_string(x='position', y='value', col='ID'), size=2)

+            if (logscale) {

+                ggplt <- ggplt + ylab('log(RPKM)')

+            } else {

+                ggplt <- ggplt + ylab('RPKM')

+            }

+            # ggplt <- ggplt + scale_color_manual(values = getDistinctColors(length(unique(df$ID))))

+            maxcols[[condition]] <- length(unique(df$ID))

         }

         ggplt <- ggplt + theme_bw() + xlab('distance from annotation in [bp]')

         if (length(region)>=2 & 'inside' %in% region) {

@@ -420,16 +514,23 @@ plotEnrichment <- function(hmm, annotation, bp.around.annotation=10000, region=c

             labels <- c(-bp.around.annotation, -bp.around.annotation/2, '0%', '50%', '100%', bp.around.annotation/2, bp.around.annotation)

             ggplt <- ggplt + scale_x_continuous(breaks=breaks, labels=labels)

         }

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

+        foldsnoinf <- setdiff(df$value, c(Inf, -Inf))

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

+        minfolds[[condition]] <- min(foldsnoinf, na.rm=TRUE)

         ggplts[[condition]] <- ggplt

     }

-    maxfolds <- unlist(maxfolds)

+    maxfold <- max(unlist(maxfolds), na.rm=TRUE)

+    minfold <- min(unlist(minfolds), na.rm=TRUE)

+    maxcol <- max(unlist(maxcols), na.rm=TRUE)

     if (statistic == 'fraction' & what %in% c('combinations','peaks')) {

         ggplts <- lapply(ggplts, function(ggplt) { ggplt + scale_y_continuous(limits=c(0,1)) })

     } else {

-        ggplts <- lapply(ggplts, function(ggplt) { ggplt + scale_y_continuous(limits=c(0,max(maxfolds, na.rm=TRUE)*1.1)) })

+        ggplts <- lapply(ggplts, function(ggplt) { ggplt + scale_y_continuous(limits=c(minfold*(1-sign(minfold)*0.1),maxfold*(1+sign(maxfold)*0.1))) })

     }

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

+    ggplts <- lapply(ggplts, function(ggplt) { ggplt + scale_color_manual(values=getDistinctColors(maxcol)) }) # Add color here like this because of weird bug

+    if (class(hmm) == class.univariate.hmm) {

+        return(ggplts[[1]])

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

         return(ggplts[[1]])

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

         return(ggplts)

@@ -448,12 +549,12 @@ plotEnrichment <- function(hmm, annotation, bp.around.annotation=10000, region=c

 #' @param annotation A \code{\link{GRanges}} object with the annotation of interest.

 #' @param bp.around.annotation An integer specifying the number of basepairs up- and downstream of the annotation for which the enrichment will be calculated.

 #' @param region A combination of \code{c('start','inside','end')} specifying the region of the annotation for which the enrichment will be calculated. Select \code{'start'} if you have a point-sized annotation like transcription start sites. Select \code{c('start','inside','end')} if you have long annotations like genes.

-#' @param what One of \code{c('combinations','peaks','counts')} specifying which statistic to calculate.

+#' @param what One of \code{c('combinations','peaks','counts')} specifying on which feature the statistic is calculated.

 #' @param num.intervals Number of intervals for enrichment 'inside' of annotation.

 #' @param statistic The statistic to calculate. Either 'fold' for fold enrichments or 'fraction' for fraction of bins falling into the annotation.

 #' @return A \code{list()} containing \code{data.frame()}s for enrichment of combinatorial states and binary states at the start, end and inside of the annotation.

 #' @importFrom S4Vectors as.factor subjectHits queryHits

-enrichmentAtAnnotation <- function(bins, info, annotation, bp.around.annotation=10000, region=c('start','inside','end'), what=c('combinations','peaks','counts'), num.intervals=21, statistic='fold') {

+enrichmentAtAnnotation <- function(bins, info, annotation, bp.around.annotation=10000, region=c('start','inside','end'), what='combinations', num.intervals=21, statistic='fold') {

     ## Check user input

     if ((!what %in% c('combinations','peaks','counts')) | length(what) > 1) {

@@ -475,7 +576,11 @@ enrichmentAtAnnotation <- function(bins, info, annotation, bp.around.annotation=

     if ('peaks' %in% what) {

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

         # Remove replicates

-        binstates <- binstates[ ,info.dedup$ID]

+        if (ncol(binstates) > 1) {

+            binstates <- binstates[ ,info.dedup$ID]

+        } else {

+            binstates <- matrix(binstates[ ,info.dedup$ID], ncol=1)

+        }

         colsums.binstates <- colSums(binstates)

     }

     if ('counts' %in% what) {

@@ -490,9 +595,9 @@ enrichmentAtAnnotation <- function(bins, info, annotation, bp.around.annotation=

         widths.annotation <- width(annotation) - 1

         annotation.1bp <- resize(annotation, 1, fix='start')

         # Initialize arrays

-        if ('peaks' %in% what) binstates.inside <- array(dim=c(num.intervals+1, length(info.dedup$ID)), dimnames=list(interval=intervals, track=info.dedup$ID))

+        if ('peaks' %in% what) binstates.inside <- array(dim=c(num.intervals+1, length(info.dedup$ID)), dimnames=list(interval=intervals, ID=info.dedup$ID))

         if ('combinations' %in% what) combinations.inside <- array(dim=c(num.intervals+1, length(levels(bins$combination))), dimnames=list(interval=intervals, combination=levels(bins$combination)))

-        if ('counts' %in% what) counts.inside <- array(dim=c(num.intervals+1, length(info$ID)), dimnames=list(interval=intervals, track=info$ID))

+        if ('counts' %in% what) counts.inside <- array(dim=c(num.intervals+1, length(info$ID)), dimnames=list(interval=intervals, ID=info$ID))

         for (interval in intervals) {

             shift <- widths.annotation * interval * c(1,-1,1)[as.integer(strand(annotation))]

@@ -506,10 +611,15 @@ enrichmentAtAnnotation <- function(bins, info, annotation, bp.around.annotation=

             index <- c(index.inside.plus, index.inside.minus)

             index <- index[index>0 & index<=length(bins)] # index could cross chromosome boundaries, but we risk it

             if ('peaks' %in% what) {

+                if (ncol(binstates) > 1) {

+                    binstates.index <- binstates[index,]

+                } else {

+                    binstates.index <- matrix(binstates[index,], ncol=1)

+                }

                 if (statistic == 'fraction') {

-                    binstates.inside[as.character(interval),] <- colSums(binstates[index,]) / length(index) # or colMeans

+                    binstates.inside[as.character(interval),] <- colSums(binstates.index) / length(index) # or colMeans

                 } else if (statistic == 'fold') {

-                    binstates.inside[as.character(interval),] <- colSums(binstates[index,]) / length(index) / colsums.binstates * length(bins)

+                    binstates.inside[as.character(interval),] <- colSums(binstates.index) / length(index) / colsums.binstates * length(bins)

                 }

             }

             if ('combinations' %in% what) {

@@ -546,17 +656,22 @@ enrichmentAtAnnotation <- function(bins, info, annotation, bp.around.annotation=

         index.start.plus <- index.start.plus[!is.na(index.start.plus)]

         index.start.minus <- index.start.minus[!is.na(index.start.minus)]

         # Occurrences at every bin position relative to feature

-        if ('peaks' %in% what) binstates.start <- array(dim=c(length(-lag:lag), length(info.dedup$ID)), dimnames=list(lag=-lag:lag, track=info.dedup$ID))