R/combinatorialStates.R

In ataudt/chromstaR: Combinatorial and Differential Chromatin State Analysis for ChIP-Seq Data

#' Get the (decimal) combinatorial states of a list of univariate HMM models
#'
#' Get the combinatorial states of a list of models generated by \code{\link{callPeaksUnivariate}}. The function returns the decimal combinatorial states for each bin (see details for an explanation of combinatorial state).
#'
#' For a given model, each genomic bin can be either called 'unmodified' or 'modified', depending on the posterior probabilities estimated by the Baum-Welch. Thus, a list of models defines a binary combinatorial state for each bin. This binary combinatorial state can be expressed as a decimal number.
#' Example: We have 4 histone modifications, and we run the univariate HMM for each of them. Then we use a false discovery rate of 0.5 to call each bin either 'unmodified' or 'modified'. The resulting binary combinatorial states can then be converted to decimal representation. The following table illustrates this:
#' 
#' \tabular{crrrrr}{
#' bin \tab modification state \tab \tab \tab \tab decimal state\cr
#'     \tab model1 \tab model2 \tab model3 \tab model4 \tab     \cr
#' 1   \tab      0 \tab      0 \tab      1 \tab      0 \tab 2   \cr
#' 2   \tab      0 \tab      0 \tab      0 \tab      0 \tab 0   \cr
#' 3   \tab      0 \tab      1 \tab      1 \tab      0 \tab 6   \cr
#' 4   \tab      0 \tab      1 \tab      1 \tab      1 \tab 7   \cr
#' }
#' @author Aaron Taudt
#' @param hmm.list A list of models generated by \code{\link{callPeaksUnivariate}}, e.g. 'list(model1,model2,...)'.
#' @param binary If \code{TRUE}, a matrix of binary instead of decimal states will be returned.
#' @return Output is a vector of integers representing the combinatorial state of each bin.
#' @seealso \code{\link{dec2bin}}, \code{\link{bin2dec}}
#' @examples
#'# Get example BAM files for 3 different marks in hypertensive rat (SHR)
#'file.path <- system.file("extdata","euratrans", package='chromstaRData')
#'files <- list.files(file.path, full.names=TRUE, pattern='SHR.*bam$')[c(1,4,6)]
#'# Bin the data
#'data(rn4_chrominfo)
#'binned.data <- list()
#'for (file in files) {
#'  binned.data[[basename(file)]] <- binReads(file, binsizes=1000, stepsizes=500,
#'                                            assembly=rn4_chrominfo, chromosomes='chr12')
#'}
#'# Obtain the univariate fits
#'models <- list()
#'for (i1 in 1:length(binned.data)) {
#'  models[[i1]] <- callPeaksUnivariate(binned.data[[i1]], max.time=60, eps=1)
#'}
#'## Get the decimal representation of the combinatorial state of this combination of models
#'states <- chromstaR:::combinatorialStates(models, binary=FALSE)
#'## Show number of each state
#'table(states)
#'
combinatorialStates = function(hmm.list, binary=FALSE) {
    
#     # Combine the input
#     models = c(list(...), hmm.list) # This will result in excessive memory usage if the list elements are large, because they will be copied to a new location
    nummod = length(hmm.list)
    numbins = length(hmm.list[[1]]$bins)

    # Get the univariate states (zero inflation = 0, unmodified = 0, modified = 1) from the hmm.list
    binary_statesmatrix = matrix(rep(NA,numbins*nummod), ncol=nummod)
    for (imod in 1:nummod) {
        binary_statesmatrix[,imod] = c(FALSE,FALSE,TRUE)[hmm.list[[imod]]$bins$state] # F,F,T corresponds to levels 'zero-inflation','unmodified','modified'
    }

    if (binary == TRUE) {
        return(binary_statesmatrix)
    } else {
    
        # Transform binary to decimal
        decimal_states = rep(0,numbins)
        for (imod in 1:nummod) {
            decimal_states = decimal_states + 2^(nummod-imod) * binary_statesmatrix[,imod]
        }

        return(decimal_states)

    }

}