@@ -1,14 +1,14 @@

 Package: AneuFinder

 Type: Package

 Title: Analysis of Copy Number Variation in Single-Cell-Sequencing Data

-Version: 0.99.1

-Date: 2015-09-19

+Version: 0.99.3

+Date: 2016-03

 Author: Aaron Taudt, Bjorn Bakker, David Porubsky

 Maintainer: Aaron Taudt <aaron.taudt@gmail.com>

 Description: This package implements functions for CNV calling, plotting, export

     and analysis from whole-genome single cell sequencing data.

 Depends:

-    R (>= 3.2.3),

+    R (>= 3.3.0),

     GenomicRanges,

     cowplot,

     AneuFinderData

@@ -91,4 +91,4 @@ importFrom(stats,stepfun)

 importFrom(stats,weighted.mean)

 importFrom(utils,read.table)

 importFrom(utils,write.table)

-useDynLib(AneuFinder)

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

@@ -25,6 +25,7 @@

 #' @param bw Bandwidth for SCE hotspot detection (see \code{\link{hotspotter}} for further details).

 #' @param pval P-value for SCE hotspot detection (see \code{\link{hotspotter}} for further details).

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

+#' @return \code{NULL}

 #' @author Aaron Taudt

 #' @import foreach

 #' @import doParallel

@@ -216,7 +216,7 @@ binReads <- function(file, format, assembly, ID=basename(file), bamindex=file, c

 		}

 		### Quality measures ###

-		qualityInfo <- list(complexity=complexity, coverage=coverage, spikyness=qc.spikyness(binned.data$counts), shannon.entropy=qc.entropy(binned.data$counts))

+		qualityInfo <- list(complexity=complexity, coverage=coverage, spikiness=qc.spikiness(binned.data$counts), shannon.entropy=qc.entropy(binned.data$counts))

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

 		attr(binned.data, 'min.mapq') <- min.mapq

@@ -63,7 +63,7 @@ correctMappability <- function(binned.data.list, same.binsize, reference, format

 		### Quality measures ###

 		## Spikyness

-		attr(binned.data, 'spikyness') <- qc.spikyness(binned.data$counts)

+		attr(binned.data, 'spikiness') <- qc.spikiness(binned.data$counts)

 		## Shannon entropy

 		attr(binned.data, 'shannon.entropy') <- qc.entropy(binned.data$counts)

@@ -200,7 +200,7 @@ correctGC <- function(binned.data.list, GC.BSgenome, same.binsize=FALSE) {

 		### Quality measures ###

 		## Spikyness

-		attr(binned.data, 'spikyness') <- qc.spikyness(binned.data$counts)

+		attr(binned.data, 'spikiness') <- qc.spikiness(binned.data$counts)

 		## Shannon entropy

 		attr(binned.data, 'shannon.entropy') <- qc.entropy(binned.data$counts)

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

 #' Use \code{exportReadCounts} to export the binned read counts from an \code{\link{aneuHMM}} object in WIGGLE format.

 #' Use \code{exportGRanges} to export a \code{\link{GRanges}} object in BED format.

 #'

+#' @return \code{NULL}

 #' @name export

 #' @author Aaron Taudt

 #'@examples

@@ -84,7 +84,7 @@ findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard", max.time=-1

 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",paste0(0:10,"-somy")), most.frequent.state="2-somy", algorithm="EM", initial.params=NULL) {

 	### Define cleanup behaviour ###

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

+	on.exit(.C("C_univariate_cleanup", PACKAGE = 'AneuFinder'))

 	## Intercept user input

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

@@ -150,7 +150,7 @@ 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,'qualityInfo')$complexity$preseqR, bhattacharyya=NA)

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

 		result$qualityInfo <- qualityInfo

 	## Convergence info

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

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

 			prob.initial[index] <- 0.5

 		}

-		hmm <- .C("R_univariate_hmm",

+		hmm <- .C("C_univariate_hmm",

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

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

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

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

 			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

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

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

+			PACKAGE = 'AneuFinder'

 		)

 		hmm$eps <- eps.try

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

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

 			message(paste0("Rerunning trial ",index2use," with eps = ",eps))

-			hmm <- .C("R_univariate_hmm",

+			hmm <- .C("C_univariate_hmm",

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

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

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

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

 				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

-				algorithm = as.integer(algorithm)

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

+				PACKAGE = 'AneuFinder'

 			)

 		}

@@ -411,7 +413,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,'qualityInfo')$complexity$preseqR, bhattacharyya=qc.bhattacharyya(result))

+			qualityInfo <- list(shannon.entropy=qc.entropy(counts), spikiness=qc.spikiness(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."))

@@ -497,7 +499,7 @@ bivariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard", m

 		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)

+		qualityInfo <- list(shannon.entropy=qc.entropy(counts), spikiness=qc.spikiness(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

@@ -712,11 +714,11 @@ bivariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard", m

 	stopTimedMessage(ptm)

 	### Define cleanup behaviour ###

-	on.exit(.C("R_multivariate_cleanup", as.integer(num.comb.states)))

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

 	### Run the multivariate HMM

 	# Call the C function

-	hmm <- .C("R_multivariate_hmm",

+	hmm <- .C("C_multivariate_hmm",

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

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

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

@@ -734,7 +736,8 @@ bivariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard", m

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

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

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

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

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

+		PACKAGE = 'AneuFinder'

 		)

 	### Check convergence ###

@@ -811,7 +814,7 @@ bivariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard", m

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

 			result$convergenceInfo <- convergenceInfo

 		## Quality info

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

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

 			result$qualityInfo <- qualityInfo

 		## Univariate infos

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

@@ -79,8 +79,8 @@ findSCEs <- function(binned.data, ID=NULL, eps=0.1, init="standard", max.time=-1

 #' @export

 #'@examples

 #'## Load an HMM

-#'file <- system.file("extdata", "primary-lung", "hmms", "AvdB150303_I_001.bam.RData",

-#'                     package="AneuFinderData")

+#'file <- list.files(system.file("extdata", "primary-lung", "hmms",

+#'                   package="AneuFinderData"), full.names=TRUE)

 #'hmm <- loadHmmsFromFiles(file)[[1]]

 #'## Check number of segments before and after filtering

 #'length(hmm$segments)

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

-#' @useDynLib AneuFinder

+#' @useDynLib AneuFinder, .registration = TRUE, .fixes = ""

 #' @import IRanges

 #' @import GenomeInfoDb

 #' @import GenomicRanges

@@ -258,7 +258,7 @@ plotBivariateHistograms <- function(bihmm) {

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

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

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

-	mask.attributes <- c(grep('complexity', names(attributes(binned.data)), value=TRUE), 'spikyness', 'shannon.entropy')

+	mask.attributes <- c(grep('complexity', names(attributes(binned.data)), value=TRUE), 'spikiness', 'shannon.entropy')

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

 	## Make fake uni.hmm and plot

@@ -369,10 +369,10 @@ plotUnivariateHistogram <- function(model, state=NULL, strand='*', chromosome=NU

 	# Quality info

 	if (is.null(model$qualityInfo$complexity)) { model$qualityInfo$complexity <- NA }

-	if (is.null(model$qualityInfo$spikyness)) { model$qualityInfo$spikyness <- NA }

+	if (is.null(model$qualityInfo$spikiness)) { model$qualityInfo$spikiness <- NA }

 	if (is.null(model$qualityInfo$shannon.entropy)) { model$qualityInfo$shannon.entropy <- NA }

 	if (is.null(model$qualityInfo$bhattacharyya)) { model$qualityInfo$bhattacharyya <- NA }

-	quality.string <- paste0('reads = ',round(sum(model$bins$counts)/1e6,2),'M, complexity = ',round(model$qualityInfo$complexity),',  spikyness = ',round(model$qualityInfo$spikyness,2),',  entropy = ',round(model$qualityInfo$shannon.entropy,2),',  bhattacharyya = ',round(model$qualityInfo$bhattacharyya,2), ', num.segments = ',length(model$segments))

+	quality.string <- paste0('reads = ',round(sum(model$bins$counts)/1e6,2),'M, complexity = ',round(model$qualityInfo$complexity),',  spikiness = ',round(model$qualityInfo$spikiness,2),',  entropy = ',round(model$qualityInfo$shannon.entropy,2),',  bhattacharyya = ',round(model$qualityInfo$bhattacharyya,2), ', num.segments = ',length(model$segments))

 	# Plot the histogram

 	ggplt <- ggplot(data.frame(counts)) + geom_histogram(aes_string(x='counts', y='..density..'), binwidth=1, color='black', fill='white') + coord_cartesian(xlim=c(0,rightxlim)) + theme_bw() + xlab("read count") + ggtitle(bquote(atop(.(model$ID), atop(.(quality.string),''))))

@@ -453,7 +453,7 @@ plotKaryogram <- function(model, both.strands=FALSE, plot.SCE=FALSE, file=NULL)

 		model <- list()

 		model$ID <- ''

 		model$bins <- binned.data

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

+		model$qualityInfo <- list(shannon.entropy=qc.entropy(binned.data$counts), spikiness=qc.spikiness(binned.data$counts), complexity=attr(binned.data, 'complexity.preseqR'), bhattacharyya=NA)

 		plot.karyogram(model, both.strands=both.strands, file=file)

 	} else if (class(model)==class.univariate.hmm) {

 		plot.karyogram(model, both.strands=both.strands, file=file)

@@ -496,10 +496,10 @@ plot.karyogram <- function(model, both.strands=FALSE, plot.SCE=FALSE, file=NULL)

 	# Quality info

 	if (is.null(model$qualityInfo$complexity)) { model$qualityInfo$complexity <- NA }

-	if (is.null(model$qualityInfo$spikyness)) { model$qualityInfo$spikyness <- NA }

+	if (is.null(model$qualityInfo$spikiness)) { model$qualityInfo$spikiness <- NA }

 	if (is.null(model$qualityInfo$shannon.entropy)) { model$qualityInfo$shannon.entropy <- NA }

 	if (is.null(model$qualityInfo$bhattacharyya)) { model$qualityInfo$bhattacharyya <- NA }

-	quality.string <- paste0('reads = ',round(sum(model$bins$counts)/1e6,2),'M, complexity = ',round(model$qualityInfo$complexity),',  spikyness = ',round(model$qualityInfo$spikyness,2),',  entropy = ',round(model$qualityInfo$shannon.entropy,2),',  bhattacharyya = ',round(model$qualityInfo$bhattacharyya,2), ', num.segments = ',length(model$segments))

+	quality.string <- paste0('reads = ',round(sum(model$bins$counts)/1e6,2),'M, complexity = ',round(model$qualityInfo$complexity),',  spikiness = ',round(model$qualityInfo$spikiness,2),',  entropy = ',round(model$qualityInfo$shannon.entropy,2),',  bhattacharyya = ',round(model$qualityInfo$bhattacharyya,2), ', num.segments = ',length(model$segments))

 	## Get SCE coordinates

 	if (plot.SCE) {

@@ -918,7 +918,7 @@ plotProfile <- function(model, both.strands=FALSE, plot.SCE=TRUE, file=NULL) {

 		model <- list()

 		model$ID <- ''

 		model$bins <- binned.data

-		model$qualityInfo <- list(shannon.entropy=qc.entropy(binned.data$counts), spikyness=qc.spikyness(binned.data$counts), complexity=attr(binned.data, 'complexity.preseqR'))

+		model$qualityInfo <- list(shannon.entropy=qc.entropy(binned.data$counts), spikiness=qc.spikiness(binned.data$counts), complexity=attr(binned.data, 'complexity.preseqR'))

 		plot.profile(model, both.strands=both.strands, plot.SCE=FALSE, file=file)

 	} else if (class(model)==class.univariate.hmm) {

 		plot.profile(model, both.strands=FALSE, plot.SCE=FALSE, file=file)

@@ -1041,10 +1041,10 @@ plot.profile <- function(model, both.strands=FALSE, plot.SCE=TRUE, file=NULL) {

 	ggplt <- ggplt + scale_x_continuous(breaks=seqlengths(model$bins)/2+cum.seqlengths.0[as.character(seqlevels(model$bins))], labels=seqlevels(model$bins))

 	# Quality info

 	if (is.null(model$qualityInfo$complexity)) { model$qualityInfo$complexity <- NA }

-	if (is.null(model$qualityInfo$spikyness)) { model$qualityInfo$spikyness <- NA }

+	if (is.null(model$qualityInfo$spikiness)) { model$qualityInfo$spikiness <- NA }

 	if (is.null(model$qualityInfo$shannon.entropy)) { model$qualityInfo$shannon.entropy <- NA }

 	if (is.null(model$qualityInfo$bhattacharyya)) { model$qualityInfo$bhattacharyya <- NA }

-	quality.string <- paste0('reads = ',round(sum(model$bins$counts)/1e6,2),'M, complexity = ',round(model$qualityInfo$complexity),',  spikyness = ',round(model$qualityInfo$spikyness,2),',  entropy = ',round(model$qualityInfo$shannon.entropy,2),',  bhattacharyya = ',round(model$qualityInfo$bhattacharyya,2), ', num.segments = ',length(model$segments))

+	quality.string <- paste0('reads = ',round(sum(model$bins$counts)/1e6,2),'M, complexity = ',round(model$qualityInfo$complexity),',  spikiness = ',round(model$qualityInfo$spikiness,2),',  entropy = ',round(model$qualityInfo$shannon.entropy,2),',  bhattacharyya = ',round(model$qualityInfo$bhattacharyya,2), ', num.segments = ',length(model$segments))

 	ggplt <- ggplt + ylab('read count') + ggtitle(bquote(atop(.(model$ID), atop(.(quality.string),''))))

 	if (!is.null(file)) {

@@ -15,12 +15,12 @@

 #' @author Aaron Taudt

 NULL

-#' @describeIn qualityControl Calculate the spikyness of a library

-qc.spikyness <- function(counts) {

+#' @describeIn qualityControl Calculate the spikiness of a library

+qc.spikiness <- function(counts) {

 	counts <- as.vector(counts)

 	sum.counts <- sum(counts)

-	spikyness <- sum(abs(diff(counts))) / sum.counts

-	return(spikyness)

+	spikiness <- sum(abs(diff(counts))) / sum.counts

+	return(spikiness)

 }

 #' @describeIn qualityControl Calculate the Shannon entropy of a library

@@ -59,7 +59,7 @@ getQC <- function(hmms) {

 			qframe[[i1]] <- data.frame(total.read.count=sum(hmm$bins$counts),

 														binsize=width(hmm$bins)[1],

 														avg.read.count=mean(hmm$bins$counts),

-														spikyness=hmm$qualityInfo$spikyness,

+														spikiness=hmm$qualityInfo$spikiness,

 														entropy=hmm$qualityInfo$shannon.entropy,

 														complexity=hmm$qualityInfo$complexity,

 														loglik=hmm$convergenceInfo$loglik,

@@ -79,7 +79,7 @@ getQC <- function(hmms) {

 #'

 #' The employed quality measures are:

 #' \itemize{

-#' \item Spikyness

+#' \item Spikiness

 #' \item Entropy

 #' \item Number of segments

 #' \item Bhattacharrya distance

@@ -89,8 +89,8 @@ getQC <- function(hmms) {

 #' @param hmms A list of \code{\link{aneuHMM}} objects or a list of files that contain such objects.

 #' @param G An integer vector specifying the number of clusters that are compared. See \code{\link[mclust:Mclust]{Mclust}} for details.

 #' @param itmax The maximum number of outer and inner iterations for the \code{\link[mclust:Mclust]{Mclust}} function. See \code{\link[mclust:emControl]{emControl}} for details.

-#' @param measures The quality measures that are used for the clustering. Supported is any combination of \code{c('spikyness','entropy','num.segments','bhattacharyya','loglik','complexity','avg.read.count','total.read.count','binsize')}. 

-#' @param orderBy The quality measure to order the clusters by. Default is \code{'spikyness'}.

+#' @param measures The quality measures that are used for the clustering. Supported is any combination of \code{c('spikiness','entropy','num.segments','bhattacharyya','loglik','complexity','avg.read.count','total.read.count','binsize')}. 

+#' @param orderBy The quality measure to order the clusters by. Default is \code{'spikiness'}.

 #' @param reverseOrder Logical indicating whether the ordering by \code{orderBy} is reversed.

 #' @return A \code{list} with the classification, parameters and the \code{\link[mclust]{Mclust}} fit.

 #' @author Aaron Taudt

@@ -108,7 +108,7 @@ getQC <- function(hmms) {

 #'## Select files from the best 2 clusters for further processing

 #'best.files <- unlist(cl$classification[1:2])

 #'

-clusterByQuality <- function(hmms, G=1:9, itmax=c(100,100), measures=c('spikyness','entropy','num.segments','bhattacharyya','loglik'), orderBy='spikyness', reverseOrder=FALSE) {

+clusterByQuality <- function(hmms, G=1:9, itmax=c(100,100), measures=c('spikiness','entropy','num.segments','bhattacharyya','loglik'), orderBy='spikiness', reverseOrder=FALSE) {

 	hmms <- loadHmmsFromFiles(hmms)

 	df <- getQC(hmms)

@@ -41,6 +41,7 @@ readConfig <- function(configfile) {

 #'

 #' @param conf A list structure with parameter values. Each entry will be written in one line.

 #' @param configfile Filename of the outputfile.

+#' @return \code{NULL}

 #' @author Aaron Taudt

 writeConfig <- function(conf, configfile) {

@@ -2,6 +2,6 @@

 .onAttach <- function(...) {

-	packageStartupMessage('Loading AneuFinder 0.99.2')

+	packageStartupMessage('Loading AneuFinder 0.99.3')

 }

@@ -85,6 +85,9 @@ Aneufinder(inputfolder, outputfolder, format, configfile = NULL, numCPU = 1,

 \item{cluster.plots}{A logical indicating whether plots should be clustered by similarity.}

 }

+\value{

+\code{NULL}

+}

 \description{

 This function is an easy-to-use wrapper to \link[AneuFinder:binning]{bin the data}, \link[AneuFinder:findCNVs]{find copy-number-variations}, \link[AneuFinder:findSCEs]{find sister-chromatid-exchange} events, plot \link[AneuFinder:heatmapGenomewide]{genomewide heatmaps}, \link[AneuFinder:plot.aneuHMM]{distributions, profiles and karyograms}.

 }

@@ -5,8 +5,8 @@

 \title{Cluster based on quality variables}

 \usage{

 clusterByQuality(hmms, G = 1:9, itmax = c(100, 100),

-  measures = c("spikyness", "entropy", "num.segments", "bhattacharyya",

-  "loglik"), orderBy = "spikyness", reverseOrder = FALSE)

+  measures = c("spikiness", "entropy", "num.segments", "bhattacharyya",

+  "loglik"), orderBy = "spikiness", reverseOrder = FALSE)

 }

 \arguments{

 \item{hmms}{A list of \code{\link{aneuHMM}} objects or a list of files that contain such objects.}

@@ -15,9 +15,9 @@ clusterByQuality(hmms, G = 1:9, itmax = c(100, 100),

 \item{itmax}{The maximum number of outer and inner iterations for the \code{\link[mclust:Mclust]{Mclust}} function. See \code{\link[mclust:emControl]{emControl}} for details.}

-\item{measures}{The quality measures that are used for the clustering. Supported is any combination of \code{c('spikyness','entropy','num.segments','bhattacharyya','loglik','complexity','avg.read.count','total.read.count','binsize')}.}

+\item{measures}{The quality measures that are used for the clustering. Supported is any combination of \code{c('spikiness','entropy','num.segments','bhattacharyya','loglik','complexity','avg.read.count','total.read.count','binsize')}.}

-\item{orderBy}{The quality measure to order the clusters by. Default is \code{'spikyness'}.}

+\item{orderBy}{The quality measure to order the clusters by. Default is \code{'spikiness'}.}

 \item{reverseOrder}{Logical indicating whether the ordering by \code{orderBy} is reversed.}

 }

@@ -30,7 +30,7 @@ This function uses the \pkg{\link{mclust}} package to cluster the input samples

 \details{

 The employed quality measures are:

 \itemize{

-\item Spikyness

+\item Spikiness

 \item Entropy

 \item Number of segments

 \item Bhattacharrya distance

@@ -40,6 +40,9 @@ exportGRanges(gr, filename, header = TRUE, trackname = NULL, score = NULL,

 \item{chromosome.format}{A character specifying the format of the chromosomes if \code{assembly} is specified. Either 'NCBI' for (1,2,3 ...) or 'UCSC' for (chr1,chr2,chr3 ...).#' @importFrom utils write.table}

 }

+\value{

+\code{NULL}

+}

 \description{

 Export copy-number-variation state or read counts as genome browser viewable file

 }

@@ -19,8 +19,8 @@ The input \code{model} with adjusted segments.

 }

 \examples{

 ## Load an HMM

-file <- system.file("extdata", "primary-lung", "hmms", "AvdB150303_I_001.bam.RData",

-                    package="AneuFinderData")

+file <- list.files(system.file("extdata", "primary-lung", "hmms",

+                  package="AneuFinderData"), full.names=TRUE)

 hmm <- loadHmmsFromFiles(file)[[1]]

 ## Check number of segments before and after filtering

 length(hmm$segments)

@@ -3,11 +3,11 @@

 \name{qualityControl}

 \alias{qc.bhattacharyya}

 \alias{qc.entropy}

-\alias{qc.spikyness}

+\alias{qc.spikiness}

 \alias{qualityControl}

 \title{Quality control measures for binned read counts}

 \usage{

-qc.spikyness(counts)

+qc.spikiness(counts)

 qc.entropy(counts)

@@ -31,7 +31,7 @@ Spikyness is defined as \eqn{K = sum(abs(diff(counts))) / sum(counts)}.

 }

 \section{Functions}{

 \itemize{

-\item \code{qc.spikyness}: Calculate the spikyness of a library

+\item \code{qc.spikiness}: Calculate the spikiness of a library

 \item \code{qc.entropy}: Calculate the Shannon entropy of a library

@@ -11,6 +11,9 @@ writeConfig(conf, configfile)

 \item{configfile}{Filename of the outputfile.}

 }

+\value{

+\code{NULL}

+}

 \description{

 Write an AneuFinder configuration file from a list structure.

 }

@@ -9,8 +9,7 @@ static double** multiD;

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

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

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

-extern "C" {

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

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

 {

 	// Define logging level

@@ -209,14 +208,12 @@ void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, d

 	delete hmm;

 	hmm = NULL; // assign NULL to defuse the additional delete in on.exit() call

 }

-} //extern

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

 // This function takes parameters from R, creates a multivariate HMM object, runs the EM 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, int* algorithm)

+void 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, int* algorithm)

 {

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

@@ -351,25 +348,20 @@ void R_multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states,

 	hmm = NULL; // assign NULL to defuse the additional delete in on.exit() call

 // 	FreeDoubleMatrix(multiD, *N);

 }

-} //extern

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

 // This function make a cleanup if anything was left over

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

-extern "C" {

-void R_univariate_cleanup()

+void univariate_cleanup()

 {

 // 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__; // This message will be shown if interrupt happens before start of C-code

 	delete hmm;

 }

-} //extern

-extern "C" {

-void R_multivariate_cleanup(int* N)

+void multivariate_cleanup(int* N)

 {

 	delete hmm;

 	FreeDoubleMatrix(multiD, *N);

 }

-} //extern

@@ -12,11 +12,15 @@

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

 // #endif

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

-// 

-// 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, int* algorithm);

-// 

-// void R_univariate_cleanup();

-// 

-// void R_multivariate_cleanup(int* N);

-// 

+extern "C"

+void univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, double* prob, int* maxiter, int* maxtime, double* eps, int* states, double* A, double* proba, double* loglik, double* weights, int* distr_type, double* initial_size, double* initial_prob, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* read_cutoff, int* algorithm);

+

+extern "C"

+void 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, int* algorithm);

+

+extern "C"

+void univariate_cleanup();

+

+extern "C"

+void multivariate_cleanup(int* N);

+

@@ -1,21 +1,22 @@

-// #include <Rinternals.h>

-// #include <R_ext/Rdynload.h>

-// #include <R_ext/Visibility.h>

-// #include "R_interface.h"

-// 

-// 

-// static const R_CMethodDef CEntries[]  = {

-//     {"R_univariate_hmm", (DL_FUNC) &R_univariate_hmm, 24},

-//     {"R_multivariate_hmm", (DL_FUNC) &R_multivariate_hmm, 18},

-//     {"R_univariate_cleanup", (DL_FUNC) &R_univariate_cleanup, 0},

-//     {"R_multivariate_cleanup", (DL_FUNC) &R_multivariate_cleanup, 1},

-//     {NULL, NULL, 0}

-// };

-// 

-// 

-// void attribute_visible R_init_aneufinder(DllInfo *dll)

-// {

-//     R_registerRoutines(dll, CEntries, NULL, NULL, NULL);

-//     R_useDynamicSymbols(dll, FALSE);

-//     R_forceSymbols(dll, TRUE);

-// }

+#include <Rinternals.h>

+#include <R_ext/Rdynload.h>

+#include "R_interface.h"

+

+

+static const R_CMethodDef CEntries[]  = {

+    {"C_univariate_hmm", (DL_FUNC) &univariate_hmm, 24, (R_NativePrimitiveArgType[24]) {INTSXP, INTSXP, INTSXP, INTSXP, REALSXP, REALSXP, INTSXP, INTSXP, REALSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, LGLSXP, INTSXP, INTSXP, INTSXP, INTSXP}},

+    {"C_multivariate_hmm", (DL_FUNC) &multivariate_hmm, 18, (R_NativePrimitiveArgType[18]) {REALSXP, INTSXP, INTSXP, INTSXP, INTSXP, INTSXP, INTSXP, REALSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, REALSXP, LGLSXP, INTSXP, INTSXP, INTSXP}},

+    {"C_univariate_cleanup", (DL_FUNC) &univariate_cleanup, 0, NULL},

+    {"C_multivariate_cleanup", (DL_FUNC) &multivariate_cleanup, 1, (R_NativePrimitiveArgType[1]) {INTSXP}},

+    {NULL, NULL, 0, NULL}

+};

+

+

+extern "C" {

+void R_init_AneuFinder(DllInfo *dll)

+{

+	R_registerRoutines(dll, CEntries, NULL, NULL, NULL);

+	R_useDynamicSymbols(dll, FALSE);

+// 	R_forceSymbols(dll, TRUE);

+}

+}

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

 /* custom error handling class */

 // extern statement to avoid 'multiple definition' errors

-extern class exception_nan: public std::exception

+const class exception_nan: public std::exception

 {

   virtual const char* what() const throw()

   {

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

 \tableofcontents

 \clearpage

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

-<<options,echo=FALSE, results='hide', message=FALSE, eval=TRUE>>=

+<<options, echo=FALSE, results='hide', message=FALSE, eval=TRUE>>=

 options(fig.width=120)

 library(AneuFinder)

 @

@@ -51,7 +51,7 @@ Aneufinder(inputfolder='folder-with-BAM', outputfolder='output-directory',

 Although in most cases the above command will produce reasonably good results, it might be worthwile to adjust the default parameters to improve performance and the quality of the results (see section \ref{sec:workflow}). You can get a description of all available parameters by typing

 \begin{scriptsize}

-<<>>==

+<<eval=TRUE>>==

 ?Aneufinder

 @

 \end{scriptsize}

@@ -90,32 +90,35 @@ heatmapGenomewide(cl$classification[[1]])

 \section{\label{sec:workflow}A detailed workflow}

-\subsection{Mappability correction}

+\subsection{\label{sec:mapcor}Mappability correction}

 The first step of your workflow should be the production of a reference file for mappability correction. Mappability correction is done via a variable-width binning approach (as compared to fixed-width bins) and requires a euploid reference. You can either simulate this reference file or take a real euploid reference. For optimal results we suggest to use a real reference, e.g. by merging BAM files of single cells from a euploid reference tissue. This can be achieved with the 'samtools merge' command (not part of R). Be careful: All CNVs that are present in the reference will lead to artifacts in the analysis later. This includes sex-chromosomes that are present in one copy only, so we advice to use a female reference and to exclude the Y-chromosome from the analysis. If you have no reference available, you can simulate one with the \Rfunction{simulateReads} command:

 \begin{scriptsize}

-<<eval=FALSE>>==

+<<eval=TRUE, message=FALSE, warning=FALSE>>==

 ## Load human genome

 library(BSgenome.Hsapiens.UCSC.hg19)

-## Get a BAM file for the estimation of quality scores

+## Get a BAM file for the estimation of quality scores (adjust this to your experiment)

 bamfile <- system.file("extdata", "BB150803_IV_074.bam", package="AneuFinderData")

 ## Simulate reads of length 51bp for human genome

-outputfile <- tempfile() # replace this with your destination file

+# We simulate reads every 5000bp for demonstration purposes, for a real

+# application you should use a much denser spacing (e.g. 500bp or less)

+simulatedReads.file <- tempfile() # replace this with your destination file

 simulateReads(BSgenome.Hsapiens.UCSC.hg19, readLength=51, bamfile=bamfile,

-              file=outputfile, every.X.bp=500)

+              file=simulatedReads.file, every.X.bp=5000)

 @

 \end{scriptsize}

 This simulated FASTQ file must then be aligned with your aligner of choice (ideally the same that you used for your other samples) and given as reference in the \Rfunction{Aneufinder} function (option \Roption{variable.width.reference}).

-\subsection{Blacklisting}

+\subsection{\label{sec:blacklist}Blacklisting}

 To further improve the quality of the results and remove artifacts caused by high mappability repeat regions, e.g. near centromers, a blacklist can be used in option \Roption{blacklist} of the \Rfunction{Aneufinder} function. All reads falling into the regions specified by the blacklist will be discarded when importing the read files.

 You can either download a blacklist from the UCSC genome browser, e.g. the ``DAC Blacklisted Regions from ENCODE/DAC(Kundaje)'' mappability track, or make your own. For optimal results, we advice to make your own blacklist from a euploid reference. The following code chunck takes a euploid reference and makes fixed-width bins of 100kb. Bins with read count above and beloow the 0.999 and 0.05 quantile are taken as blacklist:

 \begin{scriptsize}

 <<eval=TRUE, warning=FALSE, message=FALSE, fig.height=3>>==

-## Get a euploid reference

+## Get a euploid reference (adjust this to your experiment)

 bedfile <- system.file("extdata", "hg19_diploid.bam.bed.gz", package="AneuFinderData")

+

 ## Make 100kb fixed-width bins

 bins <- binReads(bedfile, format='bed', assembly='hg19',

                  binsizes=100e3, chromosomes=c(1:22,'X'))[[1]]

@@ -128,21 +131,85 @@ p <- p + geom_hline(aes(yintercept=ucutoff), color='red')

 print(p)

 blacklist <- bins[bins$counts < ucutoff & bins$counts > lcutoff]

 ## Write blacklist to file

-exportGRanges(blacklist, filename=tempfile(), header=FALSE,

+blacklist.file <- tempfile()

+exportGRanges(blacklist, filename=blacklist.file, header=FALSE,

               chromosome.format='NCBI')

 @

 \end{scriptsize}

 \subsection{Running Aneufinder}

+The function \Rfunction{Aneufinder} takes an input folder with BAM or BED files and produces an output folder with results, plots and browserfiles. The following code is an example of how to run \Rfunction{Aneufinder} with variable-width bins (see section \ref{sec:mapcor}), blacklist (see section \ref{sec:blacklist}) and GC-correction. Results will be stored in 'outputfolder/hmms' as RData objects for further processing such as quality filtering and customized plotting.

+

+\begin{scriptsize}

+<<eval=TRUE>>==

+## First, get some data and reference files (adjust this to your experiment)

+var.width.ref <- system.file("extdata", "hg19_diploid.bam.bed.gz", package="AneuFinderData")

+blacklist <- system.file("extdata", "blacklist-hg19.bed.gz", package="AneuFinderData")

+datafolder <- system.file("extdata", "B-ALL-B", package = "AneuFinderData")

+list.files(datafolder) # only 3 cells for demonstration purposes

+

+## Library for GC correction

+library(BSgenome.Hsapiens.UCSC.hg19)

+

+## Produce output files

+cnv.states <- c('zero-inflation', paste0(1:10, '-somy'))

+outputfolder <- tempdir()

+Aneufinder(inputfolder = datafolder, outputfolder = outputfolder, format = 'bed',

+           numCPU = 3, binsizes = c(5e5, 1e6), variable.width.reference = var.width.ref,

+           chromosomes = c(1:22,'X','Y'), blacklist = blacklist, states = cnv.states,

+           correction.method = 'GC', GC.BSgenome = BSgenome.Hsapiens.UCSC.hg19,

+           num.trials = 1, assembly = 'hg19')

+@

+\end{scriptsize}

 \subsection{Quality control}

+Once the function \Rfunction{Aneufinder} has completed, results will be accessible as .RData files under 'outputfolder/hmms'. Single cell sequencing is prone to noise and therefore it is a good idea to filter the results by quality to retain only high-quality cells. We found that simple filtering procedures such as cutoffs on the total number of reads etc., are insufficient to distinguish good from bad-quality libraries. Therefore, we have implemented a multivariate clustering approach that works on multiple quality metrics (see ?clusterByQuality for details) to increase robustness of the filtering. Here is an example demonstrating the usage of \Rfunction{clusterByQuality}.

+

+\begin{scriptsize}

+<<eval=TRUE, message=FALSE>>==

+## Get some pre-produced results (adjust this to your experiment)

+results <- system.file("extdata", "primary-lung", "hmms", package="AneuFinderData")

+files <- list.files(results, full.names=TRUE)

+

+## Cluster by quality

+cl <- clusterByQuality(files)

+plot(cl$Mclust, what='classification')

+print(cl$parameters)

+## Apparently, the third cluster corresponds to failed libraries

+## while the first cluster contains high-quality libraries

+

+## Select the best cluster and plot it

+selected.files <- unlist(cl$classification[[1]])

+heatmapGenomewide(selected.files)

+@

+\end{scriptsize}

 \subsection{Karyotype measures}

+This package implements two measures to quantify karyotype heterogeneity, an \emph{aneuploidy} and a \emph{heterogeneity} score. Both measures are independent of the number of cells, the length of the genome, and take into account every position in the genome. The following example compares the heterogeneity and aneuploidy between a primary lung cancer and the corresponding liver metastasis.

+\begin{scriptsize}

+<<eval=TRUE, message=FALSE>>==

+## Get some pre-produced results (adjust this to your experiment)

+results <- system.file("extdata", "primary-lung", "hmms", package="AneuFinderData")

+files.lung <- list.files(results, full.names=TRUE)

+results <- system.file("extdata", "metastasis-liver", "hmms", package="AneuFinderData")

+files.liver <- list.files(results, full.names=TRUE)

+

+## Get karyotype measures

+k.lung <- karyotypeMeasures(files.lung)

+k.liver <- karyotypeMeasures(files.liver)

+

+## Print the scores in one data.frame

+df <- rbind(lung = k.lung$genomewide, liver = k.liver$genomewide)

+print(df)

+## While the aneuploidy is similar between both cancers, the heterogeneity is

+## nearly twice as high for the primary lung cancer.

+@

+\end{scriptsize}

 \section{Session Info}

 \begin{scriptsize}

-<<>>=

+<<eval=TRUE>>=

 sessionInfo()

 warnings()

 @

new file mode 100644

Binary files /dev/null and b/vignettes/AneuFinder.pdf differ