inst/doc/additional_topics.R

In dsb: Normalize & Denoise Droplet Single Cell Protein Data (CITE-Seq)

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ---- eval=FALSE--------------------------------------------------------------
#  suppressMessages(library(SingleCellExperiment))
#  sce = SingleCellExperiment(assays = list(counts = cell.rna.raw), colData = cellmd)
#  # define the dsb normalized values as logcounts to use a common SingleCellExperiment / Bioconductor convention
#  adt = SummarizedExperiment(
#    assays = list(
#      'counts' = as.matrix(cell.adt.raw),
#      'logcounts' = as.matrix(cell.adt.dsb)
#      )
#    )
#  altExp(sce, "CITE") = adt

## ---- eval = FALSE------------------------------------------------------------
#  library(reticulate); sc = import("scanpy")
#  
#  # merge dsb-normalized protein and raw RNA data
#  combined_dat = rbind(cell.rna.raw, cell.adt.dsb)
#  s[["combined_data"]] = CreateAssayObject(data = combined_dat)
#  
#  # create Anndata Object
#  adata_seurat = sc$AnnData(
#      X   = t(GetAssayData(s,assay = "combined_data")),
#      obs = seurat@meta.data,
#      var = GetAssay(seurat)[[]]
#      )

## ---- eval = FALSE------------------------------------------------------------
#  # suggested workflow if isotype controls are not included
#  dsb_rescaled = DSBNormalizeProtein(cell_protein_matrix = cells_citeseq_mtx,
#                                     empty_drop_matrix = empty_drop_citeseq_mtx,
#                                     # do not denoise each cell's technical component
#                                     denoise.counts = FALSE)

## ---- eval = FALSE------------------------------------------------------------
#  
#  dsb_rescaled = dsb::DSBNormalizeProtein(cell_protein_matrix = cells_citeseq_mtx,
#                                          empty_drop_matrix = empty_drop_citeseq_mtx,
#                                          # denoise with background mean only
#                                          denoise.counts = TRUE,
#                                          use.isotype.control = FALSE)
#  

## ---- eval=FALSE--------------------------------------------------------------
#  # raw = Read10X see above -- path to cell ranger outs/raw_feature_bc_matrix ;
#  
#  # partial thresholding to slightly subset negative drops include all with 5 unique mRNAs
#  seurat_object = CreateSeuratObject(raw, min.genes = 5)
#  
#  # demultiplex (positive.quantile can be tuned to dataset depending on size)
#  seurat_object = HTODemux(seurat_object, assay = "HTO", positive.quantile = 0.99)
#  Idents(seurat_object) = "HTO_classification.global"
#  
#  # subset empty drop/background and cells
#  neg_object = subset(seurat_object, idents = "Negative")
#  singlet_object = subset(seurat_object, idents = "Singlet")
#  
#  ## (QC the negative object to filter out cells with high RNA content)
#  # quick example below, different crteria can be used
#  # this step depends on dataset; see main vignette for more principled filtering
#  neg_object = subset(seurat_object, idents = "Negative", nGene < 80)
#  
#  
#  # non sparse CITEseq data store more efficiently in a regular matrix
#  neg_adt_matrix = GetAssayData(neg_object, assay = "CITE", slot = 'counts') %>% as.matrix()
#  positive_adt_matrix = GetAssayData(singlet_object, assay = "CITE", slot = 'counts') %>% as.matrix()
#  
#  # normalize the data with dsb
#  dsb_norm_prot = DSBNormalizeProtein(
#                             cell_protein_matrix = cells_mtx_rawprot,
#                             empty_drop_matrix = negative_mtx_rawprot,
#                             denoise.counts = TRUE,
#                             use.isotype.control = TRUE,
#                             isotype.control.name.vec = rownames(cells_mtx_rawprot)[30:32])
#  
#  # now add the normalized dat back to the object (the singlets defined above as "object")
#  singlet_object[["CITE"]] = CreateAssayObject(data = dsb_norm_prot)
#  
#  # proceed with same tutorial workflow shown above.

## -----------------------------------------------------------------------------
library(dsb)
result.list = 
  DSBNormalizeProtein(
    cell_protein_matrix = cells_citeseq_mtx[ ,1:50], 
    empty_drop_matrix = empty_drop_citeseq_mtx, 
    denoise.counts = TRUE, 
    use.isotype.control = TRUE, 
    isotype.control.name.vec = rownames(cells_citeseq_mtx)[67:70], 
    return.stats = TRUE
  )


## -----------------------------------------------------------------------------
names(result.list)

## -----------------------------------------------------------------------------
names(result.list$protein_stats)

## -----------------------------------------------------------------------------
head(result.list$technical_stats)

## ---- eval = FALSE------------------------------------------------------------
#  dsb_norm_prot = DSBNormalizeProtein(
#                             cell_protein_matrix = cells_citeseq_mtx,
#                             empty_drop_matrix = empty_drop_citeseq_mtx,
#                             denoise.counts = TRUE,
#                             use.isotype.control = TRUE,
#                             isotype.control.name.vec = rownames(cells_citeseq_mtx)[67:70],
#                             # implement Quantile clipping
#                             quantile.clipping = TRUE
#                             # high and low otlier quantile across proteins to clip
#                             # the `quantile.clip` parameter can be adjusted:
#                             quantile.clip = c(0.001, 0.9995)
#                             )

## ---- eval = FALSE------------------------------------------------------------
#  
#  dsb_norm_prot = DSBNormalizeProtein(
#                             cell_protein_matrix = cells_citeseq_mtx,
#                             empty_drop_matrix = empty_drop_citeseq_mtx,
#                             denoise.counts = TRUE,
#                             use.isotype.control = TRUE,
#                             isotype.control.name.vec = rownames(cells_citeseq_mtx)[67:70],
#                             quantile.clipping = TRUE,
#                             scale.factor = 'mean.subtract'
#                             )
#  

## ---- eval = FALSE------------------------------------------------------------
#  # find outliers
#  pheatmap::pheatmap(apply(dsb_norm_prot, 1, function(x){
#    quantile(x,c(0.9999, 0.99, 0.98, 0.95, 0.0001, 0.01, 0.1))
#    }))
#  

## ---- eval=FALSE--------------------------------------------------------------
#  
#  dsb_object = DSBNormalizeProtein(cell_protein_matrix = dsb::cells_citeseq_mtx,
#                                   empty_drop_matrix = dsb::empty_drop_citeseq_mtx,
#                                   denoise.counts = TRUE,
#                                   isotype.control.name.vec = rownames(dsb::cells_citeseq_mtx)[67:70],
#                                   return.stats = TRUE)
#  d = as.data.frame(dsb_object$dsb_stats)
#  
#  # test correlation of background mean with the inferred dsb technical component
#  cor(d$cellwise_background_mean, d$dsb_technical_component)
#  
#  # test average isotype control value correlation with the background mean
#  isotype_names = rownames(dsb::cells_citeseq_mtx)[67:70]
#  cor(rowMeans(d[,isotype_names]), d$cellwise_background_mean)
#  

## ---- eval = FALSE------------------------------------------------------------
#  library(Seurat) # for Read10X helper function
#  
#  # path_to_reads = here("data/")
#  umi.files = list.files(path_to_reads, full.names=T, pattern = "10x" )
#  umi.list = lapply(umi.files, function(x) Read10X(data.dir = paste0(x,"/outs/raw_feature_bc_matrix/")))
#  prot = rna = list()
#  for (i in 1:length(umi.list)) {
#    prot[[i]] = umi.list[[i]]`Antibody Capture`
#    rna[[i]] = umi.list[[i]]`Gene Expression`
#    colnames(prot[[i]]) = paste0(colnames(prot[[i]]),"_", i )
#    colnames(rna[[i]]) = paste0(colnames(rna[[i]]),"_", i )
#  }
#  prot = do.call(cbind, prot)
#  rna = do.call(cbind, rna)
#  # proceed with step 1 in tutorial - define background and cell containing drops for dsb
#