Warning message:
In system("timedatectl", intern = TRUE) :
running command 'timedatectl' had status 1
library(DropletUtils)
library(LSD)
library(Matrix)
library(readxl)
library(scater)
library(scds)
library(SingleCellExperiment)
library(scran)
library(Seurat)
library(muscat)
library(ComplexHeatmap)
library(cowplot)
library(purrr)
library(tidyverse)
library(RColorBrewer)
library(viridis)
library(grid)
library(gridExtra)
library(UpSetR)
wd <- "/home/rstudio/data"
dir.inputs <- file.path(wd, "sc-rnaseq/CellRanger_outputs/outs-cellranger")
dir.objects <- file.path(wd, "sc-rnaseq/objects")
dir.create(dir.objects, recursive = T, showWarnings = F)
# increase future's maximum allowed size of objects
# to be exported from default of 500 MB to 16 GB
options(future.globals.maxSize = 8 * 2048 * 1024 ^ 2)
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
# 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)
[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"
# 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.
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
[1] 12499 91104
[1] 12499 91104
Error in file.path(dir.objects, "SCE_preprocessing.rds") :
object 'dir.objects' not found
plotColData(sce, x = "sample_id", y = "detected", colour_by = "sum") +
geom_hline(yintercept = 700, color = "red") +
coord_flip()
# 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 <- 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)
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)
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
# 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"
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)))
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
)
In consideration of the above visualizations and additional exploration with iSEE, we arrive at the following cluster annotations:
# 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)
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
)
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_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_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)
)
)
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)
# 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
)
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)
}
# 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—
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")
# 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
FindAllMarkers
astro.markers <-
Seurat::FindAllMarkers(
so.astro,
logfc.threshold = log(0.35),
min.pct = 0.25,
only.pos = TRUE,
verbose = F
)
sparse->dense coercion: allocating vector of size 1.1 GiBsparse->dense coercion: allocating vector of size 1.1 GiBsparse->dense coercion: allocating vector of size 1.1 GiBsparse->dense coercion: allocating vector of size 1.1 GiBsparse->dense coercion: allocating vector of size 1.1 GiBsparse->dense coercion: allocating vector of size 1.2 GiBsparse->dense coercion: allocating vector of size 1.2 GiBsparse->dense coercion: allocating vector of size 1.2 GiB
astro.markers.subset <- astro.markers %>%
group_by(cluster) %>%
slice_max(n = 10, order_by = avg_logFC)
astro.markers.subset
—Checkpoint—
input.spatial <- file.path(wd, "spaceranger/SpaceRanger_outputs")
spatial.samples <- list.files(input.spatial, pattern = "_S[0-9]_") %>% str_extract("GSM[0-9]+_S[0-9]+") %>% unique()
list.files(input.spatial, pattern = "_S[0-9]_")
[1] "GSM5026144_S7_barcodes.tsv.gz" "GSM5026144_S7_features.tsv.gz" "GSM5026144_S7_filtered_feature_bc_matrix.h5"
[4] "GSM5026144_S7_matrix.mtx.gz" "GSM5026144_S7_module_scores.csv" "GSM5026144_S7_tissue_hires_image.png.gz"
[7] "GSM5026145_S3_barcodes.tsv.gz" "GSM5026145_S3_features.tsv.gz" "GSM5026145_S3_filtered_feature_bc_matrix.h5"
[10] "GSM5026145_S3_matrix.mtx.gz" "GSM5026145_S3_module_scores.csv" "GSM5026145_S3_tissue_hires_image.png.gz"
[13] "GSM5026146_S4_barcodes.tsv.gz" "GSM5026146_S4_features.tsv.gz" "GSM5026146_S4_filtered_feature_bc_matrix.h5"
[16] "GSM5026146_S4_matrix.mtx.gz" "GSM5026146_S4_module_scores.csv" "GSM5026146_S4_tissue_hires_image.png.gz"
[19] "GSM5026147_S1_barcodes.tsv.gz" "GSM5026147_S1_features.tsv.gz" "GSM5026147_S1_filtered_feature_bc_matrix.h5"
[22] "GSM5026147_S1_matrix.mtx.gz" "GSM5026147_S1_module_scores.csv" "GSM5026147_S1_tissue_hires_image.png.gz"
[25] "GSM5026148_S2_barcodes.tsv.gz" "GSM5026148_S2_features.tsv.gz" "GSM5026148_S2_filtered_feature_bc_matrix.h5"
[28] "GSM5026148_S2_matrix.mtx.gz" "GSM5026148_S2_module_scores.csv" "GSM5026148_S2_tissue_hires_image.png.gz"
[31] "GSM5026149_S6_barcodes.tsv.gz" "GSM5026149_S6_features.tsv.gz" "GSM5026149_S6_filtered_feature_bc_matrix.h5"
[34] "GSM5026149_S6_matrix.mtx.gz" "GSM5026149_S6_module_scores.csv" "GSM5026149_S6_tissue_hires_image.png.gz"
so.list <- list()
for (spatial.sample in spatial.samples) {
h5.file.name <-
paste0(spatial.sample, "_filtered_feature_bc_matrix.h5")
h5.file.path <- file.path(input.spatial, h5.file.name)
sample.dgCM <- Read10X_h5(h5.file.path)
so.list[[spatial.sample]] <-
CreateSeuratObject(counts = sample.dgCM, project = "Spatial")
so.list[[spatial.sample]] <-
SCTransform(object = so.list[[spatial.sample]], return.only.var.genes = F, verbose = FALSE)
}
astro.markers <- readRDS(file.path(dir.objects, "astro.markers.rds"))
astro.markers %>%
filter(p_val_adj < 0.01) %>%
ggplot(aes(x = avg_logFC, y = cluster, color = cluster, fill = cluster)) +
geom_density_ridges(alpha = 0.2, scale = 1) +
scale_x_continuous(trans = "log2")
astro.markers.list <- astro.markers %>%
filter(p_val_adj < 0.01) %>%
select(gene, cluster) %>%
mutate(gene = word(gene, 2, sep = "\\."))
astro.markers.list <- split(astro.markers.list$gene, astro.markers.list$cluster)
for(so.name in names(so.list)){
print(so.name)
so.list[[so.name]] <- AddModuleScore(so.list[[so.name]], astro.markers.list)
module_scores <- so.list[[so.name]]@meta.data[, grep("Cluster", colnames(so.list[[so.name]]@meta.data))]
csv.filename <- paste0(so.name, "_module_scores.csv")
write.csv(module_scores, file = file.path(input.spatial, csv.filename), row.names=TRUE)
}
[1] "GSM5026144_S7"
[1] "GSM5026145_S3"
[1] "GSM5026146_S4"
[1] "GSM5026147_S1"
[1] "GSM5026148_S2"
[1] "GSM5026149_S6"