inst/doc/end_to_end_workflow.R

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

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = FALSE,
  # fig.path = "man/figures/README-",
  out.width = "100%"
)

## ---- eval = FALSE------------------------------------------------------------
#  install.packages('dsb')
#  library(dsb)
#  
#  adt_norm = 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]
#    )

## ---- eval = FALSE------------------------------------------------------------
#  library(dsb)
#  
#  # read raw data using the Seurat function "Read10X"
#  raw = Seurat::Read10X("data/raw_feature_bc_matrix/")
#  cells = Seurat::Read10X("data/filtered_feature_bc_matrix/")
#  
#  # define cell-containing barcodes and separate cells and empty drops
#  stained_cells = colnames(cells$`Gene Expression`)
#  background = setdiff(colnames(raw$`Gene Expression`), stained_cells)
#  
#  # split the data into separate matrices for RNA and ADT
#  prot = raw$`Antibody Capture`
#  rna = raw$`Gene Expression`

## ---- eval = FALSE------------------------------------------------------------
#  # create metadata of droplet QC stats used in standard scRNAseq processing
#  mtgene = grep(pattern = "^MT-", rownames(rna), value = TRUE) # used below
#  
#  md = data.frame(
#    rna.size = log10(Matrix::colSums(rna)),
#    prot.size = log10(Matrix::colSums(prot)),
#    n.gene = Matrix::colSums(rna > 0),
#    mt.prop = Matrix::colSums(rna[mtgene, ]) / Matrix::colSums(rna)
#  )
#  # add indicator for barcodes Cell Ranger called as cells
#  md$drop.class = ifelse(rownames(md) %in% stained_cells, 'cell', 'background')
#  
#  # remove barcodes with no evidence of capture in the experiment
#  md = md[md$rna.size > 0 & md$prot.size > 0, ]
#  

## ---- eval = FALSE------------------------------------------------------------
#  ggplot(md, aes(x = log10(n.gene), y = prot.size )) +
#     theme_bw() +
#     geom_bin2d(bins = 300) +
#     scale_fill_viridis_c(option = "C") +
#     facet_wrap(~drop.class)

## ---- eval = FALSE------------------------------------------------------------
#  background_drops = rownames(
#    md[ md$prot.size > 1.5 &
#        md$prot.size < 3 &
#        md$rna.size < 2.5, ]
#    )
#  background.adt.mtx = as.matrix(prot[ , background_drops])

## ---- eval = FALSE------------------------------------------------------------
#  
#  # calculate statistical thresholds for droplet filtering.
#  cellmd = md[md$drop.class == 'cell', ]
#  
#  # filter drops with + / - 3 median absolute deviations from the median library size
#  rna.mult = (3*mad(cellmd$rna.size))
#  prot.mult = (3*mad(cellmd$prot.size))
#  rna.lower = median(cellmd$rna.size) - rna.mult
#  rna.upper = median(cellmd$rna.size) + rna.mult
#  prot.lower = median(cellmd$prot.size) - prot.mult
#  prot.upper = median(cellmd$prot.size) + prot.mult
#  
#  # filter rows based on droplet qualty control metrics
#  qc_cells = rownames(
#    cellmd[cellmd$prot.size > prot.lower &
#           cellmd$prot.size < prot.upper &
#           cellmd$rna.size > rna.lower &
#           cellmd$rna.size < rna.upper &
#           cellmd$mt.prop < 0.14, ]
#    )

## ---- eval = FALSE------------------------------------------------------------
#  length(qc_cells)

## ---- eval = FALSE------------------------------------------------------------
#  cell.adt.raw = as.matrix(prot[ , qc_cells])
#  cell.rna.raw = rna[ ,qc_cells]
#  cellmd = cellmd[qc_cells, ]

## ---- eval = FALSE------------------------------------------------------------
#  # flter
#  pm = sort(apply(cell.adt.raw, 1, max))
#  pm2 = apply(background.adt.mtx, 1, max)
#  head(pm)

## ---- eval = FALSE------------------------------------------------------------
#  # remove the non staining protein
#  cell.adt.raw  = cell.adt.raw[!rownames(cell.adt.raw) == 'CD34_TotalSeqB', ]
#  background.adt.mtx = background.adt.mtx[!rownames(background.adt.mtx) == 'CD34_TotalSeqB', ]
#  

## ---- eval = FALSE------------------------------------------------------------
#  # define isotype controls
#  isotype.controls = c("IgG1_control_TotalSeqB", "IgG2a_control_TotalSeqB",
#                       "IgG2b_control_TotalSeqB")
#  
#  # normalize and denoise with dsb with
#  cells.dsb.norm = DSBNormalizeProtein(
#    cell_protein_matrix = cell.adt.raw,
#    empty_drop_matrix = background.adt.mtx,
#    denoise.counts = TRUE,
#    use.isotype.control = TRUE,
#    isotype.control.name.vec = isotype.controls
#    )
#  # note: normalization takes ~ 20 seconds
#  # system.time()
#  # user  system elapsed
#  #  20.799   0.209  21.783

## ---- eval = FALSE------------------------------------------------------------
#  # dsb with non-standard options
#  cell.adt.dsb.2 = DSBNormalizeProtein(
#    cell_protein_matrix = cell.adt.raw,
#    empty_drop_matrix = background.adt.mtx,
#    denoise.counts = TRUE,
#    use.isotype.control = TRUE,
#    isotype.control.name.vec = rownames(cell.adt.raw)[29:31],
#    define.pseudocount = TRUE,
#    pseudocount.use = 1,
#    scale.factor = 'mean.subtract',
#    quantile.clipping = TRUE,
#    quantile.clip = c(0.01, 0.99),
#    return.stats = TRUE
#    )
#  
#  

## ---- eval = FALSE------------------------------------------------------------
#  # Seurat workflow
#  library(Seurat)
#  
#  # integrating with Seurat
#  stopifnot(isTRUE(all.equal(rownames(cellmd), colnames(cell.adt.raw))))
#  stopifnot(isTRUE(all.equal(rownames(cellmd), colnames(cell.rna.raw))))
#  
#  # create Seurat object note: min.cells is a gene filter, not a cell filter
#  s = Seurat::CreateSeuratObject(counts = cell.rna.raw,
#                                 meta.data = cellmd,
#                                 assay = "RNA",
#                                 min.cells = 20)
#  
#  # add dsb normalized matrix "cell.adt.dsb" to the "CITE" data (not counts!) slot
#  s[["CITE"]] = Seurat::CreateAssayObject(data = cells.dsb.norm)
#  

## ---- eval = FALSE------------------------------------------------------------
#  # define proteins to use in clustering (non-isptype controls)
#  prots = rownames(s@assays$CITE@data)[1:28]
#  
#  # cluster and run umap
#  s = Seurat::FindNeighbors(object = s, dims = NULL,assay = 'CITE',
#                            features = prots, k.param = 30,
#                            verbose = FALSE)
#  
#  # direct graph clustering
#  s = Seurat::FindClusters(object = s, resolution = 1,
#                           algorithm = 3,
#                           graph.name = 'CITE_snn',
#                           verbose = FALSE)
#  # umap (optional)
#  # s = Seurat::RunUMAP(object = s, assay = "CITE", features = prots,
#  #                     seed.use = 1990, min.dist = 0.2, n.neighbors = 30,
#  #                     verbose = FALSE)
#  
#  # make results dataframe
#  d = cbind(s@meta.data,
#            as.data.frame(t(s@assays$CITE@data))
#            # s@reductions$umap@cell.embeddings)
#            )

## ---- eval = FALSE------------------------------------------------------------
#  library(magrittr)
#  # calculate the median protein expression separately for each cluster
#  adt_plot = d %>%
#    dplyr::group_by(CITE_snn_res.1) %>%
#    dplyr::summarize_at(.vars = prots, .funs = median) %>%
#    tibble::remove_rownames() %>%
#    tibble::column_to_rownames("CITE_snn_res.1")
#  # plot a heatmap of the average dsb normalized values for each cluster
#  pheatmap::pheatmap(t(adt_plot),
#                     color = viridis::viridis(25, option = "B"),
#                     fontsize_row = 8, border_color = NA)
#  

## -----------------------------------------------------------------------------
#  # use pearson residuals as normalized values for pca
#  DefaultAssay(s) = "RNA"
#  s = NormalizeData(s, verbose = FALSE) %>%
#    FindVariableFeatures(selection.method = 'vst', verbose = FALSE) %>%
#    ScaleData(verbose = FALSE) %>%
#    RunPCA(verbose = FALSE)
#  
#  # set up dsb values to use in WNN analysis (do not normalize with CLR, use dsb normalized values)
#  DefaultAssay(s) = "CITE"
#  VariableFeatures(s) = prots
#  s = s %>%
#    ScaleData() %>%
#    RunPCA(reduction.name = 'apca', verbose = FALSE)
#  
#  # run WNN
#  s = FindMultiModalNeighbors(
#    s, reduction.list = list("pca", "apca"),
#    dims.list = list(1:30, 1:18),
#    modality.weight.name = "RNA.weight",
#    verbose = FALSE
#  )
#  
#  # cluster
#  s <- FindClusters(s, graph.name = "wsnn",
#                    algorithm = 3,
#                    resolution = 1.5,
#                    verbose = FALSE,
#                    random.seed = 1990)

## -----------------------------------------------------------------------------
#  ## WNN with dsb values
#  # use RNA pearson residuals as normalized values for RNA pca
#  DefaultAssay(s) = "RNA"
#  s = NormalizeData(s, verbose = FALSE) %>%
#    FindVariableFeatures(selection.method = 'vst', verbose = FALSE) %>%
#    ScaleData(verbose = FALSE) %>%
#    RunPCA(verbose = FALSE)
#  
#  
#  # set up dsb values to use in WNN analysis
#  DefaultAssay(s) = "CITE"
#  VariableFeatures(s) = prots
#  
#  # run true pca to initialize dr pca slot for WNN
#  ## Not used {
#  s = ScaleData(s, assay = 'CITE', verbose = FALSE)
#  s = RunPCA(s, reduction.name = 'pdsb',features = VariableFeatures(s), verbose = FALSE)
#  # }
#  
#  # make matrix of normalized protein values to add as dr embeddings
#  pseudo = t(GetAssayData(s,slot = 'data',assay = 'CITE'))[,1:29]
#  colnames(pseudo) = paste('pseudo', 1:29, sep = "_")
#  s@reductions$pdsb@cell.embeddings = pseudo
#  
#  # run WNN directly using dsb normalized values.
#  s = FindMultiModalNeighbors(
#    object = s,
#    reduction.list = list("pca", "pdsb"),
#    weighted.nn.name = "dsb_wnn",
#    knn.graph.name = "dsb_knn",
#    modality.weight.name = "dsb_weight",
#    snn.graph.name = "dsb_snn",
#    dims.list = list(1:30, 1:29),
#    verbose = FALSE
#  )
#  
#  # cluster based on the join RNA and protein graph
#  s = FindClusters(s, graph.name = "dsb_knn",
#                   algorithm = 3,
#                   resolution = 1.5,
#                   random.seed = 1990,
#                   verbose = FALSE)

## -----------------------------------------------------------------------------
#  # create multimodal heatmap
#  vf = VariableFeatures(s,assay = "RNA")
#  
#  # find marker genes for the joint clusters
#  Idents(s) = "dsb_knn_res.1.5"
#  DefaultAssay(s)  = "RNA"
#  rnade = FindAllMarkers(s, features = vf, only.pos = TRUE, verbose = FALSE)
#  gene_plot = rnade %>%
#    dplyr::filter(avg_log2FC > 1 ) %>%
#    dplyr::group_by(cluster) %>%
#    dplyr::top_n(3) %$% gene %>% unique
#  
#  
#  cite_data = GetAssayData(s,slot = 'data',assay = 'CITE') %>% t()
#  rna_subset = GetAssayData(s,assay = 'RNA',slot = 'data')[gene_plot, ] %>%
#    as.data.frame() %>%
#    t() %>%
#    as.matrix()
#  
#  # combine into dataframe
#  d = cbind(s@meta.data, cite_data, rna_subset)
#  
#  # calculate the median protein expression per cluster
#  dat_plot = d %>%
#    dplyr::group_by(dsb_knn_res.1.5) %>%
#    dplyr::summarize_at(.vars = c(prots, gene_plot), .funs = median) %>%
#    tibble::remove_rownames() %>%
#    tibble::column_to_rownames("dsb_knn_res.1.5")
#  

## ---- fig.height=4, fig.width=3-----------------------------------------------
#  # make a combined plot
#  suppressMessages(library(ComplexHeatmap)); ht_opt$message = FALSE
#  
#  # protein heatmap
#  # protein heatmap
#  prot_col = circlize::colorRamp2(breaks = seq(-1,25, by = 1),
#                                  colors = viridis::viridis(n = 27, option = "B"))
#  p1 = Heatmap(t(dat_plot)[prots, ],
#               name = "protein",
#               col = prot_col,
#               use_raster = T,
#               row_names_gp = gpar(color = "black", fontsize = 5)
#  )
#  
#  
#  # mRNA heatmap
#  mrna = t(dat_plot)[gene_plot, ]
#  rna_col = circlize::colorRamp2(breaks = c(-2,-1,0,1,2),
#                                 colors = colorspace::diverge_hsv(n = 5))
#  p2 = Heatmap(t(scale(t(mrna))),
#               name = "mRNA",
#               col = rna_col,
#               use_raster = T,
#               clustering_method_columns = 'average',
#               column_names_gp = gpar(color = "black", fontsize = 7),
#               row_names_gp = gpar(color = "black", fontsize = 5))
#  
#  
#  # combine heatmaps
#  ht_list = p1 %v% p2
#  draw(ht_list)
#  

## ---- eval = FALSE------------------------------------------------------------
#  library(Seurat)
#  umi = Read10X(data.dir = 'data/raw_feature_bc_matrix/')
#  k = 3
#  barcode.whitelist =
#    rownames(
#      CreateSeuratObject(counts = umi,
#                         min.features = k,  # retain all barcodes with at least k raw mRNA
#                         min.cells = 800, # this just speeds up the function by removing genes.
#                         )@meta.data
#      )
#  
#  write.table(barcode.whitelist,
#  file =paste0(your_save_path,"barcode.whitelist.tsv"),
#  sep = '\t', quote = FALSE, col.names = FALSE, row.names = FALSE)
#