callPeaksMultivariate {chromstaR} | R Documentation

## Fit a Hidden Markov Model to multiple ChIP-seq samples

### Description

Fit a HMM to multiple ChIP-seq samples to determine the combinatorial state of genomic regions. Input is a list of `uniHMM`s generated by `callPeaksUnivariate`.

### Usage

    callPeaksMultivariate(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, temp.savedir = NULL)

### Arguments

`hmms` | A list of `uniHMM`s generated by `callPeaksUnivariate`, e.g. `list(hmm1,hmm2,...)` or a vector of files that contain such objects, e.g. `c("file1","file2",...)`.
`use.states` | A data.frame with combinatorial states which are used in the multivariate HMM, generated by function `stateBrewer`. If both `use.states` and `max.states` are `NULL`, the maximum possible number of combinatorial states will be used.
`max.states` | Maximum number of combinatorial states to use in the multivariate HMM. The states are ordered by occurrence as determined from the combination of univariate state calls.
`per.chrom` | If `per.chrom=TRUE` chromosomes will be treated separately. This tremendously speeds up the calculation but results might be noisier as compared to `per.chrom=FALSE`, where all chromosomes are concatenated for the HMM.
`chromosomes` | A vector specifying the chromosomes to use from the models in `hmms`. The default (`NULL`) uses all available chromosomes.
`eps` | Convergence threshold for the Baum-Welch algorithm.
`keep.posteriors` | If set to `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.
`num.threads` | Number of threads to use. Setting this to >1 may give increased performance.
`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 `NULL` is no limit.
`max.iter` | The maximum number of iterations for the Baum-Welch algorithm. The default `NULL` is no limit.
`keep.densities` | If set to `TRUE` (default=`FALSE`), densities will be available in the output. This should only be needed debugging.
`verbosity` | Verbosity level for the fitting procedure. 0 - No output, 1 - Iterations are printed.
`temp.savedir` | A directory where to store intermediate results if `per.chrom=TRUE`.

### Details

Emission distributions from the univariate HMMs are used with a Gaussian copula to generate a multivariate emission distribution for each combinatorial state. This multivariate distribution is then kept fixed and the transition probabilities are fitted with a Baum-Welch. Please refer to our manuscript at http://dx.doi.org/10.1101/038612 for a detailed description of the method.

### Value

A `multiHMM` object.

### Author(s)

Aaron Taudt, Maria Colome Tatche

### See Also

`multiHMM`, `callPeaksUnivariate`, `callPeaksReplicates`

### Examples

    # Get example BAM files for 2 different marks in hypertensive rat
    file.path <- system.file("extdata","euratrans", package='chromstaRData')
    files <- list.files(file.path, full.names=TRUE, pattern='SHR.*bam$')[c(1:2,6)]
    # Construct experiment structure
    exp <- data.frame(file=files, mark=c("H3K27me3","H3K27me3","H3K4me3"),
                     condition=rep("SHR",3), replicate=c(1:2,1), pairedEndReads=FALSE,
                     controlFiles=NA)
    states <- stateBrewer(exp, mode='combinatorial')
    # Bin the data
    data(rn4_chrominfo)
    binned.data <- list()
    for (file in files) {
     binned.data[[basename(file)]] <- binReads(file, binsizes=1000, stepsizes=500,
                                               experiment.table=exp,
                                               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)
    }
    # Call multivariate peaks
    multimodel <- callPeaksMultivariate(models, use.states=states, eps=1, max.time=60)
    # Check some plots
    heatmapTransitionProbs(multimodel)
    heatmapCountCorrelation(multimodel)

[Package chromstaR version 1.12.0 Index]