Load & reformat data

samples.df <- data.frame(
  sample.name = list.files(dir.inputs),
  stringsAsFactors = FALSE
) %>% 
  mutate(
    sample.name = str_remove(sample.name, "_(features|barcodes|matrix).+gz$")
  ) %>% 
  distinct() %>% 
  mutate(
    accession = word(sample.name, 1, sep = "_"),
    sex = word(sample.name, 2, sep = "_"),
    treatment = word(sample.name, 3, sep = "_"),
    replicate = word(sample.name, 4, sep = "_")
  ) %>% 
  mutate(
    sex = factor(sex, levels = c("Male", "Female")),
    treatment = factor(treatment, levels = c("saline", "LPS"))
  )
samples.df

# load raw counts
matrix_dirs <- list.dirs(dir.inputs, recursive = FALSE, full.names = TRUE)
names(matrix_dirs) <- basename(matrix_dirs)
sce <- read10xCounts(matrix_dirs)

# rename row/colData colnames & SCE dimnames
names(rowData(sce)) <- c("ENSEMBL", "SYMBOL", "Type")
names(colData(sce)) <- c("sample_id", "barcode")
sce$sample_id <- factor(basename(sce$sample_id), levels = samples.df$sample.name)
dimnames(sce) <- list(
    with(rowData(sce), paste(ENSEMBL, SYMBOL, sep = ".")),
    with(colData(sce), paste(barcode, sample_id, sep = ".")))

# load metadata
m <- match(sce$sample_id, samples.df$sample.name)
sce$sex <- samples.df$sex[m]
sce$treatment <- samples.df$treatment[m]
sce$replicate <- samples.df$replicate[m]

# remove undetected genes
sce <- sce[rowSums(counts(sce) > 0) > 0, ]
dim(sce)

[1]  23337 116707

Doublet removal

# split SCE by sample
cs_by_s <- split(colnames(sce), sce$sample_id)
sce_by_s <- lapply(cs_by_s, function(cs) sce[, cs])

# run 'scds'
sce_by_s <- lapply(sce_by_s, function(u) 
    cxds_bcds_hybrid(bcds(cxds(u))))

# remove doublets
sce_by_s <- lapply(sce_by_s, function(u) {
    # compute expected nb. of doublets (10x)
    n_dbl <- ceiling(0.01 * ncol(u)^2 / 1e3)
    # remove 'n_dbl' cells w/ highest doublet score
    o <- order(u$hybrid_score, decreasing = TRUE)
    u[, -o[seq_len(n_dbl)]]
})

# merge back into single SCE
sce <- do.call("cbind", sce_by_s)

Calculate QC Metrics

(mito <- grep("mt-", rownames(sce), value = TRUE))

 [1] "ENSMUSG00000064341.mt-Nd1"  "ENSMUSG00000064345.mt-Nd2"  "ENSMUSG00000064351.mt-Co1"  "ENSMUSG00000064354.mt-Co2"  "ENSMUSG00000064356.mt-Atp8"
 [6] "ENSMUSG00000064357.mt-Atp6" "ENSMUSG00000064358.mt-Co3"  "ENSMUSG00000064360.mt-Nd3"  "ENSMUSG00000065947.mt-Nd4l" "ENSMUSG00000064363.mt-Nd4" 
[11] "ENSMUSG00000064367.mt-Nd5"  "ENSMUSG00000064368.mt-Nd6"  "ENSMUSG00000064370.mt-Cytb"

sce <- addPerCellQC(sce, subsets = list(Mt = mito))
plotHighestExprs(sce, n = 20)

Filtering

# get sample-specific outliers
cols <- c("sum", "detected", "subsets_Mt_percent")
log <- c(TRUE, TRUE, FALSE)
type <- c("both", "both", "higher")

drop_cols <- paste0(cols, "_drop")
for (i in seq_along(cols))
  colData(sce)[[drop_cols[i]]] <- isOutlier(
    sce[[cols[i]]],
    nmads = 2.5,
    type = type[i],
    log = log[i],
    batch = sce$sample_id
  )

sapply(drop_cols, function(i)
  sapply(drop_cols, function(j)
    sum(sce[[i]] & sce[[j]])))

                        sum_drop detected_drop subsets_Mt_percent_drop
sum_drop                     288           102                      64
detected_drop                102          1115                     976
subsets_Mt_percent_drop       64           976                   11665

cd <- data.frame(colData(sce))
ps <- lapply(seq_along(cols), function (i) {
  p <- ggplot(cd, aes_string(x = cols[i], alpha = drop_cols[i])) +
    geom_histogram(bins = 100, show.legend = FALSE) +
    scale_alpha_manual(values = c("FALSE" = 1, "TRUE" = 0.4)) +
    facet_wrap(~ sample_id, ncol = 1, scales = "free") +
    theme_classic() + theme(strip.background = element_blank())
  if (log[i])
    p <- p + scale_x_log10()
  return(p)
})

`aes_string()` was deprecated in ggplot2 3.0.0.

plot_grid(plotlist = ps, ncol = 3)

layout(matrix(1:2, nrow = 1))
ol <- rowSums(as.matrix(colData(sce)[drop_cols])) != 0
x <- sce$sum
y <- sce$detected
heatscatter(x, y, log="xy", main = "unfiltered", 
    xlab = "Total counts", ylab = "Non-zero features")
heatscatter(x[!ol], y[!ol], log="xy", main = "filtered", 
    xlab = "Total counts", ylab = "Non-zero features")

# summary of cells kept
ns <- table(sce$sample_id)
ns_fil <- table(sce$sample_id[!ol])
print(rbind(
    unfiltered = ns, filtered = ns_fil, 
    "%" = ns_fil / ns * 100), digits = 0)

           GSM4475127_Female_saline_1 GSM4475128_Male_saline_1 GSM4475129_Female_saline_2 GSM4475130_Male_saline_2 GSM4475131_Female_saline_3
unfiltered                       4604                     5724                       5696                     3271                       7242
filtered                         4604                     5724                       5696                     3271                       7242
%                                 100                      100                        100                      100                        100
           GSM4475132_Male_saline_3 GSM4475133_Female_LPS_1 GSM4475134_Male_LPS_1 GSM4475135_Female_LPS_2 GSM4475136_Male_LPS_2 GSM4475137_Female_LPS_3
unfiltered                     5178%#24.0-1e                  7983                    6208%#22.0-1e                    9305
filtered                       5178%#24.0-1e                  7983                    6208%#22.0-1e                    9305
%                               100%#24.0-1e                   100                     100%#22.0-1e                     100
           GSM4475138_Male_LPS_3
unfiltered%#22.0-1e
filtered  %#22.0-1e
%         %#22.0-1e

# drop outlier cells
sce <- sce[, !ol]
dim(sce)

[1] 12499 91104

# require count > 1 in at least 20 cells
sce <- sce[rowSums(counts(sce) > 1) >= 20, ]
dim(sce)

[1] 12499 91104

Save SCE

saveRDS(sce, file.path(dir.objects, "SCE_preprocessing.rds"))

—Checkpoint—

Load data

sce <- readRDS(file.path(dir.objects, "SCE_preprocessing.rds"))

Error in file.path(dir.objects, "SCE_preprocessing.rds") : 
  object 'dir.objects' not found

UMI filter

plotColData(sce, x = "sample_id", y = "detected", colour_by = "sum") +
  geom_hline(yintercept = 700, color = "red") +
  coord_flip()

Integration

# create SeuratObject
so <- CreateSeuratObject(
    counts = counts(sce),
    meta.data = data.frame(colData(sce)),
    project = "LPS")
# subset by UMI count
so.UMI700 <- subset(so, subset = nCount_RNA > 700)

# split by sample
cells_by_sample <- split(colnames(sce), sce$sample_id)
so.UMI700 <- lapply(cells_by_sample, function(i) subset(so.UMI700, cells = i))

# normalize, find variable genes, and scale
so.UMI700 <- lapply(so.UMI700, NormalizeData, verbose = FALSE)
so.UMI700 <- lapply(so.UMI700, FindVariableFeatures, nfeatures = 2e3, 
    selection.method = "vst", verbose = FALSE)
so.UMI700 <- lapply(so.UMI700, ScaleData, verbose = FALSE)

# find anchors & integrate
as <- FindIntegrationAnchors(so.UMI700, verbose = FALSE)

'as(<dgCMatrix>, "dgTMatrix")' is deprecated.
Use 'as(., "TsparseMatrix")' instead.
See help("Deprecated") and help("Matrix-deprecated").
UNRELIABLE VALUE: One of the ‘future.apply’ iterations (‘future_lapply-1’) unexpectedly generated random numbers without declaring so. There is a risk that those random numbers are not statistically sound and the overall results might be invalid. To fix this, specify 'future.seed=TRUE'. This ensures that proper, parallel-safe random numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use 'future.seed = NULL', or set option 'future.rng.onMisuse' to "ignore".

so.UMI700 <- IntegrateData(anchorset = as, dims = seq_len(30), verbose = FALSE)

# scale integrated data
DefaultAssay(so.UMI700) <- "integrated"
so.UMI700 <- ScaleData(so.UMI700, verbose = FALSE)

Dimension reduction

so.UMI700 <- RunPCA(so.UMI700, npcs = 30, verbose = FALSE)

so.UMI700 <- RunTSNE(so.UMI700, reduction = "pca", dims = seq_len(20),
    seed.use = 1, do.fast = TRUE, verbose = FALSE)

so.UMI700 <- RunUMAP(so.UMI700, reduction = "pca", dims = seq_len(20),
    seed.use = 1, verbose = FALSE)

Clustering

so.UMI700 <- FindNeighbors(so.UMI700, reduction = "pca", dims = seq_len(20), verbose = FALSE)
for (res in c(0.1, 0.4, 0.6, 1))
  so.UMI700 <- FindClusters(so.UMI700, resolution = res, random.seed = 1, verbose = FALSE)

snn.colnames <-
  colnames(so.UMI700@meta.data) %>% str_subset("integrated_snn_res")

thm <- theme(aspect.ratio = 1, legend.position = "top")

ps <- lapply(snn.colnames, function(u) {
  p1 <- DimPlot(so.UMI700, reduction = "tsne", group.by = u) + thm +
    labs(title = u)
  plot_grid(p1)
})
plot_grid(plotlist = ps)

so.UMI700$group_id <- paste(so.UMI700$sex, so.UMI700$treatment, sep = "_")
so.UMI700 <- SetIdent(so.UMI700, value = "integrated_snn_res.0.6")

DR plots

thm <- theme(aspect.ratio = 1, legend.position = "none")
ps <- lapply(c("sample_id", "group_id", "ident"), function(u) {
    p1 <- DimPlot(so.UMI700, reduction = "tsne", group.by = u) + thm
    p2 <- DimPlot(so.UMI700, reduction = "umap", group.by = u)
    lgd <- get_legend(p2)
    p2 <- p2 + thm
    list(p1, p2, lgd)
    plot_grid(p1, p2, lgd, nrow = 1,
        rel_widths = c(1, 1, 0.5))
})
plot_grid(plotlist = ps, ncol = 1)

Save SeuratObject

saveRDS(so.UMI700, file.path(dir.objects, "SO_clustering.rds"))

—Checkpoint—

Load data & convert to SCE

so <- readRDS(file.path(dir.objects, "SO_clustering.rds"))
sce <- as.SingleCellExperiment(so, assay = "RNA")
colData(sce) <- as.data.frame(colData(sce)) %>%
  mutate_if(is.character, as.factor) %>%
  DataFrame(row.names = colnames(sce))

Nb. of clusters by resolution

cluster_cols <- grep("res.[0-9]", colnames(colData(sce)), value = TRUE)
sapply(colData(sce)[cluster_cols], nlevels)

integrated_snn_res.0.1 integrated_snn_res.0.4 integrated_snn_res.0.6   integrated_snn_res.1 
                    10                     18                     23                     31 

Cluster-sample counts

# Re-leveling
sce$group_id <- factor(sce$group_id, 
  levels = c("Female_saline", "Male_saline", "Female_LPS", "Male_LPS"))

# set cluster IDs to selected resolution clustering
current_snn_res <- "integrated_snn_res.1"
so <- SetIdent(so, value = current_snn_res)

reordered_clusters_levels <- levels(Idents(so)) %>% 
  length() %>% seq(0, .)
Idents(so) <- factor(Idents(so), levels = reordered_clusters_levels)

so@meta.data$cluster_id <- Idents(so)
sce$cluster_id <- Idents(so)
levels(sce$cluster_id)

 [1] "0"  "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14" "15" "16" "17" "18" "19" "20" "21" "22" "23" "24" "25" "26" "27" "28" "29" "30"

Relative cluster-abundances

n_cells <- table(sce$cluster_id, sce$sample_id)
fqs <- prop.table(n_cells, margin = 2)
mat <- as.matrix(unclass(fqs))

mat <- as.matrix(unclass(n_cells))
Heatmap(mat,
    col = rev(brewer.pal(11, "RdGy")[-6]),
    column_title = current_snn_res, 
    name = "N cells",
    cluster_rows = FALSE,
    cluster_columns = FALSE,
    row_names_side = "left",
    row_title = "cluster_id",
    rect_gp = gpar(col = "white"),
    cell_fun = function(i, j, x, y, width, height, fill)
        grid.text(mat[j, i], x = x, y = y, 
            gp = gpar(col = "white", fontsize = 8)))

Known marker genes

fs <- list(
  Astrocytes = c(
    "Aqp4",
    "Aldoc",
    "Gfap",
    "Clu",
    "Apoe",
    "Hes5",
    "Dio2",
    "Agt",
    "Sparcl1",
    "Glul",
    "Itih3"
  ),
  
  Neurons = c(
    "Rtn1",
    "Syt1",
    "Calm2",
    "Camk2b",
    "Chgb",
    "Camk2a",
    "Meg3",
    "Snhg11",
    "Slc17a7",
    "Grin2b",
    "Grin2a"
  ),
  
  Endothelial = c(
    "Cldn5",
    "Itm2a",
    "Ly6c1",
    "Bsg",
    "Flt1",
    "Pltp",
    "Klf2",
    "Slco1a4"
  ),
  
  Microglia = c(
    "C1qa",
    "C1qb",
    "Ctss",
    "Hexb",
    "P2ry12",
    "Tyrobp",
    "Trem2",
    "C1qc"
  ),
  
  Oligodendrocytes = c(
    "Mbp",
    "Mobp",
    "Plp1",
    "Kif5a",
    "Plekhb1",
    "Stmn1",
    "Qdpr",
    "Pdgfra",
    "Pllp",
    "Olig1",
    "Gjc3",
    "Neu4"
  ),
  Pericytes = c("Vtn",
                "Igfbp7",
                "Cald1",
                "Ndufa4l2",
                "Higd1b",
                "Myl9",
                "Gng11"),
  Erythrocytes = c("Alas2",
                  "Hba-a2",
                  "Hbb-bt",
                  "Hba-a1",
                  "Hbb-bs")
)

fs <- lapply(fs, sapply, function(g)
  grep(paste0(g, "$"), rownames(sce), value = TRUE))

gs <- gsub(".*\\.", "", unlist(fs))
ns <- vapply(fs, length, numeric(1))
ks <- rep.int(names(fs), ns)
labs <- sprintf("%s(%s)", gs, ks)

# split cells by cluster
cs_by_k <- split(colnames(sce), sce$cluster_id)

# compute cluster-marker means
ms_by_cluster <- lapply(fs, function(gs)
  vapply(cs_by_k, function(i)
    Matrix::rowMeans(logcounts(sce)[gs, i, drop = FALSE]),
    numeric(length(gs))))

# prep. for plotting & scale b/w 0 and 1
mat <- do.call("rbind", ms_by_cluster)
mat <- muscat:::.scale(mat)
rownames(mat) <- word(rownames(mat), 2, sep = "\\.")

cols <- muscat:::.cluster_colors[seq_along(fs)]
cols <- setNames(cols, names(fs))

row_anno <- rowAnnotation(df = data.frame(label = factor(ks, levels = names(fs))),
                          col = list(label = cols), show_annotation_name = FALSE)

Heatmap(
  mat,
  column_title = current_snn_res,
  name = "scaled avg.\nexpression",
  col = viridis(10),
  cluster_rows = FALSE,
  cluster_columns = FALSE,
  row_names_side = "left",
  column_names_rot = 0,
  column_names_centered = T,
  row_names_gp =  gpar(fontfamily = "mono", fontface = "italic"),
  left_annotation = row_anno
)

# at snn res 1
anno.res.1 <- list(
    "Astrocytes" = c(0:14, 16, 22, 23),
    "Neurons" = c(15, 18, 21, 24),
    "Endothelial" = c(17, 28, 29),
    "Microglia" = 19,
    "Oligodendrocytes" = c(20, 25, 27),
    "Pericytes" = c(26),
    "Erythrocytes" = 30
    )

anno <- anno.res.1
  
cols <- muscat:::.cluster_colors[seq_along(anno)]
cols <- setNames(cols, names(anno))

m <- match(sce$cluster_id, unlist(anno))
ns <- vapply(anno, length, numeric(1))
lab <- rep.int(names(anno), ns)[m]
sce$cell_type <- factor(lab, levels = names(anno))

# Is any cell's cell type NA? 
table(is.na(sce$cell_type))

FALSE 
86429 

# What clusters are unassigned?
assigned <- unlist(anno) %>% unname()
clusters <- 0:max(as.numeric(levels(sce$cluster_id)))
clusters[!(clusters %in% assigned)]

integer(0)

Find markers unsupervized using scran

scran_markers <- findMarkers(sce, 
    groups = sce$cluster_id, block = sce$sample_id,
    direction = "up", lfc = 2, full.stats = TRUE)

gs <- lapply(scran_markers, function(u)
  rownames(u)[u$Top == 1])
sub <- sce[unique(unlist(gs)),]
pbs <-
  aggregateData(sub,
                assay = "logcounts",
                by = "cluster_id",
                fun = "mean")
mat <- muscat:::.scale(assay(pbs))
rownames(mat) <- word(rownames(mat), 2, sep = "\\.")

col.clusters.df <-
  enframe(anno, name = "cell_type", value = "cluster") %>%
  unnest_longer(cluster) %>%
  arrange(cluster) %>%
  left_join(.,
            enframe(cols, name = "cell_type", value = "color"),
            by = "cell_type")

col.anno.clusters <- columnAnnotation("cell type" = col.clusters.df$cell_type,
                                      col = list("cell type" = cols))

Heatmap(
  mat,
  name = "scaled avg.\nexpression",
  col = viridis(10),
  cluster_rows = TRUE,
  show_row_dend = FALSE,
  cluster_columns = FALSE,
  column_names_side = "top",
  column_names_rot = 0,
  column_names_centered = TRUE,
  row_title = "cluster_id",
  column_title = current_snn_res,
  rect_gp = gpar(col = "white"),
  top_annotation = col.anno.clusters
)

Cluster-sample cell type counts

cell_types <- table(sce$cell_type, sce$sample_id)
mat <- as.matrix(unclass(cell_types))
Heatmap(mat,
    col = c("#7A92FF", "white", "#FF8492"),
    name = "N cells",
    cluster_rows = FALSE,
    cluster_columns = FALSE,
    row_names_side = "left",
    column_title = current_snn_res,
    rect_gp = gpar(col = "white"),
    cell_fun = function(i, j, x, y, width, height, fill)
        grid.text(mat[j, i], x = x, y = y, 
            gp = gpar(col = "black", fontsize = 8)))

Cell type counts by group

cell_types_per_group <- table(sce$cell_type, sce$group_id)
mat <- as.matrix(unclass(cell_types_per_group))

Heatmap(
  mat,
  col = c("#7A92FF", "white", "#FF8492"),
  name = "N cells",
  cluster_rows = FALSE,
  cluster_columns = FALSE,
  row_names_side = "left",
  column_title = current_snn_res,
  rect_gp = gpar(col = "white"),
  cell_fun = function(i, j, x, y, width, height, fill)
    grid.text(
      mat[j, i],
      x = x,
      y = y,
      gp = gpar(col = "black", fontsize = 12)
    )
)

Cell type counts by cluster

cell_types_per_group <- table(sce$cluster_id, sce$cell_type)
mat <- as.matrix(unclass(cell_types_per_group))

Heatmap(
  mat,
  col = c("#7A92FF", "white", "#FF8492"),
  name = "N cells",
  cluster_rows = FALSE,
  cluster_columns = FALSE,
  row_names_side = "left",
  column_title = current_snn_res,
  rect_gp = gpar(col = "white"),
  cell_fun = function(i, j, x, y, width, height, fill)
    grid.text(
      mat[j, i],
      x = x,
      y = y,
      gp = gpar(col = "black", fontsize = 12)
    )
)

Cell type ratios

cols <- muscat:::.cluster_colors[seq_along(anno)]
cols <- setNames(cols, names(anno))

table(sce$cell_type) %>% as.data.frame() %>% 
  mutate(
    text = Freq / sum(Freq) * 100,
    text = paste(Freq, "cells,", round(text, digits = 2), "%")
  ) %>% 
  ggplot(aes(y = Var1, x = Freq)) +
  geom_col(aes(fill = Var1)) +
  geom_text(aes(label = text, x = max(Freq) / 5)) +
  theme(legend.position = "none") +
  labs(y = NULL, x = "Cells", title = current_snn_res,) +
  scale_fill_manual(values = cols)

Heatmap of mean gene expression by cell type

# split cells by cluster
cs_by_k <- split(colnames(sce), sce$cell_type)

# compute cluster-marker means
ms_by_cluster <- lapply(fs, function(gs)
    vapply(cs_by_k, function(i)
        Matrix::rowMeans(logcounts(sce)[gs, i, drop = FALSE]), 
        numeric(length(gs))))

# prep. for plotting & scale b/w 0 and 1
mat <- do.call("rbind", ms_by_cluster)
mat <- muscat:::.scale(mat)
rownames(mat) <- word(rownames(mat), 2, sep = "\\.")

cols <- muscat:::.cluster_colors[seq_along(anno)]
cols <- setNames(cols, names(anno))

row_anno <- rowAnnotation(
  df = data.frame(label = factor(ks, levels = names(fs))),
  col = list(label = cols),
  show_legend = FALSE, show_annotation_name = FALSE
)

col_anno <- columnAnnotation(
  df = data.frame(text = factor(names(fs), levels = names(fs))),
  col = list(text = cols),
  show_legend = FALSE, show_annotation_name = FALSE
)

Heatmap(mat,
    name = "scaled avg.\nexpression",
    col = viridis::magma(10)[3:10],
    cluster_rows = FALSE,
    cluster_columns = FALSE,
    row_names_side = "left",
    column_title = NULL,
    column_title_side = "bottom",
    column_names_rot = 90,
    row_names_gp =  gpar(fontfamily = "mono", fontface = "italic"),
    left_annotation = row_anno
    )

DimRed plots

so$cell_type <- sce$cell_type

cs <- sample(colnames(so), 5e3)
.plot_dr <- function(so, dr, id)
    DimPlot(so, cells = cs, group.by = id, reduction = dr, pt.size = 0.4) +
        guides(col = guide_legend(nrow = 11, 
            override.aes = list(size = 3, alpha = 1))) +
        theme_void() + theme(aspect.ratio = 1)

ids <- c("sample_id", "group_id", "cluster_id", "cell_type")
for (id in ids) {
    p1 <- .plot_dr(so, "tsne", id)
    p2 <- .plot_dr(so, "umap", id)
    if(id == "cell_type"){
      p1 <- p1 + scale_color_manual(values = cols)
      p2 <- p2 + scale_color_manual(values = cols)
    }
    lgd <- get_legend(p1)
    p1 <- p1 + theme(legend.position = "none")
    p2 <- p2 + theme(legend.position = "none")
    ps <- plot_grid(plotlist = list(p1, p2), nrow = 1)
    p <- plot_grid(ps, lgd, nrow = 1, rel_widths = c(1, 0.2))
    
    print(p)
}

Save SCE to .rds

# normalize for visualization
sce <- logNormCounts(sce)

# reorder sample levels
m <- match(levels(sce$sample_id), sce$sample_id)
o <- order(sce$group_id[m])
sce$sample_id <- factor(sce$sample_id,
                        levels = levels(sce$sample_id)[o])

# separate ensembl IDs & gene symbols
ss <- strsplit(rownames(sce), ".", fixed = TRUE)
rowData(sce)$ensembl_id <- sapply(ss, .subset, 1)
rowData(sce)$symbol <- sapply(ss, .subset, 2)

# prep. SCE for `muscat` & write to .rds
saveRDS(
  prepSCE(
    sce,
    cluster_id = "cluster_id",
    sample_id = "sample_id",
    group_id = "group_id"
  ),
  file.path(dir.objects, "SCE_annotation.rds")
)

—Checkpoint—

Load data

# load data, DS results & source utils
sce <- readRDS(file.path(dir.objects, "SCE_annotation.rds"))

Filter and process

table(sce$cell_type)

      Astrocytes          Neurons      Endothelial        Microglia Oligodendrocytes 
           75456             3600             1274              815             1595 
       Pericytes     Erythrocytes           Remove 
             478               29             7878 

# Filter astrocytes
sce.astro <- sce[, sce$cell_type == "Astrocytes"]

# create SeuratObject
so.astro <- CreateSeuratObject(counts = counts(sce.astro),
                               meta.data = data.frame(colData(sce.astro)))
so.astro$sex <- factor(so.astro$sex, levels = c("Female", "Male"))
so.astro$treatment <-
  factor(so.astro$treatment, levels = c("saline", "LPS"))
so.astro$group_id <-
  factor(
    so.astro$group_id,
    levels = c("Female_saline", "Male_saline", "Female_LPS", "Male_LPS")
  )

# split by sample
cells_by_sample <- split(colnames(sce.astro), sce.astro$sample_id)
so.astro <-
  lapply(cells_by_sample, function(i)
    subset(so.astro, cells = i))

# normalize, find variable genes, and scale
so.astro <- lapply(so.astro, NormalizeData, verbose = FALSE)
so.astro <- lapply(
  so.astro,
  FindVariableFeatures,
  nfeatures = 2e3,
  selection.method = "vst",
  verbose = FALSE
)
so.astro <- lapply(so.astro, ScaleData, verbose = FALSE)

# find anchors & re-integrate
as <- FindIntegrationAnchors(so.astro, verbose = FALSE)
so.astro <-
  IntegrateData(anchorset = as,
                dims = seq_len(30),
                verbose = FALSE)

# scale integrated data
DefaultAssay(so.astro) <- "integrated"
so.astro <- ScaleData(so.astro, verbose = FALSE)

# Run PCA
so.astro <- RunPCA(so.astro, npcs = 50, verbose = FALSE)
ElbowPlot(so.astro, ndims = 50)

# Run t-SNE
so.astro <- RunTSNE(
  so.astro,
  reduction = "pca",
  dims = seq_len(20),
  seed.use = 1,
  do.fast = TRUE,
  verbose = FALSE
)

DimPlot(so.astro, reduction = "tsne", group.by = "group_id")

Find clusters at different resolutions

# Remove columns with previously detected clusters
so.astro@meta.data[grep("integrated_snn", names(so.astro@meta.data))] <-
  NULL

# Re-detect clusters at given res
so.astro <-
  FindNeighbors(so.astro,
                reduction = "pca",
                dims = seq_len(20),
                verbose = FALSE)
for (res in c(0.1, 0.2, 0.4))
  so.astro <-
  FindClusters(
    so.astro,
    resolution = res,
    random.seed = 1,
    verbose = FALSE
  )

clusters_colnames <-
  str_subset(names(so.astro@meta.data), "integrated_snn_res")

tSNE.byCluster.plot.list <- list()
for (clusters_colname in clusters_colnames) {
  reordered_clusters_levels <-
    so.astro@meta.data[[clusters_colname]] %>%
    levels() %>% length() %>% seq(0, .)
  
  so.astro@meta.data[[clusters_colname]] <-
    factor(so.astro@meta.data[[clusters_colname]],
           levels = reordered_clusters_levels)
  
  tSNE.byCluster.plot.list[[clusters_colname]] <- DimPlot(
    so.astro,
    reduction = "tsne",
    group.by = clusters_colname,
    split.by = "group_id"
  ) + labs(title = clusters_colname)
}
tSNE.byCluster.plot.list

$integrated_snn_res.0.1

$integrated_snn_res.0.2

$integrated_snn_res.0.4