@@ -18,7 +18,7 @@ get.reproducibility.one2many,

 get.reproducibility.many2many,

 export.hmm2bed,

-export.binned2wiggle,

+export.hmm2wiggle,

 binned2GRanges,

 hmm2GRanges,

 hmmList2GRangesList,

@@ -6,7 +6,7 @@ hmmList2GRangesList <- function(hmm.list, reduce=TRUE, numCPU=1, consensus=FALSE

 	## Transform to GRanges

 	cat('transforming to GRanges\n')

 	if (numCPU > 1) {

-		library(doParallel)

+		suppressMessages( library(doParallel) )

 		cl <- makeCluster(numCPU)

 		registerDoParallel(cl)

 		cfun <- function(...) { GRangesList(...) }

@@ -28,7 +28,7 @@ hmmList2GRangesList <- function(hmm.list, reduce=TRUE, numCPU=1, consensus=FALSE

 		}

 		if (consensus) {

 			consensus.gr <- disjoin(unlist(hmm.grl))

-			constates <- matrix(NA, ncol=length(hmm.grl), nrow=length(hmm.grl[[1]]))

+			constates <- matrix(NA, ncol=length(hmm.grl), nrow=length(consensus.gr))

 			for (i1 in 1:length(hmm.grl)) {

 				hmm.gr <- hmm.grl[[i1]]

 				splt <- split(hmm.gr, mcols(hmm.gr)$state)

@@ -38,8 +38,8 @@ hmmList2GRangesList <- function(hmm.list, reduce=TRUE, numCPU=1, consensus=FALSE

 		}

 	}

 	if (consensus) {

-		meanstates <- apply(constates, 1, mean)

-		mcols(consensus.gr) <- meanstates

+		meanstates <- apply(constates, 1, mean, na.rm=T)

+		mcols(consensus.gr)$meanstate <- meanstates

 		return(list(grl=hmm.grl, consensus=consensus.gr, constates=constates))

 	} else {

 		return(hmm.grl)

@@ -131,8 +131,8 @@ bed2GRanges <- function(bedfile, chrom.length.file, skip=1, binsize=NULL) {

 	## Inflate every range with bins

 	if (!is.null(binsize)) {

-		inflated.data <- GRangesList()

 		grl <- split(data, seqnames(data))

+		inflated.data <- GRangesList()

 		for (i1 in 1:length(grl)) {

 			rgr <- ranges(grl[[i1]])

 			widths <- width(rgr)

@@ -156,7 +156,7 @@ bed2GRanges <- function(bedfile, chrom.length.file, skip=1, binsize=NULL) {

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

 																								strand=Rle(strand('*'), sum(numbins)),

 																								state=infstates)

-			inflated.data[[i1]] <- inflated.chrom

+			suppressWarnings( inflated.data[[i1]] <- inflated.chrom )

 		}

 		inflated.data <- unlist(inflated.data)

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

@@ -1,125 +1,120 @@

-# ===============================================

-# Write color-coded tracks with univariate states

-# ===============================================

-export.hmm2bed <- function(uni.hmm.list, file="view_me_in_genome_browser", numCPU=1) {

-

-	## Intercept user input

-	if (check.univariate.modellist(uni.hmm.list)!=0) {

-		cat("Loading univariate HMMs from files ...")

-		mlist <- NULL

-		for (modelfile in uni.hmm.list) {

-			mlist[[length(mlist)+1]] <- get(load(modelfile))

-		}

-		uni.hmm.list <- mlist

-		remove(mlist)

-		cat(" done\n")

-		if (check.univariate.modellist(uni.hmm.list)!=0) stop("argument 'uni.hmm.list' expects a list of univariate hmms or a list of files that contain univariate hmms")

+# ====================================

+# Write color-coded tracks with states

+# ====================================

+export.hmm2bed <- function(hmm.list, file="view_me_in_genome_browser", numCPU=1) {

+

+	## Function definitions

+	insertchr <- function(hmm.gr) {

+		# Change chromosome names from '1' to 'chr1' if necessary

+		mask <- which(!grepl('chr', seqnames(hmm.gr)))

+		mcols(hmm.gr)$chromosome <- as.character(seqnames(hmm.gr))

+		mcols(hmm.gr)$chromosome[mask] <- sub(pattern='^', replacement='chr', mcols(hmm.gr)$chromosome[mask])

+		mcols(hmm.gr)$chromosome <- as.factor(mcols(hmm.gr)$chromosome)

+		return(hmm.gr)

 	}

+	## Load models

+	hmm.list <- loadHmmsFromFiles(hmm.list)

+

 	## Transform to GRanges

-	cat('transforming to GRanges\n')

-	if (numCPU > 1) {

-		library(doParallel)

-		cl <- makeCluster(numCPU)

-		registerDoParallel(cl)

-		cfun <- function(...) { GRangesList(...) }

-		uni.hmm.grl <- foreach (uni.hmm = uni.hmm.list, .packages='aneufinder', .combine='cfun', .multicombine=TRUE) %dopar% {

-			hmm2GRanges(uni.hmm, reduce=T)

-		}

-		consensus.gr <- disjoin(unlist(uni.hmm.grl))

-		constates <- foreach (uni.hmm.gr = uni.hmm.grl, .packages='GenomicRanges', .combine='cbind') %dopar% {

-			splt <- split(uni.hmm.gr, mcols(uni.hmm.gr)$state)

-			mind <- as.matrix(findOverlaps(consensus.gr, splt, select='first'))

-		}

-		stopCluster(cl)

-	} else {

-		uni.hmm.grl <- GRangesList()

-		for (uni.hmm in uni.hmm.list) {

-			uni.hmm.grl[[length(uni.hmm.grl)+1]] <- hmm2GRanges(uni.hmm, reduce=T)

-		}

-		consensus.gr <- disjoin(unlist(uni.hmm.grl))

-		constates <- matrix(NA, ncol=length(uni.hmm.grl), nrow=length(uni.hmm.grl[[1]]))

-		for (i1 in 1:length(uni.hmm.grl)) {

-			uni.hmm.gr <- uni.hmm.grl[[i1]]

-			splt <- split(uni.hmm.gr, mcols(uni.hmm.gr)$state)

-			mind <- as.matrix(findOverlaps(consensus.gr, splt, select='first'))

-			constates[,i1] <- mind

-		}

-	}

-	meanstates <- apply(constates, 1, mean)

-	mcols(consensus.gr) <- meanstates

+	temp <- hmmList2GRangesList(hmm.list, reduce=T, numCPU=numCPU, consensus=T)

+	consensus.gr <- insertchr(temp$consensus)

+	hmm.grl <- lapply(temp$grl, insertchr)

 	# Variables

-	nummod <- length(uni.hmm.list)

+	nummod <- length(hmm.list)

+	numstates <- hmm.list[[1]]$num.states

 	file <- paste(file,"bed", sep=".")

 	# Generate the colors

-	colors <- state.colors[state.labels]

+	colors <- state.colors[levels(hmm.grl[[1]]$state)]

 	RGBs <- t(col2rgb(colors))

 	RGBs <- apply(RGBs,1,paste,collapse=",")

+	cf <- colorRamp(colors, space='rgb', interpolate='linear')

+	consensusRGBs <- round(cf((mcols(consensus.gr)$meanstate-1)/(numstates-1)), 0)

+	itemconsensusRgb <- apply(consensusRGBs,1,paste,collapse=",")

 	# Write first line to file

-	cat("browser hide all\n", file=file)

+	cat('writing to file',file,'\n')

+# 	cat("browser hide all\n", file=file)

+	cat("", file=file)

 	### Write every model to file ###

 	for (imod in 1:nummod) {

-		uni.hmm <- uni.hmm.list[[imod]]

-		uni.hmm.gr <- uni.hmm.grl[[imod]]

-		priority <- 50 + imod

-		cat(paste("track name=",uni.hmm$ID," description=\"univariate calls for ",names(uni.hmm.list)[imod],"\" visibility=1 itemRgb=On priority=",priority,"\n", sep=""), file=file, append=TRUE)

-

-		# Change chromosome names from '1' to 'chr1' if necessary

-		mask <- which(!grepl('chr', seqnames(uni.hmm.gr)))

-		mcols(uni.hmm.gr)$chromosome <- as.character(seqnames(uni.hmm.gr))

-		mcols(uni.hmm.gr)$chromosome[mask] <- sub(pattern='^', replacement='chr', mcols(uni.hmm.gr)$chromosome[mask])

-		mcols(uni.hmm.gr)$chromosome <- as.factor(mcols(uni.hmm.gr)$chromosome)

-

-		# Collapse the calls

-		collapsed.calls <- as.data.frame(uni.hmm.gr)[c('chromosome','start','end','state')]

-

-		# Append to file

+		cat('writing hmm',imod,'/',nummod,'\r')

+		hmm <- hmm.list[[imod]]

+		hmm.gr <- hmm.grl[[imod]]

+		priority <- 51 + imod

+		cat(paste0("track name=",hmm$ID," state"," description=\"",hmm$ID," state","\" visibility=1 itemRgb=On priority=",priority,"\n"), file=file, append=TRUE)

+		collapsed.calls <- as.data.frame(hmm.gr)[c('chromosome','start','end','state')]

 		itemRgb <- RGBs[as.character(collapsed.calls$state)]

 		numsegments <- nrow(collapsed.calls)

 		df <- cbind(collapsed.calls, score=rep(0,numsegments), strand=rep(".",numsegments), thickStart=collapsed.calls$start, thickEnd=collapsed.calls$end, itemRgb=itemRgb)

 		write.table(format(df, scientific=FALSE), file=file, append=TRUE, row.names=FALSE, col.names=FALSE, quote=FALSE)

 	}

+	cat('\n')

+

+	### Write consensus model to file ###

+	cat('writing mean states to file\n')

+	priority <- 50

+	cat(paste0("track name=\"","average state","\" description=\"","average state","\" visibility=1 itemRgb=On priority=",priority,"\n"), file=file, append=TRUE)

+	collapsed.calls <- as.data.frame(consensus.gr)[c('chromosome','start','end','meanstate')]

+	numsegments <- nrow(collapsed.calls)

+	df <- cbind(collapsed.calls, score=rep(0,numsegments), strand=rep(".",numsegments), thickStart=collapsed.calls$start, thickEnd=collapsed.calls$end, itemRgb=itemconsensusRgb)

+	write.table(format(df, scientific=FALSE), file=file, append=TRUE, row.names=FALSE, col.names=FALSE, quote=FALSE)

 }

-# ====================================

-# Write signal tracks from binned data

-# ====================================

-export.binned2wiggle <- function(binned.data.list, file="view_me_in_genome_browser") {

-	

-	# Check user input

-	if (is.null(names(binned.data.list))) {

-		stop("Please name the list entries. The names will be used as track name for the Genome Browser file.")

+# =============================

+# Write signal tracks from hmms

+# =============================

+export.hmm2wiggle <- function(hmm.list, file="view_me_in_genome_browser", numCPU=1) {

+

+	## Function definitions

+	insertchr <- function(hmm.gr) {

+		# Change chromosome names from '1' to 'chr1' if necessary

+		mask <- which(!grepl('chr', seqnames(hmm.gr)))

+		mcols(hmm.gr)$chromosome <- as.character(seqnames(hmm.gr))

+		mcols(hmm.gr)$chromosome[mask] <- sub(pattern='^', replacement='chr', mcols(hmm.gr)$chromosome[mask])

+		mcols(hmm.gr)$chromosome <- as.factor(mcols(hmm.gr)$chromosome)

+		return(hmm.gr)

 	}

+	## Load models

+	hmm.list <- loadHmmsFromFiles(hmm.list)

+

+	## Transform to GRanges

+	temp <- hmmList2GRangesList(hmm.list, reduce=F, numCPU=numCPU, consensus=T)

+	hmm.grl <- temp$grl

+	consensus.gr <- insertchr(temp$consensus)

+	hmm.grl <- lapply(hmm.grl, insertchr)

+

 	# Variables

-	file <- paste(file,"wig", sep=".")

+	nummod <- length(hmm.list)

+	numstates <- hmm.list[[1]]$num.states

+	file <- paste(file,"wiggle", sep=".")

-	# Write to file

-	cat("browser hide all\n", file=file)

+	# Write first line to file

+	cat('writing to file',file,'\n')

+# 	cat("browser hide all\n", file=file)

+	cat("", file=file)

-	# Go through binned.data.list

-	for (i1 in 1:length(binned.data.list)) {

-		binned.data <- binned.data.list[[i1]]

-		names(binned.data) <- binned.data.names

-		

-		# Write track information

-		name <- names(binned.data.list)[i1]

-		binsize <- diff(binned.data$start[1:2])

-		description <- paste("binned read counts ",binsize,"bp", sep="")

-		cat(paste("track type=wiggle_0 name=\"",name,"\" description=\"",description,"\" visibility=full autoScale=on color=90,90,90 maxHeightPixels=100:50:20 graphType=bar priority=",i1,"\n", sep=""), file=file, append=TRUE)

+	### Write every model to file ###

+	for (imod in 1:nummod) {

+		cat('writing hmm',imod,'/',nummod,'\r')

+		hmm <- hmm.list[[imod]]

+		hmm.gr <- hmm.grl[[imod]]

+		priority <- 50 + imod

+		binsize <- width(hmm.gr[1])

+		cat(paste("track type=wiggle_0 name=\"",hmm$ID," read count","\" description=\"",hmm$ID," read count","\" visibility=full autoScale=on color=90,90,90 maxHeightPixels=100:50:20 graphType=bar priority=",priority,"\n", sep=""), file=file, append=TRUE)

 		# Write read data

-		for (chrom in unique(binned.data$chrom)) {

+		for (chrom in unique(hmm.gr$chromosome)) {

 			cat(paste("fixedStep chrom=",chrom," start=1 step=",binsize," span=",binsize,"\n", sep=""), file=file, append=TRUE)

-			write.table(binned.data$reads, file=file, append=TRUE, row.names=FALSE, col.names=FALSE)

+			write.table(mcols(hmm.gr[hmm.gr$chromosome==chrom])$reads, file=file, append=TRUE, row.names=FALSE, col.names=FALSE)

 		}

 	}

-

+	cat('\n')

 }

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

-findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max.time=-1, max.iter=-1, num.trials=1, eps.try=NULL, num.threads=1, output.if.not.converged=FALSE, filter.reads=TRUE) {

+findCNVs <- function(binned.data, ID, eps=0.001, init="random", max.time=-1, max.iter=-1, num.trials=1, eps.try=NULL, num.threads=1, output.if.not.converged=FALSE, read.cutoff.quantile=0.999) {

 	## Intercept user input

 	IDcheck <- ID  #trigger error if not defined

+	if (class(binned.data) != 'GRanges') {

+		binned.data <- get(load(binned.data))

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

+	}

 	if (check.positive(eps)!=0) stop("argument 'eps' expects a positive numeric")

 	if (check.integer(max.time)!=0) stop("argument 'max.time' expects an integer")

 	if (check.integer(max.iter)!=0) stop("argument 'max.iter' expects an integer")

@@ -28,28 +32,58 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max.time=-1, m

 	}

 	# Filter high reads out, makes HMM faster

-	read.cutoff <- as.integer(quantile(reads, 0.9999))

-	if (filter.reads) {

-		mask <- reads > read.cutoff

-		reads[mask] <- read.cutoff

-		numfiltered <- length(which(mask))

-		cat(paste0("Replaced read counts > ",read.cutoff," (99.99% quantile) by ",read.cutoff," in ",numfiltered," bins. Set option 'filter.reads=FALSE' to disable this filtering.\n"))

-# 		if (numfiltered > 0) {

-# 			warning(paste("There are very high read counts in your data (probably artificial). Replaced read counts > ",read.cutoff," (99.99% quantile) by ",read.cutoff," in ",numfiltered," bins. Set option 'filter.reads=FALSE' to disable this filtering.", sep=""))

-# 		}

+	read.cutoff <- as.integer(quantile(reads, read.cutoff.quantile))

+	mask <- reads > read.cutoff

+	reads[mask] <- read.cutoff

+	numfiltered <- length(which(mask))

+	if (numfiltered > 0) {

+		cat(paste0("Replaced read counts > ",read.cutoff," (",names(read.cutoff)," quantile) by ",read.cutoff," in ",numfiltered," bins. Set option 'read.cutoff.quantile=1' to disable this filtering.\n"))

 	}

-	

 	## Call univariate in a for loop to enable multiple trials

 	modellist <- list()

 	for (i_try in 1:num.trials) {

 		cat("\n\nTry ",i_try," of ",num.trials," ------------------------------\n")

+

+		## Initial parameters

+		if (init == 'random') {

+			A.initial <- matrix(runif(numstates^2), ncol=numstates)

+			A.initial <- sweep(A.initial, 1, rowSums(A.initial), "/")			

+			proba.initial <- runif(numstates)

+			size.initial <- runif(1, min=0, max=1000) * cumsum(ifelse(state.distributions=='dnbinom' | state.distributions=='dbinom', T, F))

+			prob.initial <- runif(1) * ifelse(state.distributions=='dnbinom', T, F)

+			prob.initial[state.distributions=='dgeom'] <- runif(1)

+			lambda.initial <- runif(1, min=0, max=1000) * cumsum(ifelse(state.distributions=='dpois', T, F))

+		} else if (init == 'standard') {

+			A.initial <- matrix(NA, ncol=numstates, nrow=numstates)

+			for (irow in 1:numstates) {

+				for (icol in 1:numstates) {

+					if (irow==icol) { A.initial[irow,icol] <- 0.9 }

+					else { A.initial[irow,icol] <- 0.1/(numstates-1) }

+				}

+			}

+			proba.initial <- rep(1/numstates, numstates)

+			mean.initial <- mean(reads[reads>0])/2 * cumsum(ifelse(state.distributions=='dnbinom' | state.distributions=='dbinom', T, F))

+			var.initial <- var(reads[reads>0])/2 * cumsum(ifelse(state.distributions=='dnbinom' | state.distributions=='dbinom', T, F))

+			size.initial <- rep(0,numstates)

+			prob.initial <- rep(0,numstates)

+			mask <- state.distributions=='dnbinom'

+			size.initial[mask] <- dnbinom.size(mean.initial[mask], var.initial[mask])

+			prob.initial[mask] <- dnbinom.prob(mean.initial[mask], var.initial[mask])

+			mask <- state.distributions=='dbinom'

+			size.initial[mask] <- dbinom.size(mean.initial[mask], var.initial[mask])

+			prob.initial[mask] <- dbinom.prob(mean.initial[mask], var.initial[mask])

+			prob.initial[state.distributions=='dgeom'] <- 0.9

+			lambda.initial <- mean(reads[reads>0])/2 * cumsum(ifelse(state.distributions=='dpois', T, F))

+		}

+	

 		hmm <- .C("R_univariate_hmm",

 			reads = as.integer(reads), # int* O

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

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

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

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

+			lambda = double(length=numstates), # double* lambda

 			num.iterations = as.integer(max.iter), #  int* maxiter

 			time.sec = as.integer(max.time), # double* maxtime

 			loglik.delta = as.double(eps.try), # double* eps

@@ -60,11 +94,12 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max.time=-1, m

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

 			distr.type = as.integer(state.distributions), # int* distr_type

 			ini.proc = as.integer(iniproc), # int* iniproc

-			size.initial = double(length=numstates), # double* initial_size

-			prob.initial = double(length=numstates), # double* initial_prob

-			A.initial = double(length=numstates*numstates), # double* initial_A

-			proba.initial = double(length=numstates), # double* initial_proba

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

+			size.initial = as.vector(size.initial), # double* initial_size

+			prob.initial = as.vector(prob.initial), # double* initial_prob

+			lambda.initial = as.vector(lambda.initial), # double* initial_lambda

+			A.initial = as.vector(A.initial), # double* initial_A

+			proba.initial = as.vector(proba.initial), # double* initial_proba

+			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(read.cutoff) # int* read_cutoff

@@ -95,6 +130,7 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max.time=-1, m

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

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

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

+			lambda = double(length=numstates), # double* lambda

 			num.iterations = as.integer(max.iter), #  int* maxiter

 			time.sec = as.integer(max.time), # double* maxtime

 			loglik.delta = as.double(eps), # double* eps

@@ -107,6 +143,7 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max.time=-1, m

 			ini.proc = as.integer(iniproc), # int* iniproc

 			size.initial = as.vector(hmm$size), # double* initial_size

 			prob.initial = as.vector(hmm$prob), # double* initial_prob

+			lambda.initial = as.vector(hmm$lambda), # double* initial_lambda

 			A.initial = as.vector(hmm$A), # double* initial_A

 			proba.initial = as.vector(hmm$proba), # double* initial_proba

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

@@ -117,6 +154,7 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max.time=-1, m

 	}

 	# Add useful entries

+	hmm$read.cutoff <- read.cutoff

 	hmm$ID <- ID

 	names(hmm$weights) <- state.labels

 	hmm$coordinates <- data.frame(as.character(seqnames(binned.data)), start(ranges(binned.data)), end(ranges(binned.data)))

@@ -125,32 +163,44 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max.time=-1, m

 	class(hmm) <- class.aneufinder.hmm

 	hmm$states <- factor(state.labels, levels=state.labels)[hmm$states+1]

 	hmm$eps <- eps

-	hmm$A <- matrix(hmm$A, ncol=hmm$num.states, byrow=TRUE)

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

 	rownames(hmm$A) <- state.labels

 	colnames(hmm$A) <- state.labels

-	hmm$distributions <- data.frame(type=state.distributions, size=hmm$size, prob=hmm$prob, mu=fmean(hmm$size,hmm$prob), variance=fvariance(hmm$size,hmm$prob))

-	rownames(hmm$distributions) <- state.labels

-	# Treat 'null-mixed' separately

-	if ('null-mixed' %in% state.labels) {

-		hmm$distributions['null-mixed','mu'] <- (1-hmm$distributions['null-mixed','prob'])/hmm$distributions['null-mixed','prob']

-		hmm$distributions['null-mixed','variance'] <- hmm$distributions['null-mixed','mu']/hmm$distributions['null-mixed','prob']

-		hmm$distributions['null-mixed','size'] <- NA

+	hmm$distributions <- data.frame()

+	hmm$distributions.initial <- data.frame()

+	for (idistr in 1:length(hmm$distr.type)) {

+		distr <- levels(state.distributions)[hmm$distr.type[idistr]]

+		if (distr == 'dnbinom') {

+			hmm$distributions <- rbind(hmm$distributions, data.frame(type=distr, size=hmm$size[idistr], prob=hmm$prob[idistr], lambda=NA, mu=dnbinom.mean(hmm$size[idistr],hmm$prob[idistr]), variance=dnbinom.variance(hmm$size[idistr],hmm$prob[idistr])))

+			hmm$distributions.initial <- rbind(hmm$distributions.initial, data.frame(type=distr, size=hmm$size.initial[idistr], prob=hmm$prob.initial[idistr], lambda=NA, mu=dnbinom.mean(hmm$size.initial[idistr],hmm$prob.initial[idistr]), variance=dnbinom.variance(hmm$size.initial[idistr],hmm$prob.initial[idistr])))

+		} else if (distr == 'dgeom') {

+			hmm$distributions <- rbind(hmm$distributions, data.frame(type=distr, size=NA, prob=hmm$prob[idistr], lambda=NA, mu=dgeom.mean(hmm$prob[idistr]), variance=dgeom.variance(hmm$prob[idistr])))

+			hmm$distributions.initial <- rbind(hmm$distributions.initial, data.frame(type=distr, size=NA, prob=hmm$prob.initial[idistr], lambda=NA, mu=dgeom.mean(hmm$prob.initial[idistr]), variance=dgeom.variance(hmm$prob.initial[idistr])))

+		} else if (distr == 'delta') {

+			hmm$distributions <- rbind(hmm$distributions, data.frame(type=distr, size=NA, prob=NA, lambda=NA, mu=0, variance=0))

+			hmm$distributions.initial <- rbind(hmm$distributions.initial, data.frame(type=distr, size=NA, prob=NA, lambda=NA, mu=0, variance=0))

+		} else if (distr == 'dpois') {

+			hmm$distributions <- rbind(hmm$distributions, data.frame(type=distr, size=NA, prob=NA, lambda=hmm$lambda[idistr], mu=hmm$lambda[idistr], variance=hmm$lambda[idistr]))

+			hmm$distributions.initial <- rbind(hmm$distributions.initial, data.frame(type=distr, size=NA, prob=NA, lambda=hmm$lambda.initial[idistr], mu=hmm$lambda.initial[idistr], variance=hmm$lambda.initial[idistr]))

+		} else if (distr == 'dbinom') {

+			hmm$distributions <- rbind(hmm$distributions, data.frame(type=distr, size=hmm$size[idistr], prob=hmm$prob[idistr], lambda=NA, mu=dbinom.mean(hmm$size[idistr],hmm$prob[idistr]), variance=dbinom.variance(hmm$size[idistr],hmm$prob[idistr])))

+			hmm$distributions.initial <- rbind(hmm$distributions.initial, data.frame(type=distr, size=hmm$size.initial[idistr], prob=hmm$prob.initial[idistr], lambda=NA, mu=dbinom.mean(hmm$size.initial[idistr],hmm$prob.initial[idistr]), variance=dbinom.variance(hmm$size.initial[idistr],hmm$prob.initial[idistr])))

+		}

 	}

-	hmm$A.initial <- matrix(hmm$A.initial, ncol=hmm$num.states, byrow=TRUE)

+	rownames(hmm$distributions) <- state.labels

+	rownames(hmm$distributions.initial) <- state.labels

+	hmm$A.initial <- matrix(hmm$A.initial, ncol=hmm$num.states)

 	rownames(hmm$A.initial) <- state.labels

 	colnames(hmm$A.initial) <- state.labels

-	hmm$distributions.initial <- data.frame(type=state.distributions, size=hmm$size.initial, prob=hmm$prob.initial, mu=fmean(hmm$size.initial,hmm$prob.initial), variance=fvariance(hmm$size.initial,hmm$prob.initial))

-	rownames(hmm$distributions.initial) <- state.labels

-	hmm$distributions.initial['nullsomy',2:5] <- c(0,1,0,0)

-	hmm$filter.reads <- filter.reads

 	# Delete redundant entries

 	hmm$size <- NULL

 	hmm$prob <- NULL

+	hmm$lambda <- NULL

 	hmm$size.initial <- NULL

 	hmm$prob.initial <- NULL

+	hmm$lambda.initial <- NULL

 	hmm$use.initial.params <- NULL

-	hmm$read.cutoff <- NULL

 	hmm$distr.type <- NULL

 	# Issue warnings

@@ -160,7 +210,7 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max.time=-1, m

 		}

 	}

 	if (hmm$error == 1) {

-		stop("A nan occurred during the Baum-Welch! Parameter estimation terminated prematurely. Check your read counts for very high numbers, they could be the cause for this problem.")

+		stop("A nan occurred during the Baum-Welch! Parameter estimation terminated prematurely. Check your read counts for very high numbers, they could be the cause for this problem. Try again with a lower value for 'read.cutoff.quantile'.")

 	} else if (hmm$error == 2) {

 		stop("An error occurred during the Baum-Welch! Parameter estimation terminated prematurely.")

 	}

@@ -1,8 +1,8 @@

 # =======================================================

 # Some global variables that can be used in all functions

 # =======================================================

-state.labels <- c("nullsomy","null-mixed","monosomy","disomy","trisomy","tetrasomy","multisomy")

-state.distributions <- factor(c('delta','dgeom','dnbinom','dnbinom','dnbinom','dnbinom','dnbinom'), levels=c('delta','dgeom','dnbinom'))

+state.labels <- c("nullsomy","monosomy","disomy","trisomy","tetrasomy","multisomy")

+state.distributions <- factor(c('delta','dnbinom','dnbinom','dnbinom','dnbinom','dnbinom'), levels=c('delta','dgeom','dnbinom','dpois','dbinom'))

 coordinate.names <- c("chrom","start","end")

 binned.data.names <- c(coordinate.names,"reads")

 class.aneufinder.hmm <- "aneufinder.hmm"

@@ -13,18 +13,44 @@ get.state.colors <- function() { return(state.colors[state.labels]) }

 # ============================================================================

 # Functions for a Negative Binomial to transform (mean,variance)<->(size,prob)

 # ============================================================================

-fsize <- function(mean, variance) {

+dnbinom.size <- function(mean, variance) {

 	return(mean^2 / (variance - mean))

 }

-fprob <- function(mean, variance) {

+dnbinom.prob <- function(mean, variance) {

 	return(mean/variance)

 }

-fmean <- function(size, prob) {

+dnbinom.mean <- function(size, prob) {

 	return(size/prob - size)

 }

-fvariance <- function(size, prob) {

+dnbinom.variance <- function(size, prob) {

 	return( (size - prob*size) / prob^2 )

 }

+

+dgeom.mean <- function(prob) {

+	return( (1-prob)/prob )

+}

+

+dgeom.variance <- function(prob) {

+	return( (1-prob)/prob^2 )

+}

+

+dbinom.size <- function(mean, variance) {

+	return( mean^2/(mean-variance) )

+}

+

+dbinom.prob <- function(mean, variance) {

+	return( (mean-variance)/mean )

+}

+

+dbinom.mean <- function(size, prob) {

+	return( size*prob )

+}

+

+dbinom.variance <- function(size, prob) {

+	return( size*prob * (1-prob) )

+}

+

+

@@ -80,6 +80,15 @@ plot.distribution <- function(model, state=NULL, chrom=NULL, start=NULL, end=NUL

 		} else if (model$distributions[istate,'type']=='dnbinom') {

 			# negative binomials

 			distributions[[length(distributions)+1]] <- weights[istate] * dnbinom(x, model$distributions[istate,'size'], model$distributions[istate,'prob'])

+		} else if (model$distributions[istate,'type']=='dpois') {

+			# poissons

+			distributions[[length(distributions)+1]] <- weights[istate] * dpois(x, model$distributions[istate,'lambda'])

+		} else if (model$distributions[istate,'type']=='dbinom') {

+			# binomials

+			s <- model$distributions[istate,'size']

+			p <- model$distributions[istate,'prob']

+# 			distributions[[length(distributions)+1]] <- weights[istate] * dbinom(x, model$distributions[istate,'size'], model$distributions[istate,'prob'])	# only defined for integer 'size'

+			distributions[[length(distributions)+1]] <- weights[istate] * choose(s,x) * p^x * (1-p)^(s-x)

 		}

 	}

 	distributions <- as.data.frame(distributions)

@@ -23,6 +23,7 @@ The HMM object is output of the function \code{\link{findCNVs}} and is basically

 \item{proba.initial}{Initial \option{proba} at the beginning of the Baum-Welch.}

 \item{num.threads}{Number of threads that were used in the Baum-Welch.}

 \item{error}{This is 0 unless an error occurred during the Baum-Welch.}

+\item{read.cutoff}{The read cutoff as specified by the option \option{read.cutoff.quantile}.}

 \item{eps}{Convergence threshold for the Baum-Welch.}

 \item{distributions}{Estimated parameters of the employed distributions. See \code{\link{findCNVs}} for a description of which distribution is used for which state.}

 \item{distributions.initial}{Distribution parameters at the beginning of the Baum-Welch.}

@@ -24,13 +24,13 @@ Initialization procedure for the distributions. Currently only 'standard' is imp

 \item{output.if.not.converged}{

 If set to FALSE, no output will be given in the case that the Baum-Welch did not converge to the desired \option{eps}. If set to TRUE, a model will be returned even if the Baum-Welch did not converge to the desired \option{eps}.

 }

-\item{filter.reads}{

-If set to TRUE, then read counts which are above the 99.99\% quantile of the data are set to this value. A warning is issued in this case. This filtering can increase the speed of the Baum-Welch and is unlikely to affect the results.

+\item{read.cutoff.quantile}{

+Read counts which are above this quantile are replaced by the corresponding value. This filtering can increase the speed of the Baum-Welch and is unlikely to affect the results. It may also be helpful to specify a lower cutoff if you are getting the error message 'A nan occurred during the Baum-Welch! ...'.

 }

 }

 \details{

-The Hidden Markov Model which is used to classify the bins uses 7 states: state 'nullsomy' with a delta function as emission densitiy (only zero read counts), 'null-mixed' with a geometric distribution to account, 'monosomy','disomy','trisomy','tetrasomy' and 'multisomy' with Negative Binomials as emission densities. A Baum-Welch algorithm is employed to estimate the parameters of the distributions. See our paper for a detailed description of the method.

+The Hidden Markov Model which is used to classify the bins uses 7 states: state 'nullsomy' with a delta function as emission densitiy (only zero read counts), 'null-mixed' with a geometric distribution (see \code{\link{dgeom}}), 'monosomy','disomy','trisomy','tetrasomy' and 'multisomy' with negative binomials (see \code{\link{dnbinom}}) as emission densities. A Baum-Welch algorithm is employed to estimate the parameters of the distributions. See our paper for a detailed description of the method.

 TODO: insert paper

 }

 \value{

@@ -1,13 +1,13 @@

 #include "utility.h"

 #include "scalehmm.h"

 #include "loghmm.h"

-#include <omp.h> // parallelization options

+// #include <omp.h> // parallelization options

 // ===================================================================================================================================================

 // This function takes parameters from R, creates a univariate HMM object, creates the distributions, runs the Baum-Welch and returns the result to R.

 // ===================================================================================================================================================

 extern "C" {

-void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* maxiter, int* maxtime, double* eps, int* states, double* A, double* proba, double* loglik, double* weights, int* distr_type, int* iniproc, double* initial_size, double* initial_prob, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* read_cutoff)

+void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, double* lambda, int* maxiter, int* maxtime, double* eps, int* states, double* A, double* proba, double* loglik, double* weights, int* distr_type, int* iniproc, double* initial_size, double* initial_prob, double* initial_lambda, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* read_cutoff)

 {

 	// Define logging level

@@ -16,8 +16,8 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

  	FILELog::ReportingLevel() = FILELog::FromString("ERROR");

 //  	FILELog::ReportingLevel() = FILELog::FromString("DEBUG2");

-	// Parallelization settings

-	omp_set_num_threads(*num_threads);

+// 	// Parallelization settings

+// 	omp_set_num_threads(*num_threads);

 	// Print some information

 	FILE_LOG(logINFO) << "number of states = " << *N;

@@ -75,43 +75,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 	FILE_LOG(logINFO) << "data mean = " << mean << ", data variance = " << variance;		

 	Rprintf("data mean = %g, data variance = %g\n", mean, variance);		

-	// Go through all states of the hmm and assign the density functions

-	// This loop assumes that the negative binomial states come last and are consecutive

-	double imean, ivariance;

-	for (int i_state=0; i_state<*N; i_state++)

-	{

-		if (*use_initial_params) {

-			FILE_LOG(logINFO) << "Using given parameters for size and prob";

-			Rprintf("Using given parameters for size and prob\n");

-			imean = (1-initial_prob[i_state])*initial_size[i_state] / initial_prob[i_state];

-			ivariance = imean / initial_prob[i_state];

-			FILE_LOG(logDEBUG2) << "imean = " << imean;

-			FILE_LOG(logDEBUG2) << "ivariance = " << ivariance;

-		} else {

-

-			if (*iniproc == 1)

-			{

-				// Simple initialization based on data mean, assumed to be the disomic mean

-				if (distr_type[i_state] == 1) { }

-				else if (distr_type[i_state] == 2) { }

-				else if (distr_type[i_state] == 3)

-				{

-					for (int ii_state=i_state; ii_state<*N; ii_state++)

-					{

-						imean = mean/2 * (ii_state-i_state+1);

-						ivariance = imean * 5;

-// 						ivariance = variance/2 * (i_state-1);

-						// Calculate r and p from mean and variance

-						initial_size[ii_state] = pow(imean,2)/(ivariance-imean);

-						initial_prob[ii_state] = imean/ivariance;

-					}

-					break;

-				}

-			}

-

-		}

-	}

-

+	// Create the emission densities and initialize

 	for (int i_state=0; i_state<*N; i_state++)

 	{

 		if (distr_type[i_state] == 1)

@@ -123,7 +87,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 		else if (distr_type[i_state] == 2)

 		{

 			FILE_LOG(logDEBUG1) << "Using geometric distribution for state " << i_state;

-			Geometric *d = new Geometric(O, *T, 0.9); // delete is done inside ~ScaleHMM()

+			Geometric *d = new Geometric(O, *T, initial_prob[i_state]); // delete is done inside ~ScaleHMM()

 			hmm->densityFunctions.push_back(d);

 		}

 		else if (distr_type[i_state] == 3)

@@ -132,6 +96,18 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 			NegativeBinomial *d = new NegativeBinomial(O, *T, initial_size[i_state], initial_prob[i_state]); // delete is done inside ~ScaleHMM()

 			hmm->densityFunctions.push_back(d);

 		}

+		else if (distr_type[i_state] == 4)

+		{

+			FILE_LOG(logDEBUG1) << "Using poisson for state " << i_state;

+			Poisson *d = new Poisson(O, *T, initial_lambda[i_state]); // delete is done inside ~ScaleHMM()

+			hmm->densityFunctions.push_back(d);

+		}

+		else if (distr_type[i_state] == 5)

+		{

+			FILE_LOG(logDEBUG1) << "Using binomial for state " << i_state;

+			NegativeBinomial *d = new NegativeBinomial(O, *T, initial_size[i_state], initial_prob[i_state]); // delete is done inside ~ScaleHMM()

+			hmm->densityFunctions.push_back(d);

+		}

 		else

 		{

 			FILE_LOG(logWARNING) << "Density not specified, using default negative binomial for state " << i_state;

@@ -189,7 +165,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 		proba[i] = hmm->get_proba(i);

 		for (int j=0; j<*N; j++)

 		{

-			A[i * (*N) + j] = hmm->get_A(i,j);

+			A[i * (*N) + j] = hmm->get_A(j,i);

 		}

 	}

@@ -214,6 +190,17 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 			size[i] = 0;

 			prob[i] = 1;

 		}

+		else if (hmm->densityFunctions[i]->get_name() == POISSON)

+		{

+			Poisson* d = (Poisson*)(hmm->densityFunctions[i]);

+			lambda[i] = d->get_lambda();

+		}

+		else if (hmm->densityFunctions[i]->get_name() == BINOMIAL) 

+		{

+			Binomial* d = (Binomial*)(hmm->densityFunctions[i]);

+			size[i] = d->get_size();

+			prob[i] = d->get_prob();

+		}

 	}

 	*loglik = hmm->get_logP();

 	hmm->calc_weights(weights);

@@ -91,6 +91,214 @@ void Normal::set_stdev(double stdev)

 }

+// ============================================================

+// Poisson density

+// ============================================================

+

+// Constructor and Destructor ---------------------------------

+Poisson::Poisson(int* observations, int T, double lambda)

+{

+	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	this->obs = observations;

+	this->T = T;

+	this->lambda = lambda;

+	this->lxfactorials = NULL;

+	// Precompute the lxfactorials that are used in computing the densities

+	if (this->obs != NULL)

+	{

+		this->max_obs = intMax(observations, T);

+		this->lxfactorials = (double*) calloc(max_obs+1, sizeof(double));

+		this->lxfactorials[0] = 0.0;	// Not necessary, already 0 because of calloc

+		this->lxfactorials[1] = 0.0;

+		for (int j=2; j<=max_obs; j++)

+		{

+			this->lxfactorials[j] = this->lxfactorials[j-1] + log(j);

+		}

+	}

+}

+

+Poisson::~Poisson()

+{

+	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	if (this->lxfactorials != NULL)

+	{

+		free(this->lxfactorials);

+	}

+}

+

+// Methods ----------------------------------------------------

+void Poisson::calc_logdensities(double* logdens)

+{

+	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	double logl = log(this->lambda);

+	double l = this->lambda;

+	double lxfactorial;

+	// Select strategy for computing densities

+	if (this->max_obs <= this->T)

+	{

+		FILE_LOG(logDEBUG2) << "Precomputing densities in " << __func__ << " for every obs[t], because max(O)<=T";

+		double logdens_per_read [this->max_obs+1];

+		for (int j=0; j<=this->max_obs; j++)

+		{

+			logdens_per_read[j] = j*logl - l - this->lxfactorials[j];

+		}

+		for (int t=0; t<this->T; t++)

+		{

+			logdens[t] = logdens_per_read[(int) this->obs[t]];

+			FILE_LOG(logDEBUG4) << "logdens["<<t<<"] = " << logdens[t];

+			if (isnan(logdens[t]))

+			{

+				FILE_LOG(logERROR) << __PRETTY_FUNCTION__;

+				FILE_LOG(logERROR) << "logdens["<<t<<"] = "<< logdens[t];

+				throw nan_detected;

+			}

+		}

+	}

+	else

+	{

+		FILE_LOG(logDEBUG2) << "Computing densities in " << __func__ << " for every t, because max(O)>T";

+		for (int t=0; t<this->T; t++)

+		{

+			lxfactorial = this->lxfactorials[(int) this->obs[t]];

+			logdens[t] = this->obs[t]*logl - l - lxfactorial;

+			FILE_LOG(logDEBUG4) << "logdens["<<t<<"] = " << logdens[t];

+			if (isnan(logdens[t]))

+			{

+				FILE_LOG(logERROR) << __PRETTY_FUNCTION__;

+				FILE_LOG(logERROR) << "logdens["<<t<<"] = "<< logdens[t];

+				throw nan_detected;

+			}

+		}

+	}

+} 

+

+void Poisson::calc_densities(double* dens)

+{

+	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	double logl = log(this->lambda);

+	double l = this->lambda;

+	double lxfactorial;

+	// Select strategy for computing densities

+	if (this->max_obs <= this->T)

+	{

+		FILE_LOG(logDEBUG2) << "Precomputing densities in " << __func__ << " for every obs[t], because max(O)<=T";

+		double dens_per_read [this->max_obs+1];

+		for (int j=0; j<=this->max_obs; j++)

+		{

+			dens_per_read[j] = exp( j*logl - l - this->lxfactorials[j] );

+		}

+		for (int t=0; t<this->T; t++)

+		{

+			dens[t] = dens_per_read[(int) this->obs[t]];

+			FILE_LOG(logDEBUG4) << "dens["<<t<<"] = " << dens[t];

+			if (isnan(dens[t]))

+			{

+				FILE_LOG(logERROR) << __PRETTY_FUNCTION__;

+				FILE_LOG(logERROR) << "dens["<<t<<"] = "<< dens[t];

+				throw nan_detected;

+			}

+		}

+	}

+	else

+	{

+		FILE_LOG(logDEBUG2) << "Computing densities in " << __func__ << " for every t, because max(O)>T";

+		for (int t=0; t<this->T; t++)

+		{

+			lxfactorial = this->lxfactorials[(int) this->obs[t]];

+			dens[t] = exp( this->obs[t]*logl - l - lxfactorial );

+			FILE_LOG(logDEBUG4) << "dens["<<t<<"] = " << dens[t];

+			if (isnan(dens[t]))

+			{

+				FILE_LOG(logERROR) << __PRETTY_FUNCTION__;

+				FILE_LOG(logERROR) << "dens["<<t<<"] = "<< dens[t];

+				throw nan_detected;

+			}

+		}

+	}

+} 

+

+void Poisson::update(double* weights)

+{

+	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	double numerator, denominator;

+	// Update lambda

+	numerator=denominator=0.0;

+// 	clock_t time, dtime;

+// 	time = clock();

+	for (int t=0; t<this->T; t++)

+	{

+		numerator += weights[t] * this->obs[t];

+		denominator += weights[t];

+	}

+	this->lambda = numerator/denominator; // Update of size is now done with updated lambda

+// 	dtime = clock() - time;

+// 	FILE_LOG(logDEBUG1) << "updateL(): "<<dtime<< " clicks";

+	FILE_LOG(logDEBUG1) << "l = " << this->lambda;

+

+}

+

+void Poisson::update_constrained(double** weights, int fromState, int toState)

+{

+	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	FILE_LOG(logDEBUG1) << "l = "<<this->lambda;

+	double numerator, denominator;

+	// Update lambda

+	numerator=denominator=0.0;

+// 	clock_t time, dtime;

+// 	time = clock();

+	for (int i=0; i<toState-fromState; i++)

+	{

+		for (int t=0; t<this->T; t++)

+		{

+			numerator += weights[i+fromState][t] * this->obs[t];

+			denominator += weights[i+fromState][t] * (i+1);

+		}

+	}

+	this->lambda = numerator/denominator; // Update of size is now done with old lambda

+// 	dtime = clock() - time;

+// 	FILE_LOG(logDEBUG1) << "updateL(): "<<dtime<< " clicks";

+	FILE_LOG(logDEBUG1) << "l = "<<this->lambda;

+

+}

+

+// Getter and Setter ------------------------------------------

+double Poisson::get_mean()

+{

+	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	return( this->lambda );

+}

+

+void Poisson::set_mean(double mean)

+{

+	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	this->lambda = mean;

+}

+

+double Poisson::get_variance()

+{

+	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	return( this->lambda );