@@ -89,13 +89,12 @@ Chromstar <- function(inputfolder, experiment.table, outputfolder, configfile=NU

     ## Prepare IDs and filenames

     IDs <- paste0(exp.table$mark, '-', exp.table$condition, '-rep', exp.table$replicate)

     datafiles <- file.path(inputfolder, basename(as.character(exp.table$file)))

-    names(datafiles) <- basename(datafiles)

-    names(IDs) <- basename(datafiles)

-    rownames(exp.table) <- basename(datafiles)

-    filenames <- paste0(IDs, '_binsize', binsize.string, '_stepsize', stepsize.string, '.RData')

-    names(filenames) <- basename(datafiles)

-    filenames.stepsize <- paste0(IDs, '_binsize', stepsize.string, '_stepsize', stepsize.string, '.RData') # filenames for saving counts at stepsize

-    names(filenames.stepsize) <- basename(datafiles)

+    names(datafiles) <- IDs

+    rownames(exp.table) <- IDs

+    filenames <- paste0(IDs, '_', as.character(exp.table$file), '_binsize', binsize.string, '_stepsize', stepsize.string, '.RData')

+    names(filenames) <- IDs

+    filenames.stepsize <- paste0(IDs, '_', as.character(exp.table$file), '_binsize', stepsize.string, '_stepsize', stepsize.string, '.RData') # filenames for saving counts at stepsize

+    names(filenames.stepsize) <- IDs

     ## Inputfiles

     inputfiles <- unique(unlist(strsplit(as.character(exp.table$controlFiles), '\\|')))

     inputfiles <- inputfiles[!is.na(inputfiles)]

@@ -245,19 +244,21 @@ Chromstar <- function(inputfolder, experiment.table, outputfolder, configfile=NU

     ### Count reads in bins ###

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

-    parallel.helper <- function(file, input) {

+    parallel.helper <- function(ID, input, inputfile=NULL) {

         if (!input) {

-            savename <- file.path(binpath, filenames[basename(file)])

-            savename.stepsize <- file.path(binpath, filenames.stepsize[basename(file)])

-            pairedEndReads <- exp.table[basename(file),'pairedEndReads']

+            savename <- file.path(binpath, filenames[ID])

+            savename.stepsize <- file.path(binpath, filenames.stepsize[ID])

+            pairedEndReads <- exp.table[ID,'pairedEndReads']

+            file <- datafiles[ID]

         } else {

-            savename <- file.path(binpath, inputfilenames[basename(file)])

-            savename.stepsize <- file.path(binpath, inputfilenames.stepsize[basename(file)])

-            pairedEndReads <- exp.table[grep(basename(file), exp.table$controlFiles),'pairedEndReads']

+            savename <- file.path(binpath, inputfilenames[basename(inputfile)])

+            savename.stepsize <- file.path(binpath, inputfilenames.stepsize[basename(inputfile)])

+            pairedEndReads <- exp.table[grep(basename(inputfile), exp.table$controlFiles),'pairedEndReads']

             if (any(pairedEndReads != pairedEndReads[1])) {

                 stop("Multiple definitions of 'pairedEndReads' for file ", file, ".")

             }

             pairedEndReads <- pairedEndReads[1]

+            file <- inputfile

         }

         if (!file.exists(savename)) {

             tC <- tryCatch({

@@ -265,7 +266,7 @@ Chromstar <- function(inputfolder, experiment.table, outputfolder, configfile=NU

                 if (!input) {

                     exp.table.input <- exp.table

                 }

-                binlist <- binReads(file=file, experiment.table=exp.table.input, assembly=chrom.lengths.df, pairedEndReads=pairedEndReads, binsizes=NULL, reads.per.bin=NULL, bins=list(pre.bins, pre.bins.stepsize), stepsizes=c(stepsize, stepsize), chromosomes=conf[['chromosomes']], remove.duplicate.reads=conf[['remove.duplicate.reads']], min.mapq=conf[['min.mapq']])

+                binlist <- binReads(file=file, experiment.table=exp.table.input, ID=ID, assembly=chrom.lengths.df, pairedEndReads=pairedEndReads, binsizes=NULL, reads.per.bin=NULL, bins=list(pre.bins, pre.bins.stepsize), stepsizes=c(stepsize, stepsize), chromosomes=conf[['chromosomes']], remove.duplicate.reads=conf[['remove.duplicate.reads']], min.mapq=conf[['min.mapq']])

                 ptm <- startTimedMessage("Saving to file ", savename, " ...")

                 bins <- binlist[[1]]

                 save(bins, file=savename)

@@ -280,33 +281,33 @@ Chromstar <- function(inputfolder, experiment.table, outputfolder, configfile=NU

         }

     }

     ## Unparallelized

-#     for (file in datafiles) {

-#         parallel.helper(file, input=FALSE)

+#     for (ID in IDs) {

+#         parallel.helper(ID, input=FALSE)

 #     }

 #     for (file in inputfiles) {

-#         parallel.helper(file, input=TRUE)

+#         parallel.helper(ID=NULL, input=TRUE, inputfile=file)

 #     }

     ## Parallelized

     if (numcpu > 1) {

         ptm <- startTimedMessage("Binning data ...")

-        temp <- foreach (file = datafiles, .packages=c("chromstaR")) %dopar% {

-            parallel.helper(file, input=FALSE)

+        temp <- foreach (ID = IDs, .packages=c("chromstaR")) %dopar% {

+            parallel.helper(ID, input=FALSE)

         }

         stopTimedMessage(ptm)

     } else {

-        for (file in datafiles) {

-            parallel.helper(file, input=FALSE)

+        for (ID in IDs) {

+            parallel.helper(ID, input=FALSE)

         }

     }

     if (numcpu > 1) {

         ptm <- startTimedMessage("Binning input ...")

         temp <- foreach (file = inputfiles, .packages=c("chromstaR")) %dopar% {

-            parallel.helper(file, input=TRUE)

+            parallel.helper(ID=NULL, input=TRUE, inputfile=file)

         }

         stopTimedMessage(ptm)

     } else {

         for (file in inputfiles) {

-            parallel.helper(file, input=TRUE)

+            parallel.helper(ID=NULL, input=TRUE, inputfile=file)

         }

     }

@@ -328,8 +329,6 @@ Chromstar <- function(inputfolder, experiment.table, outputfolder, configfile=NU

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

         if (!file.exists(savename)) {

             tC <- tryCatch({

-                mark <- strsplit(basename(file), '-')[[1]][1]

-                condition <- strsplit(basename(file), '-')[[1]][2]

                 input.files <- NULL

                 if (length(inputfiles)>0) {

                     input.files <- inputfiles

@@ -356,16 +355,16 @@ Chromstar <- function(inputfolder, experiment.table, outputfolder, configfile=NU

     }

     if (numcpu > 1) {

         ptm <- startTimedMessage("Univariate peak calling ...")

-        temp <- foreach (filename = names(files), .packages=c("chromstaR")) %dopar% {

-            file <- files[filename]

-            inputfiles <- inputfiles.list[[filename]]

+        temp <- foreach (ID = names(files), .packages=c("chromstaR")) %dopar% {

+            file <- files[ID]

+            inputfiles <- inputfiles.list[[ID]]

             parallel.helper(file, inputfiles)

         }

         stopTimedMessage(ptm)

     } else {

-        for (filename in names(files)) {

-            file <- files[filename]

-            inputfiles <- inputfiles.list[[filename]]

+        for (ID in names(files)) {

+            file <- files[ID]

+            inputfiles <- inputfiles.list[[ID]]

             parallel.helper(file, inputfiles)

         }

     }

@@ -9,6 +9,7 @@

 #' @aliases binning

 #' @param file A file with aligned reads. Alternatively a \code{\link{GRanges}} with aligned reads.

 #' @param experiment.table An \code{\link{experiment.table}} containing the supplied \code{file}. This is necessary to uniquely identify the file in later steps of the workflow. Set to \code{NULL} if you don't have it (not recommended).

+#' @param ID Optional ID to select a row from the \code{experiment.table}. Only necessary if the experiment table contains the same file in multiple positions in column 'file'.

 #' @inheritParams readBamFileAsGRanges

 #' @inheritParams readBedFileAsGRanges

 #' @param binsizes An integer vector specifying the bin sizes to use.

@@ -37,7 +38,7 @@

 #'                   stepsizes=500, chromosomes='chr12')

 #'print(binned)

 #'

-binReads <- function(file, experiment.table=NULL, assembly, bamindex=file, chromosomes=NULL, pairedEndReads=FALSE, min.mapq=10, remove.duplicate.reads=TRUE, max.fragment.width=1000, blacklist=NULL, binsizes=1000, stepsizes=binsizes/2, reads.per.bin=NULL, bins=NULL, variable.width.reference=NULL, use.bamsignals=TRUE, format=NULL) {

+binReads <- function(file, experiment.table=NULL, ID=NULL, assembly, bamindex=file, chromosomes=NULL, pairedEndReads=FALSE, min.mapq=10, remove.duplicate.reads=TRUE, max.fragment.width=1000, blacklist=NULL, binsizes=1000, stepsizes=binsizes/2, reads.per.bin=NULL, bins=NULL, variable.width.reference=NULL, use.bamsignals=TRUE, format=NULL) {

     ## Determine format

     if (is.null(format)) {

@@ -78,6 +79,12 @@ binReads <- function(file, experiment.table=NULL, assembly, bamindex=file, chrom

     if (!is.null(experiment.table)) {

         check.experiment.table(experiment.table)

         info <- experiment.table[basename(as.character(experiment.table$file))==basename(file),]

+        if (nrow(info) > 1) {

+            if (is.null(ID)) {

+                stop("Please specify argument 'ID' if multiple files in 'experiment.table' are identical.")

+            }

+            info <- info[ID, ]

+        }

         ID <- paste0(info$mark, '-', info$condition, '-rep', info$rep)

         info$ID <- ID

         if (pairedEndReads != info$pairedEndReads) {

@@ -331,7 +331,7 @@ runMultivariate <- function(binned.data, stepbins, info, comb.states, use.states

             get.posteriors = as.logical(get.posteriors), # bool* keep_posteriors

             densities = double(length=lenDensities), # double* densities

             keep.densities = as.logical(keep.densities), # bool* keep_densities

-            states = integer(length=length(binned.data)), # int* states

+            states = double(length=length(binned.data)), # double* states

             maxPosterior = double(length=length(binned.data)), # double* maxPosterior

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

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

@@ -292,7 +292,7 @@ callPeaksUnivariateAllChr <- function(binned.data, input.data=NULL, eps=0.01, in

                 read.counts.to.remove <- max(c(read.counts.to.remove.1, 2*minlow))

                 index.filtered <- which(counts>0 & counts<=read.counts.to.remove)

                 counts[index.filtered] <- 0

-                if (length(index.filtered)>0) {

+                if (length(index.filtered)>0 & ioffset == 1) {

                     message(paste0("Replaced read counts <= ",read.counts.to.remove," by 0. This was done because the selected bin size is considered too big for this dataset: The mean of the read counts (zeros removed) is bigger than the specified max.mean = ",max.mean,". Check the fits!"))

                 }

             }

@@ -412,6 +412,7 @@ callPeaksUnivariateAllChr <- function(binned.data, input.data=NULL, eps=0.01, in

         if (ioffset == 1) {

             ### Make return object ###

                 result <- list()

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

                 result$info <- attr(binned.data, 'info')

             ## Parameters

                 # Weights

@@ -444,8 +445,6 @@ callPeaksUnivariateAllChr <- function(binned.data, input.data=NULL, eps=0.01, in

             ## Convergence info

                 convergenceInfo <- list(eps=eps, loglik=hmm$loglik, loglik.delta=hmm$loglik.delta, num.iterations=hmm$num.iterations, time.sec=hmm$time.sec, max.mean=max.mean, read.cutoff=max(hmm$counts))

                 result$convergenceInfo <- convergenceInfo

-            ## Add class

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

         }

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

@@ -46,4 +46,4 @@ Please refer to the [vignette](https://github.com/ataudt/chromstaR/blob/master/v

 Report Errors

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

-If you encounter errors of any kind, please file an [issue here](https://github.com/ataudt/chromstaR/issues/new). I will try to react within one day.

+If you encounter errors of any kind, please file an [issue here](https://github.com/ataudt/chromstaR/issues/new). I will try to react within three working days.

@@ -5,9 +5,9 @@

 \alias{binning}

 \title{Convert aligned reads from various file formats into read counts in equidistant bins}

 \usage{

-binReads(file, experiment.table = NULL, assembly, bamindex = file,

-  chromosomes = NULL, pairedEndReads = FALSE, min.mapq = 10,

-  remove.duplicate.reads = TRUE, max.fragment.width = 1000,

+binReads(file, experiment.table = NULL, ID = NULL, assembly,

+  bamindex = file, chromosomes = NULL, pairedEndReads = FALSE,

+  min.mapq = 10, remove.duplicate.reads = TRUE, max.fragment.width = 1000,

   blacklist = NULL, binsizes = 1000, stepsizes = binsizes/2,

   reads.per.bin = NULL, bins = NULL, variable.width.reference = NULL,

   use.bamsignals = TRUE, format = NULL)

@@ -17,6 +17,8 @@ binReads(file, experiment.table = NULL, assembly, bamindex = file,

 \item{experiment.table}{An \code{\link{experiment.table}} containing the supplied \code{file}. This is necessary to uniquely identify the file in later steps of the workflow. Set to \code{NULL} if you don't have it (not recommended).}

+\item{ID}{Optional ID to select a row from the \code{experiment.table}. Only necessary if the experiment table contains the same file in multiple positions in column 'file'.}

+

 \item{assembly}{Please see \code{\link[GenomeInfoDb]{fetchExtendedChromInfoFromUCSC}} for available assemblies. Only necessary when importing BED files. BAM files are handled automatically. Alternatively a data.frame with columns 'chromosome' and 'length'.}

 \item{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.}

@@ -273,7 +273,7 @@ void univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* max

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

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

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

-void multivariate_hmm(int* O, int* T, int* N, int *Nmod, double* comb_states, double* size, double* prob, double* w, double* cor_matrix_inv, double* det, int* maxiter, int* maxtime, double* eps, double* posteriors, bool* keep_posteriors, double* densities, bool* keep_densities, int* states, double* maxPosterior, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* verbosity)

+void multivariate_hmm(int* O, int* T, int* N, int *Nmod, double* comb_states, double* size, double* prob, double* w, double* cor_matrix_inv, double* det, int* maxiter, int* maxtime, double* eps, double* posteriors, bool* keep_posteriors, double* densities, bool* keep_densities, double* states, double* maxPosterior, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* verbosity)

 {

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

@@ -20,7 +20,7 @@ extern "C"

 void univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* maxiter, int* maxtime, double* eps, double* posteriors, double* densities, bool* keep_densities, int* states, double* maxPosterior, double* A, double* proba, double* loglik, double* weights, int* iniproc, double* initial_size, double* initial_prob, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* read_cutoff, int* verbosity);

 extern "C"

-void multivariate_hmm(int* O, int* T, int* N, int *Nmod, double* comb_states, double* size, double* prob, double* w, double* cor_matrix_inv, double* det, int* maxiter, int* maxtime, double* eps, double* posteriors, bool* keep_posteriors, double* densities, bool* keep_densities, int* states, double* maxPosterior, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* verbosity);

+void multivariate_hmm(int* O, int* T, int* N, int *Nmod, double* comb_states, double* size, double* prob, double* w, double* cor_matrix_inv, double* det, int* maxiter, int* maxtime, double* eps, double* posteriors, bool* keep_posteriors, double* densities, bool* keep_densities, double* states, double* maxPosterior, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* verbosity);

 extern "C"

 void univariate_cleanup();

@@ -32,4 +32,4 @@ extern "C"

 void array3D_which_max(double* array3D, int* dim, int* ind_max);

 extern "C"

-void array2D_which_max(double* array3D, int* dim, int* ind_max, double* value_max);

\ No newline at end of file

+void array2D_which_max(double* array2D, int* dim, int* ind_max, double* value_max);

\ No newline at end of file

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

 R_NativePrimitiveArgType arg1[] = {INTSXP, INTSXP, INTSXP, REALSXP, REALSXP, INTSXP, INTSXP, REALSXP, REALSXP, REALSXP, LGLSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, REALSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, LGLSXP, INTSXP, INTSXP, INTSXP, INTSXP};

-R_NativePrimitiveArgType arg2[] = {INTSXP, INTSXP, INTSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, REALSXP, REALSXP, INTSXP, INTSXP, REALSXP, REALSXP, LGLSXP, REALSXP, LGLSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, REALSXP, REALSXP, LGLSXP, INTSXP, INTSXP, INTSXP};

+R_NativePrimitiveArgType arg2[] = {INTSXP, INTSXP, INTSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, REALSXP, REALSXP, INTSXP, INTSXP, REALSXP, REALSXP, LGLSXP, REALSXP, LGLSXP, REALSXP, REALSXP, REALSXP, REALSXP, REALSXP, REALSXP, REALSXP, LGLSXP, INTSXP, INTSXP, INTSXP};

 R_NativePrimitiveArgType arg4[] = {INTSXP};

 R_NativePrimitiveArgType arg5[] = {REALSXP, INTSXP, INTSXP};

 R_NativePrimitiveArgType arg6[] = {REALSXP, INTSXP, INTSXP, REALSXP};