@@ -7,8 +7,6 @@ S3method(plot,character)

 export(Aneufinder)

 export(bam2binned)

 export(bed2GRanges)

-export(bed2binned)

-export(bedGraph2binned)

 export(clusterByQuality)

 export(concGRangesInRange)

 export(correctGC)

@@ -9,29 +9,27 @@

 #' @param configfile A file specifying the parameters of this function (without \code{inputfolder}, \code{outputfolder} and \code{configfile}). Having the parameters in a file can be handy if many samples with the same parameter settings are to be run. If a \code{configfile} is specified, it will take priority over the command line parameters.

 #' @param numCPU The numbers of CPUs that are used. Should not be more than available on your machine.

 #' @param reuse.existing.files A logical indicating whether or not existing files in \code{outputfolder} should be reused.

-#' @param binsizes An integer vector with bin sizes. If more than one value is given, output files will be produced for each bin size.

-#' @param reads.per.bin Approximate number of desired reads per bin. The bin size will be selected accordingly. Output files are produced for each value.

-#' @param format Format of the sequencing data. One of \code{c('bam')}.

-#' @param chromosomes If only a subset of the chromosomes should be binned, specify them here.

-#' @param remove.duplicate.reads A logical indicating whether or not duplicate reads should be removed.

-#' @param min.mapq Minimum mapping quality when importing from BAM files.

+

+#' @param stepsize Fraction of the binsize that the sliding window is offset at each step. Example: If \code{stepsize=0.1} and \code{binsizes=c(200000,500000)}, the actual stepsize in basepairs is 20000 and 50000, respectively.

+#' @inheritParams align2binned

+

 #' @param correction.method Correction methods to be used for the binned read counts. Currently any combination of \code{c('GC')}.

 #' @param GC.BSgenome A \code{BSgenome} object which contains the DNA sequence that is used for the GC correction.

 #' @param callCNVs A logical indicating whether CNVs should be called.

 #' @param callSCEs A logical indicating whether SCEs should be called.

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

-#' @param max.time The maximum running time in seconds for the Baum-Welch algorithm. If this time is reached, the Baum-Welch will terminate after the current iteration finishes. \code{max.time=-1} means no limit.

-#' @param max.iter The maximum number of iterations for the Baum-Welch algorithm. \code{max.iter=-1} means no limit.

-#' @param num.trials The number of trials to find a fit where state \code{most.frequent.state} is most frequent. Each time, the HMM is seeded with different random initial values.

-#' @param states A subset or all of \code{c("zero-inflation","nullsomy","monosomy","disomy","trisomy","tetrasomy","multisomy")}. This vector defines the states that are used in the Hidden Markov Model. The order of the entries should not be changed.

+

+#' @inheritParams univariate.findCNVs

 #' @param most.frequent.state One of the states that were given in \code{states}. The specified state is assumed to be the most frequent one when calling CNVs. This can help the fitting procedure to converge into the correct fit. Default is 'disomy'.

 #' @param most.frequent.state.SCE One of the states that were given in \code{states}. The specified state is assumed to be the most frequent one when calling SCEs. This can help the fitting procedure to converge into the correct fit. Default is 'monosomy'.

+

+#' @inheritParams getSCEcoordinates

+

 #' @param cluster.plots A logical indicating whether plots should be clustered by similarity.

 #' @author Aaron Taudt

 #' @import foreach

 #' @import doParallel

 #' @export

-Aneufinder <- function(inputfolder, outputfolder, configfile=NULL, numCPU=1, reuse.existing.files=TRUE, binsizes=500000, reads.per.bin=NULL, format='bam', chromosomes=NULL, remove.duplicate.reads=TRUE, min.mapq=10, correction.method=NULL, GC.BSgenome=NULL, callCNVs=TRUE, callSCEs=FALSE, eps=0.1, max.time=60, max.iter=5000, num.trials=15, states=c('zero-inflation','nullsomy','monosomy','disomy','trisomy','tetrasomy','multisomy'), most.frequent.state='disomy', most.frequent.state.SCE='monosomy', cluster.plots=TRUE) {

+Aneufinder <- function(inputfolder, outputfolder, configfile=NULL, numCPU=1, reuse.existing.files=TRUE, binsizes=500000, reads.per.bin=NULL, stepsize=NULL, format='bam', chromosomes=NULL, remove.duplicate.reads=TRUE, min.mapq=10, correction.method=NULL, GC.BSgenome=NULL, callCNVs=TRUE, callSCEs=FALSE, eps=0.1, max.time=60, max.iter=5000, num.trials=15, states=c('zero-inflation','nullsomy','monosomy','disomy','trisomy','tetrasomy','multisomy'), most.frequent.state='disomy', most.frequent.state.SCE='monosomy', resolution=c(3,6), min.segwidth=2, min.reads=50, cluster.plots=TRUE) {

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

 ### Helper functions ###

@@ -63,7 +61,7 @@ if (class(GC.BSgenome)=='BSgenome') {

 }

 ## Put options into list and merge with conf

-params <- list(numCPU=numCPU, reuse.existing.files=reuse.existing.files, binsizes=binsizes, reads.per.bin=reads.per.bin, format=format, chromosomes=chromosomes, remove.duplicate.reads=remove.duplicate.reads, min.mapq=min.mapq, correction.method=correction.method, GC.BSgenome=GC.BSgenome, callCNVs=callCNVs, callSCEs=callSCEs, eps=eps, max.time=max.time, max.iter=max.iter, num.trials=num.trials, states=states, most.frequent.state=most.frequent.state, most.frequent.state.SCE=most.frequent.state.SCE, cluster.plots=cluster.plots)

+params <- list(numCPU=numCPU, reuse.existing.files=reuse.existing.files, binsizes=binsizes, reads.per.bin=reads.per.bin, stepsize=stepsize, format=format, chromosomes=chromosomes, remove.duplicate.reads=remove.duplicate.reads, min.mapq=min.mapq, correction.method=correction.method, GC.BSgenome=GC.BSgenome, callCNVs=callCNVs, callSCEs=callSCEs, eps=eps, max.time=max.time, max.iter=max.iter, num.trials=num.trials, states=states, most.frequent.state=most.frequent.state, most.frequent.state.SCE=most.frequent.state.SCE, resolution=resolution, min.segwidth=min.segwidth, min.reads=min.reads, cluster.plots=cluster.plots)

 conf <- c(conf, params[setdiff(names(params),names(conf))])

 ## Input checks

@@ -74,7 +72,7 @@ if (! conf[['format']] %in% c('bam')) {

 ## Helpers

 binsizes <- conf[['binsizes']]

 reads.per.bins <- conf[['reads.per.bin']]

-patterns <- c(paste0('reads.per.bin_',reads.per.bins,'_'), paste0('binsize_',format(binsizes, scientific=F, trim=T),'_'))

+patterns <- c(paste0('reads.per.bin_',reads.per.bins,'_'), paste0('binsize_',format(binsizes, scientific=T, trim=T),'_'))

 patterns <- setdiff(patterns, c('reads.per.bin__','binsize__'))

 pattern <- NULL #ease R CMD check

@@ -123,7 +121,7 @@ temp <- foreach (file = files, .packages=c('aneufinder')) %dopar% {

 	if (length(c(binsizes.todo,rpbin.todo)) > 0) {

 		tC <- tryCatch({

 			if (conf[['format']]=='bam') {

-				bam2binned(bamfile=file, binsizes=binsizes.todo, reads.per.bin=rpbin.todo, chromosomes=conf[['chromosomes']], remove.duplicate.reads=conf[['remove.duplicate.reads']], min.mapq=conf[['min.mapq']], outputfolder.binned=binpath.uncorrected, save.as.RData=TRUE, reads.store=TRUE, outputfolder.reads=readspath)

+				bam2binned(bamfile=file, binsizes=binsizes.todo, reads.per.bin=rpbin.todo, stepsize=conf[['stepsize']], chromosomes=conf[['chromosomes']], remove.duplicate.reads=conf[['remove.duplicate.reads']], min.mapq=conf[['min.mapq']], outputfolder.binned=binpath.uncorrected, save.as.RData=TRUE, reads.store=TRUE, outputfolder.reads=readspath)

 			}

 		}, error = function(err) {

 			stop(file,'\n',err)

@@ -217,7 +215,7 @@ if (conf[['callCNVs']]==TRUE) {

 	### Plotting CNV ###

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

 	if (!file.exists(CNVplotpath)) { dir.create(CNVplotpath) }

-	patterns <- c(paste0('reads.per.bin_',reads.per.bins,'_'), paste0('binsize_',format(binsizes, scientific=F, trim=T),'_'))

+	patterns <- c(paste0('reads.per.bin_',reads.per.bins,'_'), paste0('binsize_',format(binsizes, scientific=T, trim=T),'_'))

 	patterns <- setdiff(patterns, c('reads.per.bin__','binsize__'))

 	files <- list.files(CNVpath, full.names=TRUE, recursive=T, pattern='.RData$')

@@ -230,7 +228,7 @@ if (conf[['callCNVs']]==TRUE) {

 		if (!file.exists(savename)) {

 			pdf(file=savename, width=10, height=7)

 			ifiles <- list.files(CNVpath, pattern='RData$', full.names=TRUE)

-			ifiles <- grep(pattern, ifiles, value=T)

+			ifiles <- grep(gsub('\\+','\\\\+',pattern), ifiles, value=T)

 			for (ifile in ifiles) {

 				tC <- tryCatch({

 					model <- get(load(ifile))

@@ -249,7 +247,7 @@ if (conf[['callCNVs']]==TRUE) {

 	message("Plotting heatmaps ...", appendLF=F); ptm <- proc.time()

 	temp <- foreach (pattern = patterns, .packages=c('aneufinder')) %dopar% {

 		ifiles <- list.files(CNVpath, pattern='RData$', full.names=TRUE)

-		ifiles <- grep(pattern, ifiles, value=T)

+		ifiles <- grep(gsub('\\+','\\\\+',pattern), ifiles, value=T)

 		savename=file.path(CNVplotpath,paste0('genomeHeatmap_',sub('_$','',pattern),'.pdf'))

 		if (!file.exists(savename)) {

 			suppressMessages(heatmapGenomewide(ifiles, file=savename, plot.SCE=F, cluster=conf[['cluster.plots']]))

@@ -257,7 +255,7 @@ if (conf[['callCNVs']]==TRUE) {

 	}

 	temp <- foreach (pattern = patterns, .packages=c('aneufinder')) %dopar% {

 		ifiles <- list.files(CNVpath, pattern='RData$', full.names=TRUE)

-		ifiles <- grep(pattern, ifiles, value=T)

+		ifiles <- grep(gsub('\\+','\\\\+',pattern), ifiles, value=T)

 		savename=file.path(CNVplotpath,paste0('aneuploidyHeatmap_',sub('_$','',pattern),'.pdf'))

 		if (!file.exists(savename)) {

 			ggplt <- suppressMessages(heatmapAneuploidies(ifiles, cluster=conf[['cluster.plots']]))

@@ -276,7 +274,7 @@ if (conf[['callCNVs']]==TRUE) {

 		if (!file.exists(savename)) {

 			pdf(file=savename, width=20, height=5)

 			ifiles <- list.files(CNVpath, pattern='RData$', full.names=TRUE)

-			ifiles <- grep(pattern, ifiles, value=T)

+			ifiles <- grep(gsub('\\+','\\\\+',pattern), ifiles, value=T)

 			for (ifile in ifiles) {

 				tC <- tryCatch({

 					model <- get(load(ifile))

@@ -298,7 +296,7 @@ if (conf[['callCNVs']]==TRUE) {

 		if (!file.exists(savename)) {

 			pdf(file=savename, width=12*1.4, height=2*4.6)

 			ifiles <- list.files(CNVpath, pattern='RData$', full.names=TRUE)

-			ifiles <- grep(pattern, ifiles, value=T)

+			ifiles <- grep(gsub('\\+','\\\\+',pattern), ifiles, value=T)

 			for (ifile in ifiles) {

 				tC <- tryCatch({

 					model <- get(load(ifile))

@@ -320,7 +318,7 @@ if (conf[['callCNVs']]==TRUE) {

 		savename <- file.path(CNVbrowserpath,paste0(pattern))

 		if (!file.exists(paste0(savename,'.bed.gz'))) {

 			ifiles <- list.files(CNVpath, pattern='RData$', full.names=TRUE)

-			ifiles <- grep(pattern, ifiles, value=T)

+			ifiles <- grep(gsub('\\+','\\\\+',pattern), ifiles, value=T)

 			exportCNVs(ifiles, filename=savename, cluster=conf[['cluster.plots']], export.CNV=TRUE, export.SCE=FALSE)

 		}

 	}

@@ -340,7 +338,9 @@ if (conf[['callSCEs']]==TRUE) {

 		tC <- tryCatch({

 			savename <- file.path(SCEpath,basename(file))

 			if (!file.exists(savename)) {

-				model <- findSCEs(file, eps=conf[['eps']], max.time=conf[['max.time']], max.iter=conf[['max.iter']], num.trials=conf[['num.trials']], states=conf[['states']], most.frequent.state=conf[['most.frequent.state.SCE']]) 

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

+				reads.file <- list.files(readspath, full.names=T, pattern=paste0('^',attr(binned.data,'ID'),'.RData$'))

+				model <- findSCEs(file, eps=conf[['eps']], max.time=conf[['max.time']], max.iter=conf[['max.iter']], num.trials=conf[['num.trials']], states=conf[['states']], most.frequent.state=conf[['most.frequent.state.SCE']], resolution=conf[['resolution']], min.segwidth=conf[['min.segwidth']], fragments=reads.file, min.reads=conf[['min.reads']]) 

 				save(model, file=savename)

 			}

 		}, error = function(err) {

@@ -353,7 +353,7 @@ if (conf[['callSCEs']]==TRUE) {

 	### Plotting SCE ###

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

 	if (!file.exists(SCEplotpath)) { dir.create(SCEplotpath) }

-	patterns <- c(paste0('reads.per.bin_',reads.per.bins,'_'), paste0('binsize_',format(binsizes, scientific=F, trim=T),'_'))

+	patterns <- c(paste0('reads.per.bin_',reads.per.bins,'_'), paste0('binsize_',format(binsizes, scientific=T, trim=T),'_'))

 	patterns <- setdiff(patterns, c('reads.per.bin__','binsize__'))

 	files <- list.files(SCEpath, full.names=TRUE, recursive=T, pattern='.RData$')

@@ -366,7 +366,7 @@ if (conf[['callSCEs']]==TRUE) {

 		if (!file.exists(savename)) {

 			pdf(file=savename, width=10, height=7)

 			ifiles <- list.files(SCEpath, pattern='RData$', full.names=TRUE)

-			ifiles <- grep(pattern, ifiles, value=T)

+			ifiles <- grep(gsub('\\+','\\\\+',pattern), ifiles, value=T)

 			for (ifile in ifiles) {

 				tC <- tryCatch({

 					model <- get(load(ifile))

@@ -385,7 +385,7 @@ if (conf[['callSCEs']]==TRUE) {

 	message("Plotting heatmaps ...", appendLF=F); ptm <- proc.time()

 	temp <- foreach (pattern = patterns, .packages=c('aneufinder')) %dopar% {

 		ifiles <- list.files(SCEpath, pattern='RData$', full.names=TRUE)

-		ifiles <- grep(pattern, ifiles, value=T)

+		ifiles <- grep(gsub('\\+','\\\\+',pattern), ifiles, value=T)

 		savename=file.path(SCEplotpath,paste0('genomeHeatmap_',sub('_$','',pattern),'.pdf'))

 		if (!file.exists(savename)) {

 			suppressMessages(heatmapGenomewide(ifiles, file=savename, plot.SCE=T, cluster=conf[['cluster.plots']]))

@@ -393,7 +393,7 @@ if (conf[['callSCEs']]==TRUE) {

 	}

 	temp <- foreach (pattern = patterns, .packages=c('aneufinder')) %dopar% {

 		ifiles <- list.files(SCEpath, pattern='RData$', full.names=TRUE)

-		ifiles <- grep(pattern, ifiles, value=T)

+		ifiles <- grep(gsub('\\+','\\\\+',pattern), ifiles, value=T)

 		savename=file.path(SCEplotpath,paste0('aneuploidyHeatmap_',sub('_$','',pattern),'.pdf'))

 		if (!file.exists(savename)) {

 			pdf(savename, width=30, height=0.3*length(ifiles))

@@ -412,7 +412,7 @@ if (conf[['callSCEs']]==TRUE) {

 		if (!file.exists(savename)) {

 			pdf(file=savename, width=20, height=5)

 			ifiles <- list.files(SCEpath, pattern='RData$', full.names=TRUE)

-			ifiles <- grep(pattern, ifiles, value=T)

+			ifiles <- grep(gsub('\\+','\\\\+',pattern), ifiles, value=T)

 			for (ifile in ifiles) {

 				tC <- tryCatch({

 					model <- get(load(ifile))

@@ -434,7 +434,7 @@ if (conf[['callSCEs']]==TRUE) {

 		if (!file.exists(savename)) {

 			pdf(file=savename, width=12*1.4, height=2*4.6)

 			ifiles <- list.files(SCEpath, pattern='RData$', full.names=TRUE)

-			ifiles <- grep(pattern, ifiles, value=T)

+			ifiles <- grep(gsub('\\+','\\\\+',pattern), ifiles, value=T)

 			for (ifile in ifiles) {

 				tC <- tryCatch({

 					model <- get(load(ifile))

@@ -456,7 +456,7 @@ if (conf[['callSCEs']]==TRUE) {

 		savename <- file.path(SCEbrowserpath,sub('_$','',pattern))

 		if (!file.exists(paste0(savename,'.bed.gz'))) {

 			ifiles <- list.files(SCEpath, pattern='RData$', full.names=TRUE)

-			ifiles <- grep(pattern, ifiles, value=T)

+			ifiles <- grep(gsub('\\+','\\\\+',pattern), ifiles, value=T)

 			exportCNVs(ifiles, filename=savename, cluster=conf[['cluster.plots']], export.CNV=TRUE, export.SCE=TRUE)

 		}

 	}

@@ -2,10 +2,9 @@

 #' Convert aligned reads from various file formats into read counts in equidistant bins

 #'

-#' Convert aligned reads in .bam or .bed format into read counts in equidistant windows. Convert signal values in .bedGraph format to signal counts in equidistant windows.

+#' Convert aligned reads in .bam or .bed format into read counts in equidistant windows.

 #'

 #' Convert aligned reads or signal values from various file formats into read counts in equidistant windows (bins). bam2binned() and bed2binned() use 'GenomicRanges::countOverlaps()' to calculate the read counts. Therefore, with small binsizes and large read lengths, one read fragment will often overlap more than one bin.

-#' bedGraph2binned() sets the maximum signal value in a bin as value for that bin.

 #'

 #' @name binning

 #' @examples

@@ -21,13 +20,13 @@ NULL

 #' @param bamfile A file in BAM format.

 #' @param bamindex BAM index file. Can be specified without the .bai ending. If the index file does not exist it will be created and a warning is issued.

 #' @export

-bam2binned <- function(bamfile, bamindex=bamfile, ID=basename(bamfile), pairedEndReads=FALSE, outputfolder.binned="binned_data", binsizes=NULL, reads.per.bin=NULL, numbins=NULL, chromosomes=NULL, save.as.RData=FALSE, calc.complexity=TRUE, min.mapq=10, remove.duplicate.reads=TRUE, reads.store=FALSE, outputfolder.reads="data", reads.return=FALSE, reads.overwrite=FALSE, reads.only=FALSE) {

+bam2binned <- function(bamfile, bamindex=bamfile, ID=basename(bamfile), pairedEndReads=FALSE, outputfolder.binned="binned_data", binsizes=NULL, reads.per.bin=NULL, stepsize=NULL, chromosomes=NULL, save.as.RData=FALSE, calc.complexity=TRUE, min.mapq=10, remove.duplicate.reads=TRUE, reads.store=FALSE, outputfolder.reads="data", reads.return=FALSE, reads.overwrite=FALSE, reads.only=FALSE) {

 	call <- match.call()

 	underline <- paste0(rep('=',sum(nchar(call[[1]]))+3), collapse='')

 	message("\n",call[[1]],"():")

 	message(underline)

 	ptm <- proc.time()

-	binned.data <- align2binned(bamfile, format="bam", ID=ID, bamindex=bamindex, pairedEndReads=pairedEndReads, outputfolder.binned=outputfolder.binned, binsizes=binsizes, reads.per.bin=reads.per.bin, numbins=numbins, chromosomes=chromosomes, save.as.RData=save.as.RData, calc.complexity=calc.complexity, min.mapq=min.mapq, remove.duplicate.reads=remove.duplicate.reads, call=call, reads.store=reads.store, outputfolder.reads=outputfolder.reads, reads.return=reads.return, reads.overwrite=reads.overwrite, reads.only=reads.only)

+	binned.data <- align2binned(bamfile, format="bam", ID=ID, bamindex=bamindex, pairedEndReads=pairedEndReads, outputfolder.binned=outputfolder.binned, binsizes=binsizes, reads.per.bin=reads.per.bin, stepsize=stepsize, chromosomes=chromosomes, save.as.RData=save.as.RData, calc.complexity=calc.complexity, min.mapq=min.mapq, remove.duplicate.reads=remove.duplicate.reads, call=call, reads.store=reads.store, outputfolder.reads=outputfolder.reads, reads.return=reads.return, reads.overwrite=reads.overwrite, reads.only=reads.only)

 	time <- proc.time() - ptm

 	message("Time spent in ", call[[1]],"(): ",round(time[3],2),"s")

 	return(binned.data)

@@ -36,30 +35,13 @@ bam2binned <- function(bamfile, bamindex=bamfile, ID=basename(bamfile), pairedEn

 #' @describeIn binning Bin reads in BED format

 #' @inheritParams align2binned

 #' @param bedfile A file in BED format.

-#' @export

-bed2binned <- function(bedfile, chrom.length.file, ID=basename(bedfile), outputfolder.binned="binned_data", binsizes=NULL, reads.per.bin=NULL, numbins=NULL, chromosomes=NULL, save.as.RData=FALSE, calc.complexity=TRUE, min.mapq=10, remove.duplicate.reads=TRUE, reads.store=FALSE, outputfolder.reads="data", reads.return=FALSE, reads.overwrite=FALSE, reads.only=FALSE) {

+bed2binned <- function(bedfile, chrom.length.file, ID=basename(bedfile), outputfolder.binned="binned_data", binsizes=NULL, reads.per.bin=NULL, stepsize=NULL, chromosomes=NULL, save.as.RData=FALSE, calc.complexity=TRUE, min.mapq=10, remove.duplicate.reads=TRUE, reads.store=FALSE, outputfolder.reads="data", reads.return=FALSE, reads.overwrite=FALSE, reads.only=FALSE) {

 	call <- match.call()

 	underline <- paste0(rep('=',sum(nchar(call[[1]]))+3), collapse='')

 	message("\n",call[[1]],"():")

 	message(underline)

 	ptm <- proc.time()

-	binned.data <- align2binned(bedfile, format="bed", ID=ID, chrom.length.file=chrom.length.file, outputfolder.binned=outputfolder.binned, binsizes=binsizes, reads.per.bin=reads.per.bin, numbins=numbins, chromosomes=chromosomes, save.as.RData=save.as.RData, calc.complexity=calc.complexity, min.mapq=min.mapq, remove.duplicate.reads=remove.duplicate.reads, call=call, reads.store=reads.store, outputfolder.reads=outputfolder.reads, reads.return=reads.return, reads.overwrite=reads.overwrite, reads.only=reads.only)

-	time <- proc.time() - ptm

-	message("Time spent in ", call[[1]],"(): ",round(time[3],2),"s")

-	return(binned.data)

-}

-

-#' @describeIn binning Bin reads in bedGraph format

-#' @inheritParams align2binned

-#' @param bedGraphfile A file in bedGraph format.

-#' @export

-bedGraph2binned <- function(bedGraphfile, chrom.length.file, ID=basename(bedGraphfile), outputfolder.binned="binned_data", binsizes=NULL, reads.per.bin=NULL, numbins=NULL, chromosomes=NULL, save.as.RData=FALSE, calc.complexity=TRUE, min.mapq=10, remove.duplicate.reads=TRUE, reads.store=FALSE, outputfolder.reads="data", reads.return=FALSE, reads.overwrite=FALSE, reads.only=FALSE) {

-	call <- match.call()

-	underline <- paste0(rep('=',sum(nchar(call[[1]]))+3), collapse='')

-	message("\n",call[[1]],"():")

-	message(underline)

-	ptm <- proc.time()

-	binned.data <- align2binned(bedGraphfile, format="bedGraph", ID=ID, chrom.length.file=chrom.length.file, outputfolder.binned=outputfolder.binned, binsizes=binsizes, reads.per.bin=reads.per.bin, numbins=numbins, chromosomes=chromosomes, save.as.RData=save.as.RData, calc.complexity=calc.complexity, min.mapq=min.mapq, remove.duplicate.reads=remove.duplicate.reads, call=call, reads.store=reads.store, outputfolder.reads=outputfolder.reads, reads.return=reads.return, reads.overwrite=reads.overwrite, reads.only=reads.only)

+	binned.data <- align2binned(bedfile, format="bed", ID=ID, chrom.length.file=chrom.length.file, outputfolder.binned=outputfolder.binned, binsizes=binsizes, reads.per.bin=reads.per.bin, stepsize=stepsize, chromosomes=chromosomes, save.as.RData=save.as.RData, calc.complexity=calc.complexity, min.mapq=min.mapq, remove.duplicate.reads=remove.duplicate.reads, call=call, reads.store=reads.store, outputfolder.reads=outputfolder.reads, reads.return=reads.return, reads.overwrite=reads.overwrite, reads.only=reads.only)

 	time <- proc.time() - ptm

 	message("Time spent in ", call[[1]],"(): ",round(time[3],2),"s")

 	return(binned.data)

@@ -67,10 +49,10 @@ bedGraph2binned <- function(bedGraphfile, chrom.length.file, ID=basename(bedGrap

 #' Convert aligned reads from various file formats into read counts in equidistant bins

 #'

-#' Convert aligned reads in .bam or .bed format into read counts in equidistant windows. Convert signal values in .bedGraph format to signal counts in equidistant windows.

+#' Convert aligned reads in .bam or .bed format into read counts in equidistant windows.

 #'

 #' @param file A file with aligned reads.

-#' @param format One of \code{c('bam', 'bed', 'bedGraph')}.

+#' @param format One of \code{c('bam', 'bed')}.

 #' @param ID An identifier that will be used to identify the file throughout the workflow and in plotting.

 #' @param bamindex Index file if \code{format='bam'} with or without the .bai ending. If this file does not exist it will be created and a warning is issued.

 #' @param pairedEndReads Set to \code{TRUE} if you have paired-end reads in your file.

@@ -78,7 +60,7 @@ bedGraph2binned <- function(bedGraphfile, chrom.length.file, ID=basename(bedGrap

 #' @param outputfolder.binned Folder to which the binned data will be saved. If the specified folder does not exist, it will be created.

 #' @param binsizes An integer vector with bin sizes. If more than one value is given, output files will be produced for each bin size.

 #' @param reads.per.bin Approximate number of desired reads per bin. The bin size will be selected accordingly. Output files are produced for each value.

-#' @param numbins Number of bins per chromosome. Each chromosome will have a different binsize! DO NOT USE UNLESS YOU KNOW WHAT YOU ARE DOING. Output files are produced for each value.

+#' @param stepsize Fraction of the binsize that the sliding window is offset at each step. Example: If \code{stepsize=0.1} and \code{binsizes=c(200000,500000)}, the actual stepsize in basepairs is 20000 and 50000, respectively.

 #' @param chromosomes If only a subset of the chromosomes should be binned, specify them here.

 #' @param save.as.RData If set to \code{FALSE}, no output file will be written. Instead, a \link{GenomicRanges} object containing the binned data will be returned. Only the first binsize will be processed in this case.

 #' @param calc.complexity A logical indicating whether or not to estimate library complexity.

@@ -91,7 +73,7 @@ bedGraph2binned <- function(bedGraphfile, chrom.length.file, ID=basename(bedGrap

 #' @param reads.overwrite Whether or not an existing file with read fragments should be overwritten.

 #' @param reads.only If \code{TRUE} only read fragments are stored and/or returned and no binning is done.

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

-align2binned <- function(file, format, ID=basename(file), bamindex=file, pairedEndReads=FALSE, chrom.length.file, outputfolder.binned="binned_data", binsizes=200000, reads.per.bin=NULL, numbins=NULL, chromosomes=NULL, save.as.RData=FALSE, calc.complexity=TRUE, min.mapq=10, remove.duplicate.reads=TRUE, call=match.call(), reads.store=FALSE, outputfolder.reads="data", reads.return=FALSE, reads.overwrite=FALSE, reads.only=FALSE) {

+align2binned <- function(file, format, ID=basename(file), bamindex=file, pairedEndReads=FALSE, chrom.length.file, outputfolder.binned="binned_data", binsizes=200000, reads.per.bin=NULL, stepsize=NULL, chromosomes=NULL, save.as.RData=FALSE, calc.complexity=TRUE, min.mapq=10, remove.duplicate.reads=TRUE, call=match.call(), reads.store=FALSE, outputfolder.reads="data", reads.return=FALSE, reads.overwrite=FALSE, reads.only=FALSE) {

 	## Check user input

 	if (reads.return==FALSE & reads.only==FALSE) {

@@ -137,21 +119,6 @@ align2binned <- function(file, format, ID=basename(file), bamindex=file, pairedE

 			}

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

 		}

-	## BEDGraph (0-based)

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

-		# File with chromosome lengths (1-based)

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

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

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

-		# File with reads, determine classes first for faster import

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

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

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

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

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

-		# Convert to GRanges object

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

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

 	}

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

@@ -208,60 +175,44 @@ align2binned <- function(file, format, ID=basename(file), bamindex=file, pairedE

 	coverage <- list(coverage=coverage, genome.covered=genome.covered, coverage.per.chrom=coverage.per.chrom, genome.covered.per.chrom=genome.covered.per.chrom)

-	### Loop over all binsizes ###

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

 	numcountsperbp <- length(data) / sum(as.numeric(seqlengths(data)))

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

 	reads.per.bin.add <- round(binsizes * numcountsperbp, 2)

 	binsizes <- c(binsizes, binsizes.add)

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

-	if (is.null(numbins)) {

-		length.binsizes <- length(binsizes)

-	} else {

-		length.binsizes <- length(numbins)

-	}

+	length.binsizes <- length(binsizes)

+	### Loop over all binsizes ###

+	data.plus <- data[strand(data)=='+']

+	data.minus <- data[strand(data)=='-']

 	for (ibinsize in 1:length.binsizes) {

-		if (is.null(numbins)) {

-			binsize <- binsizes[ibinsize]

-			readsperbin <- reads.per.bin[ibinsize]

-			message("Binning into bin size ",binsize," with on average ",readsperbin," reads per bin")

-		} else {

-			numbin <- numbins[ibinsize]

-			message("Binning with number of bins ",numbin)

-		}

+		binsize <- binsizes[ibinsize]

+		readsperbin <- reads.per.bin[ibinsize]

+		message("Binning into bin size ",binsize," with on average ",readsperbin," reads per bin")

-		### Iterate over all chromosomes

+		### Loop over chromosomes ###

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

 		binned.data <- GenomicRanges::GRangesList()

 		for (chromosome in chroms2use) {

-			if (is.null(numbins)) {

-				## Check last incomplete bin

-				incomplete.bin <- seqlengths(data)[chromosome] %% binsize > 0

-				if (incomplete.bin) {

-					numbin <- floor(seqlengths(data)[chromosome]/binsize)	# floor: we don't want incomplete bins, ceiling: we want incomplete bins at the end

-				} else {

-					numbin <- seqlengths(data)[chromosome]/binsize

-				}

-				if (numbin == 0) {

-					warning("chromosome ",chromosome," is smaller than binsize. Skipped.")

-					next

-				}

-				## Initialize vectors

-				chroms <- rep(chromosome,numbin)

-				reads <- rep(0,numbin)

-				start <- seq(from=1, by=binsize, length.out=numbin)

-				end <- seq(from=binsize, by=binsize, length.out=numbin)

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

+			## Check last incomplete bin

+			incomplete.bin <- seqlengths(data)[chromosome] %% binsize > 0

+			if (incomplete.bin) {

+				numbin <- floor(seqlengths(data)[chromosome]/binsize)	# floor: we don't want incomplete bins, ceiling: we want incomplete bins at the end

 			} else {

-				## Initialize vectors

-				chroms <- rep(chromosome,numbin)

-				start <- round(seq(from=1, to=seqlengths(data)[chromosome], length.out=numbin+1))

-				end <- start[-1] - 1

-				end[length(end)] <- end[length(end)] + 1

-				start <- start[-length(start)]

+				numbin <- seqlengths(data)[chromosome]/binsize

 			}

+			if (numbin == 0) {

+				warning("chromosome ",chromosome," is smaller than binsize. Skipped.")

+				next

+			}

+			## Initialize vectors

+			chroms <- rep(chromosome,numbin)

+			reads <- rep(0,numbin)

+			start <- seq(from=1, by=binsize, length.out=numbin)

+			end <- seq(from=binsize, by=binsize, length.out=numbin)

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

 			## Create binned chromosome as GRanges object

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

@@ -273,44 +224,35 @@ align2binned <- function(file, format, ID=basename(file), bamindex=file, pairedE

 			suppressWarnings(

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

 			)

-		}

+

+		} ### end loop chromosomes ###

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

 		binned.data <- unlist(binned.data)

 		names(binned.data) <- NULL

-		## Count overlaps

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

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

-			mcounts <- GenomicRanges::countOverlaps(binned.data, data[strand(data)=='-'])

-			pcounts <- GenomicRanges::countOverlaps(binned.data, data[strand(data)=='+'])

-			counts <- mcounts + pcounts

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

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

-			midx <- as.matrix(findOverlaps(binned.data, data))

-			counts <- rep(0,length(binned.data))

-			signal <- mcols(data)$signal

-			rle <- rle(midx[,1])

-			read.idx <- rle$values

-			max.idx <- cumsum(rle$lengths)

-			maxvalues <- rep(NA, length(read.idx))

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

-			for (i1 in 2:length(read.idx)) {

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

-			}

-			counts[read.idx] <- maxvalues

-			mcounts <- rep(NA, length(counts))

-			pcounts <- rep(NA, length(counts))

-		}

-		if (is.null(numbins)) {

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

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

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

+		### Loop over offsets ###

+		countmatrices <- list()

+		if (is.null(stepsize)) {

+			offsets <- 0

 		} else {

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

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

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

+			offsets <- seq(from=0, to=binsize, by=as.integer(stepsize*binsize))

 		}

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

+		for (ioff in offsets) {

+			## Count overlaps

+			message(paste0("  counting overlaps for offset ",ioff," ..."), appendLF=F); ptm <- proc.time()

+			binned.data.shift <- suppressWarnings( shift(binned.data, shift=ioff) )

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

+				mcounts <- suppressWarnings( GenomicRanges::countOverlaps(binned.data.shift, data.minus) )

+				pcounts <- suppressWarnings( GenomicRanges::countOverlaps(binned.data.shift, data.plus) )

+				counts <- mcounts + pcounts

+			}

+			countmatrix <- matrix(c(counts,mcounts,pcounts), ncol=3)

+			colnames(countmatrix) <- c('counts','mcounts','pcounts')

+			countmatrices[[as.character(ioff)]] <- countmatrix

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

+		}

+		mcols(binned.data) <- as(countmatrices[['0']],'DataFrame')

+		attr(binned.data,'offset.counts') <- countmatrices

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

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

@@ -328,11 +270,7 @@ align2binned <- function(file, format, ID=basename(file), bamindex=file, pairedE

 		### Save or return the binned data ###

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

 			# Save to file

-			if (is.null(numbins)) {

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

-			} else {

-				filename <- paste0(ID,"_numbin_",format(numbin, scientific=T, trim=T),"_.RData")

-			}

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

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

 # 			attr(binned.data, 'call') <- call # do not store along with GRanges because it inflates disk usage

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

@@ -342,7 +280,7 @@ align2binned <- function(file, format, ID=basename(file), bamindex=file, pairedE

 			return(binned.data)

 		}

-	}

+	} ### end loop binsizes ###

 }

@@ -37,10 +37,15 @@ insertchr <- function(hmm.gr) {

 exportCNVs <- function(hmm.list, filename, cluster=TRUE, export.CNV=TRUE, export.SCE=TRUE) {

 	## Get segments and SCE coordinates

-	temp <- getSegments(hmm.list, cluster=cluster, getSCE=export.SCE)

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

+	temp <- getSegments(hmm.list, cluster=cluster)

 	hmm.grl <- temp$segments

+	if (cluster) {

+		hmm.list <- hmm.list[temp$clustering$order]

+	}

 	if (export.SCE) {

-		sce <- temp$sce

+		sce <- lapply(hmm.list,'[[','sce')

+		names(sce) <- lapply(hmm.list,'[[','ID')

 		sce <- sce[!unlist(lapply(sce, is.null))]

 		sce <- sce[lapply(sce, length)!=0]		

 		if (length(sce)==0) {

@@ -6,8 +6,8 @@

 #'

 #' \code{findCNVs} uses a 6-state Hidden Markov Model to classify the binned read counts: state 'nullsomy' with a delta function as emission densitiy (only zero read counts), '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

 #' @author Aaron Taudt

-#' @param method Currently only \code{c('univariate')}.

 #' @inheritParams univariate.findCNVs

+#' @param method One of \code{c('univariate','bivariate')}. In the univariate case strand information is discarded, while in the bivariate case strand information is used for the fitting.

 #' @examples

 #'## Get an example BAM file with single-cell-sequencing reads

 #'bamfile <- system.file("extdata/BB140820_I_002.bam", package="aneufinder")

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

 #'## Check the fit

 #'plot(model, type='histogram')

 #' @export

-findCNVs <- function(binned.data, ID=NULL, method='univariate', eps=0.1, init="standard", max.time=-1, max.iter=1000, num.trials=15, eps.try=10*eps, num.threads=1, count.cutoff.quantile=0.999, strand='*', states=c("zero-inflation","nullsomy","monosomy","disomy","trisomy","tetrasomy","multisomy"), most.frequent.state="disomy") {

+findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard", max.time=-1, max.iter=1000, num.trials=15, eps.try=10*eps, num.threads=1, count.cutoff.quantile=0.999, strand='*', states=c("zero-inflation","nullsomy","monosomy","disomy","trisomy","tetrasomy","multisomy"), most.frequent.state="disomy", method="univariate", algorithm="EM", allow.odd.states=TRUE, initial.params=NULL) {

 	## Intercept user input

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

@@ -37,8 +37,10 @@ findCNVs <- function(binned.data, ID=NULL, method='univariate', eps=0.1, init="s

 	ptm <- proc.time()

 	message("Find CNVs for ID = ",ID, ":")

-	if (method=='univariate') {

-		model <- univariate.findCNVs(binned.data, 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=count.cutoff.quantile, strand=strand, states=states, most.frequent.state=most.frequent.state)

+	if (method == 'univariate') {

+		model <- univariate.findCNVs(binned.data, 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=count.cutoff.quantile, strand=strand, states=states, most.frequent.state=most.frequent.state, algorithm=algorithm, initial.params=initial.params)

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

+		model <- bivariate.findCNVs(binned.data, 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=count.cutoff.quantile, allow.odd.states=allow.odd.states, states=states, most.frequent.state=most.frequent.state, initial.params=initial.params)

 	}

 	attr(model, 'call') <- call

@@ -70,8 +72,10 @@ findCNVs <- function(binned.data, ID=NULL, method='univariate', eps=0.1, init="s

 #' @param strand Run the HMM only for the specified strand. One of \code{c('+', '-', '*')}.

 #' @param states A subset or all of \code{c("zero-inflation","nullsomy","monosomy","disomy","trisomy","tetrasomy","multisomy")}. This vector defines the states that are used in the Hidden Markov Model. The order of the entries should not be changed.

 #' @param most.frequent.state One of the states that were given in \code{states} or 'none'. The specified state is assumed to be the most frequent one. This can help the fitting procedure to converge into the correct fit.

+#' @param algorithm One of \code{c('baumWelch','EM')}. The expectation maximization (\code{'EM'}) will find the most likely states and fit the best parameters to the data, the \code{'baumWelch'} will find the most likely states using the initial parameters.

+#' @param initial.params A \code{\link{aneuHMM}} object or file containing such an object from which initial starting parameters will be extracted.

 #' @return An \code{\link{aneuHMM}} object.

-univariate.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, strand='*', states=c("zero-inflation","monosomy","disomy","trisomy","tetrasomy","multisomy"), most.frequent.state="disomy") {

+univariate.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, strand='*', states=c("zero-inflation","monosomy","disomy","trisomy","tetrasomy","multisomy"), most.frequent.state="disomy", algorithm="EM", initial.params=NULL) {

 	### Define cleanup behaviour ###

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

@@ -94,6 +98,23 @@ univariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard",

 	if (check.positive.integer(num.threads)!=0) stop("argument 'num.threads' expects a positive integer")

 	if (check.strand(strand)!=0) stop("argument 'strand' expects either '+', '-' or '*'")

 	if (!most.frequent.state %in% c(states,'none')) stop("argument 'most.frequent.state' must be one of c(",paste(c(states,'none'), collapse=","),")")

+	if (!algorithm %in% c('baumWelch','EM')) {

+		stop("argument 'algorithm' expects one of c('baumWelch','EM')")

+	}

+	if (algorithm == 'baumWelch' & num.trials>1) {

+		warning("Set 'num.trials <- 1' because 'algorithm==\"baumWelch\"'.")

+		num.trials <- 1

+	}

+	initial.params <- loadHmmsFromFiles(initial.params)[[1]]

+	if (class(initial.params)!=class.univariate.hmm & !is.null(initial.params)) {

+		stop("argument 'initial.params' expects a ",class.univariate.hmm," object or file that contains such an object")

+	}

+	if (algorithm == 'baumWelch' & is.null(initial.params)) {

+		warning("'initial.params' should be specified if 'algorithm=\"baumWelch\"")

+	}

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

+		init <- 'initial.params'

+	}

 	warlist <- list()

 	if (is.null(eps.try) | num.trials==1) eps.try <- eps

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

 	temp <- initializeStates(states)

 	state.labels <- temp$labels

 	state.distributions <- temp$distributions

+	multiplicity <- temp$multiplicity

 	dependent.states.mask <- state.labels %in% c("monosomy","disomy","trisomy","tetrasomy","multisomy")

 	numstates <- length(states)

 	numbins <- length(binned.data)

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

 		select <- 'counts'

 	}

 	counts <- mcols(binned.data)[,select]

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

 	### Make return object

 		result <- list()

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

 		result$ID <- ID

 		result$bins <- binned.data

 	## Quality info

-		qualityInfo <- list(shannon.entropy=qc.entropy(counts), spikyness=qc.spikyness(counts), complexity=attr(result$bins, 'complexity.preseqR'), bhattacharyya=NA)

+		qualityInfo <- list(shannon.entropy=qc.entropy(counts), spikyness=qc.spikyness(counts), complexity=attr(result$bins,'qualityInfo')$complexity$preseqR, bhattacharyya=NA)

 		result$qualityInfo <- qualityInfo

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

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

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

 		stop(paste0("ID = ",ID,": NAs found in reads."))

 	}

@@ -156,13 +179,24 @@ univariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard",

 		message(paste0("Trial ",i_try," / ",num.trials))

 		## Initial parameters

-		if (init == 'random') {

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

+			A.initial <- initial.params$transitionProbs

+			proba.initial <- initial.params$startProbs

+			size.initial <- initial.params$distributions[,'size']

+			prob.initial <- initial.params$distributions[,'prob']

+			size.initial[is.na(size.initial)] <- 0

+			prob.initial[is.na(prob.initial)] <- 0

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

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

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

 			proba.initial <- runif(numstates)

 			# Distributions for dependent states

 			size.initial <- runif(1, min=0, max=100) * cumsum(dependent.states.mask)

 			prob.initial <- runif(1) * dependent.states.mask

+			# Assign initials for the nullsomy distribution

+			index <- which('nullsomy'==state.labels)

+			size.initial[index] <- 1

+			prob.initial[index] <- 0.5

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

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

 			for (irow in 1:numstates) {

@@ -172,8 +206,9 @@ univariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard",

 				}

 			}

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

-			mean.initial.monosomy <- mean(counts[counts>0])/2

-			var.initial.monosomy <- var(counts[counts>0])/2

+			divf <- max(multiplicity[most.frequent.state], 1)

+			mean.initial.monosomy <- mean(counts[counts>0])/divf

+			var.initial.monosomy <- var(counts[counts>0])/divf

 			if (is.na(mean.initial.monosomy)) {

 				mean.initial.monosomy <- 1

 			}

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

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

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

 			}

+			# Assign initials for the nullsomy distribution

+			index <- which('nullsomy'==state.labels)

+			size.initial[index] <- 1

+			prob.initial[index] <- 0.5

 		}

-		# Assign initials for the nullsomy distribution

-		index <- which('nullsomy'==state.labels)

-		size.initial[index] <- 1

-		prob.initial[index] <- 0.5

 		hmm <- .C("R_univariate_hmm",

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

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

 			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)

-			count.cutoff = as.integer(count.cutoff) # int* count.cutoff

+			count.cutoff = as.integer(count.cutoff), # int* count.cutoff

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

 		)

 		hmm$eps <- eps.try

@@ -237,13 +273,6 @@ univariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard",

 			# Store model in list

 			hmm$counts <- NULL

 			modellist[[as.character(i_try)]] <- hmm

-# 			# Check if most.frequent.state is more than half of actual most frequent state and stop trials

-# 			if (most.frequent.state %in% state.labels) {

-# 				imostfrequent <- which(state.labels==most.frequent.state)

-# 				if (hmm$weights[imostfrequent] / max(hmm$weights) > 0.5) {

-# 					break

-# 				}

-# 			}

 			init <- 'random'

 		} else if (num.trials == 1) {

 			if (hmm$loglik.delta > eps) {

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

 				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)

-				count.cutoff = as.integer(count.cutoff) # int* count.cutoff

+				count.cutoff = as.integer(count.cutoff), # int* count.cutoff

+				algorithm = as.integer(algorithm)

 			)

 		}

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

 			result$transitionProbs.initial <- transitionProbs.initial

 			# Initial probs

 			result$startProbs <- hmm$proba

-			names(result$startProbs) <- paste0("P(",state.labels,")")

+			names(result$startProbs) <- state.labels

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

-			names(result$startProbs.initial) <- paste0("P(",state.labels,")")

+			names(result$startProbs.initial) <- state.labels

 			# Distributions

 				distributions <- data.frame()

 				distributions.initial <- data.frame()

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

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

 			result$convergenceInfo <- convergenceInfo

 		## Quality info

-			qualityInfo <- list(shannon.entropy=qc.entropy(counts), spikyness=qc.spikyness(counts), complexity=attr(result$bins, 'complexity.preseqR'), bhattacharyya=qc.bhattacharyya(result))

+			qualityInfo <- list(shannon.entropy=qc.entropy(counts), spikyness=qc.spikyness(counts), complexity=attr(result$bins,'qualityInfo')$complexity$preseqR, bhattacharyya=qc.bhattacharyya(result))

 			result$qualityInfo <- qualityInfo

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

@@ -389,9 +419,426 @@ univariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard",

 	## Issue warnings

 	result$warnings <- warlist

-#   return(list(result, modellist))

 	## Return results

 	return(result)

 }

+#' Find copy number variations (bivariate)

+#'

+#' \code{bivariate.findCNVs} finds CNVs using read count information from both strands.

+#'

+#' @inheritParams univariate.findCNVs

+#' @param allow.odd.states If set to \code{TRUE}, all states are allowed. If \code{FALSE}, only states which have an even multiplicity (plus nullsomy-monosomy states) are allowed, e.g. disomy-disomy, monosomy-trisomy.

+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, allow.odd.states=TRUE, states=c("zero-inflation","nullsomy","monosomy","disomy","trisomy","tetrasomy","multisomy"), most.frequent.state="monosomy", algorithm='EM', initial.params=NULL) {

+

+	## Intercept user input

+	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 'counts' or a file that contains such an object")

+	}

+	if (is.null(ID)) {

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

+	}

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

+	if (check.positive.integer(num.trials)!=0) stop("argument 'num.trials' expects a positive integer")

+	if (!is.null(eps.try)) {

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

+	}

+	if (check.positive.integer(num.threads)!=0) stop("argument 'num.threads' expects a positive integer")

+	initial.params <- loadHmmsFromFiles(initial.params)[[1]]

+	if (class(initial.params)!=class.bivariate.hmm & !is.null(initial.params)) {

+		stop("argument 'initial.params' expects a ",class.bivariate.hmm," object or file that contains such an object")

+	}

+	if (algorithm == 'baumWelch' & is.null(initial.params)) {

+		warning("'initial.params' should be specified if 'algorithm=\"baumWelch\"")

+	}

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

+		init <- 'initial.params'

+	}

+

+	warlist <- list()

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

+

+	## Variables

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

+	temp <-  initializeStates(states)

+	multiplicity <- temp$multiplicity

+	state.labels <- temp$labels

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

+	}

+	

+	### Make return object

+		result <- list()

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

+		result$ID <- ID

+		result$bins <- binned.data

+	## Quality info

+		qualityInfo <- list(shannon.entropy=qc.entropy(counts), spikyness=qc.spikyness(counts), complexity=attr(result$bins,'qualityInfo')$complexity$preseqR, bhattacharyya=NA)

+		result$qualityInfo <- qualityInfo

+

+	# 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 (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=''))

+		message("Running univariate")

+		binned.data.minus <- binned.data

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

+		binned.data.minus$counts <- binned.data.minus$mcounts

+		binned.data.plus <- binned.data

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

+		binned.data.plus$counts <- binned.data.plus$pcounts

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

+

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

+		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])

+				if (allow.odd.states) {

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

+				} else {

+					state.multiplicity <- multiplicity[levels.state[i1]] + multiplicity[levels.state[i2]]

+					# Only even states and null-mono states (for sex chromosomes)

+					if (state.multiplicity %% 2 == 0 | state.multiplicity == 1) {

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