@@ -4,3 +4,5 @@

 .travis.yml

 NOTICE

 _pkgdown.yml

+^doc$

+^Meta$

@@ -33,3 +33,5 @@ etc/*

 # Celda log files with default prefix

 Celda_chain.*log.txt

 inst/doc

+doc

+Meta

@@ -17,14 +17,6 @@ calculateNativeMatrix <- function(counts, native_counts, theta, eta, phi, z, row

     .Call('_celda_calculateNativeMatrix', PACKAGE = 'celda', counts, native_counts, theta, eta, phi, z, row_index, col_index, pseudocount)

 }

-#' get row and column indices of none zero elements in the matrix

-#' 

-#' @param R_counts A matrix

-#' @return An integer matrix where each row is a row, column indices pair 

-nonzero <- function(R_counts) {

-    .Call('_celda_nonzero', PACKAGE = 'celda', R_counts)

-}

-

 cG_calcGibbsProbY_Simple <- function(counts, nGbyTS, nTSbyC, nbyTS, nbyG, y, L, index, gamma, beta, delta) {

     .Call('_celda_cG_calcGibbsProbY_Simple', PACKAGE = 'celda', counts, nGbyTS, nTSbyC, nbyTS, nbyG, y, L, index, gamma, beta, delta)

 }

@@ -77,3 +69,11 @@ fastNormPropSqrt <- function(R_counts, R_alpha) {

     .Call('_celda_fastNormPropSqrt', PACKAGE = 'celda', R_counts, R_alpha)

 }

+#' get row and column indices of none zero elements in the matrix

+#' 

+#' @param R_counts A matrix

+#' @return An integer matrix where each row is a row, column indices pair 

+nonzero <- function(R_counts) {

+    .Call('_celda_nonzero', PACKAGE = 'celda', R_counts)

+}

+

@@ -162,9 +162,11 @@ plotHeatmap <- function(counts,

     ## Select subsets of features/cells

     if (!is.null(featureIx)) {

         counts <- counts[featureIx, , drop = FALSE]

-        if (length(annotationFeature) > 1 ||

-                (length(annotationFeature) == 1 & !is.na(annotationFeature))) {

+        if (!is.null(annotationFeature) & length(annotationFeature) > 0) {

+          if(length(annotationFeature) == 1 & !is.na(annotationFeature) | 

+             length(annotationFeature) > 1) {

             annotationFeature <- annotationFeature[featureIx, , drop = FALSE]

+          }  

         }

         if (!is.null(y)) {

             y <- y[featureIx]

@@ -173,9 +175,11 @@ plotHeatmap <- function(counts,

     if (!is.null(cellIx)) {

         counts <- counts[, cellIx, drop = FALSE]

-        if (length(annotationCell) > 1 ||

-                (length(annotationCell) == 1 & !is.na(annotationCell))) {

+        if (!is.null(annotationCell) & length(annotationCell) > 0) {

+          if(length(annotationCell) == 1 & !is.na(annotationCell) | 

+             length(annotationFeature) > 1) {

             annotationCell <- annotationCell[cellIx, , drop = FALSE]

+          }  

         }

         if (!is.null(z)) {

             z <- z[cellIx]

@@ -62,6 +62,7 @@

 #'  }

 #' }

 #' @examples

+#' \dontrun{

 #' library(M3DExampleData)

 #' counts <- M3DExampleData::Mmus_example_list$data

 #' # subset 100 genes for fast clustering

@@ -81,6 +82,7 @@

 #'

 #' # Plot dendrogram

 #' plotDendro(DecTree)

+#' }

 #' @import magrittr

 #' @importFrom methods hasArg

 #' @export

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

 #' @param features A L(features) by N(samples) numeric matrix.

 #' @return A character vector of label predicitions.

 #' @examples

+#' \dontrun{

 #' library(M3DExampleData)

 #' counts <- M3DExampleData::Mmus_example_list$data

 #' # Subset 500 genes for fast clustering

@@ -26,6 +27,7 @@

 #'

 #' # Get sample estimates in training data

 #' getDecisions(DecTree$rules, features)

+#' }

 #' @export

 getDecisions <- function(rules, features) {

     features <- t(features)

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

 #' @param boxSize A numeric value. Size of rule labels. Default is 7.

 #' @param boxColor A character value. Color of rule labels. Default is `black`.

 #' @examples

+#' \dontrun{

 #' library(M3DExampleData)

 #' counts <- M3DExampleData::Mmus_example_list$data

 #' # Subset 500 genes for fast clustering

@@ -31,6 +32,7 @@

 #'

 #' # Plot dendrogram

 #' plotDendro(decTree)

+#' }

 #' @return A ggplot2 object

 #' @import ggplot2

 #' @importFrom ggdendro dendro_data ggdendrogram

@@ -145,7 +145,7 @@ decontXMarkerPlot <- function(counts, z, geneMarkers, threshold = 1, color = "re

 # geneMarkers should be a dataframe,  w/ 2 column names being `cellType` and `geneMarkers`

 # convert both `cellType` and `geneMarkers` are factors with levels being integer

-.geneMarkerProcess <- function(geneMarkers, orders = None) {

+.geneMarkerProcess <- function(geneMarkers, orders = NULL) {

   geneMarkers[, "cellName"] <- geneMarkers[, "cellType"]

   geneMarkers[, "cellType"] <- factor(geneMarkers[, "cellType"])

   levels(geneMarkers[, "cellType"]) <- 1:length(levels(geneMarkers[, "cellType"]))

deleted file mode 100644

Binary files a/celdaTsne.ori.rds and /dev/null differ

deleted file mode 100644

Binary files a/celda_res_vignettes.rds and /dev/null differ

deleted file mode 100644

Binary files a/featureViolin.png and /dev/null differ

deleted file mode 100644

Binary files a/pltMarker.png and /dev/null differ

deleted file mode 100644

Binary files a/pltTsne.ori.png and /dev/null differ

deleted file mode 100644

Binary files a/pltViolin.ori.png and /dev/null differ

deleted file mode 100644

@@ -1,41 +0,0 @@

-// [[Rcpp::depends(RcppEigen)]]

-#include <Rcpp.h>

-using namespace Rcpp;

-

-//' get row and column indices of none zero elements in the matrix

-//' 

-//' @param R_counts A matrix

-//' @return An integer matrix where each row is a row, column indices pair 

-// [[Rcpp::export]]

-SEXP nonzero(NumericMatrix R_counts) {

-

-    IntegerVector row(1);

-		IntegerVector col(1);

-		NumericVector val(1);

-

-    int nR = R_counts.nrow();

-    int nC = R_counts.ncol();

-		double x;

-

-    for (int c = 0; c < nC; c++) {

-        for (int r = 0; r < nR; r++) {

-            x = R_counts[c * nR + r];

-            if (x != 0) {

-                row.push_back(r + 1);

-                col.push_back(c + 1);

-                val.push_back(x);

-            }

-        }

-    }

-

-		row.erase(0);

-		col.erase(0);

-		val.erase(0);

-

-    List res;

-		res["row"] = row;

-		res["col"] = col;

-		res["val"] = val;

-

-		return(res);

-}

@@ -73,17 +73,6 @@ BEGIN_RCPP

     return rcpp_result_gen;

 END_RCPP

 }

-// nonzero

-SEXP nonzero(NumericMatrix R_counts);

-RcppExport SEXP _celda_nonzero(SEXP R_countsSEXP) {

-BEGIN_RCPP

-    Rcpp::RObject rcpp_result_gen;

-    Rcpp::RNGScope rcpp_rngScope_gen;

-    Rcpp::traits::input_parameter< NumericMatrix >::type R_counts(R_countsSEXP);

-    rcpp_result_gen = Rcpp::wrap(nonzero(R_counts));

-    return rcpp_result_gen;

-END_RCPP

-}

 // cG_calcGibbsProbY_Simple

 NumericVector cG_calcGibbsProbY_Simple(const IntegerMatrix counts, IntegerVector nGbyTS, IntegerMatrix nTSbyC, IntegerVector nbyTS, IntegerVector nbyG, const IntegerVector y, const int L, const int index, const double gamma, const double beta, const double delta);

 RcppExport SEXP _celda_cG_calcGibbsProbY_Simple(SEXP countsSEXP, SEXP nGbyTSSEXP, SEXP nTSbyCSEXP, SEXP nbyTSSEXP, SEXP nbyGSEXP, SEXP ySEXP, SEXP LSEXP, SEXP indexSEXP, SEXP gammaSEXP, SEXP betaSEXP, SEXP deltaSEXP) {

@@ -222,6 +211,17 @@ BEGIN_RCPP

     return rcpp_result_gen;

 END_RCPP

 }

+// nonzero

+SEXP nonzero(NumericMatrix R_counts);

+RcppExport SEXP _celda_nonzero(SEXP R_countsSEXP) {

+BEGIN_RCPP

+    Rcpp::RObject rcpp_result_gen;

+    Rcpp::RNGScope rcpp_rngScope_gen;

+    Rcpp::traits::input_parameter< NumericMatrix >::type R_counts(R_countsSEXP);

+    rcpp_result_gen = Rcpp::wrap(nonzero(R_counts));

+    return rcpp_result_gen;

+END_RCPP

+}

 RcppExport SEXP _colSumByGroup(SEXP, SEXP);

 RcppExport SEXP _colSumByGroup_numeric(SEXP, SEXP);

@@ -236,7 +236,6 @@ static const R_CallMethodDef CallEntries[] = {

     {"_celda_decontXLogLik", (DL_FUNC) &_celda_decontXLogLik, 6},

     {"_celda_decontXInitialize", (DL_FUNC) &_celda_decontXInitialize, 4},

     {"_celda_calculateNativeMatrix", (DL_FUNC) &_celda_calculateNativeMatrix, 9},

-    {"_celda_nonzero", (DL_FUNC) &_celda_nonzero, 1},

     {"_celda_cG_calcGibbsProbY_Simple", (DL_FUNC) &_celda_cG_calcGibbsProbY_Simple, 11},

     {"_celda_cG_CalcGibbsProbY_ori", (DL_FUNC) &_celda_cG_CalcGibbsProbY_ori, 13},

     {"_celda_cG_CalcGibbsProbY_fastRow", (DL_FUNC) &_celda_cG_CalcGibbsProbY_fastRow, 13},

@@ -245,6 +244,7 @@ static const R_CallMethodDef CallEntries[] = {

     {"_celda_fastNormProp", (DL_FUNC) &_celda_fastNormProp, 2},

     {"_celda_fastNormPropLog", (DL_FUNC) &_celda_fastNormPropLog, 2},

     {"_celda_fastNormPropSqrt", (DL_FUNC) &_celda_fastNormPropSqrt, 2},

+    {"_celda_nonzero", (DL_FUNC) &_celda_nonzero, 1},

     {"_colSumByGroup",         (DL_FUNC) &_colSumByGroup,         2},

     {"_colSumByGroup_numeric", (DL_FUNC) &_colSumByGroup_numeric, 2},

     {"_colSumByGroupChange",   (DL_FUNC) &_colSumByGroupChange,   4},

@@ -84,3 +84,47 @@ SEXP fastNormPropSqrt(NumericMatrix R_counts, double R_alpha) {

   return Rcpp::wrap(res);

 }

+

+

+// [[Rcpp::depends(RcppEigen)]]

+#include <Rcpp.h>

+using namespace Rcpp ;

+

+//' get row and column indices of none zero elements in the matrix

+//' 

+//' @param R_counts A matrix

+//' @return An integer matrix where each row is a row, column indices pair 

+// [[Rcpp::export]]

+SEXP nonzero(NumericMatrix R_counts) {

+  

+  IntegerVector row(1);

+  IntegerVector col(1);

+  NumericVector val(1);

+  

+  int nR = R_counts.nrow();

+  int nC = R_counts.ncol();

+  double x;

+  

+  for (int c = 0; c < nC; c++) {

+    for (int r = 0; r < nR; r++) {

+      x = R_counts[c * nR + r];

+      if (x != 0) {

+        row.push_back(r + 1);

+        col.push_back(c + 1);

+        val.push_back(x);

+      }

+    }

+  }

+  

+  row.erase(0);

+  col.erase(0);

+  val.erase(0);

+  

+  List res;

+  res["row"] = row;

+  res["col"] = col;

+  res["val"] = val;

+  

+  return(res);

+}

+

deleted file mode 100644

@@ -1,154 +0,0 @@

-title: "Decontamination of ambient RNA in single-cell RNA-seq with DecontX on 10X PBMC4k data" 

-author: "Shiyi Yang, Sean Corbett, Yusuke Koga, Zhe Wang, W. Evan Johnson, Masanao Yajima, Joshua D. Campbell"

-date: "`r Sys.Date()`"

-output:

-  BiocStyle::html_document:

-    toc: true

-vignette: >

-  %\VignetteIndexEntry{Estimate and remove cross-contamination from ambient RNA for scRNA-seq data with DecontX}

-  %\VignetteEngine{knitr::rmarkdown}

-  %\VignetteEncoding{UTF-8}

-

-

-```{r, eval=TRUE, message=FALSE}

-library(celda)

-```

-

-DecontX can take either `SingleCellExperiment` object from [SingleCellExperiment package] (https://bioconductor.org/packages/release/bioc/html/SingleCellExperiment.html) or a count matrix in the format of rows being genesand columns being cells. When using SingleCellExperiment object, DecontX is using the counts matrix extracted from the SingleCellExperiment object. 

-

-# TODO: Load PBMC4k data from 10X CellRanger output

-# Load PBMC4k data from 10X

-```{r, eval=TRUE, message=FALSE}

-library(TENxPBMCData)  # selected single cell RNAseq datasets from 10X Genomics

-pbmc <- TENxPBMCData("pbmc4k")

-colnames(pbmc) <- paste(pbmc$Sample, pbmc$Barcode, sep="_") # Rename cells

-```

-

-Take a look at the count matrix in the `pbmc` singleCellExperiment object via `assay` function from [SummarizedExperiment] (https://bioconductor.org/packages/release/bioc/html/SummarizedExperiment.html). There are total 33694 genes and 4340 cells in this PBMC data.

-```{r}

-library(SummarizedExperiment)

-pbmc_counts = SummarizedExperiment::assay(pbmc, i = "counts")

-print(pbmc_counts)

-```

-

-Among the 33694 genes, many of then are empty genes (i.e., not expressed by any cell), and many of the genes are likely to be noise due to extremely low expression among these 4340 cells.

-We filter out genes that are likely to be noise. We normally keep genes that have expressed at least 3 counts in at least 3 cells.

-After filtering, we have 4529 genes left.

-```{r}

-#pbmc_filterg = pbmc[rowSums(pbmc_counts > 2) > 2, ] 

-pbmc_filterg = pbmc[rowSums(pbmc_counts >= 3) >= 3, ] # Filtering can be applied directly on SingleCellExperiment object 

-print(pbmc_filterg)

-```

-

-

-

-# Estimate and remove contaminated counts using DecontX (2 scenarios -- with and without cell types specified)

-## Use DecontX for contamination estimation when there is cell type specified for each individual cells (Here we use celda_CG -- a Bayesian bi-clustering method on single cell RNAseq data -- to get cell clusters as well as gene modules, which enhances cell clustering. You can substitute `celda_CG` with any clustering method you have faith in. DecontX takes the expression count matrix (gene by cell count matrix) and a vector of the same length as number of cells in the count matrix specifying cell cluster. 

-```{r}

-pbmc_filterg_counts = SummarizedExperiment::assay(pbmc_filterg, "counts")

-pbmc_filterg_counts = as.matrix(pbmc_filterg_counts) # Change obejct type as basic matrix to pass into celda::celda_CG

-celda_res = celda_CG(counts = pbmc_filterg_counts, K = 19, L = 150) # specify number of cell clusters being K=19, number of gnee modules being L=150

-```

-P.s: we have analyzed this 4K PBMC data before, so we don't bother with the parameters choice for the bi-clustering method `celda_CG` here. For a detailed toturial of usage on `celda`, it is availabel [here] (celda-vignettes url) 

-

-

-

-Pass expression (gene by cell) count matrix and the cluster labels from `celda_CG` to decontX, it will return estimated results. Corresponding to the 2 types of input, output result is either in a `list`, or a SingleCellExperiment object with the result list added back to its `metadata` slot (?? Am I understanding the structure right? Is metadata a "slot" of the sce-object??). When `decontX` takes a SingleCellExperiment object, we can use function `metadata` to extract the result list and use `str` function to look at the 4 attributes ("runParams", "estimates", "contamination", "z") in it.

-```{r}

-cell_cluster = celda_res@clusters$z 

-decontx_res_sce = decontX( x = pbmc_filterg, z = cell_cluster) 

-print(decontx_res_sce)

-#decontx_res = decontX( x = pbmc_filterg_counts, z = cell_cluster) ## This works too. Difference is that the output is a list rathen than being a SingleCellExperiment object

-decontx_res = metadata(decontx_res_sce)$decontX

-str(decontx_res)

-```

-

-

-## Interpretation of DecontX results

-There are 4 attributes of DecontX results: 

-(1) runParams: parameters used for estimating contamination

-(2) estimates: parameters, decontaminated counts as well as other intermediate results estimated during contamination estimation

-(3) contamination: cell level contamination

-(4) z: cell cluster label 

-

-

-# Check DecontX decontaminated count to see if gene expressions are cleaner by looking at marker genes' expression

-## First we have to identify cell types, which is one step further from cell clustering, as we do not really know which cluster is corresponding to what cell type and multiple cell clusters can also be from the same cell type. We are looking at the 4 major cell types (T-cell, B-cell, NK-cell and monocytes) in this PBMC data.

-Combining both the clustering results from `celda_CG`, visualization on a reduced 2-dimention space using Tsne (T-stochastic Neighbor Embedding), and marker genes, we relate cell clusters (more specifically the cells in them) to cell types. The marker genes used to identify these 4 major cell types are CD3D (T-cell marker), MS4A1 and CD79A (B-cell markers), GNLY (NK-cell marker), CD14 or FCGR3A (monocytes markers). There are also other cell types (such as dendritic cells, plasma cells, etc.), we don't include these cell types in our demonstration here and we point you to [here] (celda-vignettes url on PBMC4k data) for cell type analysis using `celda`.

-```{r}

-celdaTsne.ori = celdaTsne( counts = pbmc_filterg_counts , celdaMod = celda_res ) 

-pltTsne.ori = plotDimReduceCluster( dim1 = celdaTsne.ori[,1], dim2 = celdaTsne.ori[, 2], cluster = cell_cluster, labelClusters = T)

-print(pltTsne.ori)

-```

-

-We use violin plot to look at the marker gene expression levels in each cluster identified before by `celda_CG`. Those gene names are stored in the "rowData" of the singleCellExperiment object. On the other hand, the rownames of the expression count matrix is using gene ID. We have to match marker genes' name to the gene ID to locate them in the expression count matrix.

-```{r}

-markerGenes = do.call(what = rbind, args = list(data.frame("celltype" = "Tcell", "markergenes"=c("CD3D" )), 

-		                                data.frame("celltype" = "Bcell", "markergenes"=c("MS4A1", "CD79A")),

-			                        data.frame("celltype" = "NKcell", "markergenes"=c("GNLY")),

-			                        data.frame("celltype" = "monocytes", "markergenes"=c("CD14", "FCGR3A")),

-						data.frame("celltype" = "megakaryocytes", "markergenes"=c("PPBP"))))

-print(markerGenes)

-head(rowData(pbmc_filterg))

-markerGenes_wID = merge( x = markerGenes, y = rowData(pbmc_filterg), by.x = "markergenes", by.y = "Symbol")

-markerGenes_wID[, "feature"] = paste( markerGenes_wID[, "ENSEMBL_ID"], markerGenes_wID[, "markergenes"], sep = "_")

-print(markerGenes_wID)

-

-pbmc_filterg_counts_rename = pbmc_filterg_counts

-rownames(pbmc_filterg_counts_rename) = paste(rownames(pbmc_filterg), rowData(pbmc_filterg)[, "Symbol"], sep = "_")

-pltViolin.ori = violinPlot(counts = pbmc_filterg_counts_rename, celdaMod=celda_res, features = markerGenes_wID[, "feature"])

-print(pltViolin.ori)

-```

-

-The violin plots of the marker genes' expression indicate that cluster 12, 13 are B-cells, cluster 7, 14, 15, 16, 18, 19 are T-cells, cluster 11 is NK-cells, cluster 2, 4, 5 are monocytes. For monocytes, we did not include cluster 3 as it is more likely to be megakaryocytes (PPBP+) or cluster 6 as it is likely to be doublets group (CD14+ and CD3D).

-```{r}

-cellType = list("Bcell"=c(12,13), "Tcell"=c(7, 14, 15, 16, 18, 19), "NKcell"=c(11), "monocytes"=c(2, 4, 5))

-print(cellType)

-```

-

-Now that we have identified the 4 major cell types, we can go ahead to see how DecontX makes the genes' expression clearer by comparing the decontaminated counts with the orinal counts. As we now look at properties at cell-type level than cell cluster level, again we will convert cell cluster labels into cell types. P.s: except for the aforementioned major 4 cell types, all other cell types in this PBMC data are indicated as "others" in this tutorial.

-```{r}

-# This function is to convert cluster label into cell-type label

-match.z.celltype = function( z, cellType.list   )  { 

-    z.cell = z

-    for ( i in 1:length(cellType.list) )  { 

-        z.cell[ z %in% cellType.list[[i]]  ] =  names( cellType.list[i] ) 

-    }

-    return( z.cell ) 

-}

-

-cellType_list = cellType

-cellType_list[["others"]] = c(1:19)[-unlist(cellType)]

-print(cellType_list)

-

-cell_type = match.z.celltype(z = cell_cluster, cellType.list = cellType_list)

-print(table(cell_cluster, cell_type))

-```

-

-We look at the gene expressions for these 4 major cell types by comparing the original expression count and the decontaminated counts by DecontX. Since dropout events are prevalent in single cell RNAseq data, we include more genes from now to show clearer cell-type expressions.

-```{r}

-ct_markergenes = do.call(what = rbind, args = list(data.frame("celltype" = "Tcell", "markergenes"=c("CD3D", "CD3E")),

-																									data.frame("celltype" = "Bcell", "markergenes"=c("MS4A1", "CD79A", "CD79B")),

-																									data.frame("celltype" = "NKcell", "markergenes"=c("GNLY")),

-																									data.frame("celltype" = "monocytes", "markergenes"=c("LYZ", "S100A8", "S100A9"))))

-ct_markergenes_wID = merge( x = ct_markergenes, y = rowData(pbmc_filterg), by.x = "markergenes", by.y = "Symbol")

-ct_markergenes_wID[, "feature"] = paste( ct_markergenes_wID[, "ENSEMBL_ID"], ct_markergenes_wID[, "markergenes"], sep = "_")

-print(head(ct_markergenes_wID))

-

-ct_markergenes = ct_markergenes_wID[, c("celltype", "feature")]

-colnames(ct_markergenes) = c("cellType", "geneMarkers")

-

-pltMarker = decontXMarkerPlot( counts = pbmc_filterg_counts_rename, z = cell_type, geneMarkers = ct_markergenes) 

-

-

-```

-

-

-## Use DecontX for contamination estimation when there is no cell type specified for the cells

-

-

-

-

-

deleted file mode 100644

@@ -1,99 +0,0 @@

-title: "Estimate and remove cross-contamination from ambient RNA for scRNA-seq data with DecontX"

-author: "Shiyi Yang, Sean Corbett, Yusuke Koga, Zhe Wang, W. Evan Johnson, Masanao Yajima, Joshua D. Campbell"

-date: "`r Sys.Date()`"

-output:

-  BiocStyle::html_document:

-    toc: true

-vignette: >

-  %\VignetteIndexEntry{Estimate and remove cross-contamination from ambient RNA for scRNA-seq data with DecontX}

-  %\VignetteEngine{knitr::rmarkdown}

-  %\VignetteEncoding{UTF-8}

-

-# Introduction

-

-DecontX is a Bayesian hierarchical model to estimate and remove cross-contamination from ambient RNA in single-cell RNA-seq count data generated from droplet-based sequencing devices. DecontX will take the count matrix with/without the cell labels and estimate the contamination level and deliver a decontaminted count matrix for downstream analysis. 

-

-In this vignette we will demonstrate how to use DecontX to estimate and remove contamination.  

-

-# Installation

-

-celda can be installed from Bioconductor:

-

-```{r, eval= FALSE}

-if (!requireNamespace("BiocManager", quietly=TRUE)) {

-    install.packages("BiocManager")}

-BiocManager::install("celda")

-```

-

-The package can be loaded using the `library` command.

-

-```{r, eval=TRUE, message=FALSE} 

-library(celda)

-```

-

-To see the latest updates and releases or to post a bug, see our GitHub page at https://github.com/campbio/celda. To ask questions about running celda, post a thread on Bioconductor support site at https://support.bioconductor.org/.

-

-# Reproducibility note

-

-Many functions in *celda* make use of stochastic algorithms or procedures which require the use of random number generator (RNG) for simulation or sampling. To maintain reproducibility, all these functions use a **default seed of 12345** to make sure same results are generated each time one of these functions is called. Explicitly setting the `seed` arguments is needed for greater control and randomness.

-

-# Generation of a cross-contaminated dataset 

-DecontX will take a matrix of counts (referred as observed counts) where each row is a feature, each column is a cell, and each entry in the matrix is the number of counts of each feature in each cell. To illustrate the utility of DecontX, we will apply it to a simulated dataset.

-

-In the function `simulateContaminatedMatrix`, the K parameter designates the number of cell clusters, the C parameter determines the number of cells, the G parameter determines the number of genes in the simulated dataset.

-

-```{r}

-simCounts <- simulateContaminatedMatrix(G = 300, C = 100, K = 3)

-```

-

-The `nativeCounts` is the natively expressed counts matrix, and `observedCounts` is the observed counts matrix that contains both contaminated and natively expressed transctripts. The `NByC` is the total number of observed transcripts per cell. The counts matrix which only contains contamianted transcripts can be obtained by subtracting the observed counts matrix from the observed counts matrix. 

-

-```{r}

-contamination <- simCounts$observedCounts - simCounts$nativeCounts

-```

-

-The `z` variable contains the population label for each cell.

-

-```{r}

-table(simCounts$z)

-```

-

-The `phi` and `eta` variables contain the expression distributions and contamination distributions for each population, respectively. Each column corresponds to a population, each row represents a gene. The sum of the rows equal to 1.

-

-```{r}

-colSums(simCounts$phi)

-colSums(simCounts$eta)

-```

-

-

-# Decontamination using DecontX

-DecontX uses bayesian method to estimate and remove contamination via varitaional inference.

-

-```{r, warning = FALSE, message = FALSE}

-decontxModel <- decontX(x = simCounts$observedCounts, z = simCounts$z)

-```

-

-## Check convergence

-Use log-likelihood to check convergence

-

-```{r, eval = TRUE, fig.width = 5, fig.height = 5}

-plot(decontxModel$estimates$all_cells$logLikelihood)

-```

-

-## Evaluate model performance

-

-`decontX` estimates a contamination proportion for each cell. We compare the estimated contamination proportion with the real contamination proportion.

-

-```{r, eval = TRUE, fig.width = 5, fig.height = 5}

-plot(decontxModel$contamination,

-    colSums(contamination) / simCounts$NByC, col = simCounts$z)

-abline(0, 1)

-```

-

-# Session Information

-

-```{r}

-sessionInfo()

-```

similarity index 100%

rename from vignettes/FindMarkers-analysis.Rmd

rename to vignettes/FindMarkers.Rmd

similarity index 96%

rename from vignettes/celda-analysis.Rmd

rename to vignettes/celda.Rmd

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

 ---

-title: "Analysis of single-cell 'omic count data with celda"

-author: "Sean Corbett, Yusuke Koga, Shiyi Yang, Zhe Wang, Evan Johnson, Paola Sebastiani, Masanao Yajima, Joshua Campbell"

+title: "Analysis of single-cell genomic data with celda"

+author:

+- name: Sean Corbett

+  affiliation: &id Boston University School of Medicine

+- name: Yusuke Koga

+  affiliation: *id

+- name: Shiyi Yang

+  affiliation: *id

+- name: Zhe Wang

+  affiliation: *id

+- name: Joshua Campbell

+  affiliation: *id

+  email: camp@bu.edu

 date: "`r Sys.Date()`"

 output: 

   BiocStyle::html_document:

@@ -13,9 +24,9 @@ vignette: >

 # Introduction

-**CE**llular **L**atent **D**irichlet **A**llocation (celda) is a collection of Bayesian hierarchical models to perform feature and cell bi-clustering for count data generated by single-cell platforms. This algorithm is an extension of the Latent Dirichlet Allocation (LDA) topic modeling framework that has been popular in text mining applications and has shown good performance with sparse data. celda simultaneously clusters features (i.e. gene expression) into modules based on co-expression patterns across cells and cells into subpopulations based on the probabilities of the feature modules within each cell. celda uses Dirichlet-multinomial distributions to model cells and genes so no additional normalization is required for single-cell counts. 

+**CE**llular **L**atent **D**irichlet **A**llocation (celda) is a collection of Bayesian hierarchical models to perform feature and cell bi-clustering for count data generated by single-cell platforms. This algorithm is an extension of the Latent Dirichlet Allocation (LDA) topic modeling framework that has been popular in text mining applications and has shown good performance with sparse data. celda simultaneously clusters features (i.e. gene expression) into modules based on co-expression patterns across cells and cells into subpopulations based on the probabilities of the feature modules within each cell.  

-In this vignette we will demonstrate how to use celda to perform cell and feature clustering with simulated data.

+In this vignette we will demonstrate how to use celda to perform cell and feature clustering with a simple simulated dataset.

 # Installation

@@ -307,17 +318,18 @@ celdaModel <- subsetCeldaList(celdaList = cellSplit,

 ## celdaGridSearch

-celda is able to run multiple combinations of K and L with multiple chains in parallel via the `celdaGridSearch` function. 

+celda is able to run multiple combinations of K and L with multiple chains in parallel via the `celdaGridSearch` function. Setting `verbose` to `TRUE` will print the output of each model to a text file.

 The `resamplePerplexity` function "perturbs" the original counts matrix by resampling the counts of each cell according to its normalized probability distribution. Perplexity is calculated on the resampled matrix and the procedure is repeated `resample` times. These results can be visualized with `plotGridSearchPerplexity`. The major goal is to pick the lowest K and L combination with relatively good perplexity. In general, visual inspection of the plot can be used to select the number of modules (L) or cell populations (K) where the rate of decrease in the perplexity starts to drop off.

 ```{r, eval = TRUE, message = FALSE}

 cgs <- celdaGridSearch(simCounts$counts,

     paramsTest = list(K = seq(4, 6), L = seq(9, 11)),

-    cores = 4,

+    cores = 1,

     model = "celda_CG",

     nchains = 2,

     maxIter = 100,

+    verbose = FALSE,

     bestOnly = TRUE)

 ```

@@ -347,10 +359,11 @@ If the "bestOnly" parameter is set to FALSE in the `celdaGridSearch`, then the `

 ```{r, message=FALSE}

 cgs <- celdaGridSearch(simCounts$counts,

     paramsTest = list(K = seq(4, 6), L = seq(9, 11)),

-    cores = 4,

+    cores = 1,

     model = "celda_CG",

     nchains = 2,

     maxIter = 100,

+    verbose = FALSE,

     bestOnly = FALSE)

 cgs <- resamplePerplexity(counts = simCounts$counts,

new file mode 100644

@@ -0,0 +1,58 @@

+---

+title: "Decontamination of ambient RNA in single-cell RNA-seq with DecontX on 10X PBMC4k data" 

+author: "Shiyi Yang, Joshua D. Campbell"

+date: "`r Sys.Date()`"

+output:

+  BiocStyle::html_document:

+    toc: true

+vignette: >

+  %\VignetteIndexEntry{Estimate and remove cross-contamination from ambient RNA in single-cell data with DecontX}

+  %\VignetteEngine{knitr::rmarkdown}

+  %\VignetteEncoding{UTF-8}

+---

+

+# Installation

+

+celda can be installed from Bioconductor:

+

+```{r, eval= FALSE}

+if (!requireNamespace("BiocManager", quietly=TRUE)) {

+    install.packages("BiocManager")}

+BiocManager::install("celda")

+```

+

+The package can be loaded using the `library` command.

+

+```{r, eval=TRUE, message=FALSE}

+library(celda)

+```

+

+DecontX can take either `SingleCellExperiment` object from package [SingleCellExperiment package](https://bioconductor.org/packages/release/bioc/html/SingleCellExperiment.html) or a single counts matrix as input. `decontX` will attempt to convert any input matrix to class `dgCMatrix` from package [Matrix](https://cran.r-project.org/web/packages/Matrix/index.html) before beginning any analyses.

+

+

+# Load PBMC4k data from 10X

+We will utlize the 10X PBMC 4K dataset as an example. This can be easily retrieved from the package [TENxPBMCData](http://bioconductor.org/packages/release/data/experiment/html/TENxPBMCData.html). Make sure the the column names are set before running decontX.

+

+```{r, eval=TRUE, message=FALSE}

+library(TENxPBMCData)  

+pbmc4k <- TENxPBMCData("pbmc4k")

+colnames(pbmc4k) <- paste(pbmc4k$Sample, pbmc4k$Barcode, sep="_")

+```

+

+

+To run decontX with a SingleCellExperiment object, simply use the following command. 

+

+```{r}

+pbmc4k = decontX(x = pbmc4k) 

+```

+

+

+

+# Examining decontX output

+There are 4 attributes of DecontX results: 

+(1) runParams: parameters used for estimating contamination

+(2) estimates: parameters, decontaminated counts as well as other intermediate results estimated during contamination estimation

+(3) contamination: cell level contamination

+(4) z: cell cluster label 

+

+# Plotting decontX results