new file mode 100644

@@ -0,0 +1,9 @@

+.Rproj.user

+.Rhistory

+.RData

+.Ruserdata

+singleCellTK.Rproj

+*.Rproj

+.DS_Store

+.Rprofile

+renv/

@@ -15,7 +15,7 @@ Authors@R: c(person(given="David", family="Jenkins", email="dfj@bu.edu", role=c(

                     comment = c(ORCID = "0000-0002-6247-6595")))

 Maintainer: David Jenkins <dfj@bu.edu>

 Depends:

-    R (>= 3.6.0),

+    R (>= 4.0),

     SummarizedExperiment,

     SingleCellExperiment,

     DelayedArray,

@@ -47,8 +47,8 @@ Imports:

     fields,

     ggplot2,

     ggplotify,

-    ggtree,

     ggrepel,

+    ggtree,

     gridExtra,

     GSVA (>= 1.26.0),

     GSVAdata,

@@ -118,7 +118,6 @@ Suggests:

     org.Mm.eg.db,

     stringr

 Remotes: 

-    joshua-d-campbell/DoubletFinder,

     joshua-d-campbell/harmony,

     rz2333/mnormt@1.5-6,

     joshua-d-campbell/liger,

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

 # Generated by roxygen2: do not edit by hand

-export(ComBatSCE)

 export(MAST)

 export(alignSingleCellData)

 export(calcEffectSizes)

 export(combineSCE)

+export(computeZScore)

 export(constructSCE)

 export(convertGeneIDs)

 export(convertSCEToSeurat)

@@ -15,9 +15,9 @@ export(downSampleCells)

 export(downSampleDepth)

 export(enrichRSCE)

 export(exportSCE)

+export(exportSCEtoAnnData)

 export(exportSCEtoFlatFile)

 export(featureIndex)

-export(filterSCData)

 export(findMarkerDiffExp)

 export(generateMeta)

 export(generateSimulatedData)

@@ -31,6 +31,7 @@ export(getUMAP)

 export(gsvaPlot)

 export(gsvaSCE)

 export(importAnnData)

+export(importBUStools)

 export(importCellRanger)

 export(importCellRangerV2)

 export(importCellRangerV2Sample)

@@ -43,6 +44,9 @@ export(importGeneSetsFromCollection)

 export(importGeneSetsFromGMT)

 export(importGeneSetsFromList)

 export(importGeneSetsFromMSigDB)

+export(importOptimus)

+export(importSEQC)

+export(importSTARsolo)

 export(iterateSimulations)

 export(mergeSCEColData)

 export(parseRsubreadLogs)

@@ -81,6 +85,7 @@ export(plotScrubletResults)

 export(plotTSNE)

 export(plotUMAP)

 export(qcInputProcess)

+export(readSingleCellMatrix)

 export(reportCellQC)

 export(reportDropletQC)

 export(reportQCTool)

@@ -149,6 +154,7 @@ export(subsetSCECols)

 export(subsetSCERows)

 export(summarizeSCE)

 export(thresholdGenes)

+export(trimCounts)

 export(visPlot)

 import(Biobase)

 import(DelayedArray)

deleted file mode 100644

@@ -1,79 +0,0 @@

-#' ComBatSCE

-#'

-#' Run ComBat on a SCtkExperiment object

-#'

-#' @param inSCE Input SCtkExperiment object. Required

-#' @param batch The name of a column in colData to use as the batch variable.

-#' Required

-#' @param useAssay The assay to use for ComBat. The default is "logcounts"

-#' @param par.prior TRUE indicates parametric adjustments will be used, FALSE

-#' indicates non-parametric adjustments will be used. Accepted parameters:

-#' "Parametric" or "Non-parametric"

-#' @param covariates List of other column names in colData to be added to the

-#' ComBat model as covariates

-#' @param mean.only If TRUE ComBat only corrects the mean of the batch effect

-#' @param ref.batch If given, will use the selected batch as a reference for

-#' batch adjustment.

-#'

-#' @return ComBat matrix based on inputs. You can save this matrix into the

-#' SCtkExperiment with assay()

-#' @export

-#' @examples

-#' if(requireNamespace("bladderbatch", quietly = TRUE)) {

-#'   library(bladderbatch)

-#'   data(bladderdata)

-#'

-#'   #subset for testing

-#'   dat <- bladderEset[1:50,]

-#'   dat <- as(as(dat, "SummarizedExperiment"), "SCtkExperiment")

-#'   mod <- stats::model.matrix(~as.factor(cancer), data = colData(dat))

-#'

-#'   # parametric adjustment

-#'   combat_edata1 <- ComBatSCE(inSCE = dat, useAssay = "exprs",

-#'                              batch = "batch", covariates = NULL)

-#'   assay(dat, "parametric_combat") <- combat_edata1

-#'

-#'   # non-parametric adjustment, mean-only version

-#'   combat_edata2 <- ComBatSCE(inSCE = dat, useAssay = "exprs",

-#'                              batch = "batch", par.prior = "Non-parametric",

-#'                              mean.only = TRUE, covariates = NULL)

-#'   assay(dat, "nonparametric_combat_meanonly") <- combat_edata2

-#'

-#'   # reference-batch version, with covariates

-#'   combat_edata3 <- ComBatSCE(inSCE = dat, useAssay = "exprs",

-#'                              batch = "batch", covariates = "cancer",

-#'                              ref.batch = 3)

-#'   assay(dat, "refbatch_combat_wcov") <- combat_edata3

-#'   assays(dat)

-#' }

-#'

-ComBatSCE <- function(inSCE, batch, useAssay="logcounts",

-                      par.prior="Parametric", covariates=NULL, mean.only=FALSE,

-                      ref.batch=NULL){

-

-  #prepare model matrix

-  mod <- NULL

-  if (length(covariates) > 0){

-    mod <- stats::model.matrix(

-      stats::as.formula(paste0("~", paste0(covariates, collapse = "+"))),

-      data = data.frame(SingleCellExperiment::colData(inSCE)[, covariates,

-                                                               drop = FALSE]))

-  }

-

-  #prepare parametric

-  if (par.prior == "Parametric"){

-    par.prior <- TRUE

-  } else if (par.prior == "Non-parametric") {

-    par.prior <- FALSE

-  } else {

-    stop("Invalid option given to par.prior. Accepted values are Parametric",

-         " and Non-parametric.")

-  }

-

-  resassay <-

-    sva::ComBat(dat = SummarizedExperiment::assay(inSCE, useAssay),

-                batch = SingleCellExperiment::colData(inSCE)[, batch],

-                mod = mod, par.prior = par.prior,

-                mean.only = mean.only, ref.batch = ref.batch)

-  return(resassay)

-}

new file mode 100644

@@ -0,0 +1,25 @@

+#' Compute Z-Score

+#'

+#' Computes Z-Score from an input count matrix using the formula ((x-mean(x))/sd(x))

+#' for each gene across all cells. The input count matrix can either be a base matrix,

+#' dgCMatrix or a DelayedMatrix. Computations are performed using DelayedMatrixStats

+#' package to efficiently compute the Z-Score matrix.

+#'

+#' @param counts matrix (base matrix, dgCMatrix or DelaymedMatrix)

+#'

+#' @return z-score computed counts matrix (DelayedMatrix)

+#' @export

+#'

+#' @examples

+#'

+#' data(sce_chcl, package = "scds")

+#' assay(sce_chcl, "countsZScore") <- computeZScore(assay(sce_chcl, "counts"))

+#'

+computeZScore <- function(counts) {

+    if (!methods::is(counts, "DelayedArray")) {

+        counts <- DelayedArray::DelayedArray(counts)

+    }

+    counts <- (counts - DelayedMatrixStats::rowMeans2(counts)) / DelayedMatrixStats::rowSds(counts)

+    counts[base::is.nan(counts)] <- 0

+    return(counts)

+}

new file mode 100644

@@ -0,0 +1,75 @@

+#' @title Export a \link[SingleCellExperiment]{SingleCellExperiment} R object as 

+#' Python annData object 

+#' @description Writes all assays, colData, rowData, reducedDims, and altExps objects in a

+#' \link[SingleCellExperiment]{SingleCellExperiment} to a Python annData object in the .h5ad format

+#' All parameters of Anndata.write_h5ad function (https://icb-anndata.readthedocs-hosted.com/en/stable/anndata.AnnData.write_h5ad.html)

+#' are available as parameters to this export function and set to defaults. Defaults can be

+#' overridden at function call.

+#' @param sce \link[SingleCellExperiment]{SingleCellExperiment} R object to be

+#'  exported.

+#' @param useAssay Character. The name of assay of

+#' interests that will be set as the primary matrix of the output AnnData.

+#' Default \code{"counts"}. 

+#' @param outputDir Path to the directory where .h5ad outputs will be written. Default is the current working directory.

+#' @param prefix Prefix to use for the name of the output file. Default \code{"sample"}.

+#' @param overwrite Boolean. Default \code{TRUE}.

+#' @param compression If output file compression is required, this variable accepts

+#' 'gzip' or 'lzf' as inputs. Default \code{None}.

+#' @param compressionOpts Integer. Sets the compression level

+#' @param forceDense Default \code{False} Write sparse data as a dense matrix.

+#' Refer \code{anndata.write_h5ad} documentation for details. Default \code{NULL}. 

+#' @examples

+#' \dontrun{

+#' data(sce_chcl, package = "scds")

+#' exportSCEtoAnnData(sce=sce_chcl, compression="gzip")

+#' }

+#' 

+#' @export

+exportSCEtoAnnData <- function(sce, 

+                                useAssay = 'counts',

+                                outputDir = "./",

+                                prefix = "sample",

+                                overwrite = TRUE,

+                                compression = c('None','lzf','gzip'),

+                                compressionOpts = NULL,

+                                forceDense = c('False','True')){

+  compression <- match.arg(compression)

+  forceDense <- match.arg(forceDense)

+  if (compression == 'None'){

+    compression <- NULL

+  }

+

+  if (!reticulate::py_module_available(module = "scanpy")) {

+    warning("Cannot find python module 'scanpy', please install Conda and",

+            " run sctkPythonInstallConda() or run sctkPythonInstallVirtualEnv().",

+            "If one of these have been previously run to install the modules,",

+            "make sure to run selectSCTKConda() or selectSCTKVirtualEnvironment(),",

+            " respectively, if R has been restarted since the module installation.",

+            " Alternatively, Scrublet can be installed on the local machine",

+            "with pip (e.g. pip install --user scanpy) and then the 'use_python()'",

+            " function from the 'reticulate' package can be used to select the",

+            " correct Python environment.")

+    return(sce)}

+  

+  AssayName <- SummarizedExperiment::assayNames(sce)

+  for (assay in AssayName){

+    if (!methods::is(SummarizedExperiment::assay(sce, assay), 'dgCMatrix')) {

+      SummarizedExperiment::assay(sce, assay) <- .convertToMatrix(SummarizedExperiment::assay(sce, assay))

+    }

+  }

+

+  

+  dir.create(outputDir, showWarnings = FALSE, recursive = TRUE)

+  annData <- .sce2adata(sce,useAssay)

+  fileName <- paste0(prefix,".h5ad")

+  filePath <- file.path(outputDir,fileName)

+  

+  if (file.exists(filePath) && !isTRUE(overwrite)) {

+    stop(paste0(path, " already exists. Change 'outputDir' or set 'overwrite' to TRUE."))

+    }

+

+  annData$write_h5ad(filePath,

+                     compression = compression, 

+                     compression_opts = compressionOpts,

+                     force_dense = forceDense)

+}

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

-#' @title Get runDropletQC .html report 

+#' @title Get runDropletQC .html report

 #' @description A  function to generate .html Rmarkdown report containing the visualizations of the runDropletQC function output

 #' @param inSCE A \link[SingleCellExperiment]{SingleCellExperiment} object containing

 #' the full droplet count matrix with the output from runDropletQC function

-#' @param subTitle subtitle of the QC HTML report. Default is NULL. 

-#' @param studyDesign description of the data set and experiment design. It would be shown at the top of QC HTML report. Default is NULL. 

-#' @param output_file name of the generated file. If NULL/default then the output file name will be based on the name of the Rmarkdown template 

+#' @param subTitle subtitle of the QC HTML report. Default is NULL.

+#' @param studyDesign description of the data set and experiment design. It would be shown at the top of QC HTML report. Default is NULL.

+#' @param output_file name of the generated file. If NULL/default then the output file name will be based on the name of the Rmarkdown template

 #' @param output_dir name of the output directory to save the rendered file. If NULL/default the file is stored to the current working directory

 #' @return .html file

 #' @examples

@@ -16,33 +16,33 @@

 #' @export

 reportDropletQC <- function(inSCE, output_file = NULL,

                                    output_dir = NULL,

-                                   subTitle = NULL, 

+                                   subTitle = NULL,

                                    studyDesign = NULL) {

-  

+

   if (is.null(output_dir)){

     output_dir<- getwd()

     }

-  

+

   #report_path <- tempfile(fileext = ".Rmd")

   #file.copy(system.file("rmarkdown/qc/DropletQC.Rmd", package = "singleCellTK"), report_path, overwrite = TRUE)

   ## create temp Rmd file to bypass permission issue on server

-  rmarkdown::render(system.file("rmarkdown/qc/DropletQC.Rmd", package = "singleCellTK"), 

-    params = list(object = inSCE, subTitle = subTitle, studyDesign = studyDesign), 

-    output_file = output_file, 

+  rmarkdown::render(system.file("rmarkdown/qc/DropletQC.Rmd", package = "singleCellTK"),

+    params = list(object = inSCE, subTitle = subTitle, studyDesign = studyDesign),

+    output_file = output_file,

     output_dir = output_dir,

     intermediates_dir = output_dir,

     knit_root_dir = output_dir)

  }

-#' @title Get runCellQC .html report 

+#' @title Get runCellQC .html report

 #' @description A  function to generate .html Rmarkdown report containing the visualizations of the runCellQC function output

 #' @param inSCE A \link[SingleCellExperiment]{SingleCellExperiment} object containing

 #' the filtered count matrix with the output from runCellQC function

-#' @param subTitle subtitle of the QC HTML report. Default is NULL. 

-#' @param studyDesign description of the data set and experiment design. It would be shown at the top of QC HTML report. Default is NULL. 

-#' @param output_file name of the generated file. If NULL/default then the output file name will be based on the name of the Rmarkdown template. 

+#' @param subTitle subtitle of the QC HTML report. Default is NULL.

+#' @param studyDesign description of the data set and experiment design. It would be shown at the top of QC HTML report. Default is NULL.

+#' @param output_file name of the generated file. If NULL/default then the output file name will be based on the name of the Rmarkdown template.

 #' @param output_dir name of the output directory to save the rendered file. If NULL/default the file is stored to the current working directory

 #' @return .html file

 #' @examples

@@ -55,7 +55,7 @@ reportDropletQC <- function(inSCE, output_file = NULL,

 #' @export

 reportCellQC <- function(inSCE, output_file = NULL,

                                 output_dir = NULL,

-                                subTitle = NULL, 

+                                subTitle = NULL,

                                 studyDesign = NULL) {

   if (is.null(output_dir)){

     output_dir<- getwd()

@@ -63,9 +63,9 @@ reportCellQC <- function(inSCE, output_file = NULL,

   #report_path <- tempfile(fileext = ".Rmd")

   #file.copy(system.file("rmarkdown/qc/CellQC.Rmd", package = "singleCellTK"), report_path, overwrite = TRUE)

-  rmarkdown::render(system.file("rmarkdown/qc/CellQC.Rmd", package = "singleCellTK"), 

-    params = list(object = inSCE, subTitle = subTitle, studyDesign = studyDesign), 

-    output_file = output_file, 

+  rmarkdown::render(system.file("rmarkdown/qc/CellQC.Rmd", package = "singleCellTK"),

+    params = list(object = inSCE, subTitle = subTitle, studyDesign = studyDesign),

+    output_file = output_file,

     output_dir = output_dir,

     intermediates_dir = output_dir,

     knit_root_dir = output_dir)

@@ -79,7 +79,7 @@ reportCellQC <- function(inSCE, output_file = NULL,

 #' runQCMetrics, runScrublet, runDoubletCells, runCxds, runBcds, runCxdsBcdsHybrid, runDecontX, runBarcodeRankDrops, runEmptyDrops

 #' @param algorithm Character. Specifies which QC algorithm report to generate.

 #'  Available options are "BarcodeRankDrops", "EmptyDrops", "QCMetrics", "Scrublet", "DoubletCells", "Cxds", "Bcds", "CxdsBcdsHybrid", "DoubletFinder"  and "DecontX".

-#' @param output_file name of the generated file. If NULL/default then the output file name will be based on the name of the selected QC algorithm name . 

+#' @param output_file name of the generated file. If NULL/default then the output file name will be based on the name of the selected QC algorithm name .

 #' @param output_dir name of the output directory to save the rendered file. If NULL/default the file is stored to the current working directory

 #' @return .html file

 #' @examples

@@ -87,6 +87,7 @@ reportCellQC <- function(inSCE, output_file = NULL,

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

 #' sce <- sce[, colData(sce)$type != 'EmptyDroplet']

 #' sce <- runDecontX(sce)

+#' sce <- getUMAP(sce)

 #' reportQCTool(inSCE = sce, algorithm = "DecontX")

 #' }

 #' @export

@@ -102,13 +103,13 @@ reportQCTool <- function(inSCE, algorithm=c("BarcodeRankDrops",

                                             "DecontX"),

                          output_file = NULL,

                             output_dir = NULL) {

-  

+

   algorithm <- match.arg(algorithm)

-  

+

   if (is.null(output_dir)){

     output_dir<- getwd()

   }

-  

+

   if (algorithm =="BarcodeRankDrops"){

     rmarkdown::render(system.file("rmarkdown/qc/BarcodeRankDrops.Rmd", package="singleCellTK"), params = list(

       object=inSCE), output_file = output_file, output_dir = output_dir)

@@ -152,4 +153,4 @@ reportQCTool <- function(inSCE, algorithm=c("BarcodeRankDrops",

  }

-    

+

new file mode 100644

@@ -0,0 +1,156 @@

+

+# dir <- "genecount"

+.constructSCEFromBUStoolsOutputs <- function(dir,

+    sample,

+    matrixFileName,

+    featuresFileName,

+    barcodesFileName,

+    gzipped,

+    class,

+    delayedArray) {

+

+    cb <- .readBarcodes(file.path(dir, barcodesFileName))

+    fe <- .readFeatures(file.path(dir, featuresFileName))

+    ma <- .readMatrixMM(file.path(dir, matrixFileName),

+        gzipped = gzipped,

+        class = class,

+        delayedArray = delayedArray)

+    ma <- t(ma)

+

+    coln <- paste(sample, cb[[1]], sep = "_")

+    rownames(ma) <- fe[[1]]

+

+    sce <- SingleCellExperiment::SingleCellExperiment(

+        assays = list(counts = ma))

+    SummarizedExperiment::rowData(sce) <- fe

+    SummarizedExperiment::colData(sce) <- S4Vectors::DataFrame(cb,

+        column_name = coln,

+        sample = sample,

+        row.names = coln)

+

+    return(sce)

+}

+

+

+# main function

+.importBUStools <- function(

+    BUStoolsDirs,

+    samples,

+    matrixFileNames,

+    featuresFileNames,

+    barcodesFileNames,

+    gzipped,

+    class,

+    delayedArray) {

+

+    if (length(BUStoolsDirs) != length(samples)) {

+        stop("'BUStoolsDirs' and 'samples' have unequal lengths!")

+    }

+

+    res <- vector("list", length = length(samples))

+

+    matrixFileNames <- .getVectorized(matrixFileNames, length(samples))

+    featuresFileNames <- .getVectorized(featuresFileNames, length(samples))

+    barcodesFileNames <- .getVectorized(barcodesFileNames, length(samples))

+    gzipped <- .getVectorized(gzipped, length(samples))

+

+    for (i in seq_along(samples)) {

+        dir <- file.path(BUStoolsDirs[i])

+        scei <- .constructSCEFromBUStoolsOutputs(dir,

+            sample = samples[i],

+            matrixFileName = matrixFileNames[i],

+            featuresFileName = featuresFileNames[i],

+            barcodesFileName = barcodesFileNames[i],

+            gzipped = gzipped[i],

+            class = class,

+            delayedArray = delayedArray)

+        res[[i]] <- scei

+    }

+

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

+    return(sce)

+}

+

+

+#' @name importBUStools

+#' @rdname importBUStools

+#' @title Construct SCE object from BUStools output

+#' @description Read the barcodes, features (genes), and matrix from BUStools

+#'  output. Import them

+#'  as one \link[SingleCellExperiment]{SingleCellExperiment} object. Note the

+#'  cells in the output files for BUStools 0.39.4 are not filtered.

+#' @param BUStoolsDirs A vector of paths to BUStools output files. Each sample

+#'  should have its own path. For example: \code{./genecount}.

+#'  Must have the same length as \code{samples}.

+#' @param samples A vector of user-defined sample names for the samples to be

+#'  imported. Must have the same length as \code{BUStoolsDirs}.

+#' @param matrixFileNames Filenames for the Market Exchange Format (MEX) sparse

+#'  matrix files (.mtx files). Must have length 1 or the same

+#'  length as \code{samples}.

+#' @param featuresFileNames Filenames for the feature annotation files.

+#'  Must have length 1 or the same length as \code{samples}.

+#' @param barcodesFileNames Filenames for the cell barcode list file.

+#'  Must have length 1 or the same length as \code{samples}.

+#' @param gzipped Boolean. \code{TRUE} if the BUStools output files

+#'  (barcodes.txt, genes.txt, and genes.mtx) were

+#'  gzip compressed. \code{FALSE} otherwise. This is \code{FALSE} in BUStools

+#'  0.39.4. Default \code{"auto"} which automatically detects if the

+#'  files are gzip compressed. Must have length 1 or the same length as

+#'  \code{samples}.

+#' @param class Character. The class of the expression matrix stored in the SCE

+#'  object. Can be one of "Matrix" (as returned by

+#'  \link[Matrix]{readMM} function), or "matrix" (as returned by

+#'  \link[base]{matrix} function). Default "Matrix".

+#' @param delayedArray Boolean. Whether to read the expression matrix as

+#'  \link[DelayedArray]{DelayedArray} object or not. Default \code{TRUE}.

+#' @return A \code{SingleCellExperiment} object containing the count

+#'  matrix, the gene annotation, and the cell annotation.

+#' @examples

+#' # Example #1

+#' # FASTQ files were downloaded from

+#' # https://support.10xgenomics.com/single-cell-gene-expression/datasets/3.0.0

+#' # /pbmc_1k_v3

+#' # They were concatenated as follows:

+#' # cat pbmc_1k_v3_S1_L001_R1_001.fastq.gz pbmc_1k_v3_S1_L002_R1_001.fastq.gz >

+#' # pbmc_1k_v3_R1.fastq.gz

+#' # cat pbmc_1k_v3_S1_L001_R2_001.fastq.gz pbmc_1k_v3_S1_L002_R2_001.fastq.gz >

+#' # pbmc_1k_v3_R2.fastq.gz

+#' # The following BUStools command generates the gene, cell, and

+#' # matrix files

+#'

+#' # bustools correct -w ./3M-february-2018.txt -p output.bus | \

+#' #   bustools sort -T tmp/ -t 4 -p - | \

+#' #   bustools count -o genecount/genes \

+#' #     -g ./transcripts_to_genes.txt \

+#' #     -e matrix.ec \

+#' #     -t transcripts.txt \

+#' #     --genecounts -

+#'

+#' # The top 20 genes and the first 20 cells are included in this example.

+#' sce <- importBUStools(

+#'   BUStoolsDirs = system.file("extdata/BUStools_PBMC_1k_v3_20x20/genecount/",

+#'     package = "singleCellTK"),

+#'   samples = "PBMC_1k_v3_20x20")

+#' @export

+importBUStools <- function(

+    BUStoolsDirs,

+    samples,

+    matrixFileNames = "genes.mtx",

+    featuresFileNames = "genes.genes.txt",

+    barcodesFileNames = "genes.barcodes.txt",

+    gzipped = "auto",

+    class = c("Matrix", "matrix"),

+    delayedArray = TRUE) {

+

+    class <- match.arg(class)

+

+    .importBUStools(

+        BUStoolsDirs = BUStoolsDirs,

+        samples = samples,

+        matrixFileNames = matrixFileNames,

+        featuresFileNames = featuresFileNames,

+        barcodesFileNames = barcodesFileNames,

+        gzipped = gzipped,

+        class = class,

+        delayedArray = delayedArray)

+}

new file mode 100644

@@ -0,0 +1,285 @@

+

+#' @importFrom reticulate import

+.readMatrixNpz <- function(matrixLocation,

+  colIndexLocation,

+  rowIndexLocation,

+  class,

+  delayedArray) {

+

+  ## Now importing these functions in 'reticulate_setup.R' file

+  #  sparse <- reticulate::import("scipy.sparse")

+  #  numpy <- reticulate::import("numpy")

+  if (!reticulate::py_module_available(module = "scipy.sparse")) {

+    stop("Error!", "Cannot find python module 'scipy.sparse', please install Conda and run sctkPythonInstallConda() 

+         or run sctkPythonInstallVirtualEnv(). If one of these have been previously run to install the modules,

+         make sure to run selectSCTKConda() or selectSCTKVirtualEnvironment(), respectively, if R has been

+         restarted since the module installation. Alternatively, scipy can be installed on the local machine

+         with pip (e.g. pip install scipy) and then the 'use_python()' function from the 'reticulate' package

+         can be used to select the correct Python environment.")

+  }

+  if (!reticulate::py_module_available(module = "numpy")) {

+    stop("Error!", "Cannot find python module 'numpy', please install Conda and run sctkPythonInstallConda() 

+         or run sctkPythonInstallVirtualEnv(). If one of these have been previously run to install the modules,

+         make sure to run selectSCTKConda() or selectSCTKVirtualEnvironment(), respectively, if R has been

+         restarted since the module installation. Alternatively, numpy can be installed on the local machine

+         with pip (e.g. pip install numpy) and then the 'use_python()' function from the 'reticulate' package

+         can be used to select the correct Python environment.")

+  }

+  

+  error <- try({

+    mat <- sparse$load_npz(matrixLocation)

+    colIndex <- as.vector(numpy$load(colIndexLocation, allow_pickle = TRUE))

+    rowIndex <- as.vector(numpy$load(rowIndexLocation, allow_pickle = TRUE))

+    colnames(mat) <- colIndex

+    rownames(mat) <- rowIndex

+    mat <- t(mat)

+

+    ## Convert to "dgCMatrix"

+    newM <- Matrix::Matrix(mat[,1], nrow=nrow(mat))

+    newM <- methods::as(newM, "dgCMatrix")

+    breaks <- seq(2, ncol(mat), by=1000)

+    if(length(breaks) > 2) {

+      for(i in seq(2, length(breaks))) {

+        ix <- seq(breaks[i-1], (breaks[i]-1))

+        newM <- cbind(newM, mat[,ix])

+      }

+      ix <- seq(utils::tail(breaks, n = 1), ncol(mat))

+      newM <- cbind(newM, mat[,ix])

+    } else {

+      ix <- seq(2, ncol(mat))

+      newM <- cbind(newM, mat[,ix])

+    }

+

+    colnames(newM) <- colnames(mat)

+    rownames(newM) <- rownames(mat)

+    mat <- newM

+  }, silent = TRUE)

+  

+  if(inherits(error, "try-error")) {

+    stop(paste0("importOptimus did not complete successfully. SCE could not be generated. Error given during the import process: \n\n", error))

+  }

+  

+  if (class == "matrix") {

+    mat <- as.matrix(mat)

+  }

+

+  if (isTRUE(delayedArray)) {

+    mat <- DelayedArray::DelayedArray(mat)

+  }

+  return(mat)

+}

+

+

+.readMetricsOptimus <- function(path) {

+  metrics <- data.table::fread(path)

+  return(metrics)

+}

+

+

+.readEmptyDrops <- function(path) {

+  emptyDrops <- data.table::fread(path)

+  colnames(emptyDrops) <- paste0("dropletUtils_emptyDrops_",

+    colnames(emptyDrops))

+  return(emptyDrops)

+}

+

+

+.combineColData <- function(colnames, cellMetrics, emptyDrops) {

+  cd <- data.table::data.table(CellId = colnames)

+  cd <- merge(cd,

+    cellMetrics,

+    by.x = "CellId",

+    by.y = "V1",

+    all.x = TRUE,

+    all.y = FALSE,

+    sort = FALSE)

+

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

+    cd <- merge(cd,

+      emptyDrops,

+      by.x = "CellId",

+      by.y = "dropletUtils_emptyDrops_CellId",

+      all.x = TRUE,

+      all.y = FALSE,

+      sort = FALSE)

+  }

+

+  return(cd)

+}

+

+

+.combineRowData <- function(rownames, geneMetrics) {

+  rd <- data.table::data.table(feature_ID = rownames)

+  rd <- merge(rd,

+    geneMetrics,

+    by.x = "feature_ID",

+    by.y = "V1",

+    all.x = TRUE,

+    all.y = FALSE,

+    sort = FALSE)

+  return(rd)

+}

+

+

+.constructSCEFromOptimusOutputs <- function(dir,

+  sample,

+  matrixLocation,

+  colIndexLocation,

+  rowIndexLocation,

+  cellMetricsLocation,

+  geneMetricsLocation,

+  emptyDropsLocation,

+  class,

+  delayedArray) {

+

+  mat <- .readMatrixNpz(file.path(dir, matrixLocation),

+    file.path(dir, colIndexLocation),

+    file.path(dir, rowIndexLocation),

+    class,

+    delayedArray)

+

+  cellMetrics <- .readMetricsOptimus(file.path(dir, cellMetricsLocation))

+  geneMetrics <- .readMetricsOptimus(file.path(dir, geneMetricsLocation))

+

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

+    emptyDrops <- .readEmptyDrops(file.path(dir, emptyDropsLocation))

+  }

+

+  cd <- .combineColData(colnames(mat), cellMetrics, emptyDrops)

+  rd <- .combineRowData(rownames(mat), geneMetrics)

+

+  coln <- paste(sample, colnames(mat), sep = "_")

+

+  sce <- SingleCellExperiment::SingleCellExperiment(

+    assays = list(counts = mat))

+  SummarizedExperiment::colData(sce) <- S4Vectors::DataFrame(cd,

+    column_name = coln,

+    sample = sample,

+    row.names = coln)

+  SummarizedExperiment::rowData(sce) <- S4Vectors::DataFrame(rd)

+

+  return(sce)

+}

+

+

+.checkArgsImportOptimus <- function(OptimusDirs,

+  samples) {

+

+  if (length(OptimusDirs) != length(samples)) {

+    stop("'OptimusDirs' and 'samples' have unequal lengths!")

+  }

+}

+

+

+.importOptimus <- function(OptimusDirs,

+  samples,

+  matrixLocation,

+  colIndexLocation,

+  rowIndexLocation,

+  cellMetricsLocation,

+  geneMetricsLocation,

+  emptyDropsLocation,

+  class,

+  delayedArray) {

+

+  .checkArgsImportOptimus(OptimusDirs, samples)

+

+  res <- vector("list", length = length(samples))

+

+  for (i in seq_along(samples)) {

+    scei <- .constructSCEFromOptimusOutputs(OptimusDirs[i],

+      sample = samples[i],

+      matrixLocation = matrixLocation,

+      colIndexLocation = colIndexLocation,

+      rowIndexLocation = rowIndexLocation,

+      cellMetricsLocation = cellMetricsLocation,

+      geneMetricsLocation = geneMetricsLocation,

+      emptyDropsLocation = emptyDropsLocation,

+      class = class,

+      delayedArray = delayedArray)

+    res[[i]] <- scei

+  }

+

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

+  return(sce)

+}

+

+

+#' @name importOptimus

+#' @rdname importOptimus

+#' @title Construct SCE object from Optimus output

+#' @description Read the barcodes, features (genes), and matrices from Optimus

+#'  outputs. Import them

+#'  as one \link[SingleCellExperiment]{SingleCellExperiment} object.

+#' @param OptimusDirs A vector of root directories of Optimus output files.

+#'  The paths should be something like this:

+#'  \code{/PATH/TO/bb4a2a5e-ff34-41b6-97d2-0c0c0c534530}.

+#'  Each entry in \code{OptimusDirs} is considered a sample and should have

+#'  its own path. Must have the same length as \code{samples}.

+#' @param samples A vector of user-defined sample names for the sample to be

+#'  imported. Must have the same length as \code{OptimusDirs}.

+#' @param matrixLocation Character. It is the intermediate

+#'  path to the filtered count maxtrix file saved in sparse matrix format

+#'  (\code{.npz}). Default

+#'  \code{call-MergeCountFiles/sparse_counts.npz} which works for

+#'  optimus_v1.4.0.

+#' @param colIndexLocation Character. The intermediate path to the barcode

+#'  index file. Default \code{call-MergeCountFiles/sparse_counts_col_index.npy}.

+#' @param rowIndexLocation Character. The intermediate path to the feature

+#'  (gene) index file. Default

+#'  \code{call-MergeCountFiles/sparse_counts_row_index.npy}.

+#' @param cellMetricsLocation Character. It is the intermediate

+#'  path to the cell metrics file (\code{merged-cell-metrics.csv.gz}). Default

+#'  \code{call-MergeCellMetrics/merged-cell-metrics.csv.gz} which works for

+#'  optimus_v1.4.0.

+#' @param geneMetricsLocation Character. It is the intermediate

+#'  path to the feature (gene) metrics file (\code{merged-gene-metrics.csv.gz}).

+#'  Default \code{call-MergeGeneMetrics/merged-gene-metrics.csv.gz} which works

+#'  for optimus_v1.4.0.

+#' @param emptyDropsLocation Character. It is the intermediate

+#'  path to \link[DropletUtils]{emptyDrops} metrics file

+#'  (\code{empty_drops_result.csv}).

+#'  Default \code{call-RunEmptyDrops/empty_drops_result.csv} which works for

+#'  optimus_v1.4.0.

+#' @param class Character. The class of the expression matrix stored in the SCE

+#'  object. Can be one of "Matrix" (as returned by

+#'  \link[Matrix]{readMM} function), or "matrix" (as returned by

+#'  \link[base]{matrix} function). Default "Matrix".

+#' @param delayedArray Boolean. Whether to read the expression matrix as

+#'  \link[DelayedArray]{DelayedArray} object or not. Default \code{TRUE}.

+#' @return A \link[SingleCellExperiment]{SingleCellExperiment} object

+#'  containing the count

+#'  matrix, the gene annotation, and the cell annotation.

+#' @examples

+#' \dontrun{

+#' sce <- importOptimus(OptimusDirs =

+#'   system.file("extdata/Optimus_20x1000",

+#'   package = "singleCellTK"),

+#'   samples = "Optimus_20x1000")

+#' }

+#' @export

+importOptimus <- function(OptimusDirs,

+  samples,

+  matrixLocation = "call-MergeCountFiles/sparse_counts.npz",

+  colIndexLocation = "call-MergeCountFiles/sparse_counts_col_index.npy",

+  rowIndexLocation = "call-MergeCountFiles/sparse_counts_row_index.npy",

+  cellMetricsLocation = "call-MergeCellMetrics/merged-cell-metrics.csv.gz",

+  geneMetricsLocation = "call-MergeGeneMetrics/merged-gene-metrics.csv.gz",

+  emptyDropsLocation = "call-RunEmptyDrops/empty_drops_result.csv",

+  class = c("Matrix", "matrix"),

+  delayedArray = TRUE) {

+

+  class <- match.arg(class)

+

+  .importOptimus(OptimusDirs = OptimusDirs,

+    samples = samples,

+    matrixLocation = matrixLocation,

+    colIndexLocation = colIndexLocation,

+    rowIndexLocation = rowIndexLocation,

+    cellMetricsLocation = cellMetricsLocation,

+    geneMetricsLocation = geneMetricsLocation,

+    emptyDropsLocation = emptyDropsLocation,

+    class = class,

+    delayedArray = delayedArray)

+

+}

new file mode 100644

@@ -0,0 +1,168 @@

+

+# dir <- "outs/filtered_feature_bc_matrix/"

+.constructSCEFromSTARsoloOutputs <- function(dir,

+    sample,

+    matrixFileName,

+    featuresFileName,