inst/doc/dittoSeq.R

In dittoSeq: User Friendly Single-Cell and Bulk RNA Sequencing Visualization

## ---- echo=FALSE, results="hide", message=FALSE-------------------------------
knitr::opts_chunk$set(error=FALSE, message=FALSE, warning=FALSE,
    dev="jpeg", dpi = 72, fig.width = 4.5, fig.height = 3.5)
library(BiocStyle)

## ---- eval=FALSE--------------------------------------------------------------
#  # Install BiocManager if needed
#  if (!requireNamespace("BiocManager", quietly = TRUE))
#      install.packages("BiocManager")
#  
#  # Install dittoSeq
#  BiocManager::install("dittoSeq")

## -----------------------------------------------------------------------------
library(dittoSeq)
library(scRNAseq)
library(SingleCellExperiment)
library(Seurat)
# Download data
sce <- BaronPancreasData()
# Trim to only 5 of the cell types for simplicity of vignette
sce <- sce[,meta("label",sce) %in% c(
    "acinar", "beta", "gamma", "delta", "ductal")]

## -----------------------------------------------------------------------------
# Make Seurat and grab metadata
seurat <- CreateSeuratObject(counts(sce))
seurat <- AddMetaData(seurat, sce$label, col.name = "celltype")
seurat <- AddMetaData(seurat, sce$donor, col.name = "Sample")
seurat <- AddMetaData(seurat,
                      PercentageFeatureSet(seurat, pattern = "^MT"),
                      col.name = "percent.mt")
# Basic Seurat workflow (possibly outdated, but fine for this vignette)
seurat <- NormalizeData(seurat, verbose = FALSE)
seurat <- FindVariableFeatures(object = seurat, verbose = FALSE)
seurat <- ScaleData(object = seurat, verbose = FALSE)
seurat <- RunPCA(object = seurat, verbose = FALSE)
seurat <- RunTSNE(object = seurat)
seurat <- FindNeighbors(object = seurat, verbose = FALSE)
seurat <- FindClusters(object = seurat, verbose = FALSE)

## -----------------------------------------------------------------------------
# Grab PCA, TSNE, clustering, log-norm data, and metadata to sce

# sce <- as.SingleCellExperiment(seurat)
# At the time this part of the vignette was made, the above function gave
# warnings. So... manual method
sce <- addDimReduction(
  sce, embeddings = Embeddings(seurat, reduction = "pca"), name = "PCA")
sce <- addDimReduction(
  sce, embeddings = Embeddings(seurat, reduction = "tsne"), name = "TSNE")
sce$idents <- seurat$seurat_clusters
assay(sce, "logcounts") <- GetAssayData(seurat)
sce$percent.mt <- seurat$percent.mt
sce$celltype <- seurat$celltype
sce$Sample <- seurat$Sample

## -----------------------------------------------------------------------------
dittoDimPlot(seurat, "Sample")
dittoPlot(seurat, "ENO1", group.by = "celltype")
dittoBarPlot(sce, "celltype", group.by = "Sample")

## -----------------------------------------------------------------------------
# First, let's make some mock expression and conditions data
exp <- matrix(rpois(20000, 5), ncol=20)
colnames(exp) <- paste0("sample", seq_len(ncol(exp)))
rownames(exp) <- paste0("gene", seq_len(nrow(exp)))
logexp <- logexp <- log2(exp + 1)

conditions <- factor(rep(1:4, 5))
sex <- c(rep("M", 9), rep("F", 11))

## -----------------------------------------------------------------------------
# Import
myRNA <- importDittoBulk(
    # x can be a DGEList, a DESeqDataSet, a SummarizedExperiment,
    #   or a list of data matrices
    x = list(counts = exp,
             logcounts = logexp),
    # Optional inputs:
    #   For adding metadata
    metadata = data.frame(conditions = conditions,
                          sex = sex),
    #   For adding dimensionality reductions
    reductions = list(pca = matrix(rnorm(20000), nrow=20)))

## -----------------------------------------------------------------------------
# Add metadata (metadata can alternatively be added in this way)
myRNA$conditions <- conditions
myRNA$sex <- sex

# Add dimensionality reductions (can alternatively be added this way)
# (We aren't actually calculating PCA here.)
myRNA <- addDimReduction(
    object = myRNA,
    embeddings = matrix(rnorm(20000), nrow=20),
    name = "pca",
    key = "PC")

## -----------------------------------------------------------------------------
dittoDimPlot(myRNA, "sex", size = 3, do.ellipse = TRUE)
dittoBarPlot(myRNA, "sex", group.by = "conditions")
dittoBoxPlot(myRNA, "gene1", group.by = "sex")
dittoHeatmap(myRNA, getGenes(myRNA)[1:10],
    annot.by = c("conditions", "sex"))

## -----------------------------------------------------------------------------
# Retrieve all metadata slot names
getMetas(seurat)
# Query for the presence of a metadata slot
isMeta("nCount_RNA", seurat)
# Retrieve metadata values:
meta("celltype", seurat)[1:10]
# Retrieve unique values of a metadata
metaLevels("celltype", seurat)

## -----------------------------------------------------------------------------
# Retrieve all gene names
getGenes(seurat)[1:10]
# Query for the presence of a gene(s)
isGene("CD3E", seurat)
isGene(c("CD3E","ENO1","INS","non-gene"), seurat, return.values = TRUE)
# Retrieve gene expression values:
gene("ENO1", seurat)[1:10]

## -----------------------------------------------------------------------------
# Retrieve all dimensionality reductions
getReductions(seurat)

## -----------------------------------------------------------------------------
# Getter
isBulk(sce)
isBulk(myRNA)

# Setter
mock_bulk <- setBulk(sce) # to bulk
mock_sc <- setBulk(myRNA, set = FALSE) # to single-cell

## ---- results = "hold"--------------------------------------------------------
dittoDimPlot(seurat, "celltype", reduction.use = "pca")
dittoDimPlot(sce, "ENO1")

## ---- results = "hold"--------------------------------------------------------
dittoScatterPlot(
    object = sce,
    x.var = "PPY", y.var = "INS",
    color.var = "celltype")
dittoScatterPlot(
    object = seurat,
    x.var = "nCount_RNA", y.var = "nFeature_RNA",
    color.var = "percent.mt")

## -----------------------------------------------------------------------------
dittoDimPlot(seurat, "ident",
             
             do.label = TRUE,
             labels.repel = FALSE,
             
             add.trajectory.lineages = list(
                 c("9","3"),
                 c("8","7","2","4"),
                 c("8","7","1"),
                 c("5","11","6"),
                 c("10","0")),
             trajectory.cluster.meta = "ident")

## ---- results = "hold"--------------------------------------------------------
dittoDimHex(seurat)
dittoScatterHex(seurat,
    x.var = "PPY", y.var = "INS")

## ---- results = "hold"--------------------------------------------------------
dittoDimHex(seurat, "INS")
dittoScatterHex(
    object = seurat,
    x.var = "PPY", y.var = "INS",
    color.var = "celltype",
    colors = c(1:4,7), max.density = 15)

## ---- results = "hold"--------------------------------------------------------
dittoPlot(seurat, "ENO1", group.by = "celltype",
    plots = c("vlnplot", "jitter"))
dittoRidgePlot(sce, "ENO1", group.by = "celltype")
dittoBoxPlot(seurat, "ENO1", group.by = "celltype")

## -----------------------------------------------------------------------------
dittoPlot(seurat, "ENO1", group.by = "celltype",
    plots = c("jitter", "vlnplot", "boxplot"), # <- order matters
    
    # change the color and size of jitter points
    jitter.color = "blue", jitter.size = 0.5,
    
    # change the outline color and width, and remove the fill of boxplots
    boxplot.color = "white", boxplot.width = 0.1,
    boxplot.fill = FALSE,
    
    # change how the violin plot widths are normalized across groups
    vlnplot.scaling = "count"
    )

## ---- results = "hold"--------------------------------------------------------
dittoBarPlot(seurat, "celltype", group.by = "Sample")
dittoBarPlot(seurat, "ident", group.by = "Sample",
    scale = "count")

## ---- results = "hold"--------------------------------------------------------
# Pick Genes
genes <- c("SST", "REG1A", "PPY", "INS", "CELA3A", "PRSS2", "CTRB1",
    "CPA1", "CTRB2" , "REG3A", "REG1B", "PRSS1", "GCG", "CPB1",
    "SPINK1", "CELA3B", "CLPS", "OLFM4", "ACTG1", "FTL")

# Annotating and ordering cells by some meaningful feature(s):
dittoHeatmap(seurat, genes,
    annot.by = c("celltype", "Sample"))
dittoHeatmap(seurat, genes,
    annot.by = c("celltype", "Sample"),
    order.by = "Sample")

## -----------------------------------------------------------------------------
# Add annotations
dittoHeatmap(seurat, genes,
    annot.by = c("celltype", "Sample"),
    scaled.to.max = TRUE,
    show_colnames = FALSE,
    show_rownames = FALSE)

## -----------------------------------------------------------------------------
# Highlight certain genes
dittoHeatmap(seurat, genes, annot.by = c("celltype", "Sample"),
    highlight.features = genes[1:3],
    complex = TRUE)

## -----------------------------------------------------------------------------
# Idents(seurat) <- "celltype"
# delta.marker.table <- FindMarkers(seurat, ident.1 = "delta")
# delta.genes <- rownames(delta.marker.table)[1:20]
# Idents(seurat) <- "seurat_clusters"

delta.genes <- c(
    "SST", "RBP4", "LEPR", "PAPPA2", "LY6H",
    "CBLN4", "GPX3", "BCHE", "HHEX", "DPYSL3",
    "SERPINA1", "SEC11C", "ANXA2", "CHGB", "RGS2",
    "FXYD6", "KCNIP1", "SMOC1", "RPL10", "LRFN5")

## -----------------------------------------------------------------------------
dittoDotPlot(seurat, vars = delta.genes, group.by = "celltype")
dittoDotPlot(seurat, vars = delta.genes, group.by = "celltype",
    scale = FALSE)

## -----------------------------------------------------------------------------
multi_dittoPlot(seurat, delta.genes[1:6], group.by = "celltype",
    vlnplot.lineweight = 0.2, jitter.size = 0.3)
dittoPlotVarsAcrossGroups(seurat, delta.genes, group.by = "celltype",
    main = "Delta-cell Markers")

## ---- results = "hold"--------------------------------------------------------
multi_dittoDimPlot(seurat, delta.genes[1:6])
multi_dittoDimPlotVaryCells(seurat, delta.genes[1],
    vary.cells.meta = "celltype")

## -----------------------------------------------------------------------------
# Original
dittoBarPlot(seurat, "celltype", group.by = "Sample", scale = "count")

# First 10 cells
dittoBarPlot(seurat, "celltype", group.by = "Sample", scale = "count",
    # String method
    cells.use = colnames(seurat)[1:10]
    # Index method, which would achieve the same effect
    # cells.use = 1:10
    )

# Acinar cells only
dittoBarPlot(seurat, "celltype", group.by = "Sample", scale = "count",
    # Logical method
    cells.use = meta("celltype", seurat) == "acinar")

## -----------------------------------------------------------------------------
dittoDimPlot(seurat, "celltype", split.by = "Sample")
dittoDimPlot(seurat, "ENO1", split.by = c("Sample", "celltype"))

## -----------------------------------------------------------------------------
dittoBarPlot(seurat, "celltype", group.by = "Sample",
    main = "Encounters",
    sub = "By Type",
    xlab = NULL, # NULL = remove
    ylab = "Generation 1",
    x.labels = c("Ash", "Misty", "Jessie", "James"),
    legend.title = "Types",
    var.labels.rename = c("Fire", "Water", "Grass", "Electric", "Psychic"),
    x.labels.rotate = FALSE)

## ---- results="hold"----------------------------------------------------------
# original - discrete
dittoDimPlot(seurat, "celltype")
# swapped colors
dittoDimPlot(seurat, "celltype",
    colors = 5:1)
# different colors
dittoDimPlot(seurat, "celltype",
    color.panel = c("red", "orange", "purple", "yellow", "skyblue"))

## ---- results="hold"----------------------------------------------------------
# original - expression
dittoDimPlot(seurat, "INS")
# different colors
dittoDimPlot(seurat, "INS",
    max.color = "red", min.color = "gray90")

## -----------------------------------------------------------------------------
dittoBarPlot(seurat, "celltype", group.by = "Sample",
    data.out = TRUE)

## -----------------------------------------------------------------------------
dittoHeatmap(seurat, c("SST","CPE","GPX3"), cells.use = colnames(seurat)[1:5],
    data.out = TRUE)

## ---- eval = FALSE------------------------------------------------------------
#  # These can be finicky to render in knitting, but still, example code:
#  dittoDimPlot(seurat, "INS",
#      do.hover = TRUE,
#      hover.data = c("celltype", "Sample", "ENO1", "ident", "nCount_RNA"))
#  dittoBarPlot(seurat, "celltype", group.by = "Sample",
#      do.hover = TRUE)

## -----------------------------------------------------------------------------
# Note: dpi gets re-set by the styling code of this vignette, so this is
#   just a code example, but the plot won't be quite matched.
dittoDimPlot(sce, "celltype",
    do.raster = TRUE,
    raster.dpi = 300)

## -----------------------------------------------------------------------------
dittoHeatmap(seurat, genes, scaled.to.max = TRUE,
    complex = TRUE,
    use_raster = TRUE)

## -----------------------------------------------------------------------------
sessionInfo()