callPeaksMultivariate {chromstaR} | R Documentation |
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
.
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 )
hmms |
A list of |
use.states |
A data.frame with combinatorial states which are used in the multivariate HMM, generated by function |
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 |
chromosomes |
A vector specifying the chromosomes to use from the models in |
eps |
Convergence threshold for the Baum-Welch algorithm. |
keep.posteriors |
If set to |
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 |
max.iter |
The maximum number of iterations for the Baum-Welch algorithm. The default |
keep.densities |
If set to |
verbosity |
Verbosity level for the fitting procedure. 0 - No output, 1 - Iterations are printed. |
temp.savedir |
A directory where to store intermediate results if |
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.
A multiHMM
object.
Aaron Taudt, Maria Colome Tatche
multiHMM
, callPeaksUnivariate
, callPeaksReplicates
# 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)