In maximilian-heeg/magieR: Run `magic-impute` on scRNA Seq dataset in R

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8,
  fig.height = 6
)

This vignette containes some examples on how to use the magieR wrapper to run magic on a matrix, a SingleCellExperiment and a Seurat object.

library(magieR)

Plain `magic` Example

Load the example dataset that is provided with the magic-impute package. (https://github.com/KrishnaswamyLab/MAGIC/blob/master/data/test_data.csv)

m <- utils::read.csv(system.file("extdata", "test_data.csv", package = "magieR"))
m[1:10, 1:5]

Then run magic on this matrix. It's just as easy as that. You don't have to worry about any python settings and configurations.

result <- magieR(m)
# as.data.frame makes the output look a little nicer
as.data.frame(result)[1:10, 1:5]

Use `magic` on SingleCellExperiment

You can also use magic on a SingleCellExperiment. Let's load the exemplary data from the scater vignette.

suppressPackageStartupMessages({
  library(scRNAseq)  
  library(scater)
  library(SingleCellExperiment)
})

example_sce <- ZeiselBrainData()

example_sce <- addPerCellQC(example_sce, 
    subsets=list(Mito=grep("mt-", rownames(example_sce))))

plotColData(example_sce, x = "sum", y="detected", colour_by="tissue") 

example_sce <- logNormCounts(example_sce) 

vars <- getVarianceExplained(example_sce, 
    variables=c("tissue", "total mRNA mol", "sex", "age"))
head(vars)

example_sce <- runPCA(example_sce)

set.seed(1000)
example_sce <- runTSNE(example_sce, perplexity=50, 
    dimred="PCA", n_dimred=10)
example_sce

Now we can run magieR. It uses the logcounts by default. Other assays can be used with assay=. The function returns a sce-object with a new alternative Expression. Alternative Expressions are used because if you run magic on a subset of genes with genes= the resulting dimensions of the matrix do not match the dimensions of the sce-object and therefore cannot be stored in assays.

# run magieR
example_sce <- magieR(example_sce)
example_sce

# the magic data can be accessed with
altExp(example_sce, "magic")

Let's display some genes with and without magic

plotExpression(example_sce, rownames(vars)[1:5]) + 
  ggplot2::ggtitle("Without magic")

plotExpression(altExp(example_sce, "magic"), rownames(vars)[1:5]) + 
  ggplot2::ggtitle("With magic")

Or on a TSNE

plotTSNE(example_sce, colour_by="Tspan12") + 
  ggplot2::ggtitle("Without magic")

# not sure if there is a more elegant way
# copy reduced dim to alt exp
reducedDim(altExp(example_sce, "magic"), "TSNE") <- reducedDim(example_sce, "TSNE")
plotTSNE(altExp(example_sce, "magic"), colour_by="Tspan12") + 
  ggplot2::ggtitle("With magic")

Use `magic` on a `Seurat` object

## remove all variables
rm(list=ls())

First we follow the standard Seurat workflow as described here.

suppressPackageStartupMessages({
  library(BiocFileCache)
  library(Seurat)
})

url = "https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz"

bfc <- BiocFileCache()
data <- bfcrpath(bfc, url)

tempdir <- file.path(tempdir(), "10_pbmc")
untar(data, exdir = tempdir)

pbmc.data <- Read10X(data.dir = file.path(tempdir, "filtered_gene_bc_matrices/hg19/"))

pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

pbmc <- NormalizeData(pbmc)

pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

top10 <- head(VariableFeatures(pbmc), 10)

all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)

pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
pbmc <- RunUMAP(pbmc, dims = 1:10)

pbmc

Okay, let's run magic

pbmc <- magieR(pbmc, slot="data")
pbmc

Again, make some plots.

VlnPlot(pbmc,
        features = c("CD8A", "STAT3"),
        same.y.lims = FALSE)

VlnPlot(pbmc,
        features = c("CD8A", "STAT3"),
        same.y.lims = FALSE,
        assay="magic")

FeaturePlot(pbmc, features = c("CD8A", "STAT3"))


# set the default assay to magic
Seurat::DefaultAssay(pbmc) <- "magic"
FeaturePlot(pbmc, features = c("CD8A", "STAT3"))
# set if back to default ("RNA"). You only want to use magic for visualization
# and never for statistics.
Seurat::DefaultAssay(pbmc) <- "RNA"