@@ -10,7 +10,7 @@ bedGraph2binned,

 findCNVs,

 plot.distribution,

-plot.BAIT,

+plot.chromosomes,

 plot.genome.overview,

 get.state.table,

 get.dendrogram,

@@ -34,7 +34,7 @@ align2binned <- function(file, format, index=file, chrom.length.file, outputfold

 		tab5rows <- 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]

+		classes[c(1:3,6)] <- classes.in.bed[c(1:3,6)]

 		data <- 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]+1), strand=Rle(strand("*"), nrow(data)))	# +1 to match coordinate systems

@@ -86,11 +86,15 @@ align2binned <- function(file, format, index=file, chrom.length.file, outputfold

 	chroms2use <- intersect(chroms2use, names(chrom.lengths))

 	## Determine binsize automatically

-	if (!is.null(reads.per.bin)) {

+# 	if (!is.null(reads.per.bin)) {

 		cat('Automatically determining binsizes...')

 		gr <- GenomicRanges::GRanges(seqnames=Rle(chroms2use),

 																	ranges=IRanges(start=rep(1, length(chroms2use)), end=chrom.lengths[chroms2use]))

-		autodata <- GenomicAlignments::readGAlignmentsFromBam(file, index=index, param=ScanBamParam(what=c("pos"),which=range(gr),flag=scanBamFlag(isDuplicate=F)))

+		if (format=='bam') {

+			autodata <- GenomicAlignments::readGAlignmentsFromBam(file, index=index, param=ScanBamParam(what=c("pos"),which=range(gr),flag=scanBamFlag(isDuplicate=F)))

+		} else if (format=='bed' | format=='bedGraph') {

+			autodata <- data

+		}

 		numreadsperbp <- length(autodata) / sum(as.numeric(chrom.lengths[chroms2use]))

 		## Pad binsizes and reads.per.bin with each others value

 		binsizes.add <- round(reads.per.bin / numreadsperbp, -2)

@@ -98,7 +102,7 @@ align2binned <- function(file, format, index=file, chrom.length.file, outputfold

 		binsizes <- c(binsizes, binsizes.add)

 		reads.per.bin <- c(reads.per.bin.add, reads.per.bin)

 		cat(' done\n')

-	}

+# 	}

 	### Do the loop for all binsizes

 	if (is.null(numbins)) {

@@ -167,7 +171,9 @@ align2binned <- function(file, format, index=file, chrom.length.file, outputfold

 			## Count overlaps

 			cat("counting overlaps...                     \r")

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

-				reads <- GenomicRanges::countOverlaps(i.binned.data, data[seqnames(data)==chromosome])

+				mreads <- GenomicRanges::countOverlaps(i.binned.data, data[seqnames(data)==chromosome & strand(data)=='-'])

+				preads <- GenomicRanges::countOverlaps(i.binned.data, data[seqnames(data)==chromosome & strand(data)=='+'])

+				reads <- mreads + preads

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

 				# Take the max value from all regions that fall into / overlap a given bin as read count

 				midx <- as.matrix(findOverlaps(i.binned.data, data[seqnames(data)==chromosome]))

@@ -182,14 +188,20 @@ align2binned <- function(file, format, index=file, chrom.length.file, outputfold

 					maxvalues[i1] <- max(signal[midx[(max.idx[i1-1]+1):(max.idx[i1]),2]])

 				}

 				reads[read.idx] <- maxvalues

+				mreads <- rep(NA, length(reads))

+				preads <- rep(NA, length(reads))

 			}

 			## Concatenate

 			cat("concatenate...                           \r")

 			if (is.null(numbins)) {

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

+				mcols(i.binned.data)$mreads <- mreads

+				mcols(i.binned.data)$preads <- preads

 			} else {

 				mcols(i.binned.data)$reads <- reads / (end - start)

+				mcols(i.binned.data)$mreads <- mreads / (end - start)

+				mcols(i.binned.data)$preads <- preads / (end - start)

 			}

 			if (separate.chroms==TRUE) {

@@ -197,9 +209,9 @@ align2binned <- function(file, format, index=file, chrom.length.file, outputfold

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

 					## Print to file

 					if (is.null(numbins)) {

-						filename <- paste0(basename(file),"_binsize_",binsize,"_reads.per.bin_",readsperbin,"_",chromosome,".RData")

+						filename <- paste0(basename(file),"_binsize_",format(binsize, scientific=F, trim=T),"_reads.per.bin_",readsperbin,"_",chromosome,"_.RData")

 					} else {

-						filename <- paste0(basename(file),"_numbin_",numbin,"_",chromosome,".RData")

+						filename <- paste0(basename(file),"_numbin_",format(numbin, scientific=F, trim=T),"_",chromosome,"_.RData")

 					}

 					cat("save...                                  \r")

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

@@ -217,9 +229,9 @@ align2binned <- function(file, format, index=file, chrom.length.file, outputfold

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

 				# Print to file

 				if (is.null(numbins)) {

-					filename <- paste(basename(file),"_binsize_",binsize,"_reads.per.bin_",readsperbin,".RData", sep="")

+					filename <- paste0(basename(file),"_binsize_",format(binsize, scientific=F, trim=T),"_reads.per.bin_",readsperbin,"_.RData")

 				} else {

-					filename <- paste(basename(file),"_numbin_",numbin,".RData", sep="")

+					filename <- paste0(basename(file),"_numbin_",format(numbin, scientific=F, trim=T),"_.RData")

 				}

 				cat("Saving to file ...")

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

@@ -51,12 +51,16 @@ check.integer = function(testvar) {

 check.univariate.modellist = function(modellist) {

 	if (!is(modellist,"list")) return(1)

 	for (model in modellist) {

-		if (!is(model,class.aneufinder.hmm)) return(2)

+		if (!is(model,class.univariate.hmm)) return(2)

 	}

 	return(0)

 }

 check.univariate.model = function(model) {

-	if (!is(model,class.aneufinder.hmm)) return(1)

+	if (!is(model,class.univariate.hmm)) return(1)

+	return(0)

+}

+check.multivariate.model = function(model) {

+	if (!is(model,class.multivariate.hmm)) return(1)

 	return(0)

 }

@@ -72,7 +72,7 @@ hmm2GRanges <- function(hmm, reduce=TRUE) {

 # 	library(GenomicRanges)

 	### Check user input ###

-	if (check.univariate.model(hmm)!=0) stop("argument 'hmm' expects a univariate hmm object (type ?hmm for help)")

+	if (check.univariate.model(hmm)!=0 & check.multivariate.model(hmm)!=0) stop("argument 'hmm' expects a univariate or multivariate hmm object (type ?hmm for help)")

 	if (check.logical(reduce)!=0) stop("argument 'reduce' expects TRUE or FALSE")

 	### Create GRanges ###

@@ -94,6 +94,9 @@ hmm2GRanges <- function(hmm, reduce=TRUE) {

 		for (state in ustates) {

 			red.gr <- GenomicRanges::reduce(gr[hmm$states==state])

 			mcols(red.gr)$state <- rep(factor(state, levels=levels),length(red.gr))

+			if (class(hmm)==class.multivariate.hmm) {

+				mcols(red.gr)$state.separate <- matrix(rep(strsplit(as.character(state), split=' ')[[1]], length(red.gr)), byrow=T, nrow=length(red.gr))

+			}

 			red.gr.list[[length(red.gr.list)+1]] <- red.gr

 		}

 		# Merge and sort

@@ -104,6 +107,9 @@ hmm2GRanges <- function(hmm, reduce=TRUE) {

 		mcols(gr)$reads <- hmm$reads

 		mcols(gr)$posteriors <- hmm$posteriors

 		mcols(gr)$state <- hmm$states

+		if (class(hmm)==class.multivariate.hmm) {

+			mcols(gr)$state.separate <- as.matrix(hmm$states.separate)

+		}

 		return(gr)

 	}

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

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

+findCNVs <- function(binned.data, ID, method='univariate', 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, read.cutoff.quantile=0.999) {

+

+	if (method=='univariate') {

+		model <- univariate.findCNVs(binned.data, ID, eps, init, max.time, max.iter, num.trials, eps.try, num.threads, output.if.not.converged, read.cutoff.quantile)

+	} else if (method=='bivariate') {

+		model <- bivariate.findCNVs(binned.data, ID, eps, init, max.time, max.iter, num.trials, eps.try, num.threads, output.if.not.converged, read.cutoff.quantile)

+	}

+

+	return(model)

+

+}

+

+

+univariate.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, read.cutoff.quantile=0.999) {

 	## Intercept user input

 	IDcheck <- ID  #trigger error if not defined

@@ -32,12 +45,14 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="random", max.time=-1, max

 	}

 	# Filter high reads out, makes HMM faster

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

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

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

+	read.cutoff <- as.integer(read.cutoff)

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

+		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. This filtering was done to increase the speed of the HMM and should not affect the results.\n"))

 	}

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

@@ -81,6 +96,7 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="random", max.time=-1, max

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

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

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

+			state.labels = as.integer(state.labels), # int* state_labels

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

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

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

@@ -128,6 +144,7 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="random", max.time=-1, max

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

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

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

+			state.labels = as.integer(state.labels), # int* state_labels

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

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

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

@@ -160,8 +177,9 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="random", max.time=-1, max

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

 	names(hmm$coordinates) <- coordinate.names

 	hmm$seqlengths <- seqlengths(binned.data)

-	class(hmm) <- class.aneufinder.hmm

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

+	class(hmm) <- class.univariate.hmm

+	hmm$states <- state.labels[hmm$states]

+	hmm$state.labels <- state.labels

 	hmm$eps <- eps

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

 	rownames(hmm$A) <- state.labels

@@ -224,3 +242,225 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="random", max.time=-1, max

 		return(hmm)

 	}

 }

+

+

+bivariate.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, read.cutoff.quantile=0.999) {

+

+	### Split into strands and run univariate findCNVs

+	binned.data.plus <- binned.data

+	binned.data.plus$reads <- binned.data$preads

+	plus.model <- univariate.findCNVs(binned.data.plus, ID, eps, init, max.time, max.iter, num.trials, eps.try, num.threads, output.if.not.converged, read.cutoff.quantile)

+	binned.data.minus <- binned.data

+	binned.data.minus$reads <- binned.data$mreads

+	minus.model <- univariate.findCNVs(binned.data.minus, ID, eps, init, max.time, max.iter, num.trials, eps.try, num.threads, output.if.not.converged, read.cutoff.quantile)

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

+

+	### Prepare the multivariate HMM

+		## Extract reads and other stuff

+		coordinates <- models[[1]]$coordinates

+		seqlengths <- models[[1]]$seqlengths

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

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

+		num.uni.states <- models[[1]]$num.states

+		num.models <- length(models)

+		num.bins <- models[[1]]$num.bins

+		comb.states <- paste(levels(models[[1]]$states), levels(models[[2]]$states))

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

+		comb.states.per.bin <- factor(do.call(paste, lapply(models, '[[', 'states')), levels=levels(comb.states))

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

+		reads <- matrix(c(models$plus$reads, models$minus$reads), ncol=num.models, dimnames=list(bin=1:num.bins, strand=names(models)))

+		maxreads <- max(reads)

+

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

+		cat("Computing pre z-matrix...")

+		z.per.read <- array(NA, dim=c(maxreads+1, num.models, num.uni.states))

+		xreads <- 0:maxreads

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

+			model <- models[[istrand]]

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

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

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

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

+					u <- pnbinom(xreads, size, prob)

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

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

+				}

+				qnorm_u <- qnorm(u)

+				mask <- qnorm_u==Inf

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

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

+			}

+		}

+		cat(" done\n")

+

+		## Compute the z matrix

+		cat("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(models), levels(modellist[[1]]$states)))

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

+			model <- models[[istrand]]

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

+				z.per.bin[ , istrand, istate] = z.per.read[reads[,istrand]+1, istrand, istate]

+			}

+		}

+		remove(z.per.read)

+		cat(" done\n")

+

+		## Calculate correlation matrix

+		cat("Computing inverse of correlation matrix...")

+		correlationMatrix <- array(0, dim=c(num.models,num.models,num.comb.states), dimnames=list(strand=names(models), strand=names(models), 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 (state in comb.states) {

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

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

+			# Subselect z

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

+			for (i1 in 1:length(istate)) {

+				z.temp[,i1] <- z.per.bin[mask, i1, istate[i1]]

+			}

+			temp <- tryCatch({

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

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

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

+				if (length(z.temp)>0) {

+					usestateTF[state] <- TRUE

+				} else {

+					usestateTF[state] <- FALSE

+				}

+			}, warning = function(war) {

+				usestateTF[state] <<- FALSE

+				war

+			}, error = function(err) {

+				usestateTF[state] <<- FALSE

+				err

+			})

+		}

+# 		remove(z.per.bin)

+		cat(" done\n")

+

+		# Use nullsomy state anyways

+		usestateTF[1] <- TRUE

+		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

+		cat("Calculating multivariate densities...")

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

+		for (state in comb.states) {

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

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

+			product <- 1

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

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

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

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

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

+					product <- product * dnbinom(reads[,istrand], size, prob)

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

+				}

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

+					densities[,i] <- densities[,i] * ifelse(reads[,istrand]==0, 1, 0)

+				}

+			}

+		}

+		# 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,]

+			}

+		}

+		cat(" done\n")

+		

+

+	### Run the multivariate HMM

+	# Call the C function

+	use.initial <- FALSE

+	hmm <- .C("R_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 = double(length=num.comb.states*num.comb.states), # double* initial_A

+		proba.initial = double(length=num.comb.states), # 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

+		)

+			

+	# Add useful entries

+	class(hmm) <- class.multivariate.hmm

+	hmm$coordinates <- coordinates

+	hmm$seqlengths <- seqlengths

+	hmm$densities <- densities	# reassign because of matrix layout

+	hmm$reads <- reads

+	hmm$states <- comb.states[hmm$states]

+	# Get states as factors in data.frame

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

+		toFactor <- function(column) { factor(column, levels=levels(models[[1]]$states)) }

+		hmm$states.separate <- list()

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

+			hmm$states.separate[[names(models)[i1]]] <- toFactor(matrix.states[,i1])

+		}

+		hmm$states.separate <- as.data.frame(hmm$states.separate)

+	hmm$comb.states <- comb.states

+	hmm$eps <- eps

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

+	colnames(hmm$A) <- comb.states

+	rownames(hmm$A) <- comb.states

+	names(hmm$proba) <- comb.states

+	hmm$ID <- ID

+	hmm$distributions.univariate <- distributions

+	hmm$weights.univariate <- weights

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

+	colnames(hmm$A.initial) <- comb.states

+	rownames(hmm$A.initial) <- comb.states

+	names(hmm$proba.initial) <- comb.states

+	hmm$correlation.matrix <- correlationMatrix

+	hmm$correlation.matrix.inverse <- correlationMatrixInverse

+

+	# Delete redundant entries

+	hmm$use.initial.params <- NULL

+

+	# Check convergence

+	war <- NULL

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

+		war <- warning("HMM did not converge!\n")

+	}

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

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

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

+	}

+

+	hmms <- list(bivariate=hmm, univariate.plus=models$plus, univariate.minus=models$minus)

+	class(hmms) <- class.hmm.list

+	return(hmms)

+

+}

+

+

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

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

 # Some global variables that can be used in all functions

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

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

+state.labels <- factor(c("nullsomy","monosomy","disomy","trisomy","tetrasomy","multisomy"), levels=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"

+class.univariate.hmm <- "aneufinder.univariate.hmm"

+class.multivariate.hmm <- "aneufinder.multivariate.hmm"

+class.hmm.list <- "aneufinder.hmm.list"

 state.colors <- c("mapped"="gray68","nullsomy"="gray20","null-mixed"="gray30","monosomy"="gold3","disomy"="springgreen3","trisomy"="orangered1","tetrasomy"="orangered4","multisomy"="purple3","total"="black")

 get.state.labels <- function() { return(state.labels) }

 get.state.colors <- function() { return(state.colors[state.labels]) }

@@ -48,7 +48,7 @@ plot.distribution <- function(model, state=NULL, chrom=NULL, start=NULL, end=NUL

 		states <- model$states[selectmask]

 		weights <- rep(NA, 3)

 		for (i1 in 1:length(levels(model$states))) {

-			weights[i1] <- length(which(states==state.labels[i1]))

+			weights[i1] <- length(which(states==model$state.labels[i1]))

 		}

 		weights <- weights / length(states)

 	} else {

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

 	ggplt <- ggplot(data.frame(reads)) + geom_histogram(aes(x=reads, y=..density..), binwidth=1, color='black', fill='white') + xlim(0,rightxlim) + theme_bw() + xlab("read count")

 	### Add fits to the histogram

+	c.state.labels <- as.character(model$state.labels)

 	numstates <- model$num.states

 	x <- 0:rightxlim

 	distributions <- list(x)

@@ -92,12 +93,12 @@ plot.distribution <- function(model, state=NULL, chrom=NULL, start=NULL, end=NUL

 		}

 	}

 	distributions <- as.data.frame(distributions)

-	names(distributions) <- c("x",state.labels)

+	names(distributions) <- c("x",c.state.labels)

 	# Total

 	distributions$total <- apply(distributions[-1], 1, sum)

 	# Reshape the data.frame for plotting with ggplot

-	distributions <- reshape(distributions, direction="long", varying=1+1:(numstates+1), v.names="density", timevar="state", times=c(state.labels,"total"))

+	distributions <- reshape(distributions, direction="long", varying=1+1:(numstates+1), v.names="density", timevar="state", times=c(c.state.labels,"total"))

 	### Plot the distributions

 	if (is.null(state)) {

 		ggplt <- ggplt + geom_line(aes(x=x, y=density, group=state, col=state), data=distributions)

@@ -108,9 +109,9 @@ plot.distribution <- function(model, state=NULL, chrom=NULL, start=NULL, end=NUL

 	# Make legend and colors correct

 	lmeans <- round(model$distributions[,'mu'], 2)

 	lvars <- round(model$distributions[,'variance'], 2)

-	legend <- paste(state.labels, ", mean=", lmeans, ", var=", lvars, sep='')

+	legend <- paste(c.state.labels, ", mean=", lmeans, ", var=", lvars, sep='')

 	legend <- c(legend, paste0('total, mean(data)=', round(mean(reads),2), ', var(data)=', round(var(reads),2)))

-	ggplt <- ggplt + scale_color_manual(breaks=c(state.labels, 'total'), values=state.colors[c(state.labels,'total')], labels=legend)

+	ggplt <- ggplt + scale_color_manual(breaks=c(c.state.labels, 'total'), values=state.colors[c(c.state.labels,'total')], labels=legend)

 	ggplt <- ggplt + theme(legend.position=c(1,1), legend.justification=c(1,1))

 	return(ggplt)

@@ -139,9 +140,22 @@ plot.boxplot <- function(model) {

 # ------------------------------------------------------------

-# Plot like BAIT

+# Plot state categorization for all chromosomes

 # ------------------------------------------------------------

-plot.BAIT <- function(model, file='aneufinder_BAIT_plots') {

+plot.chromosomes <- function(model, file=NULL) {

+

+	if (class(model)==class.univariate.hmm) {

+		plot.chromosomes.univariate(model, file=file)

+	} else if (class(model)==class.multivariate.hmm) {

+		plot.chromosomes.bivariate(model, file=file)

+	}

+

+}

+

+# ------------------------------------------------------------

+# Plot state categorization for all chromosomes

+# ------------------------------------------------------------

+plot.chromosomes.univariate <- function(model, file=NULL) {

 	## Convert to GRanges

 	gr <- hmm2GRanges(model, reduce=F)

@@ -157,12 +171,14 @@ plot.BAIT <- function(model, file='aneufinder_BAIT_plots') {

 	library(ggplot2)

 	nrows <- 2	# rows for plotting chromosomes

 	ncols <- ceiling(num.chroms/nrows)

-	png(file=paste0(file, '.png'), width=ncols*6, height=nrows*21, units='cm', res=150)

+	if (!is.null(file)) {

+		pdf(file=paste0(file, '.pdf'), width=ncols*2.4, height=nrows*8.3)

+	}

 	grid.newpage()

 	layout <- matrix(1:((nrows+1)*ncols), ncol=ncols, nrow=nrows+1, byrow=T)

 	pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout), heights=c(1,21,21))))

 	# Main title

-	grid.text(file, vp = viewport(layout.pos.row = 1, layout.pos.col = 1:ncols), gp=gpar(fontsize=26))

+	grid.text(model$ID, vp = viewport(layout.pos.row = 1, layout.pos.col = 1:ncols), gp=gpar(fontsize=26))

 	## Go through chromosomes and plot

 	for (i1 in 1:num.chroms) {

@@ -181,8 +197,10 @@ plot.BAIT <- function(model, file='aneufinder_BAIT_plots') {

 		# Plot the read counts

 		dfplot <- as.data.frame(grl[[i1]])

-		dfplot.points <- dfplot[dfplot$reads>=custom.xlim,]

-		dfplot$reads <- stats::runmed(dfplot$reads, 11)

+		# Set values to great for plotting to limit

+			dfplot$reads[dfplot$reads>=custom.xlim] <- custom.xlim

+			dfplot.points <- dfplot[dfplot$reads>=custom.xlim,]

+			dfplot.points$reads <- rep(custom.xlim, nrow(dfplot.points))

 		empty_theme <- theme(axis.line=element_blank(),

       axis.text.x=element_blank(),

       axis.text.y=element_blank(),

@@ -197,19 +215,111 @@ plot.BAIT <- function(model, file='aneufinder_BAIT_plots') {

       plot.background=element_blank())

 		ggplt <- ggplot(dfplot, aes(x=start, y=reads))	# data

 		ggplt <- ggplt + geom_linerange(aes(ymin=0, ymax=reads, col=state), size=0.2)	# read count

-		ggplt <- ggplt + xlim(0,maxseqlength) + ylim(-0.6*custom.xlim,custom.xlim)	# set x- and y-limits

 		ggplt <- ggplt + geom_rect(ymin=-0.05*custom.xlim-0.1*custom.xlim, ymax=-0.05*custom.xlim, xmin=0, mapping=aes(xmax=max(start)), col='white', fill='gray20')	# chromosome backbone as simple rectangle

-		ggplt <- ggplt + geom_point(data=dfplot.points, mapping=aes(x=start, y=reads, col=state))	# outliers

+		ggplt <- ggplt + geom_point(data=dfplot.points, mapping=aes(x=start, y=reads, col=state), size=5, shape=21)	# outliers

 		ggplt <- ggplt + scale_color_manual(values=state.colors, drop=F)	# do not drop levels if not present

 		ggplt <- ggplt + empty_theme	# no axes whatsoever

 		ggplt <- ggplt + ylab(paste0(seqnames(grl[[i1]])[1], "\n", pstring, "\n", pstring2))	# chromosome names

+		ggplt <- ggplt + xlim(0,maxseqlength) + ylim(-0.6*custom.xlim,custom.xlim)	# set x- and y-limits

 		ggplt <- ggplt + coord_flip()

 		suppressWarnings(

 		print(ggplt, vp = viewport(layout.pos.row = matchidx$row, layout.pos.col = matchidx$col))

 		)

 	}

-	d <- dev.off()

+	if (!is.null(file)) {

+		d <- dev.off()

+	}

+}

+

+

+

+# ------------------------------------------------------------

+# Plot state categorization for all chromosomes

+# ------------------------------------------------------------

+plot.chromosomes.bivariate <- function(model, file=NULL) {

+	

+	## Convert to GRanges

+	gr <- hmm2GRanges(model, reduce=F)

+	grl <- split(gr, seqnames(gr))

+

+	## Get some variables

+	num.chroms <- length(levels(seqnames(gr)))

+	maxseqlength <- max(seqlengths(gr))

+	custom.xlim <- model$distributions.univariate[[1]]['monosomy','mu'] * 10

+

+	## Setup page

+	library(grid)

+	library(ggplot2)

+	nrows <- 2	# rows for plotting chromosomes

+	ncols <- ceiling(num.chroms/nrows)

+	if (!is.null(file)) {

+		pdf(file=paste0(file, '.pdf'), width=ncols*2.4, height=nrows*8.3)

+	}

+	grid.newpage()

+	layout <- matrix(1:((nrows+1)*ncols), ncol=ncols, nrow=nrows+1, byrow=T)

+	pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout), heights=c(1,21,21))))

+	# Main title

+	grid.text(model$ID, vp = viewport(layout.pos.row = 1, layout.pos.col = 1:ncols), gp=gpar(fontsize=26))

+

+	## Go through chromosomes and plot

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

+		# Get the i,j matrix positions of the regions that contain this subplot

+		matchidx <- as.data.frame(which(layout == i1+ncols, arr.ind = TRUE))

+

+		# Percentage of chromosome in state

+		tstate <- table(mcols(grl[[i1]])$state)

+		pstate.all <- tstate / sum(tstate)

+		pstate <- round(pstate.all*100)[-1]	# without 'nullsomy / unmapped' state

+		pstring <- apply(pstate, 1, function(x) { paste0(": ", x, "%") })

+		pstring <- paste0(names(pstring), pstring)

+		pstring <- paste(pstring[which.max(pstate)], collapse="\n")

+		pstring2 <- round(pstate.all*100)[1]	# only 'nullsomy / unmapped'

+		pstring2 <- paste0(names(pstring2), ": ", pstring2, "%")

+

+		## Convert to data.frame for plotting and prepare the data.frame

+			dfplot <- as.data.frame(grl[[i1]])

+			dfplot$reads.minus <- - dfplot$reads.minus

+		# Set values to great for plotting to limit

+			dfplot$reads.plus[dfplot$reads.plus>=custom.xlim] <- custom.xlim

+			dfplot$reads.minus[dfplot$reads.minus<=-custom.xlim] <- -custom.xlim

+			dfplot.points.plus <- dfplot[dfplot$reads.plus>=custom.xlim,]

+			dfplot.points.plus$reads <- rep(custom.xlim, nrow(dfplot.points.plus))

+			dfplot.points.minus <- dfplot[dfplot$reads.minus<=-custom.xlim,]

+			dfplot.points.minus$reads <- rep(-custom.xlim, nrow(dfplot.points.minus))

+		# No axes and labels

+		empty_theme <- theme(axis.line=element_blank(),

+      axis.text.x=element_blank(),

+      axis.text.y=element_blank(),

+      axis.ticks=element_blank(),

+			axis.title.x=element_text(size=20),

+      axis.title.y=element_blank(),

+      legend.position="none",

+      panel.background=element_blank(),

+      panel.border=element_blank(),

+      panel.grid.major=element_blank(),

+      panel.grid.minor=element_blank(),

+      plot.background=element_blank())

+		# Plot

+		ggplt <- ggplot(dfplot, aes(x=start, y=reads.plus))	# data

+		ggplt <- ggplt + geom_linerange(aes(ymin=0, ymax=reads.plus, col=state.separate.plus), size=0.2)	# read count

+		ggplt <- ggplt + geom_linerange(aes(ymin=0, ymax=reads.minus, col=state.separate.minus), size=0.2)	# read count

+		ggplt <- ggplt + geom_rect(ymin=-0.05*custom.xlim, ymax=0.05*custom.xlim, xmin=0, mapping=aes(xmax=max(start)), col='white', fill='gray20')	# chromosome backbone as simple rectangle

+		ggplt <- ggplt + geom_point(data=dfplot.points.plus, mapping=aes(x=start, y=reads, col=state.separate.plus), size=5, shape=21)	# outliers

+		ggplt <- ggplt + geom_point(data=dfplot.points.minus, mapping=aes(x=start, y=reads, col=state.separate.minus), size=5, shape=21)	# outliers

+		ggplt <- ggplt + scale_color_manual(values=state.colors, drop=F)	# do not drop levels if not present

+		ggplt <- ggplt + empty_theme	# no axes whatsoever

+		ggplt <- ggplt + ylab(paste0(seqnames(grl[[i1]])[1], "\n", pstring, "\n", pstring2))	# chromosome names

+		ggplt <- ggplt + xlim(0,maxseqlength) + ylim(-custom.xlim,custom.xlim)	# set x- and y-limits

+		ggplt <- ggplt + coord_flip()

+		suppressWarnings(

+		print(ggplt, vp = viewport(layout.pos.row = matchidx$row, layout.pos.col = matchidx$col))

+		)

+		

+	}

+	if (!is.null(file)) {

+		d <- dev.off()

+	}

 }

@@ -241,7 +351,7 @@ plot.genome.overview <- function(modellist, file='aneufinder_genome_overview', n

 	library(ggplot2)

 	nrows <- length(uni.hmm.grl)	# rows for plotting genomes

 	ncols <- 1

-	png(file=paste0(file, '.png'), width=ncols*60, height=nrows*5, units='cm', res=150)

+	pdf(file=paste0(file, '.pdf'), width=ncols*24, height=nrows*2)

 	grid.newpage()

 	layout <- matrix(1:((nrows+1)*ncols), ncol=ncols, nrow=nrows+1, byrow=T)

 	pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout), heights=c(1,rep(10,length(uni.hmm.grl))))))

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

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

+void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, 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

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

 		{

 			posterior_per_t[iN] = hmm->get_posterior(iN, t);

 		}

-		states[t] = argMax(posterior_per_t, *N);

+		states[t] = state_labels[argMax(posterior_per_t, *N)];

 	}

 	FILE_LOG(logDEBUG1) << "Return parameters";

@@ -210,3 +210,135 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, double

 }

 } // extern C

+

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

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

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

+extern "C" {

+void R_multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states, int* maxiter, int* maxtime, double* eps, int* states, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error)

+{

+

+	// Define logging level {"ERROR", "WARNING", "INFO", "ITERATION", "DEBUG", "DEBUG1", "DEBUG2", "DEBUG3", "DEBUG4"}

+// 	FILE* pFile = fopen("chromStar.log", "w");

+// 	Output2FILE::Stream() = pFile;

+//  	FILELog::ReportingLevel() = FILELog::FromString("ITERATION");

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

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

+

+	// Parallelization settings

+// 	omp_set_num_threads(*num_threads);

+

+	// Print some information

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

+	Rprintf("number of states = %d\n", *N);

+	FILE_LOG(logINFO) << "number of bins = " << *T;

+	Rprintf("number of bins = %d\n", *T);

+	if (*maxiter < 0)

+	{

+		FILE_LOG(logINFO) << "maximum number of iterations = none";

+		Rprintf("maximum number of iterations = none\n");

+	} else {

+		FILE_LOG(logINFO) << "maximum number of iterations = " << *maxiter;

+		Rprintf("maximum number of iterations = %d\n", *maxiter);

+	}

+	if (*maxtime < 0)

+	{

+		FILE_LOG(logINFO) << "maximum running time = none";

+		Rprintf("maximum running time = none\n");

+	} else {

+		FILE_LOG(logINFO) << "maximum running time = " << *maxtime << " sec";

+		Rprintf("maximum running time = %d sec\n", *maxtime);

+	}

+	FILE_LOG(logINFO) << "epsilon = " << *eps;

+	Rprintf("epsilon = %g\n", *eps);

+	FILE_LOG(logINFO) << "number of modifications = " << *Nmod;

+	Rprintf("number of modifications = %d\n", *Nmod);

+

+	// Flush Rprintf statements to console

+	R_FlushConsole();

+

+	// Recode the densities vector to matrix representation

+// 	clock_t clocktime = clock(), dtime;

+	double** multiD = allocDoubleMatrix(*N, *T);

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

+	{

+		for (int t=0; t<*T; t++)

+		{

+			multiD[iN][t] = D[iN*(*T)+t];

+		}

+	}

+// 	dtime = clock() - clocktime;

+// 	FILE_LOG(logDEBUG1) << "recoding densities vector to matrix representation: " << dtime << " clicks";

+

+	// Create the HMM

+	FILE_LOG(logDEBUG1) << "Creating the multivariate HMM";

+	ScaleHMM* hmm = new ScaleHMM(*T, *N, *Nmod, multiD);

+	// Initialize the transition probabilities and proba

+	hmm->initialize_transition_probs(initial_A, *use_initial_params);

+	hmm->initialize_proba(initial_proba, *use_initial_params);

+	

+	// Print logproba and A

+// 	for (int iN=0; iN<*N; iN++)

+// 	{

+// 		FILE_LOG(logDEBUG) << "proba["<<iN<<"] = " <<exp(hmm->logproba[iN]);

+// 		for (int jN=0; jN<*N; jN++)

+// 		{

+// 			FILE_LOG(logDEBUG) << "A["<<iN<<"]["<<jN<<"] = " << hmm->A[iN][jN];

+// 		}

+// 	}

+

+	// Estimate the parameters

+	FILE_LOG(logDEBUG1) << "Starting Baum-Welch estimation";

+	try

+	{

+		hmm->baumWelch(maxiter, maxtime, eps);

+	}

+	catch (std::exception& e)

+	{

+		FILE_LOG(logERROR) << "Error in Baum-Welch: " << e.what();

+		Rprintf("Error in Baum-Welch: %s\n", e.what());

+		if (e.what()=="nan detected") { *error = 1; }

+		else { *error = 2; }

+	}

+	FILE_LOG(logDEBUG1) << "Finished with Baum-Welch estimation";

+	

+// 	// Compute the posteriors and save results directly to the R pointer

+// 	FILE_LOG(logDEBUG1) << "Recode posteriors into column representation";

+// 	for (int iN=0; iN<*N; iN++)

+// 	{

+// 		for (int t=0; t<*T; t++)

+// 		{

+// 			posteriors[t + iN * (*T)] = hmm->get_posterior(iN, t);

+// 		}

+// 	}

+

+	// Compute the states from posteriors

+	FILE_LOG(logDEBUG1) << "Computing states from posteriors";

+	double posterior_per_t [*N];

+	for (int t=0; t<*T; t++)

+	{

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

+		{

+			posterior_per_t[iN] = hmm->get_posterior(iN, t);

+		}

+		states[t] = comb_states[argMax(posterior_per_t, *N)];

+	}

+	

+	FILE_LOG(logDEBUG1) << "Return parameters";

+	// also return the estimated transition matrix and the initial probs

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

+	{

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

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

+		{

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

+		}

+	}

+	*loglik = hmm->get_logP();

+

+	FILE_LOG(logDEBUG1) << "Deleting the hmm";

+	delete hmm;

+	freeDoubleMatrix(multiD, *N);

+}

+} // extern C

+

@@ -1382,4 +1382,40 @@ double Geometric::get_prob()

 }

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

+// Multivariate Copula Approximation

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

+

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

+MVCopulaApproximation::MVCopulaApproximation(int** multiobservations, int T, std::vector<Density*> marginals, double* cor_matrix_inv, double cor_matrix_determinant)

+{

+	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	this->multi_obs = multiobservations;

+	this->T = T;

+	// these are the marginal distributions (we need their CDF function)

+	this->marginals = marginals;

+	this->Nmod = this->marginals.size();

+	this->cor_matrix_inv = cor_matrix_inv;

+	this->cor_matrix_determinant = cor_matrix_determinant;

+}

+

+MVCopulaApproximation::~MVCopulaApproximation()

+{

+	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	for (int imod; imod<this->Nmod; imod++)

+	{

+		delete this->marginals[imod];

+	}

+}

+

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

+

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

+DensityName MVCopulaApproximation::get_name()

+{

+	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	return(OTHER);

+}

+

+	

@@ -1,10 +1,12 @@

 #include <Rmath.h> // dnorm(), dnbinom() and digamma() etc.

+#include <vector> // storing density functions in MVCopula

 #include "utility.h" // FILE_LOG(), intMax()

 #ifndef DENSITIES_H

 #define DENSITIES_H