@@ -2,10 +2,16 @@ CHANGES IN VERSION 1.5.1

 ------------------------

 SIGNIFICANT USER-LEVEL CHANGES

+

+    o GC correction is now done with a loess-fit by default instead of the quadratic fit. This should improve accuracy.

+    

+    o A stepsize for a sliding window can be selected in addition to the binsize. This improves resolution of detected breakpoints.

+    

+    o Breakpoints are reported with confidence intervals in folder BROWSERFILES.

 BUG FIXES

-    o Can use BSgenome*NCBI* for GC correction now (worked only for BSgenome*UCSC* before).

+    o Can use BSgenome*NCBI* for GC correction (worked only for BSgenome*UCSC* before).

 CHANGES IN VERSION 1.3.4

@@ -150,30 +150,3 @@ mergeStrandseqFiles <- function(files, assembly, chromosomes=NULL, pairedEndRead

     return(reads)

 }

-

-

-createBlacklistFromReads <- function(reads, binsize=1e5, quantile.cutoff=0.99) {

-  

-    bins <- binReads(reads, binsizes = binsize)[[1]]

-    # bins.split <- split(bins, bins@seqnames)

-    # blacklist.list <- GRangesList()

-    # for (i1 in 1:length(bins.split)) {

-    #     cutoff <- quantile(bins.split[[i1]]$mcounts, quantile.cutoff)

-    #     blacklist.list[[i1]] <- bins.split[[i1]][bins.split[[i1]]$mcounts > cutoff]

-    # }

-    # blacklist <- unlist(blacklist.list, use.names=FALSE)

-    cutoff <- quantile(bins$mcounts, quantile.cutoff)

-    blacklist <- bins[bins$mcounts >= cutoff]

-    # print(as.data.frame(blacklist))

-    blacklist <- reduce(blacklist)

-    # print(as.data.frame(blacklist))

-    frac <- sum(as.numeric(width(blacklist))) / sum(as.numeric(seqlengths(blacklist)))

-    message("Fraction of genome in blacklist: ", round(frac,4))

-    

-    ggplt <- plotProfile(bins, both.strands=TRUE)

-    ggplt <- ggplt + geom_hline(data=data.frame(y=-cutoff), mapping=aes_string(yintercept='y'), col='red')

-    print(ggplt)

-    

-    return(blacklist)

-

-}

@@ -528,7 +528,7 @@ univariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard",

 		result$bins$state <- state.labels[states.step]

 		result$bins$copy.number <- multiplicity[as.character(result$bins$state)]

 	## Counts

-		result$binned.data <- binned.data.list

+		result$bincounts <- binned.data.list

 	## Segmentation

 		message("Making segmentation ...", appendLF=FALSE)

 		ptm <- proc.time()

@@ -563,10 +563,16 @@ univariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard",

 bivariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard", max.time=-1, max.iter=-1, num.trials=1, eps.try=NULL, num.threads=1, count.cutoff.quantile=0.999, states=c("zero-inflation",paste0(0:10,"-somy")), most.frequent.state="1-somy", algorithm='EM', initial.params=NULL, verbosity=1) {

 	## Intercept user input

-  binned.data <- loadFromFiles(binned.data, check.class='GRanges')[[1]]

+  binned.data <- loadFromFiles(binned.data, check.class=c('GRanges','GRangesList'))[[1]]

 	if (is.null(ID)) {

 		ID <- attr(binned.data, 'ID')

 	}

+  if (class(binned.data) == 'GRangesList') {

+    binned.data.list <- binned.data

+    binned.data <- binned.data.list[[1]]

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

+    binned.data.list <- GRangesList('0'=binned.data)

+  }

 	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")

@@ -585,392 +591,471 @@ bivariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard", m

 	if (!is.null(initial.params)) {

 		init <- 'initial.params'

 	}

-

-	warlist <- list()

-	if (is.null(eps.try)) eps.try <- eps

-

+  	

 	## Variables

-	num.bins <- length(binned.data)

 	algorithm <- factor(algorithm, levels=c('baumWelch','viterbi','EM'))

-

-	## Get counts

-	select <- 'counts'

-	counts <- matrix(c(mcols(binned.data)[,paste0('m',select)], mcols(binned.data)[,paste0('p',select)]), ncol=2, dimnames=list(bin=1:num.bins, strand=c('minus','plus')))

-	maxcounts <- max(counts)

-

-	## Filter high counts out, makes HMM faster

-	count.cutoff <- quantile(counts, count.cutoff.quantile)

-	names.count.cutoff <- names(count.cutoff)

-	count.cutoff <- ceiling(count.cutoff)

-	mask <- counts > count.cutoff

-	counts[mask] <- count.cutoff

-	numfiltered <- length(which(mask))

-	if (numfiltered > 0) {

-		message(paste0("Replaced read counts > ",count.cutoff," (",names.count.cutoff," quantile) by ",count.cutoff," in ",numfiltered," bins. Set option 'count.cutoff.quantile=1' to disable this filtering. This filtering was done to enhance performance."))

-	}

 	### Make return object

-		result <- list()

-		class(result) <- class.bivariate.hmm

-		result$ID <- ID

-		result$bins <- binned.data

+	result <- list()

+	class(result) <- class.bivariate.hmm

+	result$ID <- ID

+	result$bins <- binned.data

 	## Quality info

-		result$qualityInfo <- as.list(getQC(binned.data))

+	result$qualityInfo <- as.list(getQC(binned.data))

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

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

-		warlist[[length(warlist)+1]] <- warning(paste0("ID = ",ID,": All counts in data are zero. No HMM done."))

-		result$warnings <- warlist

-		return(result)

-	} else if (any(counts<0)) {

-		warlist[[length(warlist)+1]] <- warning(paste0("ID = ",ID,": Some counts in data are negative. No HMM done."))

-		result$warnings <- warlist

-		return(result)

-	}

+	warlist <- list()

+	if (is.null(eps.try)) eps.try <- eps

-	if (init=='initial.params') {

-		uni.transitionProbs <- initial.params$univariateParams$transitionProbs

-		uni.startProbs <- initial.params$univariateParams$startProbs

-		distributions <- initial.params$distributions

-		uni.weights <- initial.params$univariateParams$weights

-		uni.states <- names(uni.weights)

-		num.uni.states <- length(uni.states)

-		num.models <- length(distributions)

-		num.bins <- length(initial.params$bins)

-		comb.states <- factor(names(initial.params$startProbs), levels=names(initial.params$startProbs))

-		num.comb.states <- length(comb.states)

-		states.list <- list(minus=initial.params$bins$mstate, plus=initial.params$bins$pstate)

-		comb.states.per.bin <- factor(do.call(paste, states.list), levels=levels(comb.states))

-		A.initial <- initial.params$transitionProbs

-		proba.initial <- initial.params$startProbs

-		use.initial <- TRUE

-	} else {

-		### Stack the strands and run one univariate findCNVs

-		message("")

-		message(paste(rep('-',getOption('width')), collapse=''))

-		binned.data.minus <- binned.data

-		strand(binned.data.minus) <- '-'

-		binned.data.minus$counts <- counts[,'minus', drop=FALSE]

-		binned.data.plus <- binned.data

-		strand(binned.data.plus) <- '+'

-		binned.data.plus$counts <- counts[,'plus', drop=FALSE]

-		binned.data.stacked <- c(binned.data.minus, binned.data.plus)

-		mask.attributes <- c('qualityInfo', 'ID', 'min.mapq')

-		attributes(binned.data.stacked)[mask.attributes] <- attributes(binned.data)[mask.attributes]

+  ### Arrays for finding maximum posterior for each bin between offsets

+  ## Make bins with offset

+  ptm <- startTimedMessage("Making bins with offsets ...")

+  if (length(binned.data.list) > 1) {

+    stepbins <- disjoin(unlist(binned.data.list, use.names=FALSE))

+  } else {

+    stepbins <- binned.data.list[[1]]

+    mcols(stepbins) <- NULL

+  }

+  amaxPosterior.step <- array(0, dim = c(length(stepbins), 2), dimnames = list(bin=NULL, offset=c('previousOffsets', 'currentOffset'))) # to store maximum posterior for current and max-of-previous offsets

+  astates.step <- array(0, dim = c(length(stepbins), 2), dimnames = list(bin=NULL, offset=c('previousOffsets', 'currentOffset'))) # to store states for current and max-of-previous offsets

+  stopTimedMessage(ptm)

+  

+  ### Loop over offsets ###

+  for (istep in 1:length(binned.data.list)) {

+    binned.data <- binned.data.list[[istep]]

+  	num.bins <- length(binned.data)

+    if (istep > 1) {

+      ptm.offset <- startTimedMessage("Obtaining states for step = ", istep, " ...")

+      ## Run only one iteration (no updating) if we are already over istep==1

+      initial.params <- result

+      init <- 'initial.params'

+    	algorithm <- factor('baumWelch', levels=c('baumWelch','viterbi','EM'))

+      num.trials <- 1

+      verbosity <- 0

+    }

+

+  	## Get counts

+  	select <- 'counts'

+  	counts <- matrix(c(mcols(binned.data)[,paste0('m',select)], mcols(binned.data)[,paste0('p',select)]), ncol=2, dimnames=list(bin=1:num.bins, strand=c('minus','plus')))

+  	maxcounts <- max(counts)

+  

+  	## Filter high counts out, makes HMM faster

+  	count.cutoff <- quantile(counts, count.cutoff.quantile)

+  	names.count.cutoff <- names(count.cutoff)

+  	count.cutoff <- ceiling(count.cutoff)

+  	mask <- counts > count.cutoff

+  	counts[mask] <- count.cutoff

+  	numfiltered <- length(which(mask))

+  	if (numfiltered > 0 & istep == 1) {

+  		message(paste0("Replaced read counts > ",count.cutoff," (",names.count.cutoff," quantile) by ",count.cutoff," in ",numfiltered," bins. Set option 'count.cutoff.quantile=1' to disable this filtering. This filtering was done to enhance performance."))

+  	}

+  

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

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

+  		warlist[[length(warlist)+1]] <- warning(paste0("ID = ",ID,": All counts in data are zero. No HMM done."))

+  		result$warnings <- warlist

+  		return(result)

+  	} else if (any(counts<0)) {

+  		warlist[[length(warlist)+1]] <- warning(paste0("ID = ",ID,": Some counts in data are negative. No HMM done."))

+  		result$warnings <- warlist

+  		return(result)

+  	}

-		# if (method == 'HMM') {

+  	if (init=='initial.params') {

+  		uni.transitionProbs <- initial.params$univariateParams$transitionProbs

+  		uni.startProbs <- initial.params$univariateParams$startProbs

+  		distributions <- initial.params$distributions

+  		uni.weights <- initial.params$univariateParams$weights

+  		uni.states <- names(uni.weights)

+  		num.uni.states <- length(uni.states)

+  		num.models <- length(distributions)

+  		comb.states <- factor(names(initial.params$startProbs), levels=names(initial.params$startProbs))

+  		num.comb.states <- length(comb.states)

+  		states.list <- list(minus=initial.params$bins$mstate, plus=initial.params$bins$pstate)

+  		comb.states.per.bin <- factor(do.call(paste, states.list), levels=levels(comb.states))

+  		A.initial <- initial.params$transitionProbs

+  		proba.initial <- initial.params$startProbs

+  		use.initial <- TRUE

+  	} else {

+  		### Stack the strands and run one univariate findCNVs

+  		message("")

+  		message(paste(rep('-',getOption('width')), collapse=''))

+  		binned.data.minus <- binned.data

+  		strand(binned.data.minus) <- '-'

+  		binned.data.minus$counts <- counts[,'minus', drop=FALSE]

+  		binned.data.plus <- binned.data

+  		strand(binned.data.plus) <- '+'

+  		binned.data.plus$counts <- counts[,'plus', drop=FALSE]

+  		binned.data.stacked <- c(binned.data.minus, binned.data.plus)

+  		mask.attributes <- c('qualityInfo', 'ID', 'min.mapq')

+  		attributes(binned.data.stacked)[mask.attributes] <- attributes(binned.data)[mask.attributes]

+  

   		message("Running univariate HMM")

   		model.stacked <- univariate.findCNVs(binned.data.stacked, ID, eps=eps, init=init, max.time=max.time, max.iter=max.iter, num.trials=num.trials, eps.try=eps.try, num.threads=num.threads, count.cutoff.quantile=1, states=states, most.frequent.state=most.frequent.state)

   		if (is.na(model.stacked$convergenceInfo$error)) {

   		    result$warnings <- model.stacked$warnings

   		    return(result)

   		}

-# 		} else if (method == 'dnacopy') {

-#   		message("Running DNAcopy")

-#   		model.stacked <- DNAcopy.findCNVs(binned.data.stacked, ID, CNgrid.start=0.5, count.cutoff.quantile=1)

-# 		}

-		model.minus <- model.stacked

-		model.minus$bins <- model.minus$bins[strand(model.minus$bins)=='-']

-		model.minus$segments <- model.minus$segments[strand(model.minus$segments)=='-']

-		model.plus <- model.stacked

-		model.plus$bins <- model.plus$bins[strand(model.plus$bins)=='+']

-		model.plus$segments <- model.plus$segments[strand(model.plus$segments)=='+']

-

-		models <- list(minus=model.minus, plus=model.plus)

-

-		## Extract counts and other stuff

-		uni.transitionProbs <- lapply(models, '[[', 'transitionProbs')[[1]]

-		uni.startProbs <- lapply(models, '[[', 'startProbs')[[1]]

-		distributions <- lapply(models, '[[', 'distributions')

-		uni.weights <- lapply(models, '[[', 'weights')[[1]]

-		uni.states <- names(uni.weights)

-		num.uni.states <- length(uni.states)

-		num.models <- length(distributions)

-		num.bins <- length(models[[1]]$bins)

-		comb.states <- vector()

-		levels.state <- levels(models[[1]]$bins$state)

-		for (i1 in 1:length(levels.state)) {

-			for (i2 in 1:length(levels.state)) {

-				comb.state <- paste(levels.state[i1], levels.state[i2])

-				comb.states[length(comb.states)+1] <- comb.state

-			}

-		}

-		comb.states <- factor(comb.states, levels=comb.states)

-		num.comb.states <- length(comb.states)

-		states.list <- lapply(models, function(x) { x$bins$state })

-		comb.states.per.bin <- factor(do.call(paste, states.list), levels=levels(comb.states))

-		A.initial <- double(length=num.comb.states*num.comb.states)

-		proba.initial <- double(length=num.comb.states)

-		use.initial <- FALSE

-	}

-

-	### Prepare the multivariate HMM

-	message("")

-	message(paste(rep('-',getOption('width')), collapse=''))

-	message("Preparing bivariate HMM\n")

-

-	## Pre-compute z-values for each number of counts

-	ptm <- startTimedMessage("Computing pre z-matrix...")

-	z.per.count <- array(NA, dim=c(maxcounts+1, num.models, num.uni.states), dimnames=list(counts=0:maxcounts, strand=names(distributions), state=uni.states))

-	xcounts <- 0:maxcounts

-	for (istrand in 1:num.models) {

-		for (istate in 1:num.uni.states) {

-			if (distributions[[istrand]][istate,'type']=='dnbinom') {

-				size <- distributions[[istrand]][istate,'size']

-				prob <- distributions[[istrand]][istate,'prob']

-				u <- stats::pnbinom(xcounts, size, prob)

-			} else if (distributions[[istrand]][istate,'type']=='dbinom') {

-				size <- distributions[[istrand]][istate,'size']

-				prob <- distributions[[istrand]][istate,'prob']

-				u <- stats::pbinom(xcounts, size, prob)

-			} else if (distributions[[istrand]][istate,'type']=='delta') {

-				u <- rep(1, length(xcounts))

-			} else if (distributions[[istrand]][istate,'type']=='dgeom') {

-				prob <- distributions[[istrand]][istate,'prob']

-				u <- stats::pgeom(xcounts, prob)

-			}

-			qnorm_u <- stats::qnorm(u)

-			mask <- qnorm_u==Inf

-			qnorm_u[mask] <- stats::qnorm(1-1e-16)

-			z.per.count[ , istrand, istate] <- qnorm_u

-		}

-	}

-	stopTimedMessage(ptm)

-

-	## Compute the z matrix

-	ptm <- startTimedMessage("Transfering values into z-matrix...")

-	z.per.bin = array(NA, dim=c(num.bins, num.models, num.uni.states), dimnames=list(bin=1:num.bins, strand=names(distributions), state=uni.states))

-	for (istrand in 1:num.models) {

-		for (istate in 1:num.uni.states) {

-			z.per.bin[ , istrand, istate] = z.per.count[counts[,istrand]+1, istrand, istate]

+  		model.minus <- model.stacked

+  		model.minus$bins <- model.minus$bins[strand(model.minus$bins)=='-']

+  		model.minus$segments <- model.minus$segments[strand(model.minus$segments)=='-']

+  		model.plus <- model.stacked

+  		model.plus$bins <- model.plus$bins[strand(model.plus$bins)=='+']

+  		model.plus$segments <- model.plus$segments[strand(model.plus$segments)=='+']

+  

+  		models <- list(minus=model.minus, plus=model.plus)

+  

+  		## Extract counts and other stuff

+  		uni.transitionProbs <- lapply(models, '[[', 'transitionProbs')[[1]]

+  		uni.startProbs <- lapply(models, '[[', 'startProbs')[[1]]

+  		distributions <- lapply(models, '[[', 'distributions')

+  		uni.weights <- lapply(models, '[[', 'weights')[[1]]

+  		uni.states <- names(uni.weights)

+  		num.uni.states <- length(uni.states)

+  		num.models <- length(distributions)

+  		num.bins <- length(models[[1]]$bins)

+  		comb.states <- vector()

+  		levels.state <- levels(models[[1]]$bins$state)

+  		for (i1 in 1:length(levels.state)) {

+  			for (i2 in 1:length(levels.state)) {

+  				comb.state <- paste(levels.state[i1], levels.state[i2])

+  				comb.states[length(comb.states)+1] <- comb.state

+  			}

+  		}

+  		comb.states <- factor(comb.states, levels=comb.states)

+  		num.comb.states <- length(comb.states)

+  		states.list <- lapply(models, function(x) { x$bins$state })

+  		comb.states.per.bin <- factor(do.call(paste, states.list), levels=levels(comb.states))

+  		A.initial <- double(length=num.comb.states*num.comb.states)

+  		proba.initial <- double(length=num.comb.states)

+  		use.initial <- FALSE

+  	}

+  

+  	### Prepare the multivariate HMM

+  	if (verbosity >= 1 & istep == 1) {

+    	message("")

+    	message(paste(rep('-',getOption('width')), collapse=''))

+    	message("Preparing bivariate HMM\n")

+  	}

+  

+  	## Pre-compute z-values for each number of counts

+		if (verbosity >=1) {

+    	ptm <- startTimedMessage("Computing pre z-matrix...")

 		}

-	}

-	remove(z.per.count)

-	stopTimedMessage(ptm)

-

-	## Calculate correlation matrix

-	ptm <- startTimedMessage("Computing inverse of correlation matrix...")

-	correlationMatrix <- array(0, dim=c(num.models,num.models,num.comb.states), dimnames=list(strand=names(distributions), strand=names(distributions), comb.state=comb.states))

-	correlationMatrixInverse <- correlationMatrix

-	determinant <- rep(0, num.comb.states)

-	names(determinant) <- comb.states

-	usestateTF <- rep(NA, num.comb.states) # TRUE, FALSE vector for usable states

-	names(usestateTF) <- comb.states

-	for (comb.state in comb.states) {

-		state <- strsplit(comb.state, ' ')[[1]]

-		mask <- which(comb.states.per.bin==comb.state)

-		# Subselect z

-		z.temp <- matrix(NA, ncol=num.models, nrow=length(mask))

-		for (istrand in 1:num.models) {

-			z.temp[,istrand] <- z.per.bin[mask, istrand, state[istrand]]

+  	z.per.count <- array(NA, dim=c(maxcounts+1, num.models, num.uni.states), dimnames=list(counts=0:maxcounts, strand=names(distributions), state=uni.states))

+  	xcounts <- 0:maxcounts

+  	for (istrand in 1:num.models) {

+  		for (istate in 1:num.uni.states) {

+  			if (distributions[[istrand]][istate,'type']=='dnbinom') {

+  				size <- distributions[[istrand]][istate,'size']

+  				prob <- distributions[[istrand]][istate,'prob']

+  				u <- stats::pnbinom(xcounts, size, prob)

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

+  				size <- distributions[[istrand]][istate,'size']

+  				prob <- distributions[[istrand]][istate,'prob']

+  				u <- stats::pbinom(xcounts, size, prob)

+  			} else if (distributions[[istrand]][istate,'type']=='delta') {

+  				u <- rep(1, length(xcounts))

+  			} else if (distributions[[istrand]][istate,'type']=='dgeom') {

+  				prob <- distributions[[istrand]][istate,'prob']

+  				u <- stats::pgeom(xcounts, prob)

+  			}

+  			qnorm_u <- stats::qnorm(u)

+  			mask <- qnorm_u==Inf

+  			qnorm_u[mask] <- stats::qnorm(1-1e-16)

+  			z.per.count[ , istrand, istate] <- qnorm_u

+  		}

+  	}

+		if (verbosity >=1) {

+    	stopTimedMessage(ptm)

 		}

-		temp <- tryCatch({

-			if (nrow(z.temp) > 1) {

-				correlationMatrix[,,comb.state] <- cor(z.temp)

-				determinant[comb.state] <- det( correlationMatrix[,,comb.state] )

-				correlationMatrixInverse[,,comb.state] <- solve(correlationMatrix[,,comb.state])

-				usestateTF[comb.state] <- TRUE

-			} else {

-				correlationMatrix[,,comb.state] <- diag(num.models)

-				determinant[comb.state] <- det( correlationMatrix[,,comb.state] )

-				correlationMatrixInverse[,,comb.state] <- solve(correlationMatrix[,,comb.state])

-				usestateTF[comb.state] <- TRUE

-			}

-			0

-		}, warning = function(war) {

-			1

-		}, error = function(err) {

-			1

-		})

-		if (temp!=0) {

-			correlationMatrix[,,comb.state] <- diag(num.models)

-			determinant[comb.state] <- det( correlationMatrix[,,comb.state] )

-			correlationMatrixInverse[,,comb.state] <- solve(correlationMatrix[,,comb.state])

-			usestateTF[comb.state] <- TRUE

+  

+  	## Compute the z matrix

+		if (verbosity >=1) {

+    	ptm <- startTimedMessage("Transfering values into z-matrix...")

 		}

-	}

-	stopTimedMessage(ptm)

-

-	# Use only states with valid correlationMatrixInverse (all states are valid in this version)

-	correlationMatrix = correlationMatrix[,,usestateTF]

-	correlationMatrixInverse = correlationMatrixInverse[,,usestateTF]

-	comb.states = comb.states[usestateTF]

-	comb.states <- droplevels(comb.states)

-	determinant = determinant[usestateTF]

-	num.comb.states <- length(comb.states)

-

-	### Calculate multivariate densities for each state ###

-	ptm <- startTimedMessage("Calculating multivariate densities...")

-	densities <- matrix(1, ncol=num.comb.states, nrow=num.bins, dimnames=list(bin=1:num.bins, comb.state=comb.states))

-	for (comb.state in comb.states) {

-		istate <- which(comb.state==comb.states)

-		state <- strsplit(comb.state, ' ')[[1]]

-		z.temp <- matrix(NA, ncol=num.models, nrow=num.bins)

-		product <- 1

-		for (istrand in 1:num.models) {

-			z.temp[,istrand] <- z.per.bin[, istrand, state[istrand]]

-			if (distributions[[istrand]][state[istrand],'type'] == 'dnbinom') {

-				size <- distributions[[istrand]][state[istrand],'size']

-				prob <- distributions[[istrand]][state[istrand],'prob']

-				product <- product * stats::dnbinom(counts[,istrand], size, prob)

-			} else if (distributions[[istrand]][state[istrand],'type'] == 'dgeom') {

-				prob <- distributions[[istrand]][state[istrand],'prob']

-				product <- product * stats::dgeom(counts[,istrand], prob)

-			} else if (distributions[[istrand]][state[istrand],'type'] == 'delta') {

-				product <- product * ifelse(counts[,istrand]==0, 1, 0)

-			}

+  	z.per.bin = array(NA, dim=c(num.bins, num.models, num.uni.states), dimnames=list(bin=1:num.bins, strand=names(distributions), state=uni.states))

+  	for (istrand in 1:num.models) {

+  		for (istate in 1:num.uni.states) {

+  			z.per.bin[ , istrand, istate] = z.per.count[counts[,istrand]+1, istrand, istate]

+  		}

+  	}

+  	remove(z.per.count)

+		if (verbosity >=1) {

+    	stopTimedMessage(ptm)

 		}

-		exponent <- -0.5 * apply( ( z.temp %*% (correlationMatrixInverse[ , , istate] - diag(num.models)) ) * z.temp, 1, sum)

-		exponent[is.nan(exponent)] <- 0

-		densities[,istate] <- product * determinant[istate]^(-0.5) * exp( exponent )

-	}

-	# Check if densities are > 1

-# 	if (any(densities>1)) stop("Densities > 1")

-# 	if (any(densities<0)) stop("Densities < 0")

-	densities[densities>1] <- 1

-	densities[densities<0] <- 0

-	# Check if densities are 0 everywhere in some bins

-	check <- which(apply(densities, 1, sum) == 0)

-	if (length(check)>0) {

-		if (check[1]==1) {

-			densities[1,] <- rep(1e-10, ncol(densities))

-			check <- check[-1]

+  

+  	## Calculate correlation matrix

+		if (verbosity >=1) {

+    	ptm <- startTimedMessage("Computing inverse of correlation matrix...")

 		}

-		for (icheck in check) {

-			densities[icheck,] <- densities[icheck-1,]

+  	correlationMatrix <- array(0, dim=c(num.models,num.models,num.comb.states), dimnames=list(strand=names(distributions), strand=names(distributions), comb.state=comb.states))

+  	correlationMatrixInverse <- correlationMatrix

+  	determinant <- rep(0, num.comb.states)

+  	names(determinant) <- comb.states

+  	usestateTF <- rep(NA, num.comb.states) # TRUE, FALSE vector for usable states

+  	names(usestateTF) <- comb.states

+  	for (comb.state in comb.states) {

+  		state <- strsplit(comb.state, ' ')[[1]]

+  		mask <- which(comb.states.per.bin==comb.state)

+  		# Subselect z

+  		z.temp <- matrix(NA, ncol=num.models, nrow=length(mask))

+  		for (istrand in 1:num.models) {

+  			z.temp[,istrand] <- z.per.bin[mask, istrand, state[istrand]]

+  		}

+  		temp <- tryCatch({

+  			if (nrow(z.temp) > 1) {

+  				correlationMatrix[,,comb.state] <- cor(z.temp)

+  				determinant[comb.state] <- det( correlationMatrix[,,comb.state] )

+  				correlationMatrixInverse[,,comb.state] <- solve(correlationMatrix[,,comb.state])

+  				usestateTF[comb.state] <- TRUE

+  			} else {

+  				correlationMatrix[,,comb.state] <- diag(num.models)

+  				determinant[comb.state] <- det( correlationMatrix[,,comb.state] )

+  				correlationMatrixInverse[,,comb.state] <- solve(correlationMatrix[,,comb.state])

+  				usestateTF[comb.state] <- TRUE

+  			}

+  			0

+  		}, warning = function(war) {

+  			1

+  		}, error = function(err) {

+  			1

+  		})

+  		if (temp!=0) {

+  			correlationMatrix[,,comb.state] <- diag(num.models)

+  			determinant[comb.state] <- det( correlationMatrix[,,comb.state] )

+  			correlationMatrixInverse[,,comb.state] <- solve(correlationMatrix[,,comb.state])

+  			usestateTF[comb.state] <- TRUE

+  		}

+  	}

+		if (verbosity >=1) {

+    	stopTimedMessage(ptm)

 		}

-	}

-	stopTimedMessage(ptm)

-		

-	### Define cleanup behaviour ###

-	on.exit(.C("C_multivariate_cleanup", as.integer(num.comb.states), PACKAGE = 'AneuFinder'))

-

-	### Run the multivariate HMM

-	# Call the C function

-	hmm <- .C("C_multivariate_hmm",

-		densities = as.double(densities), # double* D

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

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

-		num.strands = as.integer(num.models), # int* Nmod

-		comb.states = as.integer(comb.states), # int* comb_states

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

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

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

-		states = integer(length=num.bins), # int* states

-		A = double(length=num.comb.states*num.comb.states), # double* A

-		proba = double(length=num.comb.states), # double* proba

-		loglik = double(length=1), # double* loglik

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

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

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

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

-		error = as.integer(0), # error handling

-		algorithm = as.integer(algorithm), # int* algorithm

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

-		PACKAGE = 'AneuFinder'

-		)

-			

-	### Check convergence ###

-	war <- NULL

-	if (hmm$loglik.delta > eps) {

-		war <- warning(paste0("ID = ",ID,": HMM did not converge!\n"))

-	}

-	if (hmm$error == 1) {

-		warning(paste0("ID = ",ID,": A NaN occurred during the Baum-Welch! Parameter estimation terminated prematurely. Check your library! The following factors are known to cause this error: 1) Your read counts contain very high numbers. Try again with a lower value for 'count.cutoff.quantile'. 2) Your library contains too few reads in each bin. 3) Your library contains reads for a different genome than it was aligned to."))

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

-		warning(paste0("ID = ",ID,": An error occurred during the Baum-Welch! Parameter estimation terminated prematurely. Check your library"))

-	}

-

-	### Make return object ###

-	if (hmm$error == 0) {

-		result <- list()

-		class(result) <- class.bivariate.hmm

-		result$ID <- ID

-	## Bin coordinates and states

-		result$bins <- binned.data

-		result$bins$state <- comb.states[hmm$states]

-		# Get states as factors in data.frame

-		matrix.states <- matrix(unlist(strsplit(as.character(result$bins$state), split=' ')), byrow=TRUE, ncol=num.models, dimnames=list(bin=1:num.bins, strand=names(distributions)))

-		names <- c('mstate','pstate')

-		for (i1 in 1:num.models) {

-			mcols(result$bins)[names[i1]] <- factor(matrix.states[,i1], levels=uni.states)

+  

+  	# Use only states with valid correlationMatrixInverse (all states are valid in this version)

+  	correlationMatrix = correlationMatrix[,,usestateTF]

+  	correlationMatrixInverse = correlationMatrixInverse[,,usestateTF]

+  	comb.states = comb.states[usestateTF]

+  	comb.states <- droplevels(comb.states)

+  	determinant = determinant[usestateTF]

+  	num.comb.states <- length(comb.states)

+  

+  	### Calculate multivariate densities for each state ###

+		if (verbosity >=1) {

+    	ptm <- startTimedMessage("Calculating multivariate densities...")

 		}

-	## Segmentation

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

-		suppressMessages(

-			result$segments <- as(collapseBins(as.data.frame(result$bins), column2collapseBy='state', columns2drop='width', columns2average=c('counts','mcounts','pcounts')), 'GRanges')

-		)

-		seqlevels(result$segments) <- seqlevels(result$bins) # correct order from as()

-		seqlengths(result$segments) <- seqlengths(result$bins)[names(seqlengths(result$segments))]

-		stopTimedMessage(ptm)

-	## CNV state for both strands combined

-		getnumbers <- function(x) {

-    		x <- sub("zero-inflation", "0-somy", x)

-    		x <- as.numeric(sub("-somy", "", x))

+  	densities <- matrix(1, ncol=num.comb.states, nrow=num.bins, dimnames=list(bin=1:num.bins, comb.state=comb.states))

+  	for (comb.state in comb.states) {

+  		istate <- which(comb.state==comb.states)

+  		state <- strsplit(comb.state, ' ')[[1]]

+  		z.temp <- matrix(NA, ncol=num.models, nrow=num.bins)

+  		product <- 1

+  		for (istrand in 1:num.models) {

+  			z.temp[,istrand] <- z.per.bin[, istrand, state[istrand]]

+  			if (distributions[[istrand]][state[istrand],'type'] == 'dnbinom') {

+  				size <- distributions[[istrand]][state[istrand],'size']

+  				prob <- distributions[[istrand]][state[istrand],'prob']

+  				product <- product * stats::dnbinom(counts[,istrand], size, prob)

+  			} else if (distributions[[istrand]][state[istrand],'type'] == 'dgeom') {

+  				prob <- distributions[[istrand]][state[istrand],'prob']

+  				product <- product * stats::dgeom(counts[,istrand], prob)

+  			} else if (distributions[[istrand]][state[istrand],'type'] == 'delta') {

+  				product <- product * ifelse(counts[,istrand]==0, 1, 0)

+  			}

+  		}

+  		exponent <- -0.5 * apply( ( z.temp %*% (correlationMatrixInverse[ , , istate] - diag(num.models)) ) * z.temp, 1, sum)

+  		exponent[is.nan(exponent)] <- 0

+  		densities[,istate] <- product * determinant[istate]^(-0.5) * exp( exponent )

+  	}

+  	# Check if densities are > 1

+  # 	if (any(densities>1)) stop("Densities > 1")

+  # 	if (any(densities<0)) stop("Densities < 0")

+  	densities[densities>1] <- 1

+  	densities[densities<0] <- 0

+  	# Check if densities are 0 everywhere in some bins

+  	check <- which(apply(densities, 1, sum) == 0)

+  	if (length(check)>0) {

+  		if (check[1]==1) {

+  			densities[1,] <- rep(1e-10, ncol(densities))

+  			check <- check[-1]

+  		}

+  		for (icheck in check) {

+  			densities[icheck,] <- densities[icheck-1,]

+  		}

+  	}

+		if (verbosity >=1) {

+    	stopTimedMessage(ptm)

 		}

-  	inistates <-  suppressWarnings( initializeStates(unique(c(states, paste0(sort(unique(getnumbers(result$bins$mstate) + getnumbers(result$bins$pstate))),"-somy")))) )

-  	multiplicity <- inistates$multiplicity

-  	state.labels <- inistates$states

-		# Bins

-		state <- multiplicity[as.character(result$bins$mstate)] + multiplicity[as.character(result$bins$pstate)]

-		# state[state>max(multiplicity)] <- max(multiplicity)

-		multiplicity.inverse <- names(multiplicity)

-		names(multiplicity.inverse) <- multiplicity

-		state <- multiplicity.inverse[as.character(state)]

-		state[(result$bins$mstate=='0-somy' | result$bins$pstate=='0-somy') & state=='zero-inflation'] <- '0-somy'

-    result$bins$state <- factor(state, levels=names(multiplicity))

-    result$bins$copy.number <- multiplicity[as.character(result$bins$state)]

-    result$bins$mcopy.number <- multiplicity[as.character(result$bins$mstate)]

-    result$bins$pcopy.number <- multiplicity[as.character(result$bins$pstate)]

-		# Segments

-		str <- strsplit(as.character(result$segments$state),' ')

-		result$segments$mstate <- factor(unlist(lapply(str, '[[', 1)), levels=state.labels)

-		result$segments$pstate <- factor(unlist(lapply(str, '[[', 2)), levels=state.labels)

-		state <- multiplicity[result$segments$mstate] + multiplicity[result$segments$pstate]

-		state[state>max(multiplicity)] <- max(multiplicity)

-		multiplicity.inverse <- names(multiplicity)

-		names(multiplicity.inverse) <- multiplicity

-		state <- multiplicity.inverse[as.character(state)]

-		state[(result$segments$mstate=='0-somy' | result$segments$pstate=='0-somy') & state=='zero-inflation'] <- '0-somy'

-    result$segments$state <- factor(state, levels=names(multiplicity))

-    result$segments$copy.number <- multiplicity[as.character(result$segments$state)]

-    result$segments$mcopy.number <- multiplicity[as.character(result$segments$mstate)]

-    result$segments$pcopy.number <- multiplicity[as.character(result$segments$pstate)]

-		## Parameters

-			# Weights

-			tstates <- table(result$bins$state)

-			result$weights <- tstates/sum(tstates)

-			# Transition matrices

-			result$transitionProbs <- matrix(hmm$A, ncol=num.comb.states)

-			colnames(result$transitionProbs) <- comb.states

-			rownames(result$transitionProbs) <- comb.states

-			result$transitionProbs.initial <- matrix(hmm$A.initial, ncol=num.comb.states)

-			colnames(result$transitionProbs.initial) <- comb.states

-			rownames(result$transitionProbs.initial) <- comb.states

-			# Initial probs

-			result$startProbs <- hmm$proba

-			names(result$startProbs) <- comb.states

-			result$startProbs.initial <- hmm$proba.initial

-			names(result$startProbs.initial) <- comb.states

-			# Distributions

-			result$distributions <- distributions

-		## Convergence info

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

-			result$convergenceInfo <- convergenceInfo

-		## Quality info

-  		result$qualityInfo <- as.list(getQC(result))

-		## Univariate infos

-			univariateParams <- list(transitionProbs=uni.transitionProbs, startProbs=uni.startProbs, distributions=distributions[[1]], weights=uni.weights)

-			result$univariateParams <- univariateParams

-	} else if (hmm$error == 1) {

-		warlist[[length(warlist)+1]] <- warning(paste0("ID = ",ID,": A NaN occurred during the Baum-Welch! Parameter estimation terminated prematurely. Check your library! The following factors are known to cause this error: 1) Your read counts contain very high numbers. Try again with a lower value for 'count.cutoff.quantile'. 2) Your library contains too few reads in each bin. 3) Your library contains reads for a different genome than it was aligned to."))

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

-		warlist[[length(warlist)+1]] <- warning(paste0("ID = ",ID,": An error occurred during the Baum-Welch! Parameter estimation terminated prematurely. Check your library"))

-	}

-

-	## Issue warnings

-	result$warnings <- warlist

+  		

+  	### Define cleanup behaviour ###

+  	on.exit(.C("C_multivariate_cleanup", as.integer(num.comb.states), PACKAGE = 'AneuFinder'))

+  

+  	### Run the multivariate HMM

+  	# Call the C function

+  	hmm <- .C("C_multivariate_hmm",

+  		densities = as.double(densities), # double* D

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

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

+  		num.strands = as.integer(num.models), # int* Nmod

+  		comb.states = as.integer(comb.states), # int* comb_states

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

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

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

+			maxPosterior = double(length=num.bins), # double* maxPosterior

+  		states = integer(length=num.bins), # int* states

+  		A = double(length=num.comb.states*num.comb.states), # double* A

+  		proba = double(length=num.comb.states), # double* proba

+  		loglik = double(length=1), # double* loglik

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

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

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

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

+  		error = as.integer(0), # error handling

+  		algorithm = as.integer(algorithm), # int* algorithm

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

+  		PACKAGE = 'AneuFinder'

+  		)

+  			

+  	### Check convergence ###

+  	war <- NULL

+  	if (hmm$loglik.delta > eps) {

+  		war <- warning(paste0("ID = ",ID,": HMM did not converge!\n"))

+  	}

+  	if (hmm$error == 1) {

+  		warning(paste0("ID = ",ID,": A NaN occurred during the Baum-Welch! Parameter estimation terminated prematurely. Check your library! The following factors are known to cause this error: 1) Your read counts contain very high numbers. Try again with a lower value for 'count.cutoff.quantile'. 2) Your library contains too few reads in each bin. 3) Your library contains reads for a different genome than it was aligned to."))

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

+  		warning(paste0("ID = ",ID,": An error occurred during the Baum-Welch! Parameter estimation terminated prematurely. Check your library"))

+  	}

+  	if (istep == 1) {

+    	### Make return object ###

+    	if (hmm$error == 0) {

+    		result <- list()

+    		class(result) <- class.bivariate.hmm

+    		result$ID <- ID

+      	## Bin coordinates and states

+    		result$bins <- binned.data

+    		result$bins$state <- comb.states[hmm$states]

+    		# Get states as factors in data.frame

+    		matrix.states <- matrix(unlist(strsplit(as.character(result$bins$state), split=' ')), byrow=TRUE, ncol=num.models, dimnames=list(bin=1:num.bins, strand=names(distributions)))

+    		names <- c('mstate','pstate')

+    		for (i1 in 1:num.models) {

+    			mcols(result$bins)[names[i1]] <- factor(matrix.states[,i1], levels=uni.states)

+    		}

+      	## Segmentation

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

+    		suppressMessages(

+    			result$segments <- as(collapseBins(as.data.frame(result$bins), column2collapseBy='state', columns2drop='width', columns2average=c('counts','mcounts','pcounts')), 'GRanges')

+    		)

+    		seqlevels(result$segments) <- seqlevels(result$bins) # correct order from as()

+    		seqlengths(result$segments) <- seqlengths(result$bins)[names(seqlengths(result$segments))]

+    		stopTimedMessage(ptm)

+      	## CNV state for both strands combined

+    		getnumbers <- function(x) {

+        		x <- sub("zero-inflation", "0-somy", x)

+        		x <- as.numeric(sub("-somy", "", x))

+    		}

+      	inistates <-  suppressWarnings( initializeStates(unique(c(states, paste0(sort(unique(getnumbers(result$bins$mstate) + getnumbers(result$bins$pstate))),"-somy")))) )

+      	multiplicity <- inistates$multiplicity

+      	state.labels <- inistates$states

+    		# Bins

+    		state <- multiplicity[as.character(result$bins$mstate)] + multiplicity[as.character(result$bins$pstate)]

+    		# state[state>max(multiplicity)] <- max(multiplicity)

+    		multiplicity.inverse <- names(multiplicity)

+    		names(multiplicity.inverse) <- multiplicity

+    		state <- multiplicity.inverse[as.character(state)]

+    		state[(result$bins$mstate=='0-somy' | result$bins$pstate=='0-somy') & state=='zero-inflation'] <- '0-somy'

+        result$bins$state <- factor(state, levels=names(multiplicity))

+        result$bins$copy.number <- multiplicity[as.character(result$bins$state)]

+        result$bins$mcopy.number <- multiplicity[as.character(result$bins$mstate)]

+        result$bins$pcopy.number <- multiplicity[as.character(result$bins$pstate)]

+    		# Segments

+    		str <- strsplit(as.character(result$segments$state),' ')

+    		result$segments$mstate <- factor(unlist(lapply(str, '[[', 1)), levels=state.labels)

+    		result$segments$pstate <- factor(unlist(lapply(str, '[[', 2)), levels=state.labels)

+    		state <- multiplicity[result$segments$mstate] + multiplicity[result$segments$pstate]

+    		state[state>max(multiplicity)] <- max(multiplicity)

+    		multiplicity.inverse <- names(multiplicity)

+    		names(multiplicity.inverse) <- multiplicity

+    		state <- multiplicity.inverse[as.character(state)]

+    		state[(result$segments$mstate=='0-somy' | result$segments$pstate=='0-somy') & state=='zero-inflation'] <- '0-somy'

+        result$segments$state <- factor(state, levels=names(multiplicity))

+        result$segments$copy.number <- multiplicity[as.character(result$segments$state)]

+        result$segments$mcopy.number <- multiplicity[as.character(result$segments$mstate)]

+        result$segments$pcopy.number <- multiplicity[as.character(result$segments$pstate)]

+    		## Parameters

+    			# Weights

+    			tstates <- table(result$bins$state)

+    			result$weights <- tstates/sum(tstates)

+    			# Transition matrices

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

+    			colnames(result$transitionProbs) <- comb.states

+    			rownames(result$transitionProbs) <- comb.states

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

+    			colnames(result$transitionProbs.initial) <- comb.states

+    			rownames(result$transitionProbs.initial) <- comb.states

+    			# Initial probs

+    			result$startProbs <- hmm$proba

+    			names(result$startProbs) <- comb.states

+    			result$startProbs.initial <- hmm$proba.initial

+    			names(result$startProbs.initial) <- comb.states

+    			# Distributions

+    			result$distributions <- distributions

+    		## Convergence info

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

+    			result$convergenceInfo <- convergenceInfo

+    		## Quality info

+      		result$qualityInfo <- as.list(getQC(result))

+    		## Univariate infos

+    			univariateParams <- list(transitionProbs=uni.transitionProbs, startProbs=uni.startProbs, distributions=distributions[[1]], weights=uni.weights)

+    			result$univariateParams <- univariateParams

+    	} else if (hmm$error == 1) {

+    		warlist[[length(warlist)+1]] <- warning(paste0("ID = ",ID,": A NaN occurred during the Baum-Welch! Parameter estimation terminated prematurely. Check your library! The following factors are known to cause this error: 1) Your read counts contain very high numbers. Try again with a lower value for 'count.cutoff.quantile'. 2) Your library contains too few reads in each bin. 3) Your library contains reads for a different genome than it was aligned to."))

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

+    		warlist[[length(warlist)+1]] <- warning(paste0("ID = ",ID,": An error occurred during the Baum-Welch! Parameter estimation terminated prematurely. Check your library"))

+    	}

+  	  if (hmm$error != 0) {

+      	## Issue warnings

+      	result$warnings <- warlist

+      	return(result)

+  	  }

+  	} # if (istep==1)

+    	

+    if (istep == 1) { ptm <- startTimedMessage("Collecting counts and posteriors ...") }

+    

+    ## Inflate posteriors, states, counts to new offset

+    ind <- findOverlaps(stepbins, binned.data)

+    astates.step[ind@from, 'currentOffset'] <- hmm$states[ind@to]

+    amaxPosterior.step[ind@from, 'currentOffset'] <- hmm$maxPosterior[ind@to]