@@ -255,6 +255,8 @@ importFrom(SingleCellExperiment,SingleCellExperiment)

 importFrom(SingleCellExperiment,altExp)

 importFrom(SingleCellExperiment,counts)

 importFrom(SingleCellExperiment,reducedDim)

+importFrom(SingleCellExperiment,reducedDimNames)

+importFrom(SingleCellExperiment,reducedDims)

 importFrom(SingleCellExperiment,rowSubset)

 importFrom(SummarizedExperiment,"assay<-")

 importFrom(SummarizedExperiment,"assayNames<-")

@@ -262,6 +264,7 @@ importFrom(SummarizedExperiment,"assays<-")

 importFrom(SummarizedExperiment,"colData<-")

 importFrom(SummarizedExperiment,"rowData<-")

 importFrom(SummarizedExperiment,assay)

+importFrom(SummarizedExperiment,assayNames)

 importFrom(SummarizedExperiment,assays)

 importFrom(SummarizedExperiment,colData)

 importFrom(SummarizedExperiment,rowData)

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

 .constructSCE <- function(

-  matrices,

-  features,

-  barcodes,

-  metadata,

-  reducedDims) {

-

+    matrices,

+    features,

+    barcodes,

+    metadata,

+    reducedDims) {

+  

   sce <- SingleCellExperiment::SingleCellExperiment(assays = matrices)

   SummarizedExperiment::rowData(sce) <- S4Vectors::DataFrame(features)

   SummarizedExperiment::colData(sce) <- S4Vectors::DataFrame(barcodes)

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

   } else {

     fill <- NA

   }

-

+  

   ### combine row

   if (isTRUE(combineRow) & (!is.null(row))) {

     missRow <- row[!row %in% rownames(x)]

@@ -43,14 +43,14 @@

     if (!isTRUE(sparse)) {

       missMat <- as.matrix(missMat)

     }

-

+    

     mat <- rbind(matOrigin, missMat)

     if (anyDuplicated(rownames(mat))) {

       mat <- mat[!duplicated(rownames(mat)), ]

     }

     matOrigin <- mat[row, ]

   }

-

+  

   ### combine cols

   if (isTRUE(combineCol) & (!is.null(col))) {

     missCol <- col[!col %in% colnames(x)]

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

     if (!isTRUE(sparse)) {

       missMat <- as.matrix(missMat)

     }

-

+    

     mat <- cbind(matOrigin, missMat)

     if (anyDuplicated(colnames(mat))) {

       mat <- mat[, !duplicated(colnames(mat))]

@@ -76,12 +76,12 @@

     rw[['rownames']] <- rownames(rw)

     return(rw)

   })

-

+  

   ## Get merged rowData

   by.r <- unique(c('rownames', by.r))

   unionFe <- Reduce(function(r1, r2) merge(r1, r2, by=by.r, all=TRUE), feList)

   allGenes <- unique(unlist(lapply(feList, rownames)))

-

+  

   ## rowData

   newFe <- unionFe

   if (nrow(newFe) != length(allGenes)) {

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

     cD[['rownames']] <- rownames(cD)

     return(cD)

   })

-

+  

   by.c <- unique(c("rownames", by.c))

   unionCb <- Reduce(function(c1, c2) merge(c1, c2, by=by.c, all=TRUE), cbList)

   rownames(unionCb) <- unionCb[['rownames']]

@@ -118,19 +118,19 @@

   reduceList <- lapply(sceList, SingleCellExperiment::reducedDims)

   ## get every reducedDim exists in at least one SCEs

   UnionReducedDims <- unique(unlist(lapply(sceList, SingleCellExperiment::reducedDimNames)))

-

+  

   ## for each reducedDim, get union row/cols

   reducedDims <- list()

   for (reduceDim in UnionReducedDims) {

     x <- lapply(sceList, function(x) {if (reduceDim %in% SingleCellExperiment::reducedDimNames(x)) {SingleCellExperiment::reducedDim(x, reduceDim)}})

     reducedDims[[reduceDim]] <- .getDimUnion(x)

   }

-

+  

   ## Merge reducedDim for each SCE

   redList <- list()

   for (idx in seq_along(sceList)){

     redMat <- reduceList[[idx]]

-

+    

     for (DimName in UnionReducedDims) {

       if (DimName %in% names(redMat)) {

         redMat[[DimName]] <- .getMatUnion(reducedDims[[DimName]], redMat[[DimName]],

@@ -142,10 +142,10 @@

                                           dimnames = list(colnames(sceList[[idx]]), reducedDims[[DimName]][[2]]))

       }

     }

-

+    

     redList[[idx]] <- redMat

   }

-

+  

   return(redList)

 }

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

     unique(unlist(lapply(sceList, rownames))),

     unique(unlist(lapply(sceList, colnames)))

   )

-

+  

   asList <- list()

   for (idx in seq_along(assayList)){

     assay <- assayList[[idx]]

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

     }

     asList[[idx]] <- assay

   }

-

+  

   return(asList)

 }

@@ -197,27 +197,48 @@

 # }

 .mergeMetaSCE <- function(SCE_list) {

+  # Merge highest level metadata entries (except "sctk") by sample names in

+  # SCE_list. For analysis results in "sctk", merge by SCE_list sample name if 

+  # given "all_cells" in one object, else, use existing sample names. 

+  samples <- names(SCE_list)

   sampleMeta <- lapply(SCE_list, S4Vectors::metadata)

   metaNames <- unique(unlist(lapply(sampleMeta, names)))

+  sampleSctkMeta <- lapply(SCE_list, function(x) {S4Vectors::metadata(x)$sctk})

+  sctkMetaNames <- unique(unlist(lapply(sampleMeta, 

+                                        function(x) names(x$sctk))))

+  sampleMeta$sctk <- NULL

+  metaNames <- metaNames[!metaNames %in% c("sctk")]

   NewMeta <- list()

-

+  for (meta in sctkMetaNames) {

+    for (i in seq_along(sampleSctkMeta)) {

+      # `i` is for identifying each SCE object, usually matching a sample in 

+      # the pipeline. `sampleAvail` for samples stored in metadata entry, in

+      # case users are merging merged objects. 

+      sampleAvail <- names(sampleSctkMeta[[i]][[meta]])

+      if (length(sampleAvail) == 1) {

+        NewMeta$sctk[[meta]][[samples[i]]] <- sampleSctkMeta[[i]][[meta]][[1]]

+      } else if (length(sampleAvail) > 1) {

+        names(sampleSctkMeta[[i]][[meta]]) <- 

+          paste(names(sampleSctkMeta)[i], 

+                names(sampleSctkMeta[[i]][[meta]]), sep = "_")

+        NewMeta$sctk[[meta]] <- c(NewMeta$sctk[[meta]], 

+                                  sampleSctkMeta[[i]][[meta]])

+      }

+    }

+  }

   for (meta in metaNames) {

     for (i in seq_along(sampleMeta)) {

-      NewMeta[[meta]][[i]] <- sampleMeta[[i]][[meta]]

+      NewMeta[[meta]][[samples[i]]] <- sampleMeta[[i]][[meta]]

     }

   }

-

-  if ("runBarcodeRanksMetaOutput" %in% metaNames) {

-    NewMeta[["runBarcodeRanksMetaOutput"]] <- unlist(NewMeta[["runBarcodeRanksMetaOutput"]])

-  }

-

+  

   if ("assayType" %in% metaNames) {

     assayType <- lapply(SCE_list, function(x){S4Vectors::metadata(x)$assayType})

     assayType <- BiocGenerics::Reduce(dplyr::union, assayType)

     NewMeta[["assayType"]] <- assayType

   }

-

+  

   return(NewMeta)

 }

@@ -241,7 +262,7 @@

 #' @export

 combineSCE <- function(sceList, by.r = NULL, by.c = NULL, combined = TRUE){

-  if(length(sceList) == 1){

+  if (length(sceList) == 1) {

     return(sceList[[1]])

   }

   ##  rowData

@@ -252,16 +273,21 @@ combineSCE <- function(sceList, by.r = NULL, by.c = NULL, combined = TRUE){

   redMatList <- .mergeRedimSCE(sceList)

   ## assay

   assayList <- .mergeAssaySCE(sceList)

-

+  samples <- names(sceList)

+  if (is.null(samples)) {

+    samples <- paste0("sample", seq_along(sceList))

+  }

   New_SCE <- list()

   for (i in seq(length(sceList))) {

     ## create new sce

-    New_SCE[[i]] <- .constructSCE(matrices = assayList[[i]], features = newFeList,

-                                  barcodes = newCbList[[i]],

-                                  metadata = S4Vectors::metadata(sceList[[i]]),

-                                  reducedDims = redMatList[[i]])

+    sampleName <- samples[i]

+    New_SCE[[sampleName]] <- .constructSCE(matrices = assayList[[i]], 

+                                           features = newFeList,

+                                           barcodes = newCbList[[i]],

+                                           metadata = S4Vectors::metadata(sceList[[i]]),

+                                           reducedDims = redMatList[[i]])

   }

-

+  

   if (isTRUE(combined)) {

     sce <- do.call(SingleCellExperiment::cbind, New_SCE)

     meta <- .mergeMetaSCE(New_SCE)

@@ -2,251 +2,255 @@

 ## (https://github.com/chris-mcginnis-ucsf/DoubletFinder)

 .parallel_paramSweep <- function(n, n.real.cells, real.cells, pK, pN, data,

-                                   orig.commands, PCs, sct, verbose, seed) {

-    sweep.res.list <- list()

-    list.ind <- 0

-

-    ## Make merged real-artifical data

+                                 orig.commands, PCs, sct, verbose, seed) {

+  sweep.res.list <- list()

+  list.ind <- 0

+  

+  ## Make merged real-artifical data

+  if (verbose) {

+    message(date(), " ...   Creating artificial doublets for pN = ",

+            pN[n] * 100, "%")

+  }

+  n_doublets <- ceiling(n.real.cells / (1 - pN[n]) - n.real.cells)

+  real.cells1 <- sample(real.cells, n_doublets, replace = TRUE)

+  real.cells2 <- sample(real.cells, n_doublets, replace = TRUE)

+  doublets <- (data[, real.cells1] + data[, real.cells2]) / 2

+  colnames(doublets) <- paste("X", seq(n_doublets), sep = "")

+  data_wdoublets <- cbind(data, doublets)

+  

+  ## Pre-process Seurat object

+  if (sct == FALSE) {

     if (verbose) {

-        message(date(), " ...   Creating artificial doublets for pN = ",

-                    pN[n] * 100, "%")

-    }

-    n_doublets <- ceiling(n.real.cells / (1 - pN[n]) - n.real.cells)

-    real.cells1 <- sample(real.cells, n_doublets, replace = TRUE)

-    real.cells2 <- sample(real.cells, n_doublets, replace = TRUE)

-    doublets <- (data[, real.cells1] + data[, real.cells2]) / 2

-    colnames(doublets) <- paste("X", seq(n_doublets), sep = "")

-    data_wdoublets <- cbind(data, doublets)

-

-    ## Pre-process Seurat object

-    if (sct == FALSE) {

-        if (verbose) {

-          message(date(), " ...   Creating Seurat object")

-        }

-        seu_wdoublets <- Seurat::CreateSeuratObject(counts = data_wdoublets)

-

-        if (verbose) {

-          message(date(), " ...   Normalizing Seurat object")

-        }

-        seu_wdoublets <- Seurat::NormalizeData(seu_wdoublets,

-                                       normalization.method = orig.commands$NormalizeData.RNA@params$normalization.method,

-                                       scale.factor = orig.commands$NormalizeData.RNA@params$scale.factor,

-                                       margin = orig.commands$NormalizeData.RNA@params$margin,

-                                       verbose = verbose

-        )

-

-        if (verbose) {

-          message(date(), " ...   Finding variable genes")

-        }

-        seu_wdoublets <- Seurat::FindVariableFeatures(seu_wdoublets,

-                                              selection.method = orig.commands$FindVariableFeatures.RNA$selection.method,

-                                              loess.span = orig.commands$FindVariableFeatures.RNA$loess.span,

-                                              clip.max = orig.commands$FindVariableFeatures.RNA$clip.max,

-                                              mean.function = orig.commands$FindVariableFeatures.RNA$mean.function,

-                                              dispersion.function = orig.commands$FindVariableFeatures.RNA$dispersion.function,

-                                              num.bin = orig.commands$FindVariableFeatures.RNA$num.bin,

-                                              binning.method = orig.commands$FindVariableFeatures.RNA$binning.method,

-                                              nfeatures = orig.commands$FindVariableFeatures.RNA$nfeatures,

-                                              mean.cutoff = orig.commands$FindVariableFeatures.RNA$mean.cutoff,

-                                              dispersion.cutoff = orig.commands$FindVariableFeatures.RNA$dispersion.cutoff,

-                                              verbose = verbose

-        )

-

-        if (verbose) {

-          message(date(), " ...   Scaling data")

-        }

-        seu_wdoublets <- Seurat::ScaleData(

-            seu_wdoublets,

-            features = orig.commands$ScaleData.RNA$features,

-            model.use = orig.commands$ScaleData.RNA$model.use,

-            do.scale = orig.commands$ScaleData.RNA$do.scale,

-            do.center = orig.commands$ScaleData.RNA$do.center,

-            scale.max = orig.commands$ScaleData.RNA$scale.max,

-            block.size = orig.commands$ScaleData.RNA$block.size,

-            min.cells.to.block = orig.commands$ScaleData.RNA$min.cells.to.block,

-            verbose = verbose

-        )

-

-        if (verbose) {

-          message(date(), " ...   Running PCA")

-        }

-        seu_wdoublets <- Seurat::RunPCA(

-            seu_wdoublets,

-            features = orig.commands$ScaleData.RNA$features,

-            npcs = length(PCs),

-            rev.pca = orig.commands$RunPCA.RNA$rev.pca,

-            weight.by.var = orig.commands$RunPCA.RNA$weight.by.var,

-            seed.use = seed,

-            verbose = FALSE

-        )

-    }

-

-    if (sct == TRUE) {

-        if (verbose) {

-            message(date(), " ...   Creating Seurat object")

-        }

-        seu_wdoublets <- Seurat::CreateSeuratObject(counts = data_wdoublets)

-

-        if (verbose) {

-          message(date(), " ...   Running SCTransform")

-        }

-        seu_wdoublets <- Seurat::SCTransform(seu_wdoublets)

-

-        if (verbose) {

-          message(date(), " ...   Running PCA")

-        }

-        seu_wdoublets <- Seurat::RunPCA(seu_wdoublets,

-                                npcs = length(PCs), seed.use = seed,

-                                verbose = verbose

-        )

+      message(date(), " ...   Creating Seurat object")

     }

-

-    ## Compute PC distance matrix

+    seu_wdoublets <- Seurat::CreateSeuratObject(counts = data_wdoublets)

+    

     if (verbose) {

-      message(date(), " ...   Calculating PC distance matrix")

+      message(date(), " ...   Normalizing Seurat object")

     }

-    nCells <- nrow(seu_wdoublets@meta.data)

-    pca.coord <- seu_wdoublets@reductions$pca@cell.embeddings[, PCs]

-    seu_wdoublets <- NULL

-    gc()

-    dist.mat <- fields::rdist(pca.coord)[, seq(n.real.cells)]

-

-    ## Pre-order PC distance matrix prior to iterating across pK

-    ## for pANN computations

+    seu_wdoublets <- Seurat::NormalizeData(

+      seu_wdoublets,

+      normalization.method = orig.commands$NormalizeData.RNA@params$normalization.method,

+      scale.factor = orig.commands$NormalizeData.RNA@params$scale.factor,

+      margin = orig.commands$NormalizeData.RNA@params$margin,

+      verbose = verbose

+    )

+    

     if (verbose) {

-      message(date(), " ...   Defining neighborhoods")

+      message(date(), " ...   Finding variable genes")

     }

-

-    for (i in seq(n.real.cells)) {

-        dist.mat[, i] <- order(dist.mat[, i])

+    seu_wdoublets <- Seurat::FindVariableFeatures(

+      seu_wdoublets,

+      selection.method = orig.commands$FindVariableFeatures.RNA$selection.method,

+      loess.span = orig.commands$FindVariableFeatures.RNA$loess.span,

+      clip.max = orig.commands$FindVariableFeatures.RNA$clip.max,

+      mean.function = orig.commands$FindVariableFeatures.RNA$mean.function,

+      dispersion.function = orig.commands$FindVariableFeatures.RNA$dispersion.function,

+      num.bin = orig.commands$FindVariableFeatures.RNA$num.bin,

+      binning.method = orig.commands$FindVariableFeatures.RNA$binning.method,

+      nfeatures = orig.commands$FindVariableFeatures.RNA$nfeatures,

+      mean.cutoff = orig.commands$FindVariableFeatures.RNA$mean.cutoff,

+      dispersion.cutoff = orig.commands$FindVariableFeatures.RNA$dispersion.cutoff,

+      verbose = verbose

+    )

+    

+    if (verbose) {

+      message(date(), " ...   Scaling data")

     }

-

-    ## Trim PC distance matrix for faster manipulations

-    ind <- round(nCells * max(pK)) + 5

-    dist.mat <- dist.mat[seq(ind), ]

-

-    ## Compute pANN across pK sweep

-

+    seu_wdoublets <- Seurat::ScaleData(

+      seu_wdoublets,

+      features = orig.commands$ScaleData.RNA$features,

+      model.use = orig.commands$ScaleData.RNA$model.use,

+      do.scale = orig.commands$ScaleData.RNA$do.scale,

+      do.center = orig.commands$ScaleData.RNA$do.center,

+      scale.max = orig.commands$ScaleData.RNA$scale.max,

+      block.size = orig.commands$ScaleData.RNA$block.size,

+      min.cells.to.block = orig.commands$ScaleData.RNA$min.cells.to.block,

+      verbose = verbose

+    )

+    

     if (verbose) {

-      message(date(), " ...   Computing pANN across all pK")

+      message(date(), " ...   Running PCA")

     }

-    for (k in seq_along(pK)) {

-        if (verbose) {

-            message(date(), " ...   pK = ", pK[k])

-        }

-        pk.temp <- round(nCells * pK[k])

-        pANN <- as.data.frame(matrix(0L, nrow = n.real.cells, ncol = 1))

-        colnames(pANN) <- "pANN"

-        rownames(pANN) <- real.cells

-        list.ind <- list.ind + 1

-

-        for (i in seq(n.real.cells)) {

-            neighbors <- dist.mat[2:(pk.temp + 1), i]

-            pANN$pANN[i] <- length(which(neighbors > n.real.cells)) / pk.temp

-        }

-

-        sweep.res.list[[list.ind]] <- pANN

+    seu_wdoublets <- Seurat::RunPCA(

+      seu_wdoublets,

+      features = orig.commands$ScaleData.RNA$features,

+      npcs = length(PCs),

+      rev.pca = orig.commands$RunPCA.RNA$rev.pca,

+      weight.by.var = orig.commands$RunPCA.RNA$weight.by.var,

+      seed.use = seed,

+      verbose = FALSE

+    )

+  }

+  

+  if (sct == TRUE) {

+    if (verbose) {

+      message(date(), " ...   Creating Seurat object")

     }

-

-    return(sweep.res.list)

-}

-

-.paramSweep <- function(seu, PCs = seq(10), sct = FALSE,

-                          verbose = FALSE, num.cores, seed) {

-    ## Set pN-pK param sweep ranges

-    pK <- c(0.0005, 0.001, 0.005, seq(0.01, 0.3, by = 0.01))

-    pN <- seq(0.05, 0.3, by = 0.05)

-    ## Remove pK values with too few cells

-    min.cells <- round(nrow(seu@meta.data) / (1 - 0.05) - nrow(seu@meta.data))

-    pK.test <- round(pK * min.cells)

-    pK <- pK[which(pK.test >= 1)]

-

-    ## Extract pre-processing parameters from original data analysis workflow

-    orig.commands <- seu@commands

-

-    ## Down-sample cells to 10000 (when applicable) for computational effiency

-    if (nrow(seu@meta.data) > 10000) {

-        real.cells <- rownames(seu@meta.data)[sample(seq(nrow(seu@meta.data)),

-                                                     10000, replace = FALSE)]

-        data <- seu@assays$RNA@counts[, real.cells]

-        n.real.cells <- ncol(data)

+    seu_wdoublets <- Seurat::CreateSeuratObject(counts = data_wdoublets)

+    

+    if (verbose) {

+      message(date(), " ...   Running SCTransform")

     }

-

-    if (nrow(seu@meta.data) <= 10000) {

-        real.cells <- rownames(seu@meta.data)

-        data <- seu@assays$RNA@counts

-        n.real.cells <- ncol(data)

+    seu_wdoublets <- Seurat::SCTransform(seu_wdoublets)

+    

+    if (verbose) {

+      message(date(), " ...   Running PCA")

     }

-    ## Iterate through pN, computing pANN vectors at varying pK

-    # no_cores <- detectCores()-1

-    if (is.null(num.cores)) {

-        num.cores <- 1

+    seu_wdoublets <- Seurat::RunPCA(seu_wdoublets,

+                                    npcs = length(PCs), seed.use = seed,

+                                    verbose = verbose

+    )

+  }

+  

+  ## Compute PC distance matrix

+  if (verbose) {

+    message(date(), " ...   Calculating PC distance matrix")

+  }

+  nCells <- ncol(seu_wdoublets)

+  pca.coord <- Seurat::Reductions(seu_wdoublets, "pca")@cell.embeddings[, PCs]

+  #seu_wdoublets <- NULL

+  rm(seu_wdoublets)

+  gc()

+  dist.mat <- fields::rdist(pca.coord)[, seq(n.real.cells)]

+  

+  ## Pre-order PC distance matrix prior to iterating across pK

+  ## for pANN computations

+  if (verbose) {

+    message(date(), " ...   Defining neighborhoods")

+  }

+  

+  for (i in seq(n.real.cells)) {

+    dist.mat[, i] <- order(dist.mat[, i])

+  }

+  

+  ## Trim PC distance matrix for faster manipulations

+  ind <- round(nCells * max(pK)) + 5

+  dist.mat <- dist.mat[seq(ind), ]

+  

+  ## Compute pANN across pK sweep

+  

+  if (verbose) {

+    message(date(), " ...   Computing pANN across all pK")

+  }

+  for (k in seq_along(pK)) {

+    if (verbose) {

+      message(date(), " ...   pK = ", pK[k])

     }

-

-    if (num.cores > 1) {

-        cl <- parallel::makeCluster(num.cores)

-

-        output2 <- withr::with_seed(

-            seed,

-            parallel::mclapply(as.list(seq_along(pN)),

-                     FUN = .parallel_paramSweep,

-                     n.real.cells,

-                     real.cells,

-                     pK,

-                     pN,

-                     data,

-                     orig.commands,

-                     PCs,

-                     sct, mc.cores = num.cores,

-                     verbose = verbose,

-                     seed = seed,

-                     mc.set.seed = FALSE

-            )

-        )

-        parallel::stopCluster(cl)

-    } else {

-        output2 <- lapply(as.list(seq_along(pN)),

-                          FUN = .parallel_paramSweep,

-                          n.real.cells,

-                          real.cells,

-                          pK,

-                          pN,

-                          data,

-                          orig.commands,

-                          PCs,

-                          sct,

-                          seed = seed,

-                          verbose = verbose

-        )

+    pk.temp <- round(nCells * pK[k])

+    pANN <- as.data.frame(matrix(0L, nrow = n.real.cells, ncol = 1))

+    colnames(pANN) <- "pANN"

+    rownames(pANN) <- real.cells

+    list.ind <- list.ind + 1

+    

+    for (i in seq(n.real.cells)) {

+      neighbors <- dist.mat[2:(pk.temp + 1), i]

+      pANN$pANN[i] <- length(which(neighbors > n.real.cells)) / pk.temp

     }

+    

+    sweep.res.list[[list.ind]] <- pANN

+  }

+  

+  return(sweep.res.list)

+}

-    ## Write parallelized output into list

-    sweep.res.list <- list()

-    list.ind <- 0

-    for (i in seq_along(output2)) {

-        for (j in seq_along(output2[[i]])) {

-            list.ind <- list.ind + 1

-            sweep.res.list[[list.ind]] <- output2[[i]][[j]]

-        }

-    }

-    ## Assign names to list of results

-    name.vec <- NULL

-    for (j in seq_along(pN)) {

-        name.vec <- c(name.vec, paste("pN", pN[j], "pK", pK, sep = "_"))

+.paramSweep <- function(seu, PCs = seq(10), sct = FALSE,

+                        verbose = FALSE, num.cores, seed) {

+  ## Set pN-pK param sweep ranges

+  pK <- c(0.0005, 0.001, 0.005, seq(0.01, 0.3, by = 0.01))

+  pN <- seq(0.05, 0.3, by = 0.05)

+  ## Remove pK values with too few cells

+  min.cells <- round(nrow(seu[[]]) / (1 - 0.05) - nrow(seu[[]]))

+  pK.test <- round(pK * min.cells)

+  pK <- pK[which(pK.test >= 1)]

+  

+  ## Extract pre-processing parameters from original data analysis workflow

+  orig.commands <- seu@commands

+  

+  ## Down-sample cells to 10000 (when applicable) for computational effiency

+  if (nrow(seu[[]]) > 10000) {

+    real.cells <- rownames(seu[[]])[sample(seq(nrow(seu[[]])),

+                                           10000, replace = FALSE)]

+    data <- Seurat::GetAssayData(seu, slot = "counts", 

+                                 assay = "RNA")[, real.cells]

+    n.real.cells <- ncol(data)

+  }

+  

+  if (ncol(seu) <= 10000) {

+    real.cells <- colnames(seu)

+    data <- Seurat::GetAssayData(seu, slot = "counts", assay = "RNA")

+    n.real.cells <- ncol(data)

+  }

+  ## Iterate through pN, computing pANN vectors at varying pK

+  # no_cores <- detectCores()-1

+  if (is.null(num.cores)) {

+    num.cores <- 1

+  }

+  

+  if (num.cores > 1) {

+    cl <- parallel::makeCluster(num.cores)

+    

+    output2 <- withr::with_seed(

+      seed,

+      parallel::mclapply(as.list(seq_along(pN)),

+                         FUN = .parallel_paramSweep,

+                         n.real.cells,

+                         real.cells,

+                         pK,

+                         pN,

+                         data,

+                         orig.commands,

+                         PCs,

+                         sct, mc.cores = num.cores,

+                         verbose = verbose,

+                         seed = seed,

+                         mc.set.seed = FALSE

+      )

+    )

+    parallel::stopCluster(cl)

+  } else {

+    output2 <- lapply(as.list(seq_along(pN)),

+                      FUN = .parallel_paramSweep,

+                      n.real.cells,

+                      real.cells,

+                      pK,

+                      pN,

+                      data,

+                      orig.commands,

+                      PCs,

+                      sct,

+                      seed = seed,

+                      verbose = verbose

+    )

+  }

+  

+  ## Write parallelized output into list

+  sweep.res.list <- list()

+  list.ind <- 0

+  for (i in seq_along(output2)) {

+    for (j in seq_along(output2[[i]])) {

+      list.ind <- list.ind + 1

+      sweep.res.list[[list.ind]] <- output2[[i]][[j]]

     }

-    names(sweep.res.list) <- name.vec

-    return(sweep.res.list)

-

+  }

+  ## Assign names to list of results

+  name.vec <- NULL

+  for (j in seq_along(pN)) {

+    name.vec <- c(name.vec, paste("pN", pN[j], "pK", pK, sep = "_"))

+  }

+  names(sweep.res.list) <- name.vec

+  return(sweep.res.list)

+  

 }

 .runDoubletFinder <- function(counts, seuratPcs, seuratRes, formationRate,

                               seuratNfeatures, verbose = FALSE,

                               nCores = NULL, seed = 12345) {

-

+  

   ## Convert to sparse matrix if not already in that format

   counts <- .convertToMatrix(counts)

   colnames(counts) <- gsub("_", "-", colnames(counts))

-

+  

   seurat <- suppressWarnings(Seurat::CreateSeuratObject(

     counts = counts,

     project = "seurat", min.features = 0

@@ -257,51 +261,53 @@

     verbose = verbose

   )

   seurat <- Seurat::FindVariableFeatures(seurat,

-    selection.method = "vst",

-    nfeatures = seuratNfeatures,

-    verbose = verbose

+                                         selection.method = "vst",

+                                         nfeatures = seuratNfeatures,

+                                         verbose = verbose

   )

-

+  

   allGenes <- rownames(seurat)

   seurat <- Seurat::ScaleData(seurat, features = allGenes, verbose = verbose)

-

-  numPc <- min(ncol(seurat@assays$RNA@scale.data) - 1, 50)

+  

+  numPc <- min(nrow(Seurat::GetAssayData(seurat, slot = "scale.data", 

+                                         assay = "RNA")) - 1, 

+               50)

   seurat <- Seurat::RunPCA(seurat,

-    features =

-      Seurat::VariableFeatures(object = seurat),

-    npcs = numPc, verbose = verbose, seed.use = seed

+                           features =

+                             Seurat::VariableFeatures(object = seurat),

+                           npcs = numPc, verbose = verbose, seed.use = seed

   )

   seurat <- Seurat::FindNeighbors(seurat, dims = seuratPcs, verbose = verbose)

   seurat <- Seurat::FindClusters(seurat,

-    resolution = seuratRes,

-    verbose = verbose,

-    random.seed = seed

+                                 resolution = seuratRes,

+                                 verbose = verbose,

+                                 random.seed = seed

   )

   invisible(sweepResListSeurat <- .paramSweep(seurat,

-    PCs = seuratPcs,

-    sct = FALSE,

-    num.cores = nCores,

-    verbose = verbose,

-    seed = seed

+                                              PCs = seuratPcs,

+                                              sct = FALSE,

+                                              num.cores = nCores,

+                                              verbose = verbose,

+                                              seed = seed

   ))

   sweepStatsSeurat <- .summarizeSweep(sweepResListSeurat,

-    GT = FALSE

+                                      GT = FALSE

   )

   bcmvnSeurat <- .find.pK(sweepStatsSeurat)

   pkOptimal <- as.numeric(as.matrix(bcmvnSeurat$pK[

     which.max(bcmvnSeurat$MeanBC)

   ]))

-  annotations <- seurat@meta.data$seurat_clusters

-  homotypicProp <- .modelHomotypic(annotations)

-  nExpPoi <- round(formationRate * ncol(seurat@assays$RNA))

+  #annotations <- seurat[[]]$seurat_clusters

+  #homotypicProp <- .modelHomotypic(annotations)

+  nExpPoi <- round(formationRate * ncol(seurat))

   seurat <- .doubletFinder_v3(seurat,

-    PCs = seuratPcs,

-    pN = 0.25,

-    pK = pkOptimal,

-    nExp = nExpPoi,

-    reuse.pANN = FALSE,

-    sct = FALSE,

-    verbose = FALSE

+                              PCs = seuratPcs,

+                              pN = 0.25,

+                              pK = pkOptimal,

+                              nExp = nExpPoi,

+                              reuse.pANN = FALSE,

+                              sct = FALSE,

+                              verbose = FALSE

   )

   names(seurat@meta.data)[6] <- "doubletFinderAnnScore"

   names(seurat@meta.data)[7] <- "doubletFinderLabel"

@@ -311,34 +317,39 @@

 #' @title Generates a doublet score for each cell via doubletFinder

 #' @description Uses doubletFinder to determine cells within the dataset

 #'  suspected to be doublets.

-#' @param inSCE Input SingleCellExperiment object. Must contain a counts matrix

-#' @param useAssay  Indicate which assay to use. Default "counts".

-#' @param sample Numeric vector. Each cell will be assigned a sample number.

-#' @param seed Seed for the random number generator. Default 12345.

+#' @param inSCE inSCE A \linkS4class{SingleCellExperiment} object.

+#' @param sample Character vector or colData variable name. Indicates which 

+#' sample each cell belongs to. Default \code{NULL}.

+#' @param useAssay  A string specifying which assay in the SCE to use. Default 

+#' \code{"counts"}.

+#' @param seed Seed for the random number generator, can be set to \code{NULL}. 

+#' Default \code{12345}.

 #' @param seuratNfeatures Integer. Number of highly variable genes to use.

-#'  Default 2000.

+#' Default \code{2000}.

 #' @param seuratPcs Numeric vector. The PCs used in seurat function to

-#'   determine number of clusters. Default 1:15.

-#' @param seuratRes Numeric vector. The resolution parameter used in seurat,

-#'  which adjusts the number of clusters determined via the algorithm.

-#'  Default 1.5.

-#' @param formationRate Doublet formation rate used within algorithm.

-#'  Default 0.075.

-#' @param nCores Number of cores used for running the function.

-#' @param verbose Boolean. Wheter to print messages from Seurat and

-#'  DoubletFinder. Default FALSE.

-#' @return SingleCellExperiment object containing the

-#'  'doublet_finder_doublet_score'.

+#' determine number of clusters. Default \code{1:15}.

+#' @param seuratRes Numeric vector. The resolution parameter used in Seurat,

+#' which adjusts the number of clusters determined via the algorithm. Default 

+#' \code{1.5}.

+#' @param formationRate Doublet formation rate used within algorithm. Default 

+#' \code{0.075}.

+#' @param nCores Number of cores used for running the function. Default 

+#' \code{NULL}.

+#' @param verbose Boolean. Wheter to print messages from Seurat and 

+#' DoubletFinder. Default \code{FALSE}.

+#' @return \linkS4class{SingleCellExperiment} object containing the

+#' \code{doublet_finder_doublet_score} variable in \code{colData} slot.

+#' @seealso \code{\link{runCellQC}}, \code{\link{plotDoubletFinderResults}}

 #' @examples

 #' data(scExample, package = "singleCellTK")

 #' sce <- subsetSCECols(sce, colData = "type != 'EmptyDroplet'")

 #' sce <- runDoubletFinder(sce)

 #' @export

-#' @importFrom SummarizedExperiment colData colData<-

+#' @importFrom SummarizedExperiment colData colData<- assayNames assayNames<-

 #' @importFrom SingleCellExperiment reducedDim<-

 runDoubletFinder <- function(inSCE,

-                             useAssay = "counts",

                              sample = NULL,

+                             useAssay = "counts",

                              seed = 12345,

                              seuratNfeatures = 2000,

                              seuratPcs = seq(15),

@@ -346,28 +357,25 @@ runDoubletFinder <- function(inSCE,

                              formationRate = 0.075,

                              nCores = NULL,

                              verbose = FALSE) {

-

-  # argsList <- as.list(formals(fun = sys.function(sys.parent()), envir = parent.frame()))

+  

   argsList <- mget(names(formals()), sys.frame(sys.nframe()))

-

+  argsList <- argsList[!names(argsList) %in% c("inSCE")]

+  argsList$packageVersion <- "2.0.2"

   tempSCE <- inSCE

-  SummarizedExperiment::assayNames(inSCE)[which(useAssay %in% SummarizedExperiment::assayNames(inSCE))] <- "counts"

-  useAssay <- "counts"

-