@@ -1,15 +1,37 @@

 Package: chromstaR

 Type: Package

-Title: Analysis of ChIP-Seq data for single and multiple histone modifications

-Version: 0.9.4

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

+Version: 0.98.0

 Date: 2014-05-03

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

 Maintainer: Aaron Taudt <a.s.taudt@umcg.nl>

-Description: not yet

-Depends: R (>= 3.1.0), S4Vectors, IRanges, GenomicRanges, ggplot2, reshape2

-Imports: Rsamtools, GenomicAlignments, BSgenome, BiocGenerics

-Suggests: grid, rtracklayer, mvtnorm, biomaRt, biovizBase, ggbio,

+Description: This package implements functions for combinatorial and differential analysis of ChIP-seq data. It includes uni- and multivariate peak-calling, export to genome browser viewable files, and functions for enrichment analyses.

+Depends:

+    R (>= 3.3.0),

+    GenomicRanges,

+    ggplot2,

+Imports:

+    utils,

+    grDevices,

+    graphics,

+    stats,

+    S4Vectors,

+		GenomeInfoDb,

+    IRanges,

+    reshape2,

+    Rsamtools,

+    GenomicAlignments,

+    BSgenome,

+    BiocGenerics

+Suggests:

+    grid,

+    rtracklayer,

+    biomaRt,

+    biovizBase,

+    ggbio,

     testthat

-License: GPL

+License: GPL-3

 LazyLoad: yes

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

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

+RoxygenNote: 5.0.1

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

-# Generated by roxygen2 (4.1.1): do not edit by hand

+# Generated by roxygen2: do not edit by hand

 S3method(plot,GRanges)

 S3method(plot,character)

@@ -43,17 +43,37 @@ export(summarizeDiffPeaks)

 export(summarizePeaks)

 export(transitionFrequencies)

 export(unis2pseudomulti)

-import(BiocGenerics)

+import(GenomeInfoDb)

 import(GenomicRanges)

 import(IRanges)

-import(S4Vectors)

 import(ggplot2)

 import(reshape2)

+importFrom(BiocGenerics,unlist)

 importFrom(GenomicAlignments,first)

-importFrom(GenomicAlignments,readGAlignmentPairsFromBam)

-importFrom(GenomicAlignments,readGAlignmentsFromBam)

+importFrom(GenomicAlignments,readGAlignmentPairs)

+importFrom(GenomicAlignments,readGAlignments)

 importFrom(Rsamtools,ScanBamParam)

 importFrom(Rsamtools,indexBam)

 importFrom(Rsamtools,scanBamFlag)

 importFrom(Rsamtools,scanBamHeader)

-useDynLib(chromstaR)

+importFrom(S4Vectors,as.factor)

+importFrom(S4Vectors,queryHits)

+importFrom(S4Vectors,subjectHits)

+importFrom(grDevices,col2rgb)

+importFrom(graphics,hist)

+importFrom(graphics,plot)

+importFrom(stats,aggregate)

+importFrom(stats,cutree)

+importFrom(stats,dist)

+importFrom(stats,dnbinom)

+importFrom(stats,dnorm)

+importFrom(stats,hclust)

+importFrom(stats,pnbinom)

+importFrom(stats,qnorm)

+importFrom(stats,rbinom)

+importFrom(stats,rnbinom)

+importFrom(stats,runif)

+importFrom(stats,stepfun)

+importFrom(utils,read.table)

+importFrom(utils,write.table)

+useDynLib(chromstaR, .registration = TRUE, .fixes = "")

@@ -9,6 +9,7 @@

 #' @param percentages Set to \code{TRUE} if you want to have percentages (0 to 1) instead of fold enrichments returned. Note that in this case different features are not directly comparable.

 #' @return A named array with fold enrichments. If \code{percentages=TRUE} a list with arrays of percentage (0 to 1) enrichments.

 #' @author Aaron Taudt

+#' @importFrom S4Vectors subjectHits queryHits

 #' @export

 foldEnrichment <- function(multi.hmm, featurelist, combinations=NULL, percentages=FALSE) {

@@ -35,11 +36,11 @@ foldEnrichment <- function(multi.hmm, featurelist, combinations=NULL, percentage

 			feature <- featurelist[[ifeat]]

 			ind <- findOverlaps(bins.mask, feature)

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

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

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

 			perc.combstate.in.feature[ifeat,icomb] <- sum.binsinfeature / combstate.length

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

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

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

 			perc.feature.in.combstate[ifeat,icomb] <- sum.featuresinbins / feature.lengths[[ifeat]]

@@ -105,6 +106,7 @@ expressionOverlap <- function(multi.hmm, expression, combinations=NULL, return.m

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

 #' @return A list with vectors of mean expression values per percentile for each combinatorial state. 

 #' @author Aaron Taudt

+#' @importFrom S4Vectors queryHits

 #' @export

 expressionAtPercentageOverlap <- function(multi.hmm, expression, combinations=NULL) {

@@ -122,7 +124,7 @@ expressionAtPercentageOverlap <- function(multi.hmm, expression, combinations=NU

 	for (icomb in 1:length(comb.levels)) {

 		mask <- bins$combination == comb.levels[icomb]

 		ind <- findOverlaps(expression, bins[mask])

-		rle <- rle(queryHits(ind))

+		rle <- rle(S4Vectors::queryHits(ind))

 		expression$num.bins <- 0

 		expression$num.bins[rle$values] <- rle$lengths

 		expression$genewidth <- width(expression)

@@ -206,7 +208,7 @@ transitionFrequencies <- function(multi.hmms=NULL, zero.states="", combstates=NU

 	levels.combstates <- setdiff(levels.combstates, zero.states)

 	levels.combstates <- gsub('-','_',levels.combstates)

 	for (combination in levels.combstates) {

-		mask <- intersect(grep(paste0('\\<',combination,'\\>'), freqtrans$transition), grep(paste(paste0('\\<',setdiff(levels.combstates,combination),'\\>'), collapse='|'), freqtrans$transition, invert=T))

+		mask <- intersect(grep(paste0('\\<',combination,'\\>'), freqtrans$transition), grep(paste(paste0('\\<',setdiff(levels.combstates,combination),'\\>'), collapse='|'), freqtrans$transition, invert=TRUE))

 		freqtrans$group[mask] <- paste0('stage-specific ',combination)

 		mask <- sapply(gregexpr(paste0('\\<',combination,'\\>'),freqtrans$transition), length) == num.models

 		freqtrans$group[mask] <- paste0('constant ', combination)

@@ -21,13 +21,13 @@ bin.annotation = function(annotation.file, reference.genome.file, outputfolder="

 	# Read in the data

 	cat("Reading file",basename(annotation.file),"...")

-	annotation = read.table(annotation.file, sep="\t", header=TRUE)

+	annotation = utils::read.table(annotation.file, sep="\t", header=TRUE)

 	annotation$txDiff = annotation$txEnd - annotation$txStart

 	annotation$cdsDiff = annotation$cdsEnd - annotation$cdsStart

 	cat(" done\n")

 	# Reference genome

-	reference.genome = read.table(reference.genome.file, row.names=1)

+	reference.genome = utils::read.table(reference.genome.file, row.names=1)

 	# Do the loop for all binsizes

 	for (binsize in binsizes) {

@@ -64,6 +64,7 @@ bedGraph2binned <- function(bedGraphfile, assembly, chrom.length.file=NULL, outp

 #' @param remove.duplicate.reads Logical indicating whether duplicate reads should be removed. NOTE: Duplicate removal only works if your duplicate reads are properly flagged in the BAM file.

 #' @param max.fragment.width Maximum allowed fragment length. This is to filter out erroneously wrong fragments due to mapping errors of paired end reads.

 #' @return The function produces a \code{\link{GRanges}} object with one meta data column 'reads' that contains the read count. This binned data will be either written to file (\code{save.as.RData=FALSE}) or given as return value (\code{save.as.RData=FALSE}).

+#' @importFrom utils read.table

 align2binned <- function(file, format, assembly, bamindex=file, pairedEndReads=FALSE, chrom.length.file=NULL, outputfolder="binned_data", binsizes=500, chromosomes=NULL, save.as.RData=FALSE, min.mapq=10, downsample.to.reads=NULL, remove.duplicate.reads=FALSE, max.fragment.width=1000) {

 	## Create outputfolder if not exists

@@ -72,25 +73,25 @@ align2binned <- function(file, format, assembly, bamindex=file, pairedEndReads=F

 	}

 	### Read in the data

-	message("Reading file ",basename(file)," ...", appendLF=F); ptm <- proc.time()

+	ptm <- startTimedMessage("Reading file ",basename(file)," ...")

 	## BED (0-based)

 	if (format == "bed") {

 		if (!is.null(chrom.length.file)) {

 			# File with chromosome lengths (1-based)

-			chrom.lengths.df <- read.table(chrom.length.file)

+			chrom.lengths.df <- utils::read.table(chrom.length.file)

 			chrom.lengths <- chrom.lengths.df[,2]

 			names(chrom.lengths) <- chrom.lengths.df[,1]

 		} else {

 			chrom.lengths <- getChromLengths(assembly=assembly)

 		}

 		# File with reads, determine classes first for faster import (0-based)

-		tab5rows <- read.table(file, nrows=5)

+		tab5rows <- utils::read.table(file, nrows=5)

 		classes.in.bed <- sapply(tab5rows, class)

 		classes <- rep("NULL",length(classes.in.bed))

 		classes[1:3] <- classes.in.bed[1:3]

-		data <- read.table(file, colClasses=classes)

+		data <- utils::read.table(file, colClasses=classes)

 		# Convert to GRanges object

-		data <- GenomicRanges::GRanges(seqnames=Rle(data[,1]), ranges=IRanges(start=data[,2]+1, end=data[,3]), strand=Rle(strand("*"), nrow(data)))	# start+1 to go from [0,x) -> [1,x]

+		data <- GenomicRanges::GRanges(seqnames=data[,1], ranges=IRanges(start=data[,2]+1, end=data[,3]), strand="*")	# start+1 to go from [0,x) -> [1,x]

 		seqlengths(data) <- as.integer(chrom.lengths[names(seqlengths(data))])

 	## BAM (1-based)

 	} else if (format == "bam") {

@@ -99,23 +100,23 @@ align2binned <- function(file, format, assembly, bamindex=file, pairedEndReads=F

 	} else if (format == "bedGraph") {

 		if (!is.null(chrom.length.file)) {

 			# File with chromosome lengths (1-based)

-			chrom.lengths.df <- read.table(chrom.length.file)

+			chrom.lengths.df <- utils::read.table(chrom.length.file)

 			chrom.lengths <- chrom.lengths.df[,2]

 			names(chrom.lengths) <- chrom.lengths.df[,1]

 		} else {

 			chrom.lengths <- getChromLengths(assembly=assembly)

 		}

 		# File with reads, determine classes first for faster import

-		tab5rows <- read.table(file, nrows=5)

+		tab5rows <- utils::read.table(file, nrows=5)

 		classes.in.bed <- sapply(tab5rows, class)

 		classes <- rep("NULL",length(classes.in.bed))

 		classes[1:4] <- classes.in.bed[1:4]

-		data <- read.table(file, colClasses=classes)

+		data <- utils::read.table(file, colClasses=classes)

 		# Convert to GRanges object

-		data <- GenomicRanges::GRanges(seqnames=Rle(data[,1]), ranges=IRanges(start=data[,2]+1, end=data[,3]), strand=Rle(strand("*"), nrow(data)), signal=data[,4])	# start+1 to go from [0,x) -> [1,x]

+		data <- GenomicRanges::GRanges(seqnames=data[,1], ranges=IRanges(start=data[,2]+1, end=data[,3]), strand="*", signal=data[,4])	# start+1 to go from [0,x) -> [1,x]

 		seqlengths(data) <- as.integer(chrom.lengths[names(seqlengths(data))])

 	}

-	time <- proc.time() - ptm; message(" ",round(time[3],2),"s")

+	stopTimedMessage(ptm)

 	## Select chromosomes to bin

 	if (is.null(chromosomes)) {

@@ -126,9 +127,9 @@ align2binned <- function(file, format, assembly, bamindex=file, pairedEndReads=F

 	### Downsampling

 	if (!is.null(downsample.to.reads)) {

 		if (downsample.to.reads <= length(data)) {

-			message("downsampling to ", downsample.to.reads, " reads ...", appendLF=F); ptm <- proc.time()

+			ptm <- startTimedMessage("downsampling to ", downsample.to.reads, " reads ...")

 			data <- data[sort(sample(1:length(data), size=downsample.to.reads, replace=FALSE))]

-			time <- proc.time() - ptm; message(" ",round(time[3],2),"s")

+			stopTimedMessage(ptm)

 		} else {

 			warning("No downsampling done because there are only ", length(data), " reads in the data.")

 		}

@@ -139,7 +140,7 @@ align2binned <- function(file, format, assembly, bamindex=file, pairedEndReads=F

 		message("Binning into bin size ",binsize)

 		### Iterate over all chromosomes

-		message("  binning genome ...", appendLF=F); ptm <- proc.time()

+		ptm <- startTimedMessage("  binning genome ...")

 		binned.data <- GenomicRanges::GRangesList()

 		for (chromosome in chroms2use) {

 			## Check last incomplete bin

@@ -161,9 +162,9 @@ align2binned <- function(file, format, assembly, bamindex=file, pairedEndReads=F

 # 			end[length(end)] <- seqlengths(data)[chromosome] # last ending coordinate is size of chromosome, only if incomplete bins are desired

 			## Create binned chromosome as GRanges object

-			i.binned.data <- GenomicRanges::GRanges(seqnames = Rle(chromosome, numbins),

+			i.binned.data <- GenomicRanges::GRanges(seqnames = rep(chromosome, numbins),

 							ranges = IRanges(start=start, end=end),

-							strand = Rle(strand("*"), numbins)

+							strand = "*"

 							)

 			seqlengths(i.binned.data) <- seqlengths(data)[chromosome]

@@ -171,12 +172,12 @@ align2binned <- function(file, format, assembly, bamindex=file, pairedEndReads=F

 				binned.data[[chromosome]] <- i.binned.data

 			)

 		}

-		time <- proc.time() - ptm; message(" ",round(time[3],2),"s")

+		stopTimedMessage(ptm)

 		binned.data <- unlist(binned.data)

 		names(binned.data) <- NULL

 		## Count overlaps

-		message("  counting overlaps ...", appendLF=F); ptm <- proc.time()

+		ptm <- startTimedMessage("  counting overlaps ...")

 		if (format=="bam" | format=="bed") {

 			reads <- GenomicRanges::countOverlaps(binned.data, data)

 		} else if (format=="bedGraph") {

@@ -194,8 +195,8 @@ align2binned <- function(file, format, assembly, bamindex=file, pairedEndReads=F

 			}

 			reads[read.idx] <- maxvalues

 		}

-		mcols(binned.data)$reads <- reads

-		time <- proc.time() - ptm; message(" ",round(time[3],2),"s")

+		mcols(binned.data)$counts <- reads

+		stopTimedMessage(ptm)

 		if (length(binned.data) == 0) {

 			warning(paste0("The bin size of ",binsize," is larger than any of the chromosomes."))

@@ -204,10 +205,10 @@ align2binned <- function(file, format, assembly, bamindex=file, pairedEndReads=F

 		if (save.as.RData==TRUE) {

 			# Print to file

-			filename <- paste0(basename(file),"_binsize",format(binsize, scientific=F, trim=T),".RData")

-			message("Saving to file ...", appendLF=F); ptm <- proc.time()

+			filename <- paste0(basename(file),"_binsize",format(binsize, scientific=F, trim=TRUE),".RData")

+			ptm <- startTimedMessage("Saving to file ...")

 			save(binned.data, file=file.path(outputfolder,filename) )

-			time <- proc.time() - ptm; message(" ",round(time[3],2),"s")

+			stopTimedMessage(ptm)

 		} else {

 			return(binned.data)

 		}

@@ -229,7 +230,7 @@ align2binned <- function(file, format, assembly, bamindex=file, pairedEndReads=F

 #' @param min.mapq Minimum mapping quality when importing from BAM files. Set \code{min.mapq=NULL} to keep all reads.

 #' @param max.fragment.width Maximum allowed fragment length. This is to filter out erroneously wrong fragments due to mapping errors of paired end reads.

 #' @importFrom Rsamtools indexBam scanBamHeader ScanBamParam scanBamFlag

-#' @importFrom GenomicAlignments readGAlignmentPairsFromBam readGAlignmentsFromBam first

+#' @importFrom GenomicAlignments readGAlignmentPairs readGAlignments first

 bam2GRanges <- function(file, bamindex=file, chromosomes=NULL, pairedEndReads=FALSE, keep.duplicate.reads=TRUE, min.mapq=10, max.fragment.width=1000) {

 	## Check if bamindex exists

@@ -258,25 +259,25 @@ bam2GRanges <- function(file, bamindex=file, chromosomes=NULL, pairedEndReads=FA

 		warning(paste0('Not using chromosomes ', diffs, ' because they are not in the data.'))

 	}

 	## Import the file into GRanges

-	gr <- GenomicRanges::GRanges(seqnames=Rle(chroms2use), ranges=IRanges(start=rep(1, length(chroms2use)), end=chrom.lengths[chroms2use]))

+	gr <- GenomicRanges::GRanges(seqnames=chroms2use, ranges=IRanges(start=rep(1, length(chroms2use)), end=chrom.lengths[chroms2use]))

 	if (keep.duplicate.reads) {

 		if (pairedEndReads) {

-			data.raw <- GenomicAlignments::readGAlignmentPairsFromBam(file, index=bamindex, param=Rsamtools::ScanBamParam(which=range(gr), what='mapq'))

+			data.raw <- GenomicAlignments::readGAlignmentPairs(file, index=bamindex, param=Rsamtools::ScanBamParam(which=range(gr), what='mapq'))

 		} else {

-			data.raw <- GenomicAlignments::readGAlignmentsFromBam(file, index=bamindex, param=Rsamtools::ScanBamParam(which=range(gr), what='mapq'))

+			data.raw <- GenomicAlignments::readGAlignments(file, index=bamindex, param=Rsamtools::ScanBamParam(which=range(gr), what='mapq'))

 		}

 	} else {

 		if (pairedEndReads) {

-			data.raw <- GenomicAlignments::readGAlignmentPairsFromBam(file, index=bamindex, param=Rsamtools::ScanBamParam(which=range(gr), what='mapq', flag=scanBamFlag(isDuplicate=F)))

+			data.raw <- GenomicAlignments::readGAlignmentPairs(file, index=bamindex, param=Rsamtools::ScanBamParam(which=range(gr), what='mapq', flag=scanBamFlag(isDuplicate=FALSE)))

 		} else {

-			data.raw <- GenomicAlignments::readGAlignmentsFromBam(file, index=bamindex, param=Rsamtools::ScanBamParam(which=range(gr), what='mapq', flag=scanBamFlag(isDuplicate=F)))

+			data.raw <- GenomicAlignments::readGAlignments(file, index=bamindex, param=Rsamtools::ScanBamParam(which=range(gr), what='mapq', flag=scanBamFlag(isDuplicate=FALSE)))

 		}

 	}

 	## Filter by mapping quality

 	if (pairedEndReads) {

 		message("Converting to GRanges ...", appendLF=FALSE); ptm <- proc.time()

 		data <- as(data.raw, 'GRanges') # treat as one fragment

-		time <- proc.time() - ptm; message(" ",round(time[3],2),"s")

+		stopTimedMessage(ptm)

 		message("Filtering reads ...", appendLF=FALSE); ptm <- proc.time()

 		if (!is.null(min.mapq)) {

@@ -290,11 +291,11 @@ bam2GRanges <- function(file, bamindex=file, chromosomes=NULL, pairedEndReads=FA

 			# Filter out too long fragments

 			data <- data[width(data)<=max.fragment.width]

 		}

-		time <- proc.time() - ptm; message(" ",round(time[3],2),"s")

+		stopTimedMessage(ptm)

 	} else {

 		message("Converting to GRanges ...", appendLF=FALSE); ptm <- proc.time()

 		data <- as(data.raw, 'GRanges')

-		time <- proc.time() - ptm; message(" ",round(time[3],2),"s")

+		stopTimedMessage(ptm)

 		message("Filtering reads ...", appendLF=FALSE); ptm <- proc.time()

 		if (!is.null(min.mapq)) {

@@ -304,7 +305,7 @@ bam2GRanges <- function(file, bamindex=file, chromosomes=NULL, pairedEndReads=FA

 			}

 			data <- data[mcols(data)$mapq >= min.mapq]

 		}

-		time <- proc.time() - ptm; message(" ",round(time[3],2),"s")

+		stopTimedMessage(ptm)

 	}

 	return(data)

@@ -87,7 +87,7 @@ callPeaksMultivariate <- function(modellist, use.states, num.states=NULL, chromo

 	## Prepare the HMM

 	params <- prepare.multivariate(modellist, use.states=use.states, num.states=num.states)

 	bins <- params$bins

-	reads <- params$reads

+	reads <- params$counts

 	numbins <- params$numbins

 	nummod <- params$nummod

 	comb.states2use <- params$comb.states

@@ -113,11 +113,11 @@ callPeaksMultivariate <- function(modellist, use.states, num.states=NULL, chromo

 	}

 	### Define cleanup behaviour ###

-	on.exit(.C("R_multivariate_cleanup", as.integer(nummod)))

+	on.exit(.C("C_multivariate_cleanup", as.integer(nummod)))

 	## Starting multivariate HMM

 	message("\nStarting multivariate HMM")

-	message("Using the following combinatorial states, covering ", mean(comb.states.per.bin %in% comb.states2use)*100, "% of the bins:\n", paste(comb.states2use, collapse=" "),"\n", appendLF=F)

+	message("Using the following combinatorial states, covering ", mean(comb.states.per.bin %in% comb.states2use)*100, "% of the bins:\n", paste(comb.states2use, collapse=" "),"\n", appendLF=FALSE)

 	# Prepare input for C function

 	rs <- unlist(lapply(distributions,"[",2:3,'size'))

@@ -176,7 +176,7 @@ callPeaksMultivariate <- function(modellist, use.states, num.states=NULL, chromo

 			max.time.temp <- checkpoint.after.time

 		}

 		# Call the C function

-		hmm <- .C("R_multivariate_hmm",

+		hmm <- .C("C_multivariate_hmm",

 			reads = as.integer(as.vector(reads)), # int* multiO

 			num.bins = as.integer(numbins), # int* T

 			num.states = as.integer(numstates2use), # int* N

@@ -218,9 +218,9 @@ callPeaksMultivariate <- function(modellist, use.states, num.states=NULL, chromo

 			result$IDs <- IDs

 		## Bin coordinates, posteriors and states

 			result$bins <- bins

-			result$bins$reads <- reads

+			result$bins$counts <- reads

 			if (get.posteriors) {

-				message("Transforming posteriors to `per sample` representation ...", appendLF=F); ptm <- proc.time()

+				ptm <- startTimedMessage("Transforming posteriors to `per sample` representation ...")

 				hmm$posteriors <- matrix(hmm$posteriors, ncol=hmm$num.states)

 				colnames(hmm$posteriors) <- paste0("P(",hmm$comb.states,")")

 				result$bins$score <- do.call(pmax,as.list(as.data.frame(hmm$posteriors)))

@@ -228,9 +228,9 @@ callPeaksMultivariate <- function(modellist, use.states, num.states=NULL, chromo

 				post.per.track <- hmm$posteriors %*% binstates

 				colnames(post.per.track) <- result$IDs

 				result$bins$posteriors <- post.per.track

-				time <- proc.time() - ptm; message(" ",round(time[3],2),"s")

+				stopTimedMessage(ptm)

 			}

-			message("Calculating states from posteriors ...", appendLF=F); ptm <- proc.time()

+			ptm <- startTimedMessage("Calculating states from posteriors ...")

 			if (!is.null(use.states$state)) {

 				state.levels <- levels(use.states$state)

 			} else {

@@ -241,14 +241,14 @@ callPeaksMultivariate <- function(modellist, use.states, num.states=NULL, chromo

 			} else {

 				result$bins$state <- factor(hmm$states, levels=state.levels)

 			}

-			time <- proc.time() - ptm; message(" ",round(time[3],2),"s")

+			stopTimedMessage(ptm)

 			if (keep.densities) {

 				result$bins$densities <- matrix(hmm$densities, ncol=hmm$num.states)

 				colnames(result$bins$densities) <- hmm$comb.states

 			}

 		## Segmentation

-			message("Making segmentation ...", appendLF=F); ptm <- proc.time()

+			ptm <- startTimedMessage("Making segmentation ...")

 			df <- as.data.frame(result$bins)

 			ind.readcols <- grep('^reads', names(df))

 			ind.postcols <- grep('^posteriors', names(df))

@@ -265,7 +265,7 @@ callPeaksMultivariate <- function(modellist, use.states, num.states=NULL, chromo

 			if (!keep.posteriors) {

 				result$bins$posteriors <- NULL

 			}

-			time <- proc.time() - ptm; message(" ",round(time[3],2),"s")

+			stopTimedMessage(ptm)

 		## Add combinations

 			if (!is.null(use.states)) {

 				mapping <- use.states$combination

@@ -277,7 +277,7 @@ callPeaksMultivariate <- function(modellist, use.states, num.states=NULL, chromo

 		## Parameters

 			# Weights

 			tstates <- table(hmm$states)

-			result$weights <- sort(tstates/sum(tstates), decreasing=T)

+			result$weights <- sort(tstates/sum(tstates), decreasing=TRUE)

 			result$weights.univariate <- weights

 			# Transition matrices

 			result$transitionProbs <- matrix(hmm$A, ncol=numstates2use, byrow=TRUE)

@@ -14,6 +14,7 @@

 #' @param max.distance This number is used as a cutoff to group replicates based on their distance matrix. The lower this number, the more similar replicates have to be to be grouped together.

 #' @return Output is a \code{\link{chromstaR_multivariateHMM}} object with additional entry \code{replicateInfo}. If only one \code{\link{chromstaR_univariateHMM}} was given as input, a simple list() with the \code{replicateInfo} is returned.

 #' @seealso \code{\link{chromstaR_multivariateHMM}}, \code{\link{callPeaksUnivariate}}, \code{\link{callPeaksMultivariate}}

+#' @importFrom stats hclust cutree dist

 #' @export

 callPeaksReplicates <- function(hmm.list, max.states=32, force.equal=FALSE, eps=0.01, max.iter=NULL, max.time=NULL, keep.posteriors=FALSE, num.threads=1, max.distance=0.2) {

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

 		info.df <- multimodel$replicateInfo$info

 		cor.matrix <- multimodel$replicateInfo$correlation

 		dist.matrix <- multimodel$replicateInfo$distance

-		hc <- hclust(dist.matrix)

-		info.df$group <- cutree(hc, h=max.distance)

+		hc <- stats::hclust(dist.matrix)

+		info.df$group <- stats::cutree(hc, h=max.distance)

 	### Call peaks for several replicates ###

 	} else {

@@ -39,7 +40,7 @@ callPeaksReplicates <- function(hmm.list, max.states=32, force.equal=FALSE, eps=

 		## Univariate replicateInfo

 		ids <- unlist(lapply(hmms, '[[', 'ID'))

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

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

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

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

 		rownames(info.df) <- ids

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

 			cor.matrix <- cor(binstates)

 			weight.multivariate <- apply(binstates, 2, sum) / nrow(binstates)

 			info.df$weight.multivariate <- weight.multivariate

-			dist.matrix <- dist(cor.matrix)

-			hc <- hclust(dist.matrix)

-			info.df$group <- cutree(hc, h=max.distance)

+			dist.matrix <- stats::dist(cor.matrix)

+			hc <- stats::hclust(dist.matrix)

+			info.df$group <- stats::cutree(hc, h=max.distance)

 		}

 		## Make return object

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

 #' @param eps Convergence threshold for the Baum-Welch algorithm.

 #' @param init One of the following initialization procedures:

 #' \describe{

-#' \item{\code{standard}}{The negative binomial of state 'unmodified' will be initialized with \code{mean=mean(reads)}, \code{var=var(reads)} and the negative binomial of state 'modified' with \code{mean=mean(reads)+1}, \code{var=var(reads)}. This procedure usually gives the fastest convergence.}

+#' \item{\code{standard}}{The negative binomial of state 'unmodified' will be initialized with \code{mean=mean(counts)}, \code{var=var(counts)} and the negative binomial of state 'modified' with \code{mean=mean(counts)+1}, \code{var=var(counts)}. This procedure usually gives the fastest convergence.}

 #' \item{\code{random}}{Mean and variance of the negative binomials will be initialized with random values (in certain boundaries, see source code). Try this if the \code{'standard'} procedure fails to produce a good fit.}

 #' \item{\code{empiric}}{Yet another way to initialize the Baum-Welch. Try this if the other two methods fail to produce a good fit.}

 #' }

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

 #' @param read.cutoff The default (\code{TRUE}) enables filtering of high read counts. Set \code{read.cutoff=FALSE} to disable this filtering.

 #' @param read.cutoff.quantile A quantile between 0 and 1. Should be near 1. Read counts above this quantile will be set to the read count specified by this quantile. Filtering very high read counts increases the performance of the Baum-Welch fitting procedure. However, if your data contains very few peaks they might be filtered out. If option \code{read.cutoff.absolute} is also specified, the minimum of the resulting cutoff values will be used. Set \code{read.cutoff=FALSE} to disable this filtering.

 #' @param read.cutoff.absolute Read counts above this value will be set to the read count specified by this value. Filtering very high read counts increases the performance of the Baum-Welch fitting procedure. However, if your data contains very few peaks they might be filtered out. If option \code{read.cutoff.quantile} is also specified, the minimum of the resulting cutoff values will be used. Set \code{read.cutoff=FALSE} to disable this filtering.

-#' @param max.mean If \code{mean(reads)>max.mean}, bins with low read counts will be set to 0. This is a workaround to obtain good fits in the case of large bin sizes.

+#' @param max.mean If \code{mean(counts)>max.mean}, bins with low read counts will be set to 0. This is a workaround to obtain good fits in the case of large bin sizes.

 #' @param FDR False discovery rate. code{NULL} means that the state with maximum posterior probability will be chosen, irrespective of its absolute probability (default=code{NULL}).

 #' @param control If set to \code{TRUE}, the binned data will be treated as control experiment. That means only state 'zero-inflation' and 'unmodified' will be used in the HMM.

 #' @param keep.posteriors If set to \code{TRUE} (default=\code{FALSE}), posteriors will be available in the output. This is useful to change the FDR later, but increases the necessary disk space to store the result.

@@ -75,7 +75,7 @@ callPeaksUnivariate2 <- function(binned.data, ID, prefit.on.chr, short=TRUE, eps

 #' @param eps Convergence threshold for the Baum-Welch algorithm.

 #' @param init One of the following initialization procedures:

 #' \describe{

-#' \item{\code{standard}}{The negative binomial of state 'unmodified' will be initialized with \code{mean=mean(reads)}, \code{var=var(reads)} and the negative binomial of state 'modified' with \code{mean=mean(reads)+1}, \code{var=var(reads)}. This procedure usually gives the fastest convergence.}

+#' \item{\code{standard}}{The negative binomial of state 'unmodified' will be initialized with \code{mean=mean(counts)}, \code{var=var(counts)} and the negative binomial of state 'modified' with \code{mean=mean(counts)+1}, \code{var=var(counts)}. This procedure usually gives the fastest convergence.}

 #' \item{\code{random}}{Mean and variance of the negative binomials will be initialized with random values (in certain boundaries, see source code). Try this if the \code{'standard'} procedure fails to produce a good fit.}

 #' \item{\code{empiric}}{Yet another way to initialize the Baum-Welch. Try this if the other two methods fail to produce a good fit.}

 #' }

@@ -87,7 +87,7 @@ callPeaksUnivariate2 <- function(binned.data, ID, prefit.on.chr, short=TRUE, eps

 #' @param read.cutoff The default (\code{TRUE}) enables filtering of high read counts. Set \code{read.cutoff=FALSE} to disable this filtering.

 #' @param read.cutoff.quantile A quantile between 0 and 1. Should be near 1. Read counts above this quantile will be set to the read count specified by this quantile. Filtering very high read counts increases the performance of the Baum-Welch fitting procedure. However, if your data contains very few peaks they might be filtered out. If option \code{read.cutoff.absolute} is also specified, the minimum of the resulting cutoff values will be used. Set \code{read.cutoff=FALSE} to disable this filtering.

 #' @param read.cutoff.absolute Read counts above this value will be set to the read count specified by this value. Filtering very high read counts increases the performance of the Baum-Welch fitting procedure. However, if your data contains very few peaks they might be filtered out. If option \code{read.cutoff.quantile} is also specified, the minimum of the resulting cutoff values will be used. Set \code{read.cutoff=FALSE} to disable this filtering.

-#' @param max.mean If \code{mean(reads)>max.mean}, bins with low read counts will be set to 0. This is a workaround to obtain good fits in the case of large bin sizes.

+#' @param max.mean If \code{mean(counts)>max.mean}, bins with low read counts will be set to 0. This is a workaround to obtain good fits in the case of large bin sizes.

 #' @param FDR False discovery rate. code{NULL} means that the state with maximum posterior probability will be chosen, irrespective of its absolute probability (default=code{NULL}).

 #' @param control If set to \code{TRUE}, the binned data will be treated as control experiment. That means only state 'zero-inflation' and 'unmodified' will be used in the HMM.

 #' @param keep.posteriors If set to \code{TRUE} (default=\code{FALSE}), posteriors will be available in the output. This is useful to change the FDR later, but increases the necessary disk space to store the result.

@@ -100,6 +100,7 @@ callPeaksUnivariate2 <- function(binned.data, ID, prefit.on.chr, short=TRUE, eps

 #' @param verbosity Verbosity level for the fitting procedure. 0 - No output, 1 - Iterations are printed.

 #' @author Aaron Taudt, Maria Coome Tatche

 #' @seealso \code{\link{chromstaR_univariateHMM}}, \code{\link{callPeaksMultivariate}}

+#' @importFrom stats runif

 #' @examples

 #'## Get an example BED-file with ChIP-seq reads for H3K36me3 in brain tissue

 #'bedfile <- system.file("extdata/brain/BI.Brain_Angular_Gyrus.H3K36me3.112.chr22.bed.gz",

@@ -114,7 +115,7 @@ callPeaksUnivariate2 <- function(binned.data, ID, prefit.on.chr, short=TRUE, eps

 callPeaksUnivariate <- function(binned.data, ID, 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, FDR=0.5, control=FALSE, keep.posteriors=FALSE, keep.densities=FALSE, checkpoint.after.iter=NULL, checkpoint.after.time=NULL, checkpoint.file=paste0('chromstaR_checkpoint_',ID,'.cpt'), checkpoint.overwrite=TRUE, checkpoint.use.existing=FALSE, verbosity=1) {

 	### Define cleanup behaviour ###

-	on.exit(.C("R_univariate_cleanup"))

+	on.exit(.C("C_univariate_cleanup"))

 	### Intercept user input ###

 	IDcheck <- ID  #trigger error if not defined

@@ -156,7 +157,7 @@ callPeaksUnivariate <- function(binned.data, ID, eps=0.01, init="standard", max.

 		continue.from.univariate.hmm <- TRUE

 	} else if (class(binned.data) == 'GRanges') {

 	} else {

-		stop("argument 'binned.data' expects a GRanges with meta-column 'reads' or a file that contains such an object")

+		stop("argument 'binned.data' expects a GRanges with meta-column 'counts' or a file that contains such an object")

 	}

 	### Assign variables ###

@@ -165,29 +166,29 @@ callPeaksUnivariate <- function(binned.data, ID, eps=0.01, init="standard", max.

 	}

 	numstates <- length(state.labels)

 	numbins <- length(binned.data)

-	reads <- mcols(binned.data)$reads

+	counts <- mcols(binned.data)$counts

 	iniproc <- which(init==c("standard","random","empiric")) # transform to int

 	if (keep.densities) { lenDensities <- numbins * numstates } else { lenDensities <- 1 }

-	### Check if there are reads in the data, otherwise HMM will blow up ###

-	if (!any(reads!=0)) {

-		stop("All reads in data are zero. No univariate HMM done.")

+	### Check if there are counts in the data, otherwise HMM will blow up ###

+	if (!any(counts!=0)) {

+		stop("All counts in data are zero. No univariate HMM done.")

 	}

-	### Filter high reads out, makes HMM faster ###

+	### Filter high counts out, makes HMM faster ###

 	numfiltered <- 0

 	if (!continue.from.univariate.hmm) {

 		if (read.cutoff) {

-			read.cutoff.by.quantile <- quantile(reads, read.cutoff.quantile)

+			read.cutoff.by.quantile <- quantile(counts, read.cutoff.quantile)

 			read.cutoff.by.quantile <- as.integer(read.cutoff.by.quantile)

 			read.cutoff.absolute <- min(read.cutoff.by.quantile, read.cutoff.absolute)

-			mask <- reads > read.cutoff.absolute

-			reads[mask] <- read.cutoff.absolute

+			mask <- counts > read.cutoff.absolute

+			counts[mask] <- read.cutoff.absolute

 			numfiltered <- length(which(mask))

 		}

 	} else {

-		mask <- reads > read.cutoff.absolute

-		reads[mask] <- read.cutoff.absolute

+		mask <- counts > read.cutoff.absolute

+		counts[mask] <- read.cutoff.absolute

 		numfiltered <- length(which(mask))

 	}

 	if (numfiltered > 0) {

@@ -195,29 +196,29 @@ callPeaksUnivariate <- function(binned.data, ID, eps=0.01, init="standard", max.

 	}

 	### Filter out low read counts that arise when the bin size is larger than optimal (should correct the result to near optimal again) ###

-	hist <- hist(reads[reads>0], breaks=0:max(reads), right=FALSE, plot=FALSE)

+	hist <- graphics::hist(counts[counts>0], breaks=0:max(counts), right=FALSE, plot=FALSE)

 	maxhist <- which.max(hist$counts)

 	if (maxhist-1 > max.mean) {	# -1 to get from 1-based histogram indices to (0-based) read counts

-		# Two empirical rules to remove low reads

+		# Two empirical rules to remove low counts

 		read.counts.to.remove.1 <- which(hist$counts[1:maxhist]<=hist$counts[2]) -1

 		minlow <- which.min(hist$counts[2:maxhist])

 		read.counts.to.remove <- max(c(read.counts.to.remove.1, 2*minlow))

-		index.filtered <- which(reads>0 & reads<=read.counts.to.remove)

-		reads[index.filtered] <- 0

+		index.filtered <- which(counts>0 & counts<=read.counts.to.remove)

+		counts[index.filtered] <- 0

 		if (length(index.filtered)>0) {

 			message(paste0("Replaced read counts <= ",read.counts.to.remove," by 0. This was done because the selected bin size is considered too big for this dataset: The mean of the read counts (zeros removed) is bigger than the specified max.mean = ",max.mean,". Check the fits!"))

 		}

 	}

 # 	### Initialization of emission distributions ###

-# 	reads.greater.0 <- reads[reads>0]

-# 	mean.reads <- mean(reads.greater.0)

+# 	reads.greater.0 <- counts[counts>0]

+# 	mean.counts <- mean(reads.greater.0)

 # 	var.reads <- var(reads.greater.0)

 # 	imean <- vector()

 # 	ivar <- vector()

 # 	if (init=="standard") {

-# 		imean['unmodified'] <- mean.reads

-# 		imean['modified'] <- mean.reads + 1

+# 		imean['unmodified'] <- mean.counts

+# 		imean['modified'] <- mean.counts + 1

 # 		ivar['unmodified'] <- var.reads

 # 		ivar['modified'] <- var.reads

 # 		# Make sure variance is bigger than mean for negative binomial

@@ -228,13 +229,13 @@ callPeaksUnivariate <- function(binned.data, ID, eps=0.01, init="standard", max.

 # 		}

 # 	} else if (init=="random") {

 # 		for (istate in c('unmodified','modified')) {

-# 			imean[istate] <- runif(1, min=0, max=10*mean.reads)

-# 			ivar[istate] <- runif(1, min=imean[istate], max=20*mean.reads)

+# 			imean[istate] <- stats::runif(1, min=0, max=10*mean.counts)

+# 			ivar[istate] <- stats::runif(1, min=imean[istate], max=20*mean.counts)

 # 		}

 # 	} else if (init=="empiric") {

-# 		imean['unmodified'] <- mean.reads/2

+# 		imean['unmodified'] <- mean.counts/2

 # 		ivar['unmodified'] <- imean['unmodified']*2

-# 		imean['modified'] <- mean.reads*2

+# 		imean['modified'] <- mean.counts*2

 # 		ivar['modified'] <- imean['modified']*2

 # 	} else if (init %in% seqlevels(binned.data)) {

 # 	} else if {

@@ -282,8 +283,8 @@ callPeaksUnivariate <- function(binned.data, ID, eps=0.01, init="standard", max.

 				max.time.temp <- checkpoint.after.time

 			}

 			## Call the Baum-Welch

-			hmm <- .C("R_univariate_hmm",

-				reads = as.integer(reads), # double* O

+			hmm <- .C("C_univariate_hmm",

+				counts = as.integer(counts), # double* O

 				num.bins = as.integer(numbins), # int* T

 				num.states = as.integer(numstates), # int* N

 				size = double(length=numstates), # double* size

@@ -306,7 +307,7 @@ callPeaksUnivariate <- function(binned.data, ID, eps=0.01, init="standard", max.

 				use.initial.params = as.logical(use.initial), # bool* use_initial_params

 				num.threads = as.integer(num.threads), # int* num_threads

 				error = as.integer(0), # int* error (error handling)

-				read.cutoff = as.integer(max(reads)), # int* read_cutoff

+				read.cutoff = as.integer(max(counts)), # int* read_cutoff

 				verbosity = as.integer(verbosity) # int* verbosity

 			)

 			## Adjust parameters for the next round

@@ -353,7 +354,7 @@ callPeaksUnivariate <- function(binned.data, ID, eps=0.01, init="standard", max.

 				# FDR

 				result$FDR <- FDR

 			## Convergence info

-				convergenceInfo <- list(eps=eps, loglik=hmm$loglik, loglik.delta=hmm$loglik.delta, num.iterations=hmm$num.iterations, time.sec=hmm$time.sec, read.cutoff=max(hmm$reads))

+				convergenceInfo <- list(eps=eps, loglik=hmm$loglik, loglik.delta=hmm$loglik.delta, num.iterations=hmm$num.iterations, time.sec=hmm$time.sec, read.cutoff=max(hmm$counts))

 				result$convergenceInfo <- convergenceInfo

 			## Add class

 				class(result) <- class.univariate.hmm

@@ -379,8 +380,8 @@ callPeaksUnivariate <- function(binned.data, ID, eps=0.01, init="standard", max.

 		modellist <- list()

 		for (i_try in 1:num.trials) {

 			if (verbosity>=1) message("------------------------------------ Try ",i_try," of ",num.trials," -------------------------------------")

-			hmm <- .C("R_univariate_hmm",

-				reads = as.integer(reads), # double* O

+			hmm <- .C("C_univariate_hmm",

+				counts = as.integer(counts), # double* O

 				num.bins = as.integer(numbins), # int* T

 				num.states = as.integer(numstates), # int* N

 				size = double(length=numstates), # double* size

@@ -403,7 +404,7 @@ callPeaksUnivariate <- function(binned.data, ID, eps=0.01, init="standard", max.

 				use.initial.params = as.logical(0), # bool* use_initial_params

 				num.threads = as.integer(num.threads), # int* num_threads

 				error = as.integer(0), # int* error (error handling)

-				read.cutoff = as.integer(max(reads)), # int* read_cutoff

+				read.cutoff = as.integer(max(counts)), # int* read_cutoff

 				verbosity = as.integer(verbosity) # int* verbosity

 			)

@@ -424,8 +425,8 @@ callPeaksUnivariate <- function(binned.data, ID, eps=0.01, init="standard", max.

 		# Rerun the HMM with different epsilon and initial parameters from trial run

 		if (verbosity>=1) message("------------------------- Rerunning try ",indexmax," with eps = ",eps," -------------------------")

-		hmm <- .C("R_univariate_hmm",

-			reads = as.integer(reads), # double* O

+		hmm <- .C("C_univariate_hmm",

+			counts = as.integer(counts), # double* O

 			num.bins = as.integer(numbins), # int* T

 			num.states = as.integer(numstates), # int* N

 			size = double(length=numstates), # double* size

@@ -448,7 +449,7 @@ callPeaksUnivariate <- function(binned.data, ID, eps=0.01, init="standard", max.

 			use.initial.params = as.logical(1), # bool* use_initial_params

 			num.threads = as.integer(num.threads), # int* num_threads

 			error = as.integer(0), # int* error (error handling)

-			read.cutoff = as.integer(max(reads)), # int* read_cutoff

+			read.cutoff = as.integer(max(counts)), # int* read_cutoff

 			verbosity = as.integer(verbosity) # int* verbosity

 		)

 	}

@@ -468,7 +469,7 @@ callPeaksUnivariate <- function(binned.data, ID, eps=0.01, init="standard", max.

 		result <- list()

 		result$ID <- ID

 	## Get states

-		message("Calculating states from posteriors ...", appendLF=F); ptm <- proc.time()

+		ptm <- startTimedMessage("Calculating states from posteriors ...")

 		hmm$posteriors <- matrix(hmm$posteriors, ncol=hmm$num.states)

 		colnames(hmm$posteriors) <- paste0("P(",state.labels,")")

 		threshold <- 1-FDR

@@ -480,7 +481,7 @@ callPeaksUnivariate <- function(binned.data, ID, eps=0.01, init="standard", max.

 	## Bin coordinates, posteriors and states

 		result$bins <- GRanges(seqnames=seqnames(binned.data),

 														ranges=ranges(binned.data),

-														reads=hmm$reads,

+														counts=hmm$counts,

 														state=states) 

 		result$bins$score <- apply(hmm$posteriors, 1, max)

 		result$bins$posteriors <- hmm$posteriors

@@ -488,18 +489,18 @@ callPeaksUnivariate <- function(binned.data, ID, eps=0.01, init="standard", max.

 			result$bins$densities <- matrix(hmm$densities, ncol=hmm$num.states)

 		}

 		seqlengths(result$bins) <- seqlengths(binned.data)

-		time <- proc.time() - ptm; message(" ",round(time[3],2),"s")

+		stopTimedMessage(ptm)

 	## Segmentation

-		message("Making segmentation ...", appendLF=F); ptm <- proc.time()

+		ptm <- startTimedMessage("Making segmentation ...")

 		gr <- result$bins

-		red.df <- suppressMessages(collapseBins(as.data.frame(gr), column2collapseBy='state', columns2average=c('reads','score','posteriors.P.modified.'), columns2drop=c('width','posteriors.P.zero.inflation.','posteriors.P.unmodified.')))

-		red.gr <- GRanges(seqnames=red.df[,1], ranges=IRanges(start=red.df[,2], end=red.df[,3]), strand=red.df[,4], mean.reads=red.df[,'mean.reads'], state=red.df[,'state'], score=red.df[,'mean.score'], mean.posterior.modified=red.df[,'mean.posteriors.P.modified.'])

+		red.df <- suppressMessages(collapseBins(as.data.frame(gr), column2collapseBy='state', columns2average=c('counts','score','posteriors.P.modified.'), columns2drop=c('width','posteriors.P.zero.inflation.','posteriors.P.unmodified.')))

+		red.gr <- GRanges(seqnames=red.df[,1], ranges=IRanges(start=red.df[,2], end=red.df[,3]), strand=red.df[,4], mean.counts=red.df[,'mean.counts'], state=red.df[,'state'], score=red.df[,'mean.score'], mean.posterior.modified=red.df[,'mean.posteriors.P.modified.'])

 		result$segments <- red.gr

 		seqlengths(result$segments) <- seqlengths(binned.data)

 		if (!keep.posteriors) {

 			result$bins$posteriors <- NULL

 		}

-		time <- proc.time() - ptm; message(" ",round(time[3],2),"s")

+		stopTimedMessage(ptm)

 	## Parameters

 		# Weights

 		result$weights <- hmm$weights

@@ -529,7 +530,7 @@ callPeaksUnivariate <- function(binned.data, ID, eps=0.01, init="standard", max.

 		# FDR

 		result$FDR <- FDR

 	## Convergence info

-		convergenceInfo <- list(eps=eps, loglik=hmm$loglik, loglik.delta=hmm$loglik.delta, num.iterations=hmm$num.iterations, time.sec=hmm$time.sec, max.mean=max.mean, read.cutoff=max(hmm$reads))

+		convergenceInfo <- list(eps=eps, loglik=hmm$loglik, loglik.delta=hmm$loglik.delta, num.iterations=hmm$num.iterations, time.sec=hmm$time.sec, max.mean=max.mean, read.cutoff=max(hmm$counts))

 		result$convergenceInfo <- convergenceInfo

 	## Add class

 		class(result) <- class.univariate.hmm

@@ -45,7 +45,7 @@ changeFDR <- function(model, FDR=0.5, separate.zeroinflation=TRUE, averages=TRUE

 	if (is(model,class.univariate.hmm)) {

 		if (is.null(model$bins$posteriors)) stop("Cannot recalculate states because posteriors are missing. Run 'callPeaksUnivariate' again with option 'keep.posteriors' set to TRUE.")

 		## Calculate states

-		message("Calculating states from posteriors ...", appendLF=F); ptm <- proc.time()

+		ptm <- startTimedMessage("Calculating states from posteriors ...")

 		states <- rep(NA,length(model$bins))

 		if (separate.zeroinflation) {

 			states[ model$bins$posteriors[,3]<threshold & model$bins$posteriors[,2]<=model$bins$posteriors[,1] ] <- 1

@@ -57,9 +57,9 @@ changeFDR <- function(model, FDR=0.5, separate.zeroinflation=TRUE, averages=TRUE

 			states <- state.labels[2:3][states]

 		}

 		model$bins$state <- states

-		time <- proc.time() - ptm; message(" ",round(time[3],2),"s")

+		stopTimedMessage(ptm)

 		## Redo segmentation

-		message("Making segmentation ...", appendLF=F); ptm <- proc.time()

+		ptm <- startTimedMessage("Making segmentation ...")

 		gr <- model$bins

 		df <- as.data.frame(gr)

 		if (averages==TRUE) {

@@ -71,7 +71,7 @@ changeFDR <- function(model, FDR=0.5, separate.zeroinflation=TRUE, averages=TRUE

 		}	

 		model$segments <- red.gr

 		seqlengths(model$segments) <- seqlengths(model$bins)

-		time <- proc.time() - ptm; message(" ",round(time[3],2),"s")

+		stopTimedMessage(ptm)

 # 		## Redo weights

 # 		model$weights <- table(model$bins$state) / length(model$bins)

@@ -79,7 +79,7 @@ changeFDR <- function(model, FDR=0.5, separate.zeroinflation=TRUE, averages=TRUE

 	} else if (is(model,class.multivariate.hmm)) {

 		if (is.null(model$bins$posteriors)) stop("Cannot recalculate states because posteriors are missing. Run 'callPeaksMultivariate' again with option 'keep.posteriors' set to TRUE.")

 		## Calculate states

-		message("Calculating states from posteriors ...", appendLF=F); ptm <- proc.time()

+		ptm <- startTimedMessage("Calculating states from posteriors ...")

 		states <- factor(bin2dec(model$bins$posteriors >= threshold), levels=levels(model$bins$state))

 		model$bins$state <- states

 		## Combinations

@@ -89,9 +89,9 @@ changeFDR <- function(model, FDR=0.5, separate.zeroinflation=TRUE, averages=TRUE

 # 			names(mapping) <- levels(model$bins$state)

 			model$bins$combination <- factor(mapping[as.character(model$bins$state)], levels=mapping[as.character(levels(model$bins$state))])

 		}

-		time <- proc.time() - ptm; message(" ",round(time[3],2),"s")

+		stopTimedMessage(ptm)

 		## Redo segmentation

-		message("Making segmentation ...", appendLF=F); ptm <- proc.time()

+		ptm <- startTimedMessage("Making segmentation ...")

 		df <- as.data.frame(model$bins)

 		ind.readcols <- grep('^reads', names(df))

 		ind.postcols <- grep('^posteriors', names(df))

@@ -116,7 +116,7 @@ changeFDR <- function(model, FDR=0.5, separate.zeroinflation=TRUE, averages=TRUE

 		}

 		model$segments <- red.gr

 		seqlengths(model$segments) <- seqlengths(model$bins)

-		time <- proc.time() - ptm; message(" ",round(time[3],2),"s")

+		stopTimedMessage(ptm)

 # 		## Redo weights

 # 		model$weights <- table(model$bins$state) / length(model$bins)

 	} else {

@@ -1,6 +1,6 @@

-#' Combinatorial and differential ChIP-seq analysis

+#' Combinatorial and differential chromatin state analysis for ChIP-seq data

 #'

 #' 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 read our publication for a detailed description of the method. TODO: insert publication.

 #'

@@ -60,7 +60,7 @@ collapseBins = function(data, column2collapseBy=NULL, columns2sumUp=NULL, column

 	ind_maxcols <- columns2getMax

 	## Make the comparison vector

-	message('Making comparison vector ...', appendLF=F); ptm <- proc.time()

+	ptm <- startTimedMessage('Making comparison vector ...')

 	if (is.null(column2collapseBy)) {

 		c <- data$start

 		cShift1 <- rep(NA,length(c))

@@ -85,13 +85,13 @@ collapseBins = function(data, column2collapseBy=NULL, columns2sumUp=NULL, column

 	compare[1] <- TRUE

 	numcollapsedbins <- length(which(compare==TRUE))

 	numbins <- nrow(data)

-	time <- proc.time() - ptm; message(" ",round(time[3],2),"s")

+	stopTimedMessage(ptm)

 	if (any(is.na(compare))) {