Package 'chromstaR'

Description

This package implements functions for the combinatorial and differential analysis of ChIP-seq data. It was developed for histone modifications with a broad profile but is also suitable for the analysis of transcription factor binding data. A Hidden Markov Model with a mixture of Negative Binomials as emission densities is used to call peaks. Please refer to our manuscript at http://dx.doi.org/10.1101/038612 for a detailed description of the method.

Details

The main function of this package is Chromstar. For a detailed introduction type browseVignettes("chromstaR") and read the vignette. Here is an overview of all plotting functions.

Author(s)

Aaron Taudt, Maria Colome-Tatche, Matthias Heinig, Minh Anh Nguyen

Description

A GRanges-class object which contains binned read counts as meta data column counts. It is output of the binReads function.

Description

Convert aligned reads in .bam or .bed(.gz) format into read counts in equidistant windows.

Usage

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

Arguments

Details

Convert aligned reads from .bam or .bed(.gz) files into read counts in equidistant windows (bins). This function uses GenomicRanges::countOverlaps to calculate the read counts, or alternatively bamsignals::bamProfile if option use.bamsignals is set (only effective for .bam files).

Value

If only one bin size was specified for option binsizes, the function returns a single GRanges-class object with meta data column 'counts' that contains the read count. If multiple binsizes were specified , the function returns a named list() of GRanges-class objects.

Examples

## Get an example BAM file with ChIP-seq reads
file <- system.file("extdata", "euratrans",
                   "lv-H3K27me3-BN-male-bio2-tech1.bam",
                    package="chromstaRData")
## Bin the file into bin size 1000bp
data(rn4_chrominfo)
data(experiment_table)
binned <- binReads(file, experiment.table=experiment_table,
                  assembly=rn4_chrominfo, binsizes=1000,
                  stepsizes=500, chromosomes='chr12')
print(binned)

Description

Fit a HMM to multiple ChIP-seq samples to determine the combinatorial state of genomic regions. Input is a list of uniHMMs 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

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)

Description

Fit an HMM to multiple ChIP-seq replicates and derive correlation measures. Input is a list of uniHMMs generated by callPeaksUnivariate.

Usage

callPeaksReplicates(
  hmm.list,
  max.states = 32,
  force.equal = FALSE,
  eps = 0.01,
  max.iter = NULL,
  max.time = NULL,
  keep.posteriors = TRUE,
  num.threads = 1,
  max.distance = 0.2,
  per.chrom = TRUE
)

Arguments

Value

Output is a multiHMM object with additional entry replicateInfo. If only one uniHMM was given as input, a simple list() with the replicateInfo is returned.

Author(s)

Aaron Taudt

See Also

multiHMM, callPeaksUnivariate, callPeaksMultivariate

Examples

# Let's get some example data with 3 replicates
file.path <- system.file("extdata","euratrans", package='chromstaRData')
files <- list.files(file.path, pattern="H3K27me3.*SHR.*bam$", full.names=TRUE)[1:3]
# Obtain chromosome lengths. This is only necessary for BED files. BAM files are
# handled automatically.
data(rn4_chrominfo)
# Define experiment structure
exp <- data.frame(file=files, mark='H3K27me3', condition='SHR', replicate=1:3,
                 pairedEndReads=FALSE, controlFiles=NA)
# We use bin size 1000bp and chromosome 12 to keep the example quick
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')
}
# The univariate fit is obtained for each replicate
models <- list()
for (i1 in 1:length(binned.data)) {
 models[[i1]] <- callPeaksUnivariate(binned.data[[i1]], max.time=60, eps=1)
}
# Obtain peak calls considering information from all replicates
multi.model <- callPeaksReplicates(models, force.equal=TRUE, max.time=60, eps=1)

Description

Fit a HMM to a ChIP-seq sample to determine the modification state of genomic regions, e.g. call peaks in the sample.

Usage

callPeaksUnivariate(
  binned.data,
  control.data = NULL,
  prefit.on.chr = NULL,
  short = TRUE,
  eps = 0.1,
  init = "standard",
  max.time = NULL,
  max.iter = 5000,
  num.trials = 1,
  eps.try = NULL,
  num.threads = 1,
  read.cutoff = TRUE,
  read.cutoff.quantile = 1,
  read.cutoff.absolute = 500,
  max.mean = Inf,
  post.cutoff = 0.5,
  control = FALSE,
  keep.posteriors = FALSE,
  keep.densities = FALSE,
  verbosity = 1
)

Arguments

Details

This function is similar to callPeaksUnivariateAllChr but allows to pre-fit on a single chromosome instead of the whole genome. This gives a significant performance increase and can help to converge into a better fit in case of unsteady quality for some chromosomes.

Value

A uniHMM object.

Author(s)

Aaron Taudt, Maria Colome Tatche

See Also

uniHMM, callPeaksMultivariate

Examples

## Get an example BAM file with ChIP-seq reads
file <- system.file("extdata", "euratrans",
                      "lv-H3K27me3-BN-male-bio2-tech1.bam",
                       package="chromstaRData")
## Bin the BED file into bin size 1000bp
data(rn4_chrominfo)
data(experiment_table)
binned <- binReads(file, experiment.table=experiment_table,
                  assembly=rn4_chrominfo, binsizes=1000,
                  stepsizes=500, chromosomes='chr12')
## Fit the univariate Hidden Markov Model
hmm <- callPeaksUnivariate(binned, max.time=60, eps=1)
## Check if the fit is ok
plotHistogram(hmm)

Description

Fit a HMM to a ChIP-seq sample to determine the modification state of genomic regions, e.g. call peaks in the sample.

Usage

callPeaksUnivariateAllChr(
  binned.data,
  control.data = NULL,
  eps = 0.01,
  init = "standard",
  max.time = NULL,
  max.iter = NULL,
  num.trials = 1,
  eps.try = NULL,
  num.threads = 1,
  read.cutoff = TRUE,
  read.cutoff.quantile = 1,
  read.cutoff.absolute = 500,
  max.mean = Inf,
  post.cutoff = 0.5,
  control = FALSE,
  keep.posteriors = FALSE,
  keep.densities = FALSE,
  verbosity = 1
)

Arguments

Details

The Hidden Markov Model which is used to classify the bins uses 3 states: state 'zero-inflation' with a delta function as emission densitiy (only zero read counts), 'unmodified' and 'modified' with Negative Binomials as emission densities. A Baum-Welch algorithm is employed to estimate the parameters of the distributions. Please refer to our manuscript at http://dx.doi.org/10.1101/038612 for a detailed description of the method.

Value

A uniHMM object.

Author(s)

Aaron Taudt, Maria Coome Tatche

See Also

uniHMM, callPeaksMultivariate

Description

Adjusts the peak calls of a uniHMM, multiHMM or combinedMultiHMM object with a cutoff on the maximum-posterior within each peak. Higher values of maxPost.cutoff mean less sensitive and more precise peak calls. Remaining peaks are kept intact, as opposed to function changePostCutoff, where broad peaks are fragmented. This function was formerly called 'changeFDR' and is still available for backwards compatibiltiy.

Usage

changeMaxPostCutoff(model, maxPost.cutoff = 0.99, invert = FALSE)

changeFDR(model, fdr = 0.01, invert = FALSE)

Arguments

Details

Each peak has a maximum-posterior (maxPostInPeak, between 0 and 1) associated. The sensitivity is adjusted with a simple cutoff on maxPostInPeak, e.g. for maxPost.cutoff = 0.99 only peaks with maxPostInPeak >= 0.99 will be selected.

Value

The input object is returned with adjusted peak calls.

Functions

• changeFDR: This function was renamed to 'changeMaxPostCutoff' in chromstaR 1.5.1 but it still available for backwards compatibility.

Author(s)

Aaron Taudt

See Also

changePostCutoff

Examples

## Get an example uniHMM ##
file <- system.file("data","H3K27me3-BN-rep1.RData", package="chromstaR")
model <- get(load(file))
## Compare fits with different fdrs
plotHistogram(model) + ylim(0,0.25) + ylim(0,0.3)
plotHistogram(changeMaxPostCutoff(model, maxPost.cutoff=0.99)) + ylim(0,0.3)
plotHistogram(changeMaxPostCutoff(model, maxPost.cutoff=1-1e-12)) + ylim(0,0.3)

## Get an example multiHMM ##
file <- system.file("data","multivariate_mode-combinatorial_condition-SHR.RData",
                    package="chromstaR")
model <- get(load(file))
genomicFrequencies(model)
model.new <- changeMaxPostCutoff(model, maxPost.cutoff=0.9999, invert=FALSE)
genomicFrequencies(model.new)

## Get an example combinedMultiHMM ##
file <- system.file("data","combined_mode-differential.RData",
                    package="chromstaR")
model <- get(load(file))
genomicFrequencies(model)
model.new <- changeMaxPostCutoff(model, maxPost.cutoff=0.9999, invert=FALSE)
genomicFrequencies(model.new)

Description

Adjusts the peak calls of a uniHMM, multiHMM or combinedMultiHMM object with the given posterior cutoff.

Usage

changePostCutoff(model, post.cutoff = 0.5)

Arguments

Details

Posterior probabilities are between 0 and 1. Peaks are called if the posteriors for a state (univariate) or sample (multivariate) are >= post.cutoff.

Value

The input object is returned with adjusted peak calls.

Author(s)

Aaron Taudt

See Also

changeMaxPostCutoff

Examples

## Get an example BAM file with ChIP-seq reads
file <- system.file("extdata", "euratrans",
                      "lv-H3K27me3-BN-male-bio2-tech1.bam",
                       package="chromstaRData")
## Bin the BED file into bin size 1000bp
data(rn4_chrominfo)
data(experiment_table)
binned <- binReads(file, experiment.table=experiment_table,
                  assembly=rn4_chrominfo, binsizes=1000,
                  stepsizes=500, chromosomes='chr12')
plotHistogram(binned)
## Fit HMM
model <- callPeaksUnivariate(binned, keep.posteriors=TRUE, verbosity=0)
## Compare fits with different post.cutoffs
plotHistogram(changePostCutoff(model, post.cutoff=0.01)) + ylim(0,0.25)
plotHistogram(model) + ylim(0,0.25)
plotHistogram(changePostCutoff(model, post.cutoff=0.99)) + ylim(0,0.25)

## Get an example multiHMM ##
file <- system.file("data","multivariate_mode-combinatorial_condition-SHR.RData",
                    package="chromstaR")
model <- get(load(file))
genomicFrequencies(model)
model.new <- changePostCutoff(model, post.cutoff=0.9999)
genomicFrequencies(model.new)

## Get an example combinedMultiHMM ##
file <- system.file("data","combined_mode-differential.RData",
                    package="chromstaR")
model <- get(load(file))
genomicFrequencies(model)
model.new <- changePostCutoff(model, post.cutoff=0.9999)
genomicFrequencies(model.new)

Wrapper function for the chromstaR package

Description

This function performs binning, univariate peak calling and multivariate peak calling from a list of input files.

Usage

Chromstar(
  inputfolder,
  experiment.table,
  outputfolder,
  configfile = NULL,
  numCPU = 1,
  binsize = 1000,
  stepsize = binsize/2,
  assembly = NULL,
  chromosomes = NULL,
  remove.duplicate.reads = TRUE,
  min.mapq = 10,
  format = NULL,
  prefit.on.chr = NULL,
  eps.univariate = 0.1,
  max.time = NULL,
  max.iter = 5000,
  read.cutoff.absolute = 500,
  keep.posteriors = TRUE,
  mode = "differential",
  max.states = 128,
  per.chrom = TRUE,
  eps.multivariate = 0.01,
  exclusive.table = NULL
)

Arguments

Value

NULL

Examples

## Prepare the file paths. Exchange this with your input and output directories.
inputfolder <- system.file("extdata","euratrans", package="chromstaRData")
outputfolder <- file.path(tempdir(), 'SHR-example')
## Define experiment structure
data(experiment_table_SHR)
## Define assembly
# This is only necessary if you have BED files, BAM files are handled automatically.
# For common assemblies you can also specify them as 'hg19' for example.
data(rn4_chrominfo)
## Run ChromstaR
Chromstar(inputfolder, experiment.table=experiment_table_SHR,
         outputfolder=outputfolder, numCPU=4, binsize=1000, assembly=rn4_chrominfo,
         prefit.on.chr='chr12', chromosomes='chr12', mode='combinatorial', eps.univariate=1,
         eps.multivariate=1)

Description

chromstaR defines several objects.

• uniHMM: Returned by callPeaksUnivariate.

• multiHMM: Returned by callPeaksMultivariate and callPeaksReplicates.

• combinedMultiHMM: Returned by combineMultivariates.

Description

The function will collapse consecutive bins which have, for example, the same combinatorial state.

Usage

collapseBins(
  data,
  column2collapseBy = NULL,
  columns2sumUp = NULL,
  columns2average = NULL,
  columns2getMax = NULL,
  columns2drop = NULL
)

Arguments

Details

The following tables illustrate the principle of the collapsing:

Input data:

Output data:

Value

A data.frame.

Author(s)

Aaron Taudt

Examples

## Load example data
## Get an example multiHMM
file <- system.file("data","multivariate_mode-combinatorial_condition-SHR.RData",
                    package="chromstaR")
model <- get(load(file))
df <- as.data.frame(model$bins)
shortdf <- collapseBins(df, column2collapseBy='state', columns2sumUp='width', columns2average=6:9)

Description

Get the combinatorial states of a list of models generated by callPeaksUnivariate. The function returns the decimal combinatorial states for each bin (see details for an explanation of combinatorial state).

Usage

combinatorialStates(hmm.list, binary = FALSE)

Arguments

Details

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:

Value

Output is a vector of integers representing the combinatorial state of each bin.

Author(s)

Aaron Taudt

See Also

dec2bin, 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)

Description

The combined multivariate HMM object is output of the function combineMultivariates and is a list() with various entries. The class() attribute of this list was set to "combinedMultiHMM". For a given hmm, the entries can be accessed with the list operators 'hmm[[]]' or 'hmm$'.

Value

A list() with the following entries:

See Also

combineMultivariates, uniHMM, multiHMM

Description

Combine combinatorial states from several multiHMM objects. Combinatorial states can be combined for objects containing multiple marks (mode='combinatorial') or multiple conditions (mode='differential').

Usage

combineMultivariates(hmms, mode)

Arguments

Value

A combinedMultiHMM objects with combinatorial states for each condition.

Author(s)

Aaron Taudt

Examples

### Multivariate peak calling for spontaneous hypertensive rat (SHR) ###
# Get example BAM files for 2 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:2,4:5)]
# Construct experiment structure
exp <- data.frame(file=files, mark=c("H3K27me3","H3K27me3","H3K4me3","H3K4me3"),
                 condition=rep("SHR",4), replicate=c(1:2,1:2), 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.SHR <- callPeaksMultivariate(models, use.states=states, eps=1, max.time=60)

#'### Multivariate peak calling for brown norway (BN) rat ###
# Get example BAM files for 2 different marks in brown norway rat
file.path <- system.file("extdata","euratrans", package='chromstaRData')
files <- list.files(file.path, full.names=TRUE, pattern='BN.*bam$')[c(1:2,3:4)]
# Construct experiment structure
exp <- data.frame(file=files, mark=c("H3K27me3","H3K27me3","H3K4me3","H3K4me3"),
                 condition=rep("BN",4), replicate=c(1:2,1:2), 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.BN <- callPeaksMultivariate(models, use.states=states, eps=1, max.time=60)

### Combine multivariates ###
hmms <- list(multimodel.SHR, multimodel.BN)
comb.model <- combineMultivariates(hmms, mode='combinatorial')

Description

Convert combinatorial states in decimal representation to combinatorial states in binary representation and vice versa.

Usage

dec2bin(dec, colnames = NULL, ndigits = NULL)

bin2dec(bin)

Arguments

Details

chromstaR uses decimal numbers to represent combinatorial states of peaks. These functions serve as a convenient way to get from the efficient decimal representation to a more human-readable binary representation.

Value

A vector of integers for bin2dec and a matrix of logicals with one state per row for dec2bin.

Functions

• dec2bin: Decimal to binary conversion.

• bin2dec: Binary to decimal conversion.

Author(s)

Aaron Taudt

Examples

decimal.states <- c(0:31)
binary.states <- dec2bin(decimal.states, colnames=paste0('mark',1:5))
control.decimal.states <- bin2dec(binary.states)

Description

Plotting functions for enrichment analysis of multiHMM or combinedMultiHMM objects with any annotation of interest, specified as a GRanges-class object.

Usage

plotFoldEnrichHeatmap(
  hmm,
  annotations,
  what = "combinations",
  combinations = NULL,
  marks = NULL,
  plot = TRUE,
  logscale = TRUE
)

plotEnrichCountHeatmap(
  hmm,
  annotation,
  bp.around.annotation = 10000,
  max.rows = 1000,
  combinations = NULL,
  colorByCombinations = sortByCombinations,
  sortByCombinations = is.null(sortByColumns),
  sortByColumns = NULL
)

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

Arguments

Value

A ggplot object containing the plot or a list() with ggplot objects if several plots are returned. For plotFoldEnrichHeatmap a named array with fold enrichments if plot=FALSE.

Functions

• plotFoldEnrichHeatmap: Compute the fold enrichment of combinatorial states for multiple annotations.

• plotEnrichCountHeatmap: Plot read counts around annotation as heatmap.

• plotEnrichment: Plot fold enrichment of combinatorial states around and inside of annotation.

Author(s)

Aaron Taudt

See Also

plotting

Examples

### Get an example multiHMM ###
file <- system.file("data","multivariate_mode-combinatorial_condition-SHR.RData",
                    package="chromstaR")
model <- get(load(file))

### Obtain gene coordinates for rat from biomaRt ###
library(biomaRt)
ensembl <- useEnsembl(biomart='ENSEMBL_MART_ENSEMBL', dataset='rnorvegicus_gene_ensembl')
genes <- getBM(attributes=c('ensembl_gene_id', 'chromosome_name', 'start_position',
                           'end_position', 'strand', 'external_gene_name',
                           'gene_biotype'),
              mart=ensembl)
# Transform to GRanges for easier handling
genes <- GRanges(seqnames=paste0('chr',genes$chromosome_name),
                ranges=IRanges(start=genes$start, end=genes$end),
                strand=genes$strand,
                name=genes$external_gene_name, biotype=genes$gene_biotype)
# Rename chrMT to chrM
seqlevels(genes)[seqlevels(genes)=='chrMT'] <- 'chrM'
print(genes)

### Make the enrichment plots ###
# We expect promoter [H3K4me3] and bivalent-promoter signatures [H3K4me3+H3K27me3]
# to be enriched at transcription start sites.
   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.
   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
   plotEnrichment(model, genes, region = 'start') +
   facet_wrap(~ combination)
 
# Another form of visualization that shows every TSS in a heatmap
# If transparency is not supported try to plot to pdf() instead.
   tss <- resize(genes, width = 3, fix = 'start')
   plotEnrichCountHeatmap(model, tss) +
   theme(strip.text.x = element_text(size=6))
 
# 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()

Description

The function calculates the enrichment of a genomic feature with peaks or combinatorial states. Input is a multiHMM object (containing the peak calls and combinatorial states) and a GRanges-class object containing the annotation of interest (e.g. transcription start sites or genes).

Usage

enrichmentAtAnnotation(
  bins,
  info,
  annotation,
  bp.around.annotation = 10000,
  region = c("start", "inside", "end"),
  what = "combinations",
  num.intervals = 21,
  statistic = "fold"
)

Arguments

Value

A list() containing data.frame()s for enrichment of combinatorial states and binary states at the start, end and inside of the annotation.

Author(s)

Aaron Taudt

Description

A data.frame specifying the structure of the experiment.

Format

A data.frame with columns 'file', 'mark', 'condition', 'replicate', 'pairedEndReads' and 'controlFiles'. Avoid the use of special characters like '-' or '+' as this will confuse the internal file management.

Examples

data(experiment_table)
print(experiment_table)

Description

These functions allow to export chromstaR-objects as files which can be uploaded to a genome browser. Peak calls are exported in BED format (.bed.gz), read counts in wiggle format (.wig.gz) as RPKM values, and combinatorial states are exported in BED format (.bed.gz).

Usage

exportPeaks(
  model,
  filename,
  header = TRUE,
  separate.files = TRUE,
  trackname = NULL
)

exportCounts(
  model,
  filename,
  header = TRUE,
  separate.files = TRUE,
  trackname = NULL
)

exportCombinations(
  model,
  filename,
  header = TRUE,
  separate.files = TRUE,
  trackname = NULL,
  exclude.states = "[]",
  include.states = NULL
)

Arguments

Value

NULL

Functions

• exportPeaks: Export peak calls in BED format.

• exportCounts: Export read counts as RPKM values in wiggle format.

• exportCombinations: Export combinatorial states in BED format.

Examples

## Get an example multiHMM
file <- system.file("data","combined_mode-differential.RData",
                    package="chromstaR")
model <- get(load(file))
## Export peak calls and combinatorial states
exportPeaks(model, filename=tempfile())
exportCombinations(model, filename=tempfile())

Description

Export GRanges as genome browser viewable file

Usage

exportGRangesAsBedFile(
  gr,
  trackname,
  filename,
  namecol = "combination",
  scorecol = "score",
  colorcol = NULL,
  colors = NULL,
  header = TRUE,
  append = FALSE
)

Arguments

Details

Export regions from GRanges-class as a file which can be uploaded into a genome browser. Regions are exported in BED format (.bed.gz).

Value

NULL

Author(s)

Aaron Taudt

See Also

exportPeaks, exportCounts, exportCombinations

Examples

### Export regions with read counts above 20 ###
# Get an example BAM file with ChIP-seq reads
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,
                  stepsizes=500, chromosomes='chr12')
plotHistogram(binned)
# Export regions with read count above 20
exportGRangesAsBedFile(binned[binned$counts[,1] > 20], filename=tempfile(),
             trackname='read counts above 20')

Description

Make fixed-width bins based on given bin size.

Usage

fixedWidthBins(
  bamfile = NULL,
  assembly = NULL,
  chrom.lengths = NULL,
  chromosome.format,
  binsizes = 1e+06,
  chromosomes = NULL
)

Arguments

Value

A list() of GRanges-class objects with fixed-width bins.

Author(s)

Aaron Taudt

Examples

## Make fixed-width bins of size 500kb and 1Mb
data(rn4_chrominfo)
chrom.lengths <- rn4_chrominfo$length
names(chrom.lengths) <- rn4_chrominfo$chromosome
bins <- fixedWidthBins(chrom.lengths=chrom.lengths, binsizes=c(5e5,1e6))
bins

## Make bins using NCBI server (requires internet connection)
# bins <- fixedWidthBins(assembly='mm10', chromosome.format='NCBI', binsizes=c(5e5,1e6))

Description

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

Format

A data.frame.

Examples

data(genes_rn4)
head(genes_rn4)

Description

Get the genomewide frequency of each combinatorial state.

Usage

genomicFrequencies(multi.hmm, combinations = NULL, per.mark = FALSE)

Arguments

Value

A table with frequencies of each combinatorial state.

Author(s)

Aaron Taudt

Examples

## Get an example multiHMM
file <- system.file("data","multivariate_mode-combinatorial_condition-SHR.RData",
                    package="chromstaR")
model <- get(load(file))
genomicFrequencies(model)

Description

Get a DataFrame with combinations from a GRanges-class object.

Usage

getCombinations(gr)

Arguments

Value

A DataFrame.

Examples

### Get an example multiHMM ###
file <- system.file("data","multivariate_mode-combinatorial_condition-SHR.RData",
                    package="chromstaR")
model <- get(load(file))
### Get the combinations
bin.combs <- getCombinations(model$bins)
print(bin.combs)
seg.combs <- getCombinations(model$segments)
print(seg.combs)

Description

Get a set of distinct colors selected from colors.

Usage

getDistinctColors(
  n,
  start.color = "blue4",
  exclude.colors = c("white", "black", "gray", "grey", "\\<yellow\\>", "yellow1",
    "lemonchiffon"),
  exclude.brightness.above = 1,
  exclude.rgb.above = 210
)

Arguments

Details

The function computes the euclidian distance between all colors and iteratively selects those that have the furthest closes distance to the set of already selected colors.

Value

A character vector with colors.

Author(s)

Aaron Taudt

Examples

cols <- getDistinctColors(5)
pie(rep(1,5), labels=cols, col=cols)

Description

Get the colors that are used for plotting.

Usage

getStateColors(labels = NULL)

Arguments

Value

A character vector with colors.

See Also

plotting

Examples

cols <- getStateColors()
pie(1:length(cols), col=cols, labels=names(cols))

Description

Plot a heatmap that shows the binary presence/absence of marks for the different combinations.

Usage

heatmapCombinations(model = NULL, marks = NULL, emissionProbs = NULL)

Arguments

Value

A ggplot object.

Author(s)

Aaron Taudt

See Also

plotting

Examples

marks <- c('H3K4me3','H3K27me3','H4K20me1')
heatmapCombinations(marks=marks)

file <- system.file("data","multivariate_mode-combinatorial_condition-SHR.RData",
                    package="chromstaR")
heatmapCombinations(file)

Description

Heatmap of read count correlations (see cor).

Usage

heatmapCountCorrelation(model, cluster = TRUE)

Arguments

Value

A ggplot object.

See Also

plotting

Examples

## Get an example multiHMM ##
file <- system.file("data","multivariate_mode-combinatorial_condition-SHR.RData",
                    package="chromstaR")
model <- get(load(file))
## Plot count correlations as heatmap
heatmapCountCorrelation(model)

Description

Plot a heatmap of transition probabilities for a multiHMM model.

Usage

heatmapTransitionProbs(
  model = NULL,
  reorder.states = TRUE,
  transitionProbs = NULL
)

Arguments

Value

A ggplot object.

See Also

plotting

Examples

## Get an example multiHMM ##
file <- system.file("data","multivariate_mode-combinatorial_condition-SHR.RData",
                    package="chromstaR")
model <- get(load(file))
## Plot transition probabilites as heatmap
heatmapTransitionProbs(model, reorder.states=TRUE)

Description

Wrapper to load chromstaR objects from file and check the class of the loaded objects.

Usage

loadHmmsFromFiles(
  files,
  check.class = c("GRanges", "uniHMM", "multiHMM", "combinedMultiHMM")
)

Arguments

Value

A list of chromstaR-object.

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,
                  stepsizes=500, chromosomes='chr12')
## Fit the univariate Hidden Markov Model
hmm <- callPeaksUnivariate(binned, max.time=60, eps=1)
temp.file <- tempfile()
save(hmm, file=temp.file)
loaded.hmm <- loadHmmsFromFiles(temp.file)[[1]]
class(loaded.hmm)

Merge several multiHMMs into one object

Description

Merge several multiHMMs into one object. This can be done to merge fits for separate chromosomes into one object for easier handling. Merging will only be done if all models have the same IDs.

Usage

mergeChroms(multi.hmm.list, filename = NULL)

Arguments

Value

A multiHMM object or NULL, depending on option filename.

Author(s)

Aaron Taudt

Description

A combinedMultiHMM object for demonstration purposes in examples of package chromstaR.

Format

A combinedMultiHMM object.

Examples

## Get an example combinedMultiHMM
file <- system.file("data","combined_mode-differential.RData",
                    package="chromstaR")
model <- get(load(file))

Description

A multiHMM object for demonstration purposes in examples of package chromstaR.

Format

A multiHMM object.

Examples

## Get an example multiHMM
file <- system.file("data","multivariate_mode-combinatorial_condition-SHR.RData",
                    package="chromstaR")
model <- get(load(file))

Description

A uniHMM object for demonstration purposes in examples of package chromstaR.

Format

A uniHMM object.

Examples

## Get an example uniHMM
file <- system.file("data","H3K27me3-BN-rep1.RData", package="chromstaR")
model <- get(load(file))

Description

The multivariate HMM object is output of the function callPeaksMultivariate and is a list() with various entries. The class() attribute of this list was set to "multiHMM". For a given hmm, the entries can be accessed with the list operators 'hmm[[]]' or 'hmm$'.

Value

A list() with the following entries:

See Also

callPeaksMultivariate, uniHMM, combinedMultiHMM

Examples

## Get an example multiHMM
file <- system.file("data","multivariate_mode-combinatorial_condition-SHR.RData",
                    package="chromstaR")
model <- get(load(file))

Description

Make segmentation from bins for a multiHMM object.

Usage

multivariateSegmentation(bins, column2collapseBy = "state")

Arguments

Value

A GRanges-class with segmented regions.

Description

Get the expression values that overlap with each combinatorial state.

Usage

plotExpression(hmm, expression, combinations = NULL, return.marks = FALSE)

Arguments

Value

A ggplot2 object if a multiHMM was given or a named list with ggplot2 objects if a combinedMultiHMM was given.

Author(s)

Aaron Taudt

See Also

plotting

Examples

## Load an example multiHMM
file <- system.file("data","multivariate_mode-combinatorial_condition-SHR.RData",
                    package="chromstaR")
model <- get(load(file))

## Obtain expression data
data(expression_lv)
head(expression_lv)

## We need to get coordinates for each of the genes
library(biomaRt)
ensembl <- useEnsembl(biomart='ENSEMBL_MART_ENSEMBL', dataset='rnorvegicus_gene_ensembl')
genes <- getBM(attributes=c('ensembl_gene_id', 'chromosome_name', 'start_position',
                           'end_position', 'strand', 'external_gene_name',
                           'gene_biotype'),
              mart=ensembl)
expr <- merge(genes, expression_lv, by='ensembl_gene_id')
# Transform to GRanges
expression.SHR <- GRanges(seqnames=paste0('chr',expr$chromosome_name),
                         ranges=IRanges(start=expr$start, end=expr$end),
                         strand=expr$strand, name=expr$external_gene_name,
                         biotype=expr$gene_biotype,
                         expression=expr$expression_SHR)
# We apply an asinh transformation to reduce the effect of outliers
expression.SHR$expression <- asinh(expression.SHR$expression)

## Plot
plotExpression(model, expression.SHR) +
 theme(axis.text.x=element_text(angle=0, hjust=0.5)) +
 ggtitle('Expression of genes overlapping combinatorial states')
plotExpression(model, expression.SHR, return.marks=TRUE) +
 ggtitle('Expression of marks overlapping combinatorial states')

#' Plot a genome browser view #' #' Plot a simple genome browser view. This is useful for scripted genome browser snapshots. #' #' @param counts A GRanges-class object with meta-data column 'counts'. #' @param peaklist A named list() of GRanges-class objects containing peak coordinates. #' @param chr,start,end Chromosome, start and end coordinates for the plot. #' @param countcol A character giving the color for the counts. #' @param peakcols A character vector with colors for the peaks in peaklist. #' @param style One of c('peaks', 'density'). #' @param peakTrackHeight Relative height of the tracks given in peaklist compared to the counts. #' @return A ggplot object. #' @examples #'## Get an example multiHMM ## #'file <- system.file("data","multivariate_mode-combinatorial_condition-SHR.RData", #' package="chromstaR") #'model <- get(load(file)) #'## Plot genome browser snapshot #'bins <- model$bins #'bins$counts <- model$bins$counts.rpkm[,1] #'plotGenomeBrowser(counts=bins, peaklist=model$peaks, #' chr='chr12', start=1, end=1e6) #' plotGenomeBrowser2 <- function(counts, peaklist=NULL, chr, start, end, countcol='black', peakcols=NULL, style='peaks', peakTrackHeight=5) ## Select ranges to plot ranges2plot <- reduce(counts[counts@seqnames == chr & start(counts) >= start & start(counts) <= end]) ## Counts counts <- subsetByOverlaps(counts, ranges2plot) if (style == 'peaks') df <- data.frame(x=(start(counts)+end(counts))/2, counts=counts$counts) # plot triangles centered at middle of the bin ggplt <- ggplot(df) + geom_area(aes_string(x='x', y='counts')) + theme(panel.grid = element_blank(), panel.background = element_blank(), axis.text.x = element_blank(), axis.title = element_blank(), axis.ticks.x = element_blank(), axis.line = element_blank()) maxcounts <- max(counts$counts) ggplt <- ggplt + scale_y_continuous(breaks=c(0, maxcounts)) else if (style == 'density') df <- data.frame(xmin=start(counts), xmax=end(counts), counts=counts$counts) ggplt <- ggplot(df) + geom_rect(aes_string(xmin='xmin', xmax='xmax', ymin=0, ymax=4, alpha='counts')) + theme(panel.grid = element_blank(), panel.background = element_blank(), axis.text = element_blank(), axis.title = element_blank(), axis.ticks = element_blank(), axis.line = element_blank()) else stop("Unknown value '", style, "' for parameter 'style'. Must be one of c('peaks', 'density').") ## Peaks if (!is.null(peaklist)) if (is.null(peakcols)) peakcols <- getDistinctColors(length(peaklist)) for (i1 in 1:length(peaklist)) p <- peakTrackHeight peaks <- subsetByOverlaps(peaklist[[i1]], ranges2plot) if (length(peaks) > 0) df <- data.frame(start=start(peaks), end=end(peaks), ymin=-p*i1, ymax=-p*i1+0.9*p) ggplt <- ggplt + geom_rect(data=df, mapping=aes_string(xmin='start', xmax='end', ymin='ymin', ymax='ymax'), col=peakcols[i1], fill=peakcols[i1]) trackname <- names(peaklist)[i1] df <- data.frame(x=start(counts)[1], y=-p*i1+0.5*p, label=trackname) ggplt <- ggplt + geom_text(data=df, mapping=aes_string(x='x', y='y', label='label'), vjust=0.5, hjust=0.5, col=peakcols[i1]) return(ggplt) Plot a genome browser view

Description

Plot a simple genome browser view of chromstaR-objects. This is useful for scripted genome browser snapshots.

Usage

plotGenomeBrowser(
  model,
  chr,
  start,
  end,
  style = "peaks",
  peakHeight = 0.2,
  peakColor = "blue",
  same.yaxis = TRUE
)

Arguments

Value

A list() of ggplot objects.

Examples

## Get an example uniHMM ##
file <- system.file("data","H3K27me3-BN-rep1.RData", package="chromstaR")
model <- get(load(file))
plotGenomeBrowser(model, chr='chr12', start=1, end=1e6, style='peaks',
                 peakHeight=0.1)
                 
## Get an example multiHMM ##
file <- system.file("data","multivariate_mode-combinatorial_condition-SHR.RData",
                    package="chromstaR")
model <- get(load(file))
plotGenomeBrowser(model, chr='chr12', start=1, end=1e6, style='peaks',
                 peakHeight=0.1)
                 
## Get an example combinedMultiHMM ##
file <- system.file("data","combined_mode-differential.RData",
                    package="chromstaR")
model <- get(load(file))
plotlist <- plotGenomeBrowser(model, chr='chr12', start=1, end=1e6, style='peaks',
                 peakHeight=0.1)

Description

Plot a histogram of binned read counts with fitted mixture distributions from a uniHMM object.

Usage

plotHistogram(
  model,
  state = NULL,
  chromosomes = NULL,
  start = NULL,
  end = NULL,
  linewidth = 1
)

Arguments

Value

A ggplot object.

See Also

plotting

Examples

## Get an example BAM file with ChIP-seq reads
file <- system.file("extdata", "euratrans",
                      "lv-H3K27me3-BN-male-bio2-tech1.bam",
                       package="chromstaRData")
## Bin the BED file into bin size 1000bp
data(rn4_chrominfo)
data(experiment_table)
binned <- binReads(file, experiment.table=experiment_table,
                  assembly=rn4_chrominfo, binsizes=1000,
                  stepsizes=500, chromosomes='chr12')
plotHistogram(binned)
## Fit the univariate Hidden Markov Model
hmm <- callPeaksUnivariate(binned, max.time=60, eps=1)
## Check if the fit is ok
plotHistogram(hmm)

Description

Plot histograms of binned read counts with fitted mixture distributions from a multiHMM object.

Usage

plotHistograms(model, ...)

Arguments

Value

A ggplot object.

See Also

plotting

Description

This page provides an overview of all chromstaR plotting functions.

Details

Plotting functions that work on uniHMM objects:

plotHistogram

Read count histogram with fitted mixture distributions.

Plotting functions that work on multiHMM objects:

heatmapCountCorrelation

Heatmap of read count correlations.

heatmapTransitionProbs

Heatmap of transition probabilities of the Hidden Markov Model.

heatmapCombinations

Binary presence/absence pattern of combinatorial states.

plotExpression

Boxplot of expression values that overlap combinatorial states.

Plotting functions that work on multiHMM and combinedMultiHMM objects:

heatmapCountCorrelation

Heatmap of read count correlations.

plotEnrichCountHeatmap

Heatmap of read counts around annotation.

plotEnrichment

Enrichment of combinatorial states around annotation.

plotFoldEnrichHeatmap

Enrichment of combinatorial states at multiple annotations.

plotExpression

Boxplot of expression values that overlap combinatorial states.

Other plotting functions:

heatmapCombinations

Binary presence/absence pattern of combinatorial states.

Description

Print combinedMultiHMM object

Usage

## S3 method for class 'combinedMultiHMM'
print(x, ...)

Arguments

Value

An invisible NULL.

Description

Print multiHMM object

Usage

## S3 method for class 'multiHMM'
print(x, ...)

Arguments

Value

An invisible NULL.

Description

Print uniHMM object

Usage

## S3 method for class 'uniHMM'
print(x, ...)

Arguments

Value

An invisible NULL.

Description

Import aligned reads from a BAM file into a GRanges-class object.

Usage

readBamFileAsGRanges(
  bamfile,
  bamindex = bamfile,
  chromosomes = NULL,
  pairedEndReads = FALSE,
  remove.duplicate.reads = FALSE,
  min.mapq = 10,
  max.fragment.width = 1000,
  blacklist = NULL,
  what = "mapq"
)

Arguments

Value

A GRanges-class object containing the reads.

Examples

## Get an example BAM file with ChIP-seq reads
bamfile <- system.file("extdata", "euratrans", "lv-H3K4me3-BN-female-bio1-tech1.bam",
                      package="chromstaRData")
## Read the file into a GRanges object
reads <- readBamFileAsGRanges(bamfile, chromosomes='chr12', pairedEndReads=FALSE,
                    min.mapq=10, remove.duplicate.reads=TRUE)
print(reads)

Description

Import aligned reads from a BED file into a GRanges-class object.

Usage

readBedFileAsGRanges(
  bedfile,
  assembly,
  chromosomes = NULL,
  remove.duplicate.reads = FALSE,
  min.mapq = 10,
  max.fragment.width = 1000,
  blacklist = NULL
)

Arguments

Value

A GRanges-class object containing the reads.

Examples

## Get an example BED file with single-cell-sequencing reads
bedfile <- system.file("extdata", "liver-H3K4me3-BN-male-bio2-tech1.bed.gz",
                       package="chromstaRData")
## Read the file into a GRanges object
data(rn4_chrominfo)
reads <- readBedFileAsGRanges(bedfile, assembly=rn4_chrominfo, chromosomes='chr12',
                    min.mapq=10, remove.duplicate.reads=TRUE)
print(reads)

Description

Read a chromstaR configuration file into a list structure. The configuration file has to be specified in INI format. R expressions can be used and will be evaluated.

Usage

readConfig(configfile)

Arguments

Value

A list with one entry for each element in configfile.

Author(s)

Aaron Taudt

Description

This is a simple convenience function to read a bed(.gz)-file into a GRanges-class object. The bed-file is expected to have the following fields: chromosome, start, end, name, score, strand.

Usage

readCustomBedFile(
  bedfile,
  col.names = c("chromosome", "start", "end", "name", "score", "strand"),
  col.classes = NULL,
  skip = 0,
  chromosome.format = "NCBI",
  sep = ""
)

Arguments

Value

A GRanges-class object with the contents of the bed-file.

Author(s)

Aaron Taudt

Examples

## Get an example BED file
bedfile <- system.file("extdata", "liver-H3K4me3-BN-male-bio2-tech1.bed.gz",
                       package="chromstaRData")
## Import the file and skip the first 10 lines
data <- readCustomBedFile(bedfile, skip=10)

Description

Remove a condition from a combinedMultiHMM object.

Usage

removeCondition(model, conditions)

Arguments

Value

The input combinedMultiHMM object with specified conditions removed.

Examples

## Get an example HMM
file <- system.file("data","combined_mode-differential.RData",
                    package="chromstaR")
model <- get(load(file))

## Print available conditions
print(unique(model$info$condition))

## Remove condition SHR
new.model <- removeCondition(model, conditions='SHR')

Description

Use simulations to find the best bin size among a set of input files. There is no guarantee that the bin size will be the best for your data, since it is only "best" in terms of fewest miscalls for simulated data. However, it can give you a hint what bin size to choose.

Usage

scanBinsizes(
  files.binned,
  outputfolder,
  chromosomes = "chr10",
  eps = 0.01,
  max.iter = 100,
  max.time = 300,
  repetitions = 3,
  plot.progress = FALSE
)

Arguments

Details

The function first runs callPeaksUnivariate on the given binned.data files. From the estimated parameters it generates simulated data and calls the peaks on this simulated data. Because the data is simulated, the fraction of miscalls can be precisely calculated.

Value

A ggplot object with a bar plot of the number of miscalls dependent on the bin size.

Author(s)

Aaron Taudt

Description

Various scores used in chromstaR.

Usage

differentialScoreMax(mat, info, FUN = "-")

differentialScoreSum(mat, info, FUN = "-")

Arguments

Value

A numeric vector.

Functions

• differentialScoreMax: Maximum differential score. Values are between 0 and 1. A value of 1 means that at least one mark is maximally different between conditions.

• differentialScoreSum: Additive differential score. Values are between 0 and N, where N is the number of marks. A value around 1 means that approximately 1 mark is different, a value of 2 means that 2 marks are different etc.

Author(s)

Aaron Taudt

Description

Simulate known states, read counts and read coordinates using a multivariate Hidden Markov Model.

Usage

simulateMultivariate(
  bins,
  transition,
  emissions,
  weights,
  correlationMatrices,
  combstates,
  IDs,
  fragLen = 50
)

Arguments

Value

A list() with entries $bins containing the simulated states and read count, $reads with simulated read coordinates.

Description

Simulate read coordinates using read counts as input.

Usage

simulateReadsFromCounts(bins, fragLen = 50)

Arguments

Value

A GRanges-class with read coordinates.

Description

Simulate known states, read counts and read coordinates using a univariate Hidden Markov Model with three states ("zero-inflation", "unmodified" and "modified").

Usage

simulateUnivariate(bins, transition, emission, fragLen = 50)

Arguments

Value

A list with entries $bins containing the simulated states and read count, $reads with simulated read coordinates and $transition and $emission.

Description

This function returns all combinatorial (decimal) states that are consistent with a given abstract specification.

Usage

state.brewer(
  replicates = NULL,
  differential.states = FALSE,
  min.diff = 1,
  common.states = FALSE,
  conditions = NULL,
  tracks2compare = NULL,
  sep = "+",
  statespec = NULL,
  diffstatespec = NULL,
  exclusive.table = NULL,
  binary.matrix = NULL
)

Arguments

Details

The binary modification state (unmodified=0 or modified=1) of multiple ChIP-seq samples defines a (decimal) combinatorial state such as:

Value

A data.frame with combinations and their corresponding (decimal) combinatorial states.

Author(s)

Aaron Taudt, David Widmann

Examples

# Get all combinatorial states where sample1=0, sample2=1, sample3=(0 or 1),
#  sample4=sample5
chromstaR:::state.brewer(statespec=c('0.A','1.B','x.C','r.D','r.D'))

# Get all combinatorial states where sample1=sample2=sample3, sample4=sample5
chromstaR:::state.brewer(statespec=c('r.A','r.A','r.A','r.B','r.B'))

# Get all combinatorial states where sample1=sample5, sample2=sample3=1,
#  sample4=(0 or 1)
chromstaR:::state.brewer(statespec=c('r.A','1.B','1.C','x.D','r.A'))

Description

This function computes combinatorial states from an experiment.table.

Usage

stateBrewer(
  experiment.table,
  mode,
  differential.states = FALSE,
  common.states = FALSE,
  exclusive.table = NULL,
  binary.matrix = NULL
)

Arguments

Details

The binary modification state (unmodified=0 or modified=1) of multiple ChIP-seq samples defines a (decimal) combinatorial state such as:

Value

A data.frame with combinations and their corresponding (decimal) combinatorial states.

Author(s)

Aaron Taudt

Examples

## Construct an experiment table
data(experiment_table)
print(experiment_table)
## Construct combinatorial states
stateBrewer(experiment_table, mode='combinatorial')
stateBrewer(experiment_table, mode='differential')
stateBrewer(experiment_table, mode='full', common.states=TRUE)

## Exclude states with exclusive.table
excl <- data.frame(mark=c('H3K4me3','H3K27me3'),
                             group=c(1,1))
stateBrewer(experiment_table, mode='full', exclusive.table=excl)

Description

Normalize read counts to a given read depth. Reads counts are randomly removed from the input to match the specified read depth.

Usage

subsample(binned.data, sample.reads)

Arguments

Value

A GRanges-class object with downsampled read counts.

Author(s)

Aaron Taudt

Description

Get a table of transition frequencies between combinatorial states of different multiHMMs.

Usage

transitionFrequencies(
  multi.hmms = NULL,
  combined.hmm = NULL,
  zero.states = "[]",
  combstates = NULL
)

Arguments

Value

A data.frame with transition frequencies.

Author(s)

Aaron Taudt

Examples

## Get an example combinedMultiHMM
file <- system.file("data","combined_mode-differential.RData",
                    package="chromstaR")
model <- get(load(file))
freqs <- transitionFrequencies(combined.hmm=model)
freqs$table

Description

The univariate HMM object is output of the function callPeaksUnivariate and is a list() with various entries. The class() attribute of this list was set to "uniHMM". For a given hmm, the entries can be accessed with the list operators 'hmm[[]]' or 'hmm$'.

Value

A list() with the following entries:

See Also

callPeaksUnivariate, multiHMM, combinedMultiHMM

Description

Combine multiple uniHMMs to a multiHMM without running callPeaksMultivariate. This should only be done for comparison purposes.

Usage

unis2pseudomulti(hmms)

Arguments

Details

Use this function if you want to combine ChIP-seq samples without actually running a multivariate Hidden Markov Model. The resulting object will be of class multiHMM but will not be truly multivariate.

Value

A multiHMM object.

Author(s)

Aaron Taudt

Examples

# Get example BAM files for 2 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)]
# 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)
}
## Combine the univariate HMMs without fitting a multivariate HMM
names(models) <- c('H3K27me3','H3K4me3')
pseudo.multi.HMM <- unis2pseudomulti(models)
## Compare frequencies with real multivariate HMM
exp <- data.frame(file=files, mark=c("H3K27me3","H3K4me3"),
                 condition=rep("SHR",2), replicate=c(1,1), pairedEndReads=FALSE,
                 controlFiles=NA)
states <- stateBrewer(exp, mode='combinatorial')
real.multi.HMM <- callPeaksMultivariate(models, use.states=states, eps=1, max.time=60)
genomicFrequencies(real.multi.HMM)
genomicFrequencies(pseudo.multi.HMM)

Description

Make variable-width bins based on a reference BAM file. This can be a simulated file (produced by TODO: insert link and aligned with your favourite aligner) or a real reference.

Usage

variableWidthBins(reads, binsizes, chromosomes = NULL)

Arguments

Details

Variable-width bins are produced by first binning the reference BAM file with fixed-width bins and selecting the desired number of reads per bin as the (non-zero) maximum of the histogram. A new set of bins is then generated such that every bin contains the desired number of reads.

Value

A list() of GRanges-class objects with variable-width bins.

Author(s)

Aaron Taudt

Examples

## Get an example BAM file with ChIP-seq reads
bamfile <- system.file("extdata", "euratrans", "lv-H3K4me3-BN-female-bio1-tech1.bam",
                      package="chromstaRData")
## Read the file into a GRanges object
reads <- readBamFileAsGRanges(bamfile, chromosomes='chr12', pairedEndReads=FALSE,
                    min.mapq=10, remove.duplicate.reads=TRUE)
## Make variable-width bins of size 1000bp
bins <- variableWidthBins(reads, binsizes=1000)
## Plot the distribution of binsizes
hist(width(bins[['1000']]), breaks=50)

Description

Write a chromstaR configuration file from a list structure.

Usage

writeConfig(conf, configfile)

Arguments

Value

NULL

Author(s)

Aaron Taudt

Description

Density, distribution function, quantile function and random generation for the zero-inflated negative binomial distribution with parameters w, size and prob.

Usage

dzinbinom(x, w, size, prob, mu)

pzinbinom(q, w, size, prob, mu, lower.tail = TRUE)

qzinbinom(p, w, size, prob, mu, lower.tail = TRUE)

rzinbinom(n, w, size, prob, mu)

Arguments

Details

The zero-inflated negative binomial distribution with size $= n$ and prob $= p$ has density

$p(x) = w + (1-w) \frac{\Gamma(x+n)}{\Gamma(n) x!} p^n (1-p)^x$

for $x = 0$, $n > 0$, $0 < p \le 1$ and $0 \le w \le 1$.

$p(x) = (1-w) \frac{\Gamma(x+n)}{\Gamma(n) x!} p^n (1-p)^x$

for $x = 1, 2, \ldots$, $n > 0$, $0 < p \le 1$ and $0 \le w \le 1$.

Value

dzinbinom gives the density, pzinbinom gives the distribution function, qzinbinom gives the quantile function, and rzinbinom generates random deviates.

Functions

• dzinbinom: gives the density

• pzinbinom: gives the cumulative distribution function

• qzinbinom: gives the quantile function

• rzinbinom: random number generation

Author(s)

Matthias Heinig, Aaron Taudt

See Also

Distributions for standard distributions, including dbinom for the binomial, dnbinom for the negative binomial, dpois for the Poisson and dgeom for the geometric distribution, which is a special case of the negative binomial.