Mouse Muscle Stem Cell Project

Tabula Muris analysis on stem cells

Génome Québec Innovation Centre, McGill University, Montreal, Canada

Prepare analysis workflow

Set filepaths and parameters

set.seed(42)
knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
options(
  readr.show_progress = FALSE,
  digits = 2
)

Load packages

suppressPackageStartupMessages({
  library(scater)
  library(scran)
  library(Seurat)
  library(tidyverse)
  library(BiocParallel)
  library(biomaRt)
  theme_set(theme_bw())
})

Helper functions

get_counts_stats <- function(sce10x_stem, cluster_to_remove) {
  gene_counts_by_cluster <-
    t(assay(sce10x_stem)) %>%
    as.matrix() %>%
    as_tibble() %>%
    group_by(cluster = colData(sce10x_stem)$cluster_condition) %>%
    summarise_all(list(~ mean(.)))

  umi_counts_max_1 <-
    gene_counts_by_cluster[c(2, 4, 5), -1] %>%
    summarize_all(list(~ max(.))) %>%
    mutate(dummy = "dummy") %>%
    pivot_longer(-dummy, names_to = "gene_name", values_to = "max_yng") %>%
    select(-dummy)

  umi_counts_max_2 <-
    gene_counts_by_cluster[-5, -1] %>%
    summarize_all(list(~ max(.))) %>%
    mutate(dummy = "dummy") %>%
    pivot_longer(-dummy, names_to = "gene_name", values_to = "max_stem12") %>%
    select(-dummy, -gene_name)

  stats <-
    perFeatureQCMetrics(sce10x_stem,
      exprs_values = c("counts"),
      flatten = TRUE
    ) %>%
    as_tibble() %>%
    bind_cols(., rowData(sce10x_stem) %>%
      as_tibble()) %>%
    bind_cols(., umi_counts_max_1, umi_counts_max_2)

  return(stats)
}

get_lfc <- function(contrast, name, dds) {
  lfcShrink(dds,
    contrast = contrast,
    type = "ashr",
    svalue = T
  ) %>%
    as_tibble() %>%
    dplyr::select(-baseMean) %>%
    mutate(svalue = abs(svalue)) %>%
    set_names(paste(c("lfc", "lfc_se", "svalue"), name, sep = "_"))
}

plot_expression <- function(gene, data) {
  plotExpression(data,
    features = gene,
    x = "sample",
    exprs_values = "logcounts",
    colour_by = "condition",
    point_size = 1
  ) +
    facet_wrap(~ colData(data)$clusters, ncol = 1)
}

Define file paths

data_dir <- "./data"
figures_dir <- file.path("./figures")

Analysis and plots

Experiment Data

load muscle satellite cells

sce10x <-
  readRDS(file.path(
    data_dir,
    "preprocessed",
    "sce10x_filtered_final.rds"
  ))

sce10x_stem <- sce10x[, colData(sce10x)$cell_type == "stem"]
colData(sce10x_stem) <- colData(sce10x_stem)[, c("sum", "detected", "sample", "clusters", "cluster_condition", "condition")]
assay(sce10x_stem, 2) <- NULL
assay(sce10x_stem, 2) <- NULL
assay(sce10x_stem, 3) <- NULL
assay(sce10x_stem, 3) <- NULL
reducedDims(sce10x_stem) <- NULL
reducedDims(sce10x_stem) <- NULL
colData(sce10x_stem)$cluster <- "muscs"
colData(sce10x_stem) <- colData(sce10x_stem)[, c(1, 2, 3, 7, 5, 6, 4)]
rm(sce10x)

table(colData(sce10x_stem)$clusters, colData(sce10x_stem)$sample)

       
        yng1 yng2 yng3 aged1 aged2 aged3 aged4
  stem1  563 1043  383   155   300   154   340
  stem2  386  650  540   366   382   172   399
  stem3  250  220  126     0     0     0     0

Discard genes not expressed in at least 10 cells

n_exprs_genes <-
  nexprs(sce10x_stem,
    detection_limit = 0,
    byrow = TRUE
  )
keep <- n_exprs_genes >= 10
table(keep)

keep
FALSE  TRUE 
12928 20546 

sce10x_stem <- sce10x_stem[keep, ]

Compute the mean log normalized counts for each cluster

gene_counts_by_cluster <-
  t(assay(sce10x_stem, "logcounts")) %>%
  as.matrix() %>%
  as_tibble() %>%
  group_by(cluster = colData(sce10x_stem)$cluster_condition) %>%
  summarise_all(list(~ mean(.)))
gene_counts_long <-
  t(gene_counts_by_cluster[, -1]) %>%
  as_tibble(rownames = "gene") %>%
  set_names(c("gene", gene_counts_by_cluster[, 1]$cluster))

The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
Using compatibility `.name_repair`.

gene_counts_long

ggplot(
  gene_counts_long,
  aes(x = stem1_yng, y = stem2_yng)
) +
  geom_point(
    size = .5,
    alpha = 0.5
  ) +
  geom_text(aes(
    label = gene
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  ) +
  geom_smooth(alpha = 0.2, colour = "grey") +
  coord_fixed(ratio = 1)

ggplot(
  gene_counts_long,
  aes(x = stem1_aged, y = stem2_aged)
) +
  geom_point(
    size = .5,
    alpha = 0.5
  ) +
  geom_text(aes(
    label = gene
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  ) +
  geom_smooth(alpha = 0.2, colour = "grey") +
  coord_fixed(ratio = 1)

ggplot(
  gene_counts_long,
  aes(x = stem1_yng, y = stem1_aged)
) +
  geom_point(
    size = .5,
    alpha = 0.5
  ) +
  geom_text(aes(
    label = gene
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  ) +
  geom_smooth(alpha = 0.2, colour = "grey") +
  coord_fixed(ratio = 1)

ggplot(
  gene_counts_long,
  aes(x = stem2_yng, y = stem2_aged)
) +
  geom_point(
    size = .5,
    alpha = 0.5
  ) +
  geom_text(aes(
    label = gene
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  ) +
  geom_smooth(alpha = 0.2, colour = "grey") +
  coord_fixed(ratio = 1)

ggplot(
  gene_counts_long,
  aes(x = stem1_yng, y = stem3_yng)
) +
  geom_point(
    size = .5,
    alpha = 0.5
  ) +
  geom_text(aes(
    label = gene
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  ) +
  geom_smooth(alpha = 0.2, colour = "grey") +
  coord_fixed(ratio = 1)

ggplot(
  gene_counts_long,
  aes(x = stem1_yng, y = stem3_yng)
) +
  geom_point(
    size = .5,
    alpha = 0.5
  ) +
  geom_text(aes(
    label = gene
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  ) +
  geom_smooth(alpha = 0.2, colour = "grey") +
  coord_fixed(ratio = 1)

ggplot(
  gene_counts_long,
  aes(x = stem2_yng, y = stem3_yng)
) +
  geom_point(
    size = .5,
    alpha = 0.5
  ) +
  geom_text(aes(
    label = gene
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  ) +
  geom_smooth(alpha = 0.2, colour = "grey") +
  coord_fixed(ratio = 1)

Tabula Muris Data

load Tabula Muris data

load(file.path(
  data_dir,
  "tabula_muris",
  "facs_Limb_Muscle_seurat_tiss.Robj"
))
facs_limb_muscle_tabula <- tiss
rm(tiss)
facs_limb_muscle_tabula <- UpdateSeuratObject(facs_limb_muscle_tabula)

Updating from v2.X to v3.X
Validating object structure
Updating object slots
Ensuring keys are in the proper strucutre
Ensuring feature names don't have underscores or pipes
Object representation is consistent with the most current Seurat version

facs_limb_muscle_tabula <- as.SingleCellExperiment(facs_limb_muscle_tabula)
facs_limb_muscle_tabula

class: SingleCellExperiment 
dim: 23341 1090 
metadata(0):
assays(3): counts logcounts scaledata
rownames(23341): 0610005C13Rik 0610007C21Rik ... l7Rn6 zsGreen-transgene
rowData names(0):
colnames(1090): A1.B002765.3_38_F.1.1 A1.D042103.3_11_M.1 ... P9.B002769.3_39_F.1.1 P9.D042186.3_8_M.1.1
colData names(18): nReads orig.ident ... nFeature_RNA ident
reducedDimNames(2): PCA TSNE
mainExpName: RNA
altExpNames(0):

load(file.path(
  data_dir,
  "tabula_muris",
  "droplet_Limb_Muscle_seurat_tiss.Robj"
))

droplet_limb_muscle_tabula <- tiss
rm(tiss)
droplet_limb_muscle_tabula <- UpdateSeuratObject(droplet_limb_muscle_tabula)

Updating from v2.X to v3.X
Validating object structure
Updating object slots
Ensuring keys are in the proper strucutre
Ensuring feature names don't have underscores or pipes
Object representation is consistent with the most current Seurat version

droplet_limb_muscle_tabula <- as.SingleCellExperiment(droplet_limb_muscle_tabula)
droplet_limb_muscle_tabula

class: SingleCellExperiment 
dim: 23341 4536 
metadata(0):
assays(3): counts logcounts scaledata
rownames(23341): Xkr4 Rp1 ... Tdtom-transgene zsGreen-transgene
rowData names(0):
colnames(4536): 10X_P7_14_AAACCTGCAGGTCCAC 10X_P7_14_AAACCTGCATCAGTAC ... 10X_P7_15_TTTGTCAGTCTCCACT 10X_P7_15_TTTGTCAGTTGTGGAG
colData names(16): orig.ident channel ... nFeature_RNA ident
reducedDimNames(2): PCA TSNE
mainExpName: RNA
altExpNames(0):

Keep muscle satellite and macrophage cells

# target_cell_types <- c("mesenchymal stem cell", "skeletal muscle satellite cell" ,"endothelial cell", "macrophage")
target_cell_types <- c("skeletal muscle satellite cell", "macrophage")
facs_limb_muscle_tabula <- facs_limb_muscle_tabula[, colData(facs_limb_muscle_tabula)$cell_ontology_class %in% target_cell_types]
droplet_limb_muscle_tabula <- droplet_limb_muscle_tabula[, colData(droplet_limb_muscle_tabula)$cell_ontology_class %in% target_cell_types]

Process and combine data from facs and droplet experiments

droplet_coldata <- colData(droplet_limb_muscle_tabula)[c("nCount_RNA", "nFeature_RNA", "mouse.id", "cell_ontology_class", "cluster.ids")]
droplet_coldata$condition <- "mid-aged"
colnames(droplet_coldata) <- c("sum", "detected", "sample", "cluster", "cluster_condition", "condition")
colData(droplet_limb_muscle_tabula) <- droplet_coldata
colData(droplet_limb_muscle_tabula)

DataFrame with 662 rows and 6 columns
                                 sum  detected   sample                cluster cluster_condition   condition
                           <numeric> <integer> <factor>            <character>       <character> <character>
10X_P7_14_AAACGGGAGTTTCCTT     19328      3305   3-F-56             macrophage                 7    mid-aged
10X_P7_14_AAATGCCCACTCAGGC      5321      1204   3-F-56             macrophage                 7    mid-aged
10X_P7_14_AAATGCCGTCAACATC      2069       892   3-F-56             macrophage                 7    mid-aged
10X_P7_14_AAATGCCTCGATCCCT      4580      1415   3-F-56 skeletal muscle sate..                 5    mid-aged
10X_P7_14_AACACGTAGTAGGTGC      3978       950   3-F-56             macrophage                 7    mid-aged
...                              ...       ...      ...                    ...               ...         ...
10X_P7_15_TTTGCGCTCGGAAACG      3186      1120   3-F-57             macrophage                 7    mid-aged
10X_P7_15_TTTGGTTAGAGCTGGT      3357      1094   3-F-57 skeletal muscle sate..                 5    mid-aged
10X_P7_15_TTTGGTTAGCGATGAC      5408      1809   3-F-57             macrophage                 7    mid-aged
10X_P7_15_TTTGGTTAGGTGATAT      8579      1929   3-F-57             macrophage                 7    mid-aged
10X_P7_15_TTTGTCAGTTGTGGAG      3463      1131   3-F-57 skeletal muscle sate..                 5    mid-aged

facs_coldata <- colData(facs_limb_muscle_tabula)[c("nCount_RNA", "nFeature_RNA", "mouse.id", "cell_ontology_class", "cluster.ids")]
facs_coldata$condition <- "mid-aged"
colnames(facs_coldata) <- c("sum", "detected", "sample", "cluster", "cluster_condition", "condition")
colData(facs_limb_muscle_tabula) <- facs_coldata
colData(facs_limb_muscle_tabula)

DataFrame with 585 rows and 6 columns
                             sum  detected   sample                cluster cluster_condition   condition
                       <numeric> <integer> <factor>            <character>       <character> <character>
A1.B002765.3_38_F.1.1     878318      2149   3_38_F skeletal muscle sate..                 0    mid-aged
A1.D042103.3_11_M.1       774558      1194   3_11_M skeletal muscle sate..                 0    mid-aged
A1.D042186.3_8_M.1.1      914257      2751   3_8_M  skeletal muscle sate..                 2    mid-aged
A10.B002765.3_38_F.1.1   1543007      2333   3_38_F skeletal muscle sate..                 0    mid-aged
A10.B002769.3_39_F.1.1   2521710      3556   3_39_F skeletal muscle sate..                 2    mid-aged
...                          ...       ...      ...                    ...               ...         ...
P7.B002765.3_38_F.1.1    1134974      3213   3_38_F             macrophage                 5    mid-aged
P7.D042473.3_10_M.1.1     963668      3584   3_10_M             macrophage                 5    mid-aged
P8.B002765.3_38_F.1.1     949789      2342   3_38_F             macrophage                 5    mid-aged
P8.B002769.3_39_F.1.1    2217452      3040   3_39_F             macrophage                 5    mid-aged
P9.B002765.3_38_F.1.1    1446574      2749   3_38_F             macrophage                 5    mid-aged

reducedDim(droplet_limb_muscle_tabula) <- NULL
reducedDim(facs_limb_muscle_tabula) <- NULL
reducedDim(droplet_limb_muscle_tabula) <- NULL
reducedDim(facs_limb_muscle_tabula) <- NULL
assay(droplet_limb_muscle_tabula, 3) <- NULL
assay(facs_limb_muscle_tabula, 3) <- NULL

colData(droplet_limb_muscle_tabula)$cluster <-
  case_when(
    colData(droplet_limb_muscle_tabula)$cluster == "macrophage" ~ "tabula_droplet_macrophage",
    colData(droplet_limb_muscle_tabula)$cluster == "skeletal muscle satellite cell" ~ "tabula_droplet_muscs"
  )

colData(facs_limb_muscle_tabula)$cluster <-
  case_when(
    colData(facs_limb_muscle_tabula)$cluster == "macrophage" ~ "tabula_facs_macrophage",
    colData(facs_limb_muscle_tabula)$cluster == "skeletal muscle satellite cell" ~ "tabula_facs_muscs"
  )

colData(droplet_limb_muscle_tabula)$cluster_condition <- 
  paste(colData(droplet_limb_muscle_tabula)$cluster, 
        colData(droplet_limb_muscle_tabula)$cluster_condition, sep = "_")
colData(facs_limb_muscle_tabula)$cluster_condition <- 
  paste(colData(facs_limb_muscle_tabula)$cluster, 
        colData(facs_limb_muscle_tabula)$cluster_condition, sep = "_")
colData(facs_limb_muscle_tabula)$clusters <- "facs_tabula"
colData(droplet_limb_muscle_tabula)$clusters <- "droplet_tabula"

droplet_limb_muscle_tabula <- droplet_limb_muscle_tabula[rownames(facs_limb_muscle_tabula), ]
tabula_muris <- cbind(droplet_limb_muscle_tabula, facs_limb_muscle_tabula)
# tabula_muris <- droplet_limb_muscle_tabula
tabula_muris

class: SingleCellExperiment 
dim: 23341 1247 
metadata(0):
assays(2): counts logcounts
rownames(23341): 0610005C13Rik 0610007C21Rik ... l7Rn6 zsGreen-transgene
rowData names(0):
colnames(1247): 10X_P7_14_AAACGGGAGTTTCCTT 10X_P7_14_AAATGCCCACTCAGGC ... P8.B002769.3_39_F.1.1 P9.B002765.3_38_F.1.1
colData names(7): sum detected ... condition clusters
reducedDimNames(0):
mainExpName: RNA
altExpNames(0):

Discard genes not expressed in at least 10 cells

n_exprs_genes <-
  nexprs(tabula_muris,
    detection_limit = 0,
    byrow = TRUE
  )
keep <- n_exprs_genes >= 10
table(keep)

keep
FALSE  TRUE 
10714 12627 

tabula_muris <- tabula_muris[keep, ]

Map gene names to ENSEMBL gene IDS

ensembl93 <- useEnsembl(
  biomart = "genes",
  dataset = "mmusculus_gene_ensembl",
  version = 93
)
mart <-
  useDataset(
    "mmusculus_gene_ensembl",
    ensembl93
  ) # from 2018

gene_metadata <-
  getBM(
    filters = "external_gene_name",
    attributes = c(
      "ensembl_gene_id_version",
      "external_gene_name",
      "gene_biotype",
      "chromosome_name",
      "description"
    ),
    values = rownames(tabula_muris),
    mart = mart
  )

chr_to_keep <- c(seq(1, 19), "MT", "X", "Y")
gene_metadata <-
  gene_metadata %>%
  filter(chromosome_name %in% chr_to_keep) %>%
  filter(!duplicated(external_gene_name))

gene_metadata <-
  gene_metadata %>%
  left_join(tibble(external_gene_name = rownames(tabula_muris)),
    .,
    by = "external_gene_name"
  ) %>%
  dplyr::rename(chr = "chromosome_name")
rowData(tabula_muris) <- gene_metadata

Discard genes not in standard chromosomes and those with duplicate ids

tabula_muris_filtered <- tabula_muris[!is.na(rowData(tabula_muris)$ensembl_gene_id_version), ]
tabula_muris_filtered

class: SingleCellExperiment 
dim: 11528 1247 
metadata(0):
assays(2): counts logcounts
rownames(11528): 0610009B22Rik 0610009L18Rik ... Zzef1 Zzz3
rowData names(5): external_gene_name ensembl_gene_id_version gene_biotype chr description
colnames(1247): 10X_P7_14_AAACGGGAGTTTCCTT 10X_P7_14_AAATGCCCACTCAGGC ... P8.B002769.3_39_F.1.1 P9.B002765.3_38_F.1.1
colData names(7): sum detected ... condition clusters
reducedDimNames(0):
mainExpName: RNA
altExpNames(0):

table(colData(tabula_muris_filtered)$cluster)

tabula_droplet_macrophage      tabula_droplet_muscs    tabula_facs_macrophage         tabula_facs_muscs 
                      308                       354                        45                       540 

# gene_counts_by_cluster_tabula_muris <-
#   t(assay(tabula_muris_filtered, "logcounts")) %>%
#   as.matrix() %>%
#   as_tibble() %>%
#   group_by(cluster = colData(tabula_muris_filtered)$cluster) %>%
#   summarise_all(list(~ mean(.)))
# gene_counts_long_tabula_muris <-
#   t(gene_counts_by_cluster_tabula_muris[, -1]) %>%
#   as_tibble(rownames = "gene") %>%
#   set_names(c("gene", gene_counts_by_cluster_tabula_muris[, 1]$cluster))
# gene_counts_long_tabula_muris

Combined Data

Join data using common genes only

commmon_genes <- 
  intersect(rowData(sce10x_stem)$ensembl_gene_id_version, 
            rowData(tabula_muris_filtered)$ensembl_gene_id_version)
sce10x_stem_filtered <- sce10x_stem
rowData(sce10x_stem_filtered) <- rowData(sce10x_stem_filtered)[, c(2, 1, 3, 4, 5)]
rownames(sce10x_stem_filtered) <- rowData(sce10x_stem)$ensembl_gene_id_version
rownames(tabula_muris_filtered) <- rowData(tabula_muris_filtered)$ensembl_gene_id_version
sce10x_stem_filtered <- sce10x_stem_filtered[commmon_genes, ]
tabula_muris_filtered_2 <- tabula_muris_filtered[commmon_genes, ]
rowData(sce10x_stem_filtered) <- rowData(tabula_muris_filtered_2)

rownames(sce10x_stem_filtered) <- rowData(sce10x_stem_filtered)$external_gene_name
rownames(tabula_muris_filtered_2) <- rowData(tabula_muris_filtered_2)$external_gene_name
combined <- cbind(sce10x_stem_filtered, tabula_muris_filtered_2)
combined

class: SingleCellExperiment 
dim: 9978 7676 
metadata(1134): tximetaInfo quantInfo ... txomeInfo txdbInfo
assays(2): counts logcounts
rownames(9978): Lypla1 Tcea1 ... Uty Ddx3y
rowData names(5): external_gene_name ensembl_gene_id_version gene_biotype chr description
colnames(7676): CTACAGAGTGGCCCAT_d1 CTAACCCCACACCTGG_d1 ... P8.B002769.3_39_F.1.1 P9.B002765.3_38_F.1.1
colData names(7): sum detected ... condition clusters
reducedDimNames(0):
mainExpName: RNA
altExpNames(0):

Compute the mean log normalized counts for each cluster

gene_counts_by_cluster_combined <-
  t(assay(combined, "logcounts")) %>%
  as.matrix() %>%
  as_tibble() %>%
  group_by(cluster = colData(combined)$cluster) %>%
  summarise_all(list(~ mean(.)))
gene_counts_long_combined <-
  t(gene_counts_by_cluster_combined[, -1]) %>%
  as_tibble(rownames = "gene") %>%
  set_names(c("gene", gene_counts_by_cluster_combined[, 1]$cluster))
gene_counts_long_combined

Compute correlations

cor_dt <- cor(gene_counts_long_combined[ , -1])
cor_dt

                          muscs tabula_droplet_macrophage tabula_droplet_muscs tabula_facs_macrophage tabula_facs_muscs
muscs                      1.00                      0.70                 0.85                   0.52              0.78
tabula_droplet_macrophage  0.70                      1.00                 0.84                   0.83              0.69
tabula_droplet_muscs       0.85                      0.84                 1.00                   0.62              0.83
tabula_facs_macrophage     0.52                      0.83                 0.62                   1.00              0.71
tabula_facs_muscs          0.78                      0.69                 0.83                   0.71              1.00

corrplot::corrplot(cor_dt, 
                   order = 'AOE', 
                   type = 'upper', 
                   tl.pos = 'd',
                   tl.cex=.7)
corrplot::corrplot(cor_dt, 
                   add = TRUE, 
                   type = 'lower',
                   method = 'number',
                   order = 'AOE',
                   diag = FALSE, 
                   tl.pos = 'n', 
                   cl.pos = 'n')

comparison across the three datasets

MuSC Tabula Muris droplet vs MuSC experiment droplet

ggplot(
  gene_counts_long_combined,
  aes( x= tabula_droplet_muscs, y = muscs)
) +
  geom_point(
    size = .5,
    alpha = 0.5
  ) +
  geom_text(aes(
    label = gene
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  ) +
  geom_smooth(alpha = 0.2, colour = "grey") +
  coord_fixed(ratio = 1)

MuSCs Tabula Muris droplet vs MuSCs Tabula Muris FACS

ggplot(
  gene_counts_long_combined,
  aes(x = tabula_facs_muscs , y=tabula_droplet_muscs)
) +
  geom_point(
    size = .5,
    alpha = 0.5
  ) +
  geom_text(aes(
    label = gene
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  ) +
  geom_smooth(alpha = 0.2, colour = "grey") +
  coord_fixed(ratio = 1)

Macrophage Tabula Muris droplet vs MuSCs experiment droplet

ggplot(
  gene_counts_long_combined,
  aes(x = tabula_droplet_macrophage , y=muscs)
) +
  geom_point(
    size = .5,
    alpha = 0.5
  ) +
  geom_text(aes(
    label = gene
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  ) +
  geom_smooth(alpha = 0.2, colour = "grey") +
  coord_fixed(ratio = 1)

Macrophage Tabula Muris droplet vs MuSC experiment droplet

ggplot(
  gene_counts_long_combined,
  aes(x = tabula_droplet_macrophage, y = tabula_droplet_muscs)
) +
  geom_point(
    size = .5,
    alpha = 0.5
  ) +
  geom_text(aes(
    label = gene
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  ) +
  geom_smooth(alpha = 0.2, colour = "grey") +
  coord_fixed(ratio = 1)

Macrophage Tabula Muris facs vs MuSC experiment droplet

ggplot(
  gene_counts_long_combined,
  aes(x = tabula_facs_macrophage, y = muscs)
) +
  geom_point(
    size = .5,
    alpha = 0.5
  ) +
  geom_text(aes(
    label = gene
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  ) +
  geom_smooth(alpha = 0.2, colour = "grey") +
  coord_fixed(ratio = 1)

Plot gene expression

dec.combined <- modelGeneVar(combined[, !(colData(combined)$clusters =="facs_tabula")])

# Visualizing the fit:
fit.combined <- metadata(dec.combined)
plot(fit.combined$mean, fit.combined$var,
  xlab = "Mean of log-expression",
  ylab = "Variance of log-expression"
)
curve(fit.combined$trend(x), col = "dodgerblue", add = TRUE, lwd = 2)

hvg.combined.var <- getTopHVGs(dec.combined, n = 300)

plotHeatmap(combined[, !(colData(combined)$clusters =="facs_tabula")],
  color = colorRampPalette(RColorBrewer::brewer.pal(n = 11, name = "RdBu"))(11),
  features = hvg.combined.var,
  #  cluster_cols=FALSE,
  cluster_rows = TRUE,
  center = TRUE,
  # zlim= c(-3,3),
  order_columns_by = c("cluster", "sum"),
  #  gaps_row = seq(0,5*5,by=5)
)

# plotExpression(combined,
#   features = "Map4k4",
#   x = "cluster",
#   exprs_values = "logcounts",
#   colour_by = "clusters",
#   point_size = 1
# ) +
#   facet_wrap(~ colData(combined)$clusters, ncol = 1)

sessionInfo()

R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] biomaRt_2.48.3              BiocParallel_1.26.2         forcats_0.5.1               stringr_1.4.0               dplyr_1.0.7                 purrr_0.3.4                
 [7] readr_2.0.2                 tidyr_1.1.4                 tibble_3.1.6                tidyverse_1.3.1             SeuratObject_4.0.3          Seurat_4.0.5               
[13] scran_1.20.1                scater_1.20.1               ggplot2_3.3.5               scuttle_1.2.1               SingleCellExperiment_1.14.1 SummarizedExperiment_1.22.0
[19] Biobase_2.52.0              GenomicRanges_1.44.0        GenomeInfoDb_1.28.4         IRanges_2.26.0              S4Vectors_0.30.2            BiocGenerics_0.38.0        
[25] MatrixGenerics_1.4.3        matrixStats_0.61.0         

loaded via a namespace (and not attached):
  [1] utf8_1.2.2                reticulate_1.22           tidyselect_1.1.1          RSQLite_2.2.8             AnnotationDbi_1.54.1      htmlwidgets_1.5.4        
  [7] grid_4.1.2                Rtsne_0.15                munsell_0.5.0             ScaledMatrix_1.0.0        codetools_0.2-18          ica_1.0-2                
 [13] statmod_1.4.36            future_1.23.0             miniUI_0.1.1.1            withr_2.4.2               colorspace_2.0-2          filelock_1.0.2           
 [19] knitr_1.36                rstudioapi_0.13           ROCR_1.0-11               tensor_1.5                listenv_0.8.0             labeling_0.4.2           
 [25] GenomeInfoDbData_1.2.6    polyclip_1.10-0           pheatmap_1.0.12           farver_2.1.0              bit64_4.0.5               rprojroot_2.0.2          
 [31] parallelly_1.28.1         vctrs_0.3.8               generics_0.1.1            xfun_0.28                 BiocFileCache_2.0.0       R6_2.5.1                 
 [37] ggbeeswarm_0.6.0          rsvd_1.0.5                locfit_1.5-9.4            cachem_1.0.6              bitops_1.0-7              spatstat.utils_2.2-0     
 [43] DelayedArray_0.18.0       assertthat_0.2.1          promises_1.2.0.1          scales_1.1.1              beeswarm_0.4.0            gtable_0.3.0             
 [49] beachmat_2.8.1            globals_0.14.0            goftest_1.2-3             rlang_0.4.12              splines_4.1.2             lazyeval_0.2.2           
 [55] spatstat.geom_2.3-0       broom_0.7.10              yaml_2.2.1                reshape2_1.4.4            abind_1.4-5               modelr_0.1.8             
 [61] backports_1.3.0           httpuv_1.6.3              tools_4.1.2               ellipsis_0.3.2            spatstat.core_2.3-1       jquerylib_0.1.4          
 [67] RColorBrewer_1.1-2        ggridges_0.5.3            Rcpp_1.0.7                plyr_1.8.6                progress_1.2.2            sparseMatrixStats_1.4.2  
 [73] zlibbioc_1.38.0           RCurl_1.98-1.5            prettyunits_1.1.1         rpart_4.1-15              deldir_1.0-6              pbapply_1.5-0            
 [79] viridis_0.6.2             cowplot_1.1.1             zoo_1.8-9                 haven_2.4.3               ggrepel_0.9.1             cluster_2.1.2            
 [85] fs_1.5.0                  magrittr_2.0.1            data.table_1.14.2         scattermore_0.7           lmtest_0.9-39             reprex_2.0.1             
 [91] RANN_2.6.1                fitdistrplus_1.1-6        hms_1.1.1                 patchwork_1.1.1           mime_0.12                 evaluate_0.14            
 [97] xtable_1.8-4              XML_3.99-0.8              readxl_1.3.1              gridExtra_2.3             compiler_4.1.2            KernSmooth_2.23-20       
[103] crayon_1.4.2              htmltools_0.5.2           mgcv_1.8-38               later_1.3.0               tzdb_0.2.0                lubridate_1.8.0          
[109] DBI_1.1.1                 corrplot_0.92             dbplyr_2.1.1              rappdirs_0.3.3            MASS_7.3-54               Matrix_1.3-4             
[115] cli_3.1.0                 metapod_1.0.0             igraph_1.2.8              pkgconfig_2.0.3           plotly_4.10.0             spatstat.sparse_2.0-0    
[121] xml2_1.3.2                vipor_0.4.5               bslib_0.3.1               dqrng_0.3.0               XVector_0.32.0            rvest_1.0.2              
[127] digest_0.6.28             sctransform_0.3.2         RcppAnnoy_0.0.19          Biostrings_2.60.2         spatstat.data_2.1-0       rmarkdown_2.11           
[133] cellranger_1.1.0          leiden_0.3.9              uwot_0.1.10               edgeR_3.34.1              DelayedMatrixStats_1.14.3 curl_4.3.2               
[139] shiny_1.7.1               lifecycle_1.0.1           nlme_3.1-152              jsonlite_1.7.2            BiocNeighbors_1.10.0      viridisLite_0.4.0        
[145] limma_3.48.3              fansi_0.5.0               pillar_1.6.4              lattice_0.20-45           KEGGREST_1.32.0           fastmap_1.1.0            
[151] httr_1.4.2                survival_3.2-13           glue_1.5.0                png_0.1-7                 bit_4.0.4                 bluster_1.2.1            
[157] stringi_1.7.5             sass_0.4.0                blob_1.2.2                BiocSingular_1.8.1        memoise_2.0.0             irlba_2.3.3              
[163] future.apply_1.8.1