Focal white matter lesions drive grey matter inflammation and synapse loss.

Code ↔ Paper

10 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 10 matches

[1] § Methods › DBiT-seq bioinformatics analysis ↔ DBiTseq_analysis/ST_Analysis_IO_20250929_final.Rmd, lines 702–789 · score 0.76 · get.knn, dilated, edgelist, node, radius, distances
[2] § Methods › DBiT-seq bioinformatics analysis ↔ DBiTseq_analysis/ST_Analysis_Lesion_20250929.Rmd, lines 702–789 · score 0.76 · get.knn, dilated, edgelist, node, radius, distances
[3] § Methods › Bulk RNA-seq deconvolution and bioinformatics analysis ↔ Bulk_deconvolution/Bulk_RNor2Mm_deconv_2025_EA_1.Rmd, lines 203–224 · score 0.73 · SampleID, SingleCellExperiment, SubClass, medulla, metadata, Bulk
[4] § Distinct transcriptional changes within the circuit › Focal white matter lesions evoke transcriptional changes in the IO ↔ DBiTseq_analysis/ST_Analysis_IO_20250929_final.Rmd, lines 3314–3341 · score 0.70 · Ajap1, Calm1, Cdh11, Cox2, Cox3, Gria1
[5] § Methods › DBiT-seq bioinformatics analysis ↔ DBiTseq_analysis/ST_Analysis_Lesion_20250929.Rmd, lines 1033–1060 · score 0.70 · FindClusters, FindNeighbors, elbow, vst, log, variable
[6] § Methods › DBiT-seq bioinformatics analysis ↔ DBiTseq_analysis/ST_analysis_Healthy.rmd, lines 1046–1073 · score 0.69 · FindClusters, FindNeighbors, elbow, vst, log, variable
[7] § Distinct transcriptional changes within the circuit › Focal white matter lesions evoke transcriptional changes in the IO ↔ DBiTseq_analysis/ST_Analysis_IO_20250929_final.Rmd, lines 2691–2741 · score 0.69 · retrograde endocannabinoid, oxidative phosphorylation, Parkinson, Alzheimer, disease, enrichment
[8] § Methods › DBiT-seq bioinformatics analysis ↔ DBiTseq_analysis/ST_Analysis_IO_20250929_final.Rmd, lines 3342–3365 · score 0.68 · AddModuleScore, Csf1r, nbin, Maf, Spi1, Aif1
[9] § Methods › DBiT-seq bioinformatics analysis ↔ Bulk_deconvolution/Bulk_RNor2Mm_deconv_2025_EA_1.Rmd, lines 146–183 · score 0.57 · FindClusters, FindNeighbors, Neighbour, UMAP, resolution, clustering
[10] § Methods › DBiT-seq bioinformatics analysis ↔ DBiTseq_analysis/ST_Analysis_Lesion_20250929.Rmd, lines 1033–1060 · score 0.54 · FindClusters, FindNeighbors, PCs, Neighbour, resolution, clustering

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

R Markdown · 5,962 lines · 195 KB · no license · 4 matches

---
title: "ST_Analysis_SeuratOnly_IO"
authors: "Bastien Herve, Joseph Wong"
output: html_document
date: "2025-09-29"
---
```{r}
reticulate::use_condaenv("r-reticulate", required = TRUE)
```
```{r}
library(Seurat)
library(Signac)
library(hdf5r)
library(biovizBase)
library(GenomicRanges)
library(stringr)
library(ggplot2)
library(sf)
library(scales)
library(grid)
library(gridExtra)
library(patchwork)
library(plyr)
library(dplyr)
library(Rmagic)
library(clustree)
library(pheatmap)
library(png)
library(imager)
library(RColorBrewer)
library(future)
library(pbapply)
library(cluster)
library(corrplot)
library(reticulate)
library(presto)
library(circlize)
library(ggcorrplot)
library(scatterpie)
library(glmGamPoi)
library(muscat)
library(DESeq2)
library(purrr)
library(Matrix.utils)
library(magrittr)
library(scales)
library(EnhancedVolcano)
library(UpSetR)
library(ReactomePA)
library(enrichplot)
library(biomaRt)
library(org.Mm.eg.db)
library(enrichR)
library(FNN)
library(igraph)
# additional
# library(MASS)
# library(data.table)
```
```{r}
sessionInfo()
```
```{r}
# set up parameters
seed <- 424242
options(repr.plot.width=16, repr.plot.height=12)
options(future.globals.maxSize = 110000 * 1024^2)
options(repr.matrix.max.rows=100, repr.matrix.max.cols=100)
pixel_size_plot = 2.5
empty_lanes_cutoff_mean <- 500
nb_pixels_dbit = 10000
names_Ctrl <- c("Ctrlrep1", "Ctrlrep2", "Ctrlrep3", "Ctrlrep4")
names_D7 <- c("D7rep1", "D7rep2", "D7rep3", "D7rep4")
names_D14 <- c("D14rep1", "D14rep2", "D14rep3", "D14rep4")
names_D28 <- c("D28rep1", "D28rep2", "D28rep3")
sample_names <- c("P35102_1001","P35102_1002","P35102_1003","P35102_1004","P35102_1005","P35102_1006","P35102_1007","P35102_1008","P35102_1009","P35102_1010","P35102_1011","P35102_1012","P35102_1013","P35102_1014","P35102_1015")
names(sample_names) <- c(names_Ctrl,names_D7,names_D14,names_D28)
```
```{r}
orig.ident_colors <- c("#536553", "#657B65", "#8FA38F", "#BCC8BC",
"#B88100", "#F5AB00", "#FFD470", "#FFE7AD",
"#E92A0C", "#F5563D", "#F88877", "#FDDDD8",
"#97615E", "#C8A9A7", "#E0CECD")
names(orig.ident_colors) <- c("Ctrlrep1", "Ctrlrep2", "Ctrlrep3", "Ctrlrep4",
"D7rep1", "D7rep2", "D7rep3", "D7rep4",
"D14rep1", "D14rep2", "D14rep3", "D14rep4",
"D28rep1", "D28rep2", "D28rep3")
show_col(orig.ident_colors)
orig.ident_colors_contrast <- c("#97615E", "#E92A0C", "#F5AB00", "#6E876E",
"#97615E", "#E92A0C", "#F5AB00", "#6E876E",
"#97615E", "#E92A0C", "#F5AB00", "#6E876E",
"#97615E", "#E92A0C", "#F5AB00")
names(orig.ident_colors_contrast) <- c("Ctrlrep1", "Ctrlrep2", "Ctrlrep3", "Ctrlrep4",
"D7rep1", "D7rep2", "D7rep3", "D7rep4",
"D14rep1", "D14rep2", "D14rep3", "D14rep4",
"D28rep1", "D28rep2", "D28rep3")
show_col(orig.ident_colors_contrast)
orig.ident_merge_colors <- c("#6E876E", "#F5AB00", "#E92A0C", "#97615E")
names(orig.ident_merge_colors) <- c("Ctrl", "D7", "D14", "D28")
show_col(orig.ident_merge_colors)
model_colors <- c("#6E876E", "#E92A0C")
names(model_colors) <- c("Ctrl", "Lesion")
show_col(model_colors)
duo_colors <- c("#E92A0C", "#2A75CB")
show_col(duo_colors)
repair_colors <- duo_colors
names(repair_colors) <- c("repaired", "non repaired")
wheeler_colors <- c("#3772FF", "#160F29", "#A9FFCB", "#FDCA40", "#A77E58", "#DE6449", "#C03221", "#791E94", "#EFC3E6", "#0CCA4A", "#D78521", "#000080", "#C3C9E9", "#FFFF00")
names(wheeler_colors) <- c("Astro", "DC", "Endo", "Epen", "Epith", "MicroMacro", "MonoMacro", "Neurons", "Neutro", "OL", "Peri", "SMC", "Stromal", "Tcells")
show_col(wheeler_colors)
zeisel_colors <- c("#3772FF", "#FDCA40", "#A77E58", "#DE6449", "#C03221", "#791E94", "#0CCA4A", "#D78521", "#000080")
names(zeisel_colors) <- c("Astrocyte", "Ependymal", "Blood", "Immune", "Bergmann-glia", "Neurons", "Oligos", "Ttr", "Vascular")
show_col(zeisel_colors)
```
```{r}
PrctCellExpringGene <- function(object, genes, group.by = "all"){
if(group.by == "all"){
prct = unlist(lapply(genes,calc_helper, object=object))
result = data.frame(Markers = genes, Cell_proportion = prct)
return(result)
}
else{
list = SplitObject(object, group.by)
factors = names(list)
results = lapply(list, PrctCellExpringGene, genes=genes)
for(i in 1:length(factors)){
results[[i]]$Feature = factors[i]
}
combined = do.call("rbind", results)
return(combined)
}
}
calc_helper <- function(object,genes){
counts = object[['RNA']]@counts
ncells = ncol(counts)
if(genes %in% row.names(counts)){
sum(counts[genes,]>0)/ncells
}else{return(NA)}
}
NbCellExpringGene <- function(object, genes, group.by = "all"){
if(group.by == "all"){
nb = unlist(lapply(genes,calc_helper_nb, object=object))
result = data.frame(Markers = genes, Cell_nb = nb)
return(result)
}
else{
list = SplitObject(object, group.by)
factors = names(list)
results = lapply(list, NbCellExpringGene, genes=genes)
for(i in 1:length(factors)){
results[[i]]$Feature = factors[i]
}
combined = do.call("rbind", results)
return(combined)
}
}
calc_helper_nb <- function(object,genes){
counts = object[['RNA']]@counts
if(genes %in% row.names(counts)){
length(counts[genes,][counts[genes,] > 0])
}else{return(NA)}
}
DoMultiBarHeatmap <- function (object, features = NULL, cells = NULL, group.by = "ident", additional.group.by = NULL, additional.group.sort.by = NULL, cols.use = NULL, group.bar = TRUE, disp.min = -2.5, disp.max = NULL, layer = "scale.data", assay = NULL, label = TRUE, size = 5.5, hjust = 0, angle = 45, raster = TRUE, draw.lines = TRUE, lines.width = NULL, group.bar.height = 0.02, combine = TRUE){
cells <- cells %||% colnames(x = object)
if (is.numeric(x = cells)) {
cells <- colnames(x = object)[cells]
}
assay <- assay %||% DefaultAssay(object = object)
DefaultAssay(object = object) <- assay
features <- features %||% VariableFeatures(object = object)
## Why reverse???
features <- rev(x = unique(x = features))
disp.max <- disp.max %||% ifelse(test = layer == "scale.data", yes = 2.5, no = 6)
possible.features <- rownames(x = LayerData(object = object, assay = assay, layer = layer))
if (any(!features %in% possible.features)) {
bad.features <- features[!features %in% possible.features]
features <- features[features %in% possible.features]
if (length(x = features) == 0) {
stop("No requested features found in the ", layer,
" slot for the ", assay, " assay.")
}
warning("The following features were omitted as they were not found in the ", layer, " slot for the ", assay, " assay: ", paste(bad.features, collapse = ", "))
}
if (!is.null(additional.group.sort.by)) {
if (any(!additional.group.sort.by %in% additional.group.by)) {
bad.sorts <- additional.group.sort.by[!additional.group.sort.by %in% additional.group.by]
additional.group.sort.by <- additional.group.sort.by[additional.group.sort.by %in% additional.group.by]
if (length(x = bad.sorts) > 0) {
warning("The following additional sorts were omitted as they were not a subset of additional.group.by : ",
paste(bad.sorts, collapse = ", "))
}
}
}
data <- as.data.frame(x = as.matrix(x = t(x = LayerData(object = object, assay = assay, layer = layer)[features, cells, drop = FALSE])))
object <- suppressMessages(expr = StashIdent(object = object, save.name = "ident"))
group.by <- group.by %||% "ident"
groups.use <- object[[c(group.by, additional.group.by[!additional.group.by %in% group.by])]][cells, , drop = FALSE]
plots <- list()
for (i in group.by) {
data.group <- data
if (!is.null(additional.group.by)) {
additional.group.use <- additional.group.by[additional.group.by!=i]
if (!is.null(additional.group.sort.by)){
additional.sort.use = additional.group.sort.by[additional.group.sort.by != i]
} else {
additional.sort.use = NULL
}
} else {
additional.group.use = NULL
additional.sort.use = NULL
}
group.use <- groups.use[, c(i, additional.group.use), drop = FALSE]
for(colname in colnames(group.use)){
if (!is.factor(x = group.use[[colname]])) {
group.use[[colname]] <- factor(x = group.use[[colname]])
}
}
if (draw.lines) {
lines.width <- lines.width %||% ceiling(x = nrow(x = data.group) * 0.0025)
placeholder.cells <- sapply(X = 1:(length(x = levels(x = group.use[[i]])) * lines.width), FUN = function(x) {
return(Seurat:::RandomName(length = 20))
})
placeholder.groups <- data.frame(rep(x = levels(x = group.use[[i]]), times = lines.width))
group.levels <- list()
group.levels[[i]] = levels(x = group.use[[i]])
for (j in additional.group.use) {
group.levels[[j]] <- levels(x = group.use[[j]])
placeholder.groups[[j]] = NA
}
colnames(placeholder.groups) <- colnames(group.use)
rownames(placeholder.groups) <- placeholder.cells
group.use <- sapply(group.use, as.vector)
rownames(x = group.use) <- cells
group.use <- rbind(group.use, placeholder.groups)
for (j in names(group.levels)) {
group.use[[j]] <- factor(x = group.use[[j]], levels = group.levels[[j]])
}
na.data.group <- matrix(data = NA, nrow = length(x = placeholder.cells), ncol = ncol(x = data.group), dimnames = list(placeholder.cells, colnames(x = data.group)))
data.group <- rbind(data.group, na.data.group)
}
order_expr <- paste0('order(', paste(c(i, additional.sort.use), collapse=','), ')')
group.use = with(group.use, group.use[eval(parse(text=order_expr)), , drop=F])
plot <- Seurat:::SingleRasterMap(data = data.group, raster = raster, disp.min = disp.min, disp.max = disp.max, feature.order = features, cell.order = rownames(x = group.use), group.by = group.use[[i]])
if (group.bar) {
pbuild <- ggplot_build(plot = plot)
group.use2 <- group.use
cols <- list()
na.group <- Seurat:::RandomName(length = 20)
for (colname in rev(x = colnames(group.use2))) {
if (colname == i) {
colid = paste0('Identity (', colname, ')')
} else {
colid = colname
}
# Default
cols[[colname]] <- c(scales::hue_pal()(length(x = levels(x = group.use[[colname]]))))
#Overwrite if better value is provided
if (!is.null(cols.use[[colname]])) {
req_length = length(x = levels(group.use))
if (length(cols.use[[colname]]) < req_length){
warning("Cannot use provided colors for ", colname, " since there aren't enough colors.")
} else {
if (!is.null(names(cols.use[[colname]]))) {
if (all(levels(group.use[[colname]]) %in% names(cols.use[[colname]]))) {
cols[[colname]] <- as.vector(cols.use[[colname]][levels(group.use[[colname]])])
} else {
warning("Cannot use provided colors for ", colname, " since all levels (", paste(levels(group.use[[colname]]), collapse=","), ") are not represented.")
}
} else {
cols[[colname]] <- as.vector(cols.use[[colname]])[c(1:length(x = levels(x = group.use[[colname]])))]
}
}
}
# Add white if there's lines
if (draw.lines) {
levels(x = group.use2[[colname]]) <- c(levels(x = group.use2[[colname]]), na.group)
group.use2[placeholder.cells, colname] <- na.group
cols[[colname]] <- c(cols[[colname]], "#FFFFFF")
}
names(x = cols[[colname]]) <- levels(x = group.use2[[colname]])
y.range <- diff(x = pbuild$layout$panel_params[[1]]$y.range)
y.pos <- max(pbuild$layout$panel_params[[1]]$y.range) + y.range * 0.015
y.max <- y.pos + group.bar.height * y.range
pbuild$layout$panel_params[[1]]$y.range <- c(pbuild$layout$panel_params[[1]]$y.range[1], y.max)
plot <- suppressMessages(plot +
annotation_raster(raster = t(x = cols[[colname]][group.use2[[colname]]]), xmin = -Inf, xmax = Inf, ymin = y.pos, ymax = y.max) +
annotation_custom(grob = grid::textGrob(label = colid, hjust = 0, gp = gpar(cex = 0.75)), ymin = mean(c(y.pos, y.max)), ymax = mean(c(y.pos, y.max)), xmin = Inf, xmax = Inf) +
coord_cartesian(ylim = c(0, y.max), clip = "off")
)
if ((colname == i) && label) {
x.max <- max(pbuild$layout$panel_params[[1]]$x.range)
x.divs <- pbuild$layout$panel_params[[1]]$x.major %||% pbuild$layout$panel_params[[1]]$x$break_positions()
group.use$x <- x.divs
label.x.pos <- tapply(X = group.use$x, INDEX = group.use[[colname]], FUN = median) * x.max
label.x.pos <- data.frame(group = names(x = label.x.pos), label.x.pos)
plot <- plot + geom_text(stat = "identity", data = label.x.pos, aes_string(label = "group", x = "label.x.pos"), y = y.max + y.max * 0.03 * 0.5, angle = angle, hjust = hjust, size = size)
plot <- suppressMessages(plot + coord_cartesian(ylim = c(0, y.max + y.max * 0.002 * max(nchar(x = levels(x = group.use[[colname]]))) * size), clip = "off"))
}
}
}
plot <- plot + theme(line = element_blank())
plots[[i]] <- plot
}
if (combine) {
plots <- CombinePlots(plots = plots)
}
return(plots)
}
plot_count_on_tissue <- function(object=NULL, img_data=NULL, gene=NULL, assay=NULL, slot=NULL, image=NULL, readPNG_dim=NULL, transp=NULL, limit_min=NULL, limit_max=NULL, pixel_size=NULL, palette_col=NULL){
DefaultAssay(object) <- assay
img_data[,"value"] <- GetAssayData(object = object, slot = slot)[gene,img_data[,"cell_names"]]
img_data_pixel <- as.data.frame(img_data) %>%
st_as_sf(coords = c("x1", "y1")) %>%
group_by(x,y,var,value) %>%
dplyr::summarise(do_union = F) %>%
st_cast("POLYGON") %>%
st_cast("MULTIPOLYGON") %>%
mutate(var = as.factor(var))
ggplot(img_data_pixel) +
annotation_custom(image) +
geom_sf(aes(fill = value),shape=15, size = pixel_size, alpha = transp) +
scale_fill_gradientn(colours = palette_col, limits = c(limit_min,limit_max), oob=squish) +
scale_x_continuous(expand = c(0, 0), limits = c(0, readPNG_dim[2])) +
scale_y_continuous(expand = c(0, 0), limits = c(0, readPNG_dim[1])) +
theme_void()
}
scale01 <- function(x){(x-min(x))/(max(x)-min(x))}
plot_diff_count_on_tissue <- function(object=NULL, gene=NULL, assay=NULL, slot=NULL, image=NULL, rescale=NULL, transp=NULL, limit=NULL, pixel_size=NULL, palette_col=NULL){
if(length(gene) > 1 & length(assay) > 1){
print("please only compare 2 assays or 2 genes not both")
} else{
if (length(gene) > 1){
DefaultAssay(object) <- assay
df <- data.frame(
x = object$x,
y = object$y,
value1 = GetAssayData(object = object, slot = slot)[gene[1],],
value2 = GetAssayData(object = object, slot = slot)[gene[2],]
)
gene_title <- paste0(gene[1],"/",gene[2])
assay_title <- assay
}else if (length(assay) > 1){
df <- data.frame(
x = object$x,
y = object$y
)
DefaultAssay(object) <- assay[1]
df$value1 = GetAssayData(object = object, slot = slot)[gene,]
DefaultAssay(object) <- assay[2]
df$value2 = GetAssayData(object = object, slot = slot)[gene,]
gene_title <- gene
assay_title <- paste0(assay[1],"/",assay[2])
}
df$value1 <- scale01(df$value1)
df$value2 <- scale01(df$value2)
df$value <- df$value1-df$value2
df$value1 <- NULL
df$value2 <- NULL
if(!is.null(rescale)){
df$x=rescale(df$x, to=c(rescale[1],rescale[2]))
df$y=rescale(df$y, to=c(rescale[3],rescale[4]))
}
ggplot(df, aes(x, y)) +
annotation_custom(image) +
geom_point(aes(colour = value),shape=15, size = pixel_size, alpha = transp) +
scale_colour_gradientn(colours = palette_col, limits = c(-limit,limit), oob=squish) +
scale_x_continuous(expand=c(-0.005,0), lim=c(0,101)) +
scale_y_continuous(expand=c(-0.006,0), lim=c(0,101)) +
theme(axis.line=element_blank(),axis.text.x=element_blank(),
axis.text.y=element_blank(),axis.ticks=element_blank(),
axis.title.x=element_blank(),axis.title.y=element_blank(),
panel.background=element_blank(),panel.border=element_blank(),panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),plot.background=element_blank()) +
ggtitle(paste0(gene_title,"_",assay_title))
}
}
ExportGroupBW <- function(
object,
assay = NULL,
group.by = NULL,
idents = NULL,
normMethod = "RC",
tileSize = 100,
minCells = 5,
cutoff = NULL,
chromosome = NULL,
outdir = NULL,
verbose=TRUE
) {
# Check if temporary directory exist
if (!dir.exists(outdir)){
dir.create(outdir)
}
if (!requireNamespace("rtracklayer", quietly = TRUE)) {
message("Please install rtracklayer. http://www.bioconductor.org/packages/rtracklayer/")
return(NULL)
}
assay <- SetIfNull(x = assay, y = DefaultAssay(object = object))
DefaultAssay(object = object) <- assay
group.by <- SetIfNull(x = group.by, y = 'ident')
Idents(object = object) <- group.by
idents <- SetIfNull(x = idents, y = levels(x = object))
GroupsNames <- names(x = table(object[[group.by]])[table(object[[group.by]]) > minCells])
GroupsNames <- GroupsNames[GroupsNames %in% idents]
# Check if output files already exist
lapply(X = GroupsNames, FUN = function(x) {
fn <- paste0(outdir, .Platform$file.sep, x, ".bed")
if (file.exists(fn)) {
message(sprintf("The group \"%s\" is already present in the destination folder and will be overwritten !",x))
file.remove(fn)
}
})
# Splitting fragments file for each idents in group.by
SplitFragments(
object = object,
assay = assay,
group.by = group.by,
idents = idents,
outdir = outdir,
file.suffix = "",
append = TRUE,
buffer_length = 256L,
verbose = verbose
)
# Column to normalized by
if(!is.null(x = normMethod)) {
if (tolower(x = normMethod) %in% c('rc', 'ncells', 'none')){
normBy <- normMethod
} else{
normBy <- object[[normMethod, drop = FALSE]]
}
}
# Get chromosome information
if(!is.null(x = chromosome)){
seqlevels(object) <- chromosome
}
availableChr <- names(x = seqlengths(object))
chromLengths <- seqlengths(object)
chromSizes <- GRanges(
seqnames = availableChr,
ranges = IRanges(
start = rep(1, length(x = availableChr)),
end = as.numeric(x = chromLengths)
)
)
if (verbose) {
message("Creating tiles")
}
# Create tiles for each chromosome, from GenomicRanges
tiles <- unlist(
x = slidingWindows(x = chromSizes, width = tileSize, step = tileSize)
)
if (verbose) {
message("Creating bigwig files at ", outdir)
}
# Run the creation of bigwig for each cellgroups
if (nbrOfWorkers() > 1) {
mylapply <- future_lapply
} else {
mylapply <- ifelse(test = verbose, yes = pblapply, no = lapply)
}
covFiles <- mylapply(
GroupsNames,
FUN = CreateBWGroup,
availableChr,
chromLengths,
tiles,
normBy,
tileSize,
normMethod,
cutoff,
outdir
)
return(covFiles)
}
CreateBWGroup <- function(
groupNamei,
availableChr,
chromLengths,
tiles,
normBy,
tileSize,
normMethod,
cutoff,
outdir
) {
if (!requireNamespace("rtracklayer", quietly = TRUE)) {
message("Please install rtracklayer. http://www.bioconductor.org/packages/rtracklayer/")
return(NULL)
}
normMethod <- tolower(x = normMethod)
# Read the fragments file associated to the group
fragi <- rtracklayer::import(
paste0(outdir, .Platform$file.sep, groupNamei, ".bed"), format = "bed"
)
cellGroupi <- unique(x = fragi$name)
# Open the writing bigwig file
covFile <- file.path(
outdir,
paste0(groupNamei, "-TileSize-",tileSize,"-normMethod-",normMethod,".bw")
)
covList <- lapply(X = seq_along(availableChr), FUN = function(k) {
fragik <- fragi[seqnames(fragi) == availableChr[k],]
tilesk <- tiles[BiocGenerics::which(S4Vectors::match(seqnames(tiles), availableChr[k], nomatch = 0) > 0)]
if (length(x = fragik) == 0) {
tilesk$reads <- 0
# If fragments
} else {
# N Tiles
nTiles <- chromLengths[availableChr[k]] / tileSize
# Add one tile if there is extra bases
if (nTiles%%1 != 0) {
nTiles <- trunc(x = nTiles) + 1
}
# Create Sparse Matrix
matchID <- S4Vectors::match(mcols(fragik)$name, cellGroupi)
# For each tiles of this chromosome, create start tile and end tile row,
# set the associated counts matching with the fragments
mat <- sparseMatrix(
i = c(trunc(x = start(x = fragik) / tileSize),
trunc(x = end(x = fragik) / tileSize)) + 1,
j = as.vector(x = c(matchID, matchID)),
x = rep(1, 2*length(x = fragik)),
dims = c(nTiles, length(x = cellGroupi))
)
# Max count for a cells in a tile is set to cutoff
if (!is.null(x = cutoff)){
mat@x[mat@x > cutoff] <- cutoff
}
# Sums the cells
mat <- rowSums(x = mat)
tilesk$reads <- mat
# Normalization
if (!is.null(x = normMethod)) {
if (normMethod == "rc") {
tilesk$reads <- tilesk$reads * 10^4 / length(fragi$name)
} else if (normMethod == "ncells") {
tilesk$reads <- tilesk$reads / length(cellGroupi)
} else if (normMethod == "none") {
} else {
if (!is.null(x = normBy)){
tilesk$reads <- tilesk$reads * 10^4 / sum(normBy[cellGroupi, 1])
}
}
}
}
tilesk <- coverage(tilesk, weight = tilesk$reads)[[availableChr[k]]]
tilesk
})
names(covList) <- availableChr
covList <- as(object = covList, Class = "RleList")
rtracklayer::export.bw(object = covList, con = covFile)
return(covFile)
}
SetIfNull <- function(x, y) {
if (is.null(x = x)) {
return(y)
} else {
return(x)
}
}
repair_lanes <- function(seurat_obj, assay = "RNA", idents=NULL, to_fix_x, to_fix_y) {
if(!is.null(idents)){
seurat_tmp <- subset(seurat_obj, subset=orig.ident==idents)
}else{
seurat_tmp <- colnames(seurat_obj)
}
mat <- seurat_tmp[["RNA"]]$counts
coord <- [email hidden][, c("x","y")]
lengths_x <- sum(lengths(lapply(to_fix_x, function(x){sort(unique(coord[coord$x == x, ]$y))})))
lengths_y <- sum(lengths(lapply(to_fix_y, function(y){sort(unique(coord[coord$y == y, ]$x))})))
pb = txtProgressBar(min = 0, max = sum(lengths_x, lengths_y), style = 3)
pb_counter=0
for(to_fix in to_fix_x){
all_y_list <- sort(unique(coord[coord$x == to_fix, ]$y))
for (i in all_y_list) {
pb_counter=pb_counter+1
setTxtProgressBar(pb,pb_counter)
if (i==1){
idx_to_average <- c(which(coord$x == to_fix + 1 & coord$y == i))
}else if(i==sqrt(nb_pixels_dbit)){
idx_to_average <- c(which(coord$x == to_fix - 1 & coord$y == i))
}else{
idx_to_average <- c(which(coord$x == to_fix - 1 & coord$y == i), which(coord$x == to_fix + 1 & coord$y == i))
}
idx_to_fix <- which(coord$x == to_fix & coord$y == i)
if (length(idx_to_average) > 0) {
mat[, idx_to_fix] <- round(rowMeans(mat[, idx_to_average, drop = FALSE]))
}
}
}
for(to_fix in to_fix_y){
all_x_list <- sort(unique(coord[coord$y == to_fix, ]$x))
for (i in all_x_list) {
pb_counter=pb_counter+1
setTxtProgressBar(pb,pb_counter)
if (i==1){
idx_to_average <- c(which(coord$y == to_fix + 1 & coord$x == i))
}else if(i==sqrt(nb_pixels_dbit)){
idx_to_average <- c(which(coord$y == to_fix - 1 & coord$x == i))
}else{
idx_to_average <- c(which(coord$y == to_fix - 1 & coord$x == i), which(coord$y == to_fix + 1 & coord$x == i))
}
idx_to_fix <- which(coord$y == to_fix & coord$x == i)
if (length(idx_to_average) > 0) {
mat[, idx_to_fix] <- round(rowMeans(mat[, idx_to_average, drop = FALSE]))
}
}
}
close(pb)
return(mat)
}
save_pheatmap_pdf <- function(x, filename, width=7, height=7) {
stopifnot(!missing(x))
stopifnot(!missing(filename))
pdf(filename, width=width, height=height)
grid::grid.newpage()
grid::grid.draw(x$gtable)
dev.off()
}
order_heatmap_genes <- function(x){
if(ncol(x) == 1 | nrow(x) == 1) { return(rownames(x))
} else {
xs <- split(x, apply(x, 1, function(i) colnames(x)[which.max(i)]))
if(ncol(x) == length(orig.ident_merge_colors)-1) {
for(tp in names(xs)){
xs[[tp]] <- xs[[tp]][order(xs[[tp]][,tp], decreasing=TRUE),]
}
}
sapply(names(xs), function(i){
order_heatmap_genes(xs[[i]][, colnames(x) != i, drop = FALSE ])
}, simplify = FALSE)
}
}
dilate <- function(mat, radius = 1) {
padded <- matrix(0, nrow = nrow(mat) + 2 * radius, ncol = ncol(mat) + 2 * radius)
padded[(radius + 1):(nrow(padded) - radius), (radius + 1):(ncol(padded) - radius)] <- mat
result <- matrix(0, nrow = nrow(mat), ncol = ncol(mat))
for (i in 1:nrow(mat)) {
for (j in 1:ncol(mat)) {
neighborhood <- padded[i:(i + 2*radius), j:(j + 2*radius)]
if (sum(neighborhood) > 0) {
result[i, j] <- 1
}
}
}
return(result)
}
cluster_compartment <- function(coords, k, dilation_radius, percentile_thresh, extra_clust_ratio, grid_size){
knn_res <- get.knn(coords, k = k)
edges <- do.call(rbind, lapply(1:nrow(coords), function(i) {
cbind(i, knn_res$nn.index[i, ])
}))
g <- graph_from_edgelist(edges, directed = FALSE)
clusters <- components(g)
main_cluster_id <- which.max(clusters$csize)
potential_extra_cluster <- which(clusters$csize>max(clusters$csize)*extra_clust_ratio)
potential_extra_cluster <- unique(c(potential_extra_cluster,main_cluster_id))
main_cluster_nodes <- which(clusters$membership %in% c(potential_extra_cluster))
main_coords <- coords[main_cluster_nodes, ]
avg_dists <- rowMeans(knn_res$nn.dist)
threshold <- quantile(avg_dists, percentile_thresh)
core_indices <- which(avg_dists <= threshold)
core_coords <- main_coords[core_indices, ]
avg_dists <- rowMeans(knn_res$nn.dist[main_cluster_nodes,])
threshold <- quantile(avg_dists, percentile_thresh)
keep_indices <- which(avg_dists <= threshold)
filtered_coords <- main_coords[keep_indices, ]
grid <- matrix(0, nrow = grid_size, ncol = grid_size)
for (i in 1:nrow(filtered_coords)) {
x <- filtered_coords[i, "x"]
y <- filtered_coords[i, "y"]
if (x >= 1 && x <= grid_size && y >= 1 && y <= grid_size) {
grid[y, x] <- 1
}
}
grid_filled <- dilate(grid, radius = dilation_radius)
my_mat <- t(grid_filled)
d1 <- expand.grid(x = 1:grid_size, y = 1:grid_size)
out <- transform(d1, value = my_mat[as.matrix(d1)])
return(out[out$value==1,])
}
barplot_with_error_and_dots <- function(data, x, y, error = c("se", "sd"), title) {
error <- match.arg(error)
x <- enquo(x)
y <- enquo(y)
summary_data <- data %>%
group_by(!!x) %>%
summarise(
mean = mean(!!y, na.rm = TRUE),
sd = sd(!!y, na.rm = TRUE),
se = sd / sqrt(n()),
.groups = 'drop'
) %>%
mutate(error_value = ifelse(error == "sd", sd, se))
# Base plot
p <- ggplot(data, aes(x = !!x, y = !!y, fill = !!x)) +
geom_bar(data = summary_data,
aes(y = mean, fill = !!x),
stat = "identity",
width = 0.6,
color = "black") +
scale_fill_manual(values = orig.ident_merge_colors) +
geom_errorbar(data = summary_data,
aes(y = mean,
ymin = mean - error_value,
ymax = mean + error_value),
width = 0.2) +
geom_jitter(color = "black",
width = 0.15, height = 0,
size = 2, alpha = 0.8) +
theme_minimal() +
labs(y = title, x = "") +
theme(legend.position = "none")
return(p)
}
```
```{r}
options(repr.plot.width=5, repr.plot.height=5)
color_palette_30 <- c(
"#ffff00",
"#8b4513",
"#483d8b",
"#008000",
"#008b8b",
"#4682b4",
"#000080",
"#daa520",
"#7f007f",
"#8fbc8f",
"#b03060",
"#ff4500",
"#00ff00",
"#deb887",
"#556b2f",
"#9400d3",
"#00ff7f",
"#dc143c",
"#00ffff",
"#0000ff",
"#adff2f",
"#1e90ff",
"#fa8072",
"#90ee90",
"#add8e6",
"#ff1493",
"#7b68ee",
"#ee82ee",
"#ffb6c1")
seurat_clusters_color <- list()
n_clusters <- 22
seurat_clusters_color <- color_palette_30[1:n_clusters]
names(seurat_clusters_color) <- seq(1:n_clusters)-1
show_col(seurat_clusters_color)
```
```{r}
seurat_obj <- readRDS("/rds/project/rds-qqstKNzWvy8/Omar_Spatial/P35102_seurat_obj_clustering_DIET_20250819.rds")
```
```{r}
DefaultAssay(seurat_obj) <- "RNA"
seurat_obj <- JoinLayers(seurat_obj)
seurat_obj <- NormalizeData(seurat_obj, normalization.method = "LogNormalize", scale.factor = 10000)
DefaultAssay(seurat_obj) <- "RNA_repaired"
seurat_obj <- NormalizeData(seurat_obj, normalization.method = "LogNormalize", scale.factor = 10000)
```
```{r}
for (sample_name in names(orig.ident_colors)){
annot_plot <- [email hidden][seurat_obj$orig.ident==sample_name,]
options(repr.plot.width=5.3, repr.plot.height=4.6)
df_tmp <- data.frame(
x = annot_plot$x,
y = annot_plot$y,
value = annot_plot$seurat_clusters
)
print(ggplot(df_tmp, aes(x, y)) +
geom_point(aes(colour = value),shape=15, size = 1)+
scale_colour_manual(values=seurat_clusters_color) +
scale_x_continuous(expand=c(-0.005,0), lim=c(0,101)) +
scale_y_continuous(expand=c(-0.006,0), lim=c(0,101)) +
theme_void() +
ggtitle(sample_name))
}
```
```{r}
seurat_obj$InfOlive <- "NA"
for (sample_name in names(orig.ident_colors)){
message(sample_name)
annot_plot_sample <- [email hidden][seurat_obj$orig.ident==sample_name,]
annot_plot_cluster <- annot_plot_sample[annot_plot_sample$seurat_clusters==5,]
#dilation of 5 for the core+pheriphery
mat_peri_compartment <- cluster_compartment(annot_plot_cluster[,c("x","y")],k=5,dilation_radius=5,percentile_thresh=0.5,extra_clust_ratio=0.5,grid_size=sqrt(nb_pixels_dbit))
#dilation of 2 for the core
mat_core_compartment <- cluster_compartment(annot_plot_cluster[,c("x","y")],k=5,dilation_radius=2,percentile_thresh=0.5,extra_clust_ratio=0.5,grid_size=sqrt(nb_pixels_dbit))
# Record periphery first
for(line in 1:nrow(mat_peri_compartment)){
annot_plot_sample[annot_plot_sample$x==mat_peri_compartment[line,"x"] & annot_plot_sample$y==mat_peri_compartment[line,"y"], "InfOlive"] <- "periphery"
}
# Record core second
for(line in 1:nrow(mat_core_compartment)){
annot_plot_sample[annot_plot_sample$x==mat_core_compartment[line,"x"] & annot_plot_sample$y==mat_core_compartment[line,"y"], "InfOlive"] <- "core"
}
[email hidden][rownames(annot_plot_sample),"InfOlive"] <- annot_plot_sample[,"InfOlive", drop=FALSE]
}
```
```{r}
for (sample_name in names(orig.ident_colors)){
annot_plot <- [email hidden][seurat_obj$orig.ident==sample_name,]
options(repr.plot.width=5.3, repr.plot.height=4.6)
df_tmp <- data.frame(
x = annot_plot$x,
y = annot_plot$y,
value = annot_plot$InfOlive
)
print(ggplot(df_tmp, aes(x, y)) +
geom_point(aes(colour = value),shape=15, size = 1)+
scale_x_continuous(expand=c(-0.005,0), lim=c(0,101)) +
scale_y_continuous(expand=c(-0.006,0), lim=c(0,101)) +
theme_void() +
ggtitle(sample_name))
}
```
# Process the inferior olive
```{r}
# read the rds
seurat_obj_IO <- readRDS("~/rds/rds-rds-karadottir-qqstKNzWvy8/Omar_Spatial/P35102_seurat_obj_IO_DIET.rds")
orig.ident_colors <- orig.ident_colors[!names(orig.ident_colors) %like% "D28"]
orig.ident_merge_colors <- orig.ident_merge_colors[names(orig.ident_merge_colors) != "D28"]
seurat_obj_IO <- subset(seurat_obj_IO, subset=orig.ident_merge!="D28")
seurat_obj_IO$orig.ident <- factor(x = seurat_obj_IO$orig.ident, levels = names(orig.ident_colors))
seurat_obj_IO$orig.ident_merge <- factor(x = seurat_obj_IO$orig.ident_merge, levels = names(orig.ident_merge_colors))
```
```{r}
DefaultAssay(seurat_obj_IO) <- "RNA"
seurat_obj_IO <- NormalizeData(seurat_obj_IO, normalization.method = "LogNormalize", scale.factor = 10000)
DefaultAssay(seurat_obj_IO) <- "RNA_repaired"
seurat_obj_IO <- NormalizeData(seurat_obj_IO, normalization.method = "LogNormalize", scale.factor = 10000)
seurat_obj_IO <- FindVariableFeatures(seurat_obj_IO, selection.method = "vst", nfeatures = 1000)
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
VariableFeaturePlot(object = seurat_obj_IO, selection.method = "vst", log = TRUE)
```
```{r}
seurat_obj_IO <- ScaleData(seurat_obj_IO, features = rownames(seurat_obj_IO))
```
```{r}
seurat_obj_IO <- RunPCA(seurat_obj_IO)
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
ElbowPlot(seurat_obj_IO, reduction = "pca", ndims = 50)
```
```{r}
nb_pcs = 20
seurat_obj_IO <- FindNeighbors(seurat_obj_IO, reduction = "pca", dims = 1:nb_pcs, prune.SNN = 0)
```
```{r}
cluster_res = 1
DefaultAssay(seurat_obj_IO) <- "RNA_repaired"
seurat_obj_IO <- FindClusters(seurat_obj_IO, resolution = cluster_res)
```
```{r}
table(seurat_obj_IO$seurat_clusters)
```
```{r}
seurat_clusters_IO_color <- list()
n_clusters <- 9
seurat_clusters_IO_color <- color_palette_30[1:n_clusters]
names(seurat_clusters_IO_color) <- seq(1:n_clusters)-1
seurat_obj_IO <- RunUMAP(seurat_obj_IO, reduction = "pca", dims = 1:nb_pcs)
options(repr.plot.width=25, repr.plot.height=8)
DimPlot(seurat_obj_IO, group.by="orig.ident", split.by="orig.ident_merge", cols=orig.ident_colors, shuffle=TRUE, seed=seed, ncol=4, pt.size=1, label=FALSE, raster=FALSE)
```
```{r}
options(repr.plot.width=11, repr.plot.height=9)
for(timepoint in names(orig.ident_merge_colors)){
highlight_cells <- rownames([email hidden][[email hidden]$orig.ident_merge==timepoint,])
print(DimPlot(seurat_obj_IO, cells.highlight = highlight_cells, pt.size=0.5, cols.highlight = orig.ident_merge_colors[timepoint], raster=FALSE))
}
```
```{r}
options(repr.plot.width=11, repr.plot.height=9)
DimPlot(seurat_obj_IO, group.by="seurat_clusters", cols=seurat_clusters_IO_color, shuffle=TRUE, seed=seed, pt.size=1, label=FALSE, raster=FALSE)
```
```{r}
options(repr.plot.width=25, repr.plot.height=8)
DimPlot(seurat_obj_IO, group.by="seurat_clusters", split.by="orig.ident_merge", cols=seurat_clusters_IO_color, shuffle=TRUE, seed=seed, ncol=4, pt.size=1, label=FALSE, raster=FALSE)
```
```{r}
DefaultAssay(seurat_obj_IO) <- "RNA_repaired"
options(repr.plot.width=25, repr.plot.height=60)
FeaturePlot(seurat_obj_IO, feature=c("Calb1"), pt.size=0.5, label=FALSE, raster=FALSE)
```
```{r}
options(repr.plot.width=25, repr.plot.height=60)
p <- FeaturePlot(seurat_obj_IO, feature = "Calb1", pt.size = 1, label = FALSE, raster = FALSE, order = TRUE) +
NoAxes()
# Change the color using a viridis palette
p1 <- p + viridis::scale_color_viridis(option = "magma", direction = -1)
print(p1)
ggsave("IO_calb1Ex_umap_20250915.svg", p1)
```
```{r}
DefaultAssay(seurat_obj_IO) <- "RNA_repaired"
options(repr.plot.width=25, repr.plot.height=60)
FeaturePlot(seurat_obj_IO, feature=c("Calb1", "Aif1"), split.by="orig.ident_merge", ncol=4, pt.size=0.1, label=FALSE, raster=FALSE)
```
```{r}
options(repr.plot.width=15, repr.plot.height=15)
set.seed(115)
mat = as.matrix(table(seurat_obj_IO$orig.ident_merge, seurat_obj_IO$seurat_clusters))
for (row in 1:nrow(mat)) mat[row,] <- as.integer(mat[row,]/(sum(mat[row,])/min(rowSums(mat))))
circos.clear()
par(cex = 0.8)
chordDiagram(mat, annotationTrack = c("grid", "axis"), preAllocateTracks = list(track.height = max(strwidth(unlist(dimnames(mat))))/3), big.gap = 20, small.gap = 2, order = c(colnames(mat), names(orig.ident_merge_colors)), grid.col=c(rev(seurat_clusters_IO_color), rev(orig.ident_merge_colors)))
circos.track(track.index = 1, panel.fun = function(x, y) {circos.text(CELL_META$xcenter, CELL_META$ylim[1], CELL_META$sector.index, facing = "clockwise", niceFacing = TRUE, adj = c(-0.3, 0.5))}, bg.border = NA)
```
```{r}
pdf("Chord_Plot_IO_20250910.pdf", width=8, height=8)
options(repr.plot.width=15, repr.plot.height=15)
set.seed(115)
mat = as.matrix(table(seurat_obj_IO$orig.ident_merge, seurat_obj_IO$seurat_clusters))
for (row in 1:nrow(mat)) mat[row,] <- as.integer(mat[row,]/(sum(mat[row,])/min(rowSums(mat))))
circos.clear()
par(cex = 0.8)
chordDiagram(mat, annotationTrack = c("grid", "axis"), preAllocateTracks = list(track.height = max(strwidth(unlist(dimnames(mat))))/3), big.gap = 20, small.gap = 2, order = c(colnames(mat), names(orig.ident_merge_colors)), grid.col=c(rev(seurat_clusters_IO_color), rev(orig.ident_merge_colors)))
circos.track(track.index = 1, panel.fun = function(x, y) {circos.text(CELL_META$xcenter, CELL_META$ylim[1], CELL_META$sector.index, facing = "clockwise", niceFacing = TRUE, adj = c(-0.3, 0.5))}, bg.border = NA)
dev.off()
```
```{r}
t <- table(Cluster=seurat_obj_IO$seurat_clusters, Batch=[email hidden][["orig.ident_merge"]])
t <- t[,rev(names(orig.ident_merge_colors))]
options(repr.plot.width=10, repr.plot.height=10)
barplot(t(t), xlab = "Cluster", ylab="Number of cells", legend = TRUE, ylim = c(0, round_any(as.integer(max(rowSums(t)))+500, 20, f = ceiling)), col = rev(orig.ident_merge_colors), args.legend = list(bty = "n", x = "top", ncol = 3))
```
```{r}
t <- table(Cluster=seurat_obj_IO$seurat_clusters, Batch=[email hidden][["orig.ident"]])
t <- t[,rev(names(orig.ident_colors))]
options(repr.plot.width=10, repr.plot.height=10)
barplot(t(t), xlab = "Cluster", ylab="Number of cells", legend = TRUE, ylim = c(0, round_any(as.integer(max(rowSums(t)))+1000, 20, f = ceiling)), col = rev(orig.ident_colors), args.legend = list(bty = "n", x = "top", ncol = 3))
```
```{r}
DefaultAssay(seurat_obj_IO) <- "RNA"
Idents(object = seurat_obj_IO) <- "seurat_clusters"
```
```{r}
seurat_obj_IO.markers.seurat_clusters <- FindAllMarkers(seurat_obj_IO, slot = "data", min.pct = 0.1, logfc.threshold = 0.5, only.pos = TRUE)
```
```{r}
seurat_obj_IO.markers.seurat_clusters_pval <- seurat_obj_IO.markers.seurat_clusters[seurat_obj_IO.markers.seurat_clusters$p_val_adj < 0.05,]
seurat_obj_IO.markers.seurat_clusters_best <- seurat_obj_IO.markers.seurat_clusters_pval %>%
filter(avg_log2FC > 0.5) %>%
filter(pct.1 > 0.3) %>%
arrange(cluster, desc(avg_log2FC)) %>%
group_by(cluster)
markers_heatmap <- seurat_obj_IO.markers.seurat_clusters_best %>% top_n(n = 10, wt = avg_log2FC)
topDiffCluster <- seurat_obj_IO.markers.seurat_clusters_best %>% top_n(n = 50, wt = avg_log2FC)
sample_df <- as.data.frame(lapply(split(topDiffCluster, topDiffCluster$cluster), function(x) c(x$gene,rep("None", 50-length(x$gene)))))
colnames(sample_df) <- levels(topDiffCluster$cluster)
sample_df
```
```{r}
# export the cluster marker for all pixels
write.csv(sample_df, "IO_allPixel_clusterMarkers_20250915.csv", row.names = TRUE)
```
```{r}
options(repr.plot.width=30, repr.plot.height=20)
DoMultiBarHeatmap(
subset(seurat_obj_IO, downsample = 200),
features = unique(markers_heatmap$gene),
cells = NULL,
group.by = "seurat_clusters",
additional.group.by = c("orig.ident_merge"),
additional.group.sort.by = c("orig.ident_merge"),
cols.use = list(seurat_clusters=seurat_clusters_IO_color, orig.ident_merge=orig.ident_merge_colors),
group.bar = TRUE,
disp.min = -2.5,
disp.max = NULL,
layer = "scale.data",
assay = "RNA_repaired",
label = TRUE,
size = 5.5,
hjust = 0,
angle = 45,
raster = TRUE,
draw.lines = TRUE,
lines.width = NULL,
group.bar.height = 0.02,
combine = TRUE
) + theme(text = element_text(size = 8))
```
```{r}
for (sample_name in names(orig.ident_colors)){
annot_plot <- [email hidden][seurat_obj_IO$orig.ident==sample_name,]
options(repr.plot.width=5.3, repr.plot.height=4.6)
df_tmp <- data.frame(
x = annot_plot$x,
y = annot_plot$y,
value = annot_plot$seurat_clusters
)
print(ggplot(df_tmp, aes(x, y)) +
geom_point(aes(colour = value),shape=15, size = 1)+
scale_color_manual(values=seurat_clusters_IO_color) +
scale_x_continuous(expand=c(-0.005,0), lim=c(0,101)) +
scale_y_continuous(expand=c(-0.006,0), lim=c(0,101)) +
theme_void() +
ggtitle(sample_name))
}
```
```{r}
# load the given .rds file instead of using subset
# seurat_obj_IO_calb1 <- subset(seurat_obj_IO, subset = Calb1 > 0, slot = "counts")
```
# 7. Calbindin neurons in the IO
```{r}
# load the given .rds file instead of running
seurat_obj_IO_calb1 <- readRDS("~/rds/rds-rds-karadottir-qqstKNzWvy8/Omar_Spatial/seurat_obj_IO_calb1_DIET.rds")
```
```{r}
options(repr.plot.width=11, repr.plot.height=9)
DimPlot(seurat_obj_IO, cells.highlight = colnames(seurat_obj_IO_calb1), pt.size=0.5,cols.highlight = c("red", "grey50"), raster=FALSE)
```
```{r}
percent_micro <- (as.vector(table(seurat_obj_IO_calb1$orig.ident))[1:12]*100)/as.vector(table(seurat_obj_IO$orig.ident))[1:12]
mat <- tibble(
group = c("Ctrl","Ctrl","Ctrl","Ctrl","D7","D7","D7","D7","D14","D14","D14","D14"),
measurement = percent_micro
)
mat$group <- factor(mat$group, levels = c("Ctrl", "D7", "D14"))
options(repr.plot.width=4, repr.plot.height=6)
barplot_with_error_and_dots(mat, x = group, y = measurement, error = "se", title="Calb1+ pixels (% of IO)")
```
```{r}
DefaultAssay(seurat_obj_IO_calb1) <- "RNA"
seurat_obj_IO_calb1 <- NormalizeData(seurat_obj_IO_calb1, normalization.method = "LogNormalize", scale.factor = 10000)
DefaultAssay(seurat_obj_IO_calb1) <- "RNA_repaired"
seurat_obj_IO_calb1 <- NormalizeData(seurat_obj_IO_calb1, normalization.method = "LogNormalize", scale.factor = 10000)
seurat_obj_IO_calb1 <- FindVariableFeatures(seurat_obj_IO_calb1, selection.method = "vst", nfeatures = 1000)
options(repr.plot.width=10, repr.plot.height=7)
VariableFeaturePlot(object = seurat_obj_IO_calb1, selection.method = "vst", log = TRUE)
```
```{r}
seurat_obj_IO_calb1 <- ScaleData(seurat_obj_IO_calb1, features = rownames(seurat_obj_IO_calb1))
```
```{r}
seurat_obj_IO_calb1 <- RunPCA(seurat_obj_IO_calb1)
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
ElbowPlot(seurat_obj_IO_calb1, reduction = "pca", ndims = 50)
```
```{r}
nb_pcs = 20
seurat_obj_IO_calb1 <- FindNeighbors(seurat_obj_IO_calb1, reduction = "pca", dims = 1:nb_pcs, prune.SNN = 0)
#
cluster_res = 0.8
seurat_obj_IO_calb1 <- FindClusters(seurat_obj_IO_calb1, resolution = cluster_res)
```
```{r}
table(seurat_obj_IO_calb1$seurat_clusters)
```
```{r}
seurat_clusters_IO_calb1_color <- list()
n_clusters <- 5
seurat_clusters_IO_calb1_color <- color_palette_30[1:n_clusters]
names(seurat_clusters_IO_calb1_color) <- seq(1:n_clusters)-1
seurat_obj_IO_calb1 <- RunUMAP(seurat_obj_IO_calb1, reduction = "pca", dims = 1:nb_pcs)
```
```{r}
options(repr.plot.width=25, repr.plot.height=8)
DimPlot(seurat_obj_IO_calb1, group.by="orig.ident", split.by="orig.ident_merge", cols=orig.ident_colors, shuffle=TRUE, seed=seed, ncol=4, pt.size=1, label=FALSE, raster=FALSE)
```
```{r}
options(repr.plot.width=11, repr.plot.height=9)
for(timepoint in names(orig.ident_merge_colors)){
highlight_cells <- rownames([email hidden][[email hidden]$orig.ident_merge==timepoint,])
print(DimPlot(seurat_obj_IO_calb1, cells.highlight = highlight_cells, pt.size=0.5, cols.highlight = orig.ident_merge_colors[timepoint], raster=FALSE))
}
```
```{r}
options(repr.plot.width=11, repr.plot.height=9)
DimPlot(seurat_obj_IO_calb1, group.by="seurat_clusters", cols=seurat_clusters_IO_calb1_color, shuffle=TRUE, seed=seed, pt.size=1, label=FALSE, raster=FALSE) + NoAxes()
```
```{r}
# export the umap in .svg file
svg(file="IO_calb1_neuron_20250910.svg", width=4, height=4)
DimPlot(seurat_obj_IO_calb1, group.by="seurat_clusters", cols=seurat_clusters_IO_calb1_color, shuffle=TRUE, seed=seed, pt.size=1, label=FALSE, raster=FALSE)+ NoAxes()
dev.off()
```
```{r}
options(repr.plot.width=11, repr.plot.height=9)
DimPlot(seurat_obj_IO_calb1, group.by="orig.ident_merge", cols=orig.ident_merge_colors, shuffle=TRUE, seed=seed, pt.size=1, label=FALSE, raster=FALSE)
```
```{r}
# export the umap in .svg file
svg(file="IO_calb1_neuron_days_20250903.svg", width=11, height=9)
DimPlot(seurat_obj_IO_calb1, group.by="orig.ident_merge", cols=orig.ident_merge_colors, shuffle=TRUE, seed=seed, pt.size=1, label=FALSE, raster=FALSE)
dev.off()
```
```{r}
options(repr.plot.width=25, repr.plot.height=8)
DimPlot(seurat_obj_IO_calb1, group.by="seurat_clusters", split.by="orig.ident_merge", cols=seurat_clusters_IO_calb1_color, shuffle=TRUE, seed=seed, ncol=4, pt.size=1, label=FALSE, raster=FALSE)
```
```{r}
# Extract metadata
df <- [email hidden] %>%
dplyr::select(seurat_clusters, orig.ident_merge)
# Count cells per cluster and per day
df_summary <- df %>%
group_by(orig.ident_merge, seurat_clusters) %>%
summarise(n = n(), .groups = "drop") %>%
group_by(orig.ident_merge) %>%
mutate(freq = n / sum(n)) # normalize to fractions
# Plot stacked barplot
ggplot(df_summary, aes(x = orig.ident_merge, y = freq, fill = seurat_clusters)) +
geom_bar(stat = "identity", position = "fill") +
ylab("Fraction of Cells") +
xlab("Day") +
scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggplot(df_summary, aes(x = orig.ident_merge, y = freq, fill = seurat_clusters)) +
geom_bar(stat = "identity", position = "fill") +
geom_text(aes(label = scales::percent(freq, accuracy = 1)),
position = position_stack(vjust = 0.5), size = 3) +
ylab("Fraction of Cells") +
xlab("Day") +
scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))+
scale_fill_manual(values = seurat_clusters_IO_calb1_color)
```
```{r}
# export the bar plot
barPlot <- ggplot(df_summary, aes(x = orig.ident_merge, y = freq, fill = seurat_clusters)) +
geom_bar(stat = "identity", position = "fill") +
geom_text(aes(label = scales::percent(freq, accuracy = 1)),
position = position_stack(vjust = 0.5), size = 3) +
ylab("Fraction of Cells") +
xlab("Day") +
scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))+
scale_fill_manual(values = seurat_clusters_IO_calb1_color)
ggsave(filename = "IO_calb1Neuron_Cluster_barplot_20250903.svg", plot = barPlot, width = 10, height = 8, units = "in")
```
```{r}
options(repr.plot.width=15, repr.plot.height=15)
set.seed(115)
mat = as.matrix(table(seurat_obj_IO_calb1$orig.ident_merge, seurat_obj_IO_calb1$seurat_clusters))
for (row in 1:nrow(mat)) mat[row,] <- as.integer(mat[row,]/(sum(mat[row,])/min(rowSums(mat))))
circos.clear()
par(cex = 0.8)
chordDiagram(mat, annotationTrack = c("grid", "axis"), preAllocateTracks = list(track.height = max(strwidth(unlist(dimnames(mat))))/3), big.gap = 20, small.gap = 2, order = c(colnames(mat), names(orig.ident_merge_colors)), grid.col=c(rev(seurat_clusters_IO_calb1_color), rev(orig.ident_merge_colors)))
circos.track(track.index = 1, panel.fun = function(x, y) {circos.text(CELL_META$xcenter, CELL_META$ylim[1], CELL_META$sector.index, facing = "clockwise", niceFacing = TRUE, adj = c(-0.3, 0.5))}, bg.border = NA)
```
```{r}
pdf("Chord_plot_IO_NeuronClus.pdf", width=8, height=8)
options(repr.plot.width=15, repr.plot.height=15)
set.seed(115)
mat = as.matrix(table(seurat_obj_IO_calb1$orig.ident_merge, seurat_obj_IO_calb1$seurat_clusters))
for (row in 1:nrow(mat)) mat[row,] <- as.integer(mat[row,]/(sum(mat[row,])/min(rowSums(mat))))
circos.clear()
par(cex = 0.8)
chordDiagram(mat, annotationTrack = c("grid", "axis"), preAllocateTracks = list(track.height = max(strwidth(unlist(dimnames(mat))))/3), big.gap = 20, small.gap = 2, order = c(colnames(mat), names(orig.ident_merge_colors)), grid.col=c(rev(seurat_clusters_IO_calb1_color), rev(orig.ident_merge_colors)))
circos.track(track.index = 1, panel.fun = function(x, y) {circos.text(CELL_META$xcenter, CELL_META$ylim[1], CELL_META$sector.index, facing = "clockwise", niceFacing = TRUE, adj = c(-0.3, 0.5))}, bg.border = NA)
dev.off()
```
```{r}
df_table <- as.data.frame(table(seurat_obj_IO_calb1$orig.ident, seurat_obj_IO_calb1$seurat_clusters))
colnames(df_table) <- c("group", "cluster", "count")
t <- table(seurat_obj_IO_calb1$orig.ident, seurat_obj_IO_calb1$seurat_clusters)
t_prop <- t/as.vector(table(seurat_obj_IO$orig.ident))*100
rownames(t_prop) <- c("Ctrl","Ctrl","Ctrl","Ctrl","D7","D7","D7","D7","D14","D14","D14","D14")
df_table <- as.data.frame(t_prop)
colnames(df_table) <- c("group", "cluster", "count")
df_summary <- df_table %>%
group_by(group, cluster) %>%
summarise(
mean = mean(count),
se = sd(count) / sqrt(n()),
.groups = "drop"
)
options(repr.plot.width=10, repr.plot.height=10)
ggplot(df_summary, aes(x = cluster, y = mean, fill = group)) +
# Bars
geom_bar(stat = "identity", position = position_dodge(width = 0.9), color = "black") +
# Error bars
geom_errorbar(aes(ymin = mean - se, ymax = mean + se),
position = position_dodge(width = 0.9),
width = 0.2) +
# Manual fill colors
scale_fill_manual(values = orig.ident_merge_colors) +
# Points dodged and jittered
geom_point(
data = df_table,
aes(x = cluster, y = count, fill = group), # <- fill included so dodge works
position = position_jitterdodge(jitter.width = 0.4, dodge.width = 0.9),
size = 2, shape = 21, stroke = 0.5,
color = "black", # border
show.legend = FALSE
) +
theme_minimal() +
labs(x = "", y = "% of pixels in the IO") +
ggtitle("Number of calbindin neurons pixels divided by each IO for each sample") +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank())
```
```{r}
t <- table(Cluster=seurat_obj_IO_calb1$seurat_clusters, Batch=[email hidden][["orig.ident_merge"]])
t <- t[,rev(names(orig.ident_merge_colors))]
options(repr.plot.width=10, repr.plot.height=10)
barplot(t(t), xlab = "Cluster", ylab="Number of cells", legend = TRUE, ylim = c(0, round_any(as.integer(max(rowSums(t)))+100, 20, f = ceiling)), col = rev(orig.ident_merge_colors), args.legend = list(bty = "n", x = "top", ncol = 3))
```
```{r}
t <- table(Cluster=seurat_obj_IO_calb1$seurat_clusters, Batch=[email hidden][["orig.ident"]])
t <- t[,rev(names(orig.ident_colors))]
options(repr.plot.width=10, repr.plot.height=10)
barplot(t(t), xlab = "Cluster", ylab="Number of cells", legend = TRUE, ylim = c(0, round_any(as.integer(max(rowSums(t)))+150, 20, f = ceiling)), col = rev(orig.ident_colors), args.legend = list(bty = "n", x = "top", ncol = 3))
```
```{r}
DefaultAssay(seurat_obj_IO_calb1) <- "RNA"
Idents(object = seurat_obj_IO_calb1) <- "seurat_clusters"
seurat_obj_IO_calb1.markers.seurat_clusters <- FindAllMarkers(seurat_obj_IO_calb1, slot = "data", min.pct = 0.1, logfc.threshold = 0, only.pos = TRUE)
```
```{r}
write.csv( seurat_obj_IO_calb1.markers.seurat_clusters, "IO_calb1_neuron_clusterMarkers_20250916.csv")
```
```{r}
seurat_obj_IO_calb1.markers.seurat_clusters_pval <- seurat_obj_IO_calb1.markers.seurat_clusters[seurat_obj_IO_calb1.markers.seurat_clusters$p_val_adj < 0.05,]
seurat_obj_IO_calb1.markers.seurat_clusters_best <- seurat_obj_IO_calb1.markers.seurat_clusters_pval %>%
filter(avg_log2FC > 0.5) %>%
filter(pct.1 > 0.3) %>%
arrange(cluster, desc(avg_log2FC)) %>%
group_by(cluster)
markers_heatmap <- seurat_obj_IO_calb1.markers.seurat_clusters_best %>% top_n(n = 10, wt = avg_log2FC)
topDiffCluster <- seurat_obj_IO_calb1.markers.seurat_clusters_best %>% top_n(n = 200, wt = avg_log2FC)
sample_df <- as.data.frame(lapply(split(topDiffCluster, topDiffCluster$cluster), function(x) c(x$gene,rep("None", 200-length(x$gene)))))
colnames(sample_df) <- levels(topDiffCluster$cluster)
sample_df
```
```{r}
options(repr.plot.width=30, repr.plot.height=20)
DoMultiBarHeatmap(
subset(seurat_obj_IO_calb1, downsample = 200),
features = unique(markers_heatmap$gene),
cells = NULL,
group.by = "seurat_clusters",
additional.group.by = c("orig.ident_merge"),
additional.group.sort.by = c("orig.ident_merge"),
cols.use = list(seurat_clusters=seurat_clusters_IO_calb1_color, orig.ident_merge=orig.ident_merge_colors),
group.bar = TRUE,
disp.min = -2.5,
disp.max = NULL,
layer = "scale.data",
assay = "RNA_repaired",
label = TRUE,
size = 5.5,
hjust = 0,
angle = 45,
raster = TRUE,
draw.lines = TRUE,
lines.width = NULL,
group.bar.height = 0.02,
combine = TRUE
) + theme(text = element_text(size = 8))
```
```{r}
for (sample_name in names(orig.ident_colors)){
annot_plot <- [email hidden][seurat_obj_IO_calb1$orig.ident==sample_name,]
annot_plot_plus <- [email hidden][(seurat_obj_IO$orig.ident==sample_name & seurat_obj_IO$InfOlive=="core"),]
annot_plot_plus <- annot_plot_plus[!(rownames(annot_plot_plus) %in% rownames(annot_plot)),]
annot_plot_plus$seurat_clusters <- 99
annot_plot <- annot_plot[,c("x","y","seurat_clusters")]
annot_plot_plus <- annot_plot_plus[,c("x","y","seurat_clusters")]
annot_plot$seurat_clusters <- as.character(annot_plot$seurat_clusters)
annot_plot_plus$seurat_clusters <- as.character(annot_plot_plus$seurat_clusters)
annot_plot <- rbind(annot_plot_plus,annot_plot)
annot_plot$seurat_clusters <- factor(x = annot_plot$seurat_clusters, levels = c(names(seurat_clusters_IO_calb1_color),99))
options(repr.plot.width=5.3, repr.plot.height=4.6)
df_tmp <- data.frame(
x = annot_plot$x,
y = annot_plot$y,
value = annot_plot$seurat_clusters
)
tmp_colors <- c(seurat_clusters_IO_calb1_color,"lightgrey")
names(tmp_colors) <- c(names(seurat_clusters_IO_calb1_color),99)
print(ggplot(df_tmp, aes(x, y)) +
geom_point(aes(colour = value),shape=15, size = 1)+
scale_color_manual(values=tmp_colors) +
scale_x_continuous(expand=c(-0.005,0), lim=c(0,101)) +
scale_y_continuous(expand=c(-0.006,0), lim=c(0,101)) +
theme_void() +
ggtitle(sample_name))
}
```
## Run GO on the neuron cluster markers
```{r}
# extract gene lists from sample_df
Clus0_markers <- sample_df[,1]
Clus1_markers <- sample_df[,2]
Clus2_markers <- sample_df[,3]
Clus3_markers <- sample_df[,4]
Clus4_markers <- sample_df[,5]
# dbs set up
dbs <- c("GO_Molecular_Function_2023","GO_Cellular_Component_2023","GO_Biological_Process_2023","KEGG_2019_Mouse") # ,"HDSigDB_Mouse_2021","Tabula_Muris","WikiPathways_2024_Mouse")
#
enrichedR.list <- list()
enrichedR.list[["Clus0_markers"]] <- enrichr(Clus0_markers, dbs)
enrichedR.list[["Clus1_markers"]] <- enrichr(Clus1_markers, dbs)
enrichedR.list[["Clus2_markers"]] <- enrichr(Clus2_markers, dbs)
enrichedR.list[["Clus3_markers"]] <- enrichr(Clus3_markers, dbs)
enrichedR.list[["Clus4_markers"]] <- enrichr(Clus4_markers, dbs)
```
```{r}
library(openxlsx)
# Create an empty list to store the data frames for export
export_list <- list()
# Loop through your enrichedR.list and extract a specific database for each cluster
for (cluster_name in names(enrichedR.list)) {
# For example, extract the "GO_Biological_Process_2023" results
export_list[[cluster_name]] <- enrichedR.list[[cluster_name]]$GO_Biological_Process_2023
}
# Export the prepared list to a single Excel file
write.xlsx(export_list, file = "enrichment_IO_calb1_neuron_GO_KEGG.xlsx")
```
```{r}
pdf("enrichment_plot_IO_NeuronClus_200gene_20250915.pdf", width=12, height=8)
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
p <- plotEnrich(
enr_df,
showTerms = 50,
numChar = 100,
y = "Count",
orderBy = "FDR",
title = paste0(dge," in ", go_db),
xlab = ""
)
print(
p + theme(
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 6),
plot.title = element_text(size = 10),
plot.margin = margin(10, 10, 10, 10)
)
)
}
}
dev.off()
```
## DGE on IO neurons Calb1+
```{r}
condition <- "IO"
HVGs.list <- list()
res.list <- list()
HVGs.list[[condition]] <- list()
res.list[[condition]] <- list()
```
```{r}
DefaultAssay(seurat_obj_IO_calb1) <- 'RNA'
table(seurat_obj_IO_calb1$orig.ident_merge)
```
```{r}
seurat_obj_sce_calb1 <- as.SingleCellExperiment(seurat_obj_IO_calb1, assay="RNA")
seurat_obj_sce_calb1 <- prepSCE(seurat_obj_sce_calb1, kid = "InfOlive", gid = "orig.ident_merge", sid = "orig.ident", drop = FALSE)
kids <- purrr::set_names(levels(seurat_obj_sce_calb1$cluster_id))
# Total number of clusters
nk <- length(kids)
# Named vector of sample names
sids <- purrr::set_names(levels(seurat_obj_sce_calb1$sample_id))
# Total number of samples
ns <- length(sids)
kids
sids
```
```{r}
## Determine the number of cells per sample
table(seurat_obj_sce_calb1$sample_id)
## Turn named vector into a numeric vector of number of cells per sample
n_cells <- as.numeric(table(seurat_obj_sce_calb1$sample_id))
## Determine how to reoder the samples (rows) of the metadata to match the order of sample names in sids vector
m <- match(sids, seurat_obj_sce_calb1$sample_id)
## Create the sample level metadata by combining the reordered metadata with the number of cells corresponding to each sample.
ei <- data.frame(colData(seurat_obj_sce_calb1)[m, ], n_cells, row.names = NULL) %>% dplyr::select(-"cluster_id")
ei
```
```{r}
# Aggregate the counts per sample_id and cluster_id
# Subset metadata to only include the cluster and sample IDs to aggregate across
groups <- colData(seurat_obj_sce_calb1)[, c("cluster_id", "sample_id")]
# Aggregate across cluster-sample groups
pb <- aggregate.Matrix(t(counts(seurat_obj_sce_calb1)), groupings = groups, fun = "sum")
splitf <- sapply(stringr::str_split(rownames(pb), pattern = "_", n = 2), `[`, 1)
pb <- split.data.frame(pb, factor(splitf)) %>% lapply(function(u) magrittr::set_colnames(t(u), stringr::str_extract(rownames(u), "(?<=_)[:alnum:]+")))
class(pb)
# Explore the different components of list
str(pb)
```
```{r}
options(width = 100)
table(seurat_obj_sce_calb1$cluster_id, seurat_obj_sce_calb1$group_id)
```
```{r}
options(width = 100)
table(seurat_obj_sce_calb1$cluster_id, seurat_obj_sce_calb1$seurat_clusters)
```
```{r}
# prep. data.frame for plotting
get_sample_ids <- function(x){pb[[x]] %>% colnames()}
de_samples <- purrr::map(1:length(kids), get_sample_ids) %>% unlist()
samples_list <- purrr::map(1:length(kids), get_sample_ids)
get_cluster_ids <- function(x){rep(names(pb)[x],each = length(samples_list[[x]]))}
de_cluster_ids <- purrr::map(1:length(kids), get_cluster_ids) %>% unlist()
gg_df <- data.frame(cluster_id = de_cluster_ids, sample_id = de_samples)
gg_df <- left_join(gg_df, ei[, c("sample_id", "group_id")])
metadata <- gg_df %>% dplyr::select(cluster_id, sample_id, group_id)
levels(metadata$cluster_id) <- unique(metadata$cluster_id)
# Generate vector of cluster IDs
clusters <- levels(metadata$cluster_id)
clusters
cluster_choice <- "core"
table_sample_low <- table(seurat_obj_sce_calb1$cluster_id, seurat_obj_sce_calb1$sample_id)[cluster_choice,][table(seurat_obj_sce_calb1$cluster_id, seurat_obj_sce_calb1$sample_id)[cluster_choice,] <= 5]
table_sample_low
sample_keep <- names(orig.ident_colors)[!names(orig.ident_colors) %in% names(table_sample_low)]
table_group_low <-table(seurat_obj_sce_calb1$cluster_id, seurat_obj_sce_calb1$group_id)[cluster_choice,][table(seurat_obj_sce_calb1$cluster_id, seurat_obj_sce_calb1$group_id)[cluster_choice,] <= 30]
table_group_low
group_keep <- names(orig.ident_merge_colors)[!names(orig.ident_merge_colors) %in% names(table_group_low)]
cluster_metadata <- metadata[which(metadata$cluster_id == cluster_choice), ]
# Assign the rownames of the metadata to be the sample IDs
rownames(cluster_metadata) <- cluster_metadata$sample_id
#Remove sample and timepoint with low counts (sample<=5, timepoint<=30)
cluster_metadata <- cluster_metadata[(cluster_metadata$sample_id %in% sample_keep & cluster_metadata$group_id %in% group_keep),]
new_group_levels <- levels(cluster_metadata$group_id)[levels(cluster_metadata$group_id) %in% group_keep]
cluster_metadata$group_id <- droplevels(cluster_metadata$group_id)
levels(cluster_metadata$group_id) <- new_group_levels
# Subset the counts to only this cluster
counts <- pb[[cluster_choice]]
cluster_counts <- data.frame(counts[, which(colnames(counts) %in% rownames(cluster_metadata))])
# Check that all of the row names of the metadata are the same and in the same order as the column names of the counts in order to use as input to DESeq2
all(rownames(cluster_metadata) == colnames(cluster_counts))
```
```{r}
dds <- DESeqDataSetFromMatrix(cluster_counts, colData = cluster_metadata, design = ~ group_id)
```
```{r}
# Transform counts for data visualization
rld <- rlog(dds, blind=TRUE)
```
```{r}
# Plot PCA
options(repr.plot.width=8, repr.plot.height=5)
DESeq2::plotPCA(rld, intgroup = "group_id")
```
```{r}
# Extract the rlog matrix from the object and compute pairwise correlation values
rld_mat <- assay(rld)
rld_cor <- cor(rld_mat)
```
```{r}
# Plot heatmap
options(repr.plot.width=10, repr.plot.height=8)
pheatmap(rld_cor, annotation = cluster_metadata[, "group_id", drop=F], annotation_colors = list(group_id=orig.ident_merge_colors))
```
```{r}
dds <- DESeq(dds)
```
```{r}
libsize = data.frame(x=sizeFactors(dds), y=colSums(assay(dds)))
options(repr.plot.width=7, repr.plot.height=7)
ggplot(data=libsize, aes(x=x, y=y)) + geom_point() + geom_smooth(method="lm") + xlab("Estimated size factor") + ylab("Library size")
```
```{r}
data.frame(colData(dds))
```
```{r}
# Plot dispersion estimates
options(repr.plot.width=12, repr.plot.height=12)
plotDispEsts(dds)
```
```{r}
HVGs.list[[condition]][["upregulated"]] <- list()
HVGs.list[[condition]][["downregulated"]] <- list()
```
## 9.1 Contrasts
```{r}
condition1 = "D7"
condition2 = "Ctrl"
contrast <- c("group_id", condition1, condition2)
# resultsNames(dds)
res <- results(dds, contrast = contrast, alpha = 0.05)
summary(res)
options(repr.plot.width=10, repr.plot.height=7)
hist(res$pvalue, breaks=0:20/20, col="grey50", border="white", xlim=c(0,1), main="Histogram of p-values", xlab="p-value")
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
hist(res$padj, breaks=0:20/20, col="grey50", border="white", xlim=c(0,1), main="Histogram of adj p-values", xlab="adj p-value")
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
qs <- c(0, quantile(res$baseMean[res$baseMean > 0], 0:7/7))
# cut the genes into the bins
bins <- cut(res$baseMean, qs)
# rename the levels of the bins using the middle point
levels(bins) <- paste0("~",round(.5*qs[-1] + .5*qs[-length(qs)]))
# calculate the ratio of $p$ values less than .01 for each bin
ratios <- tapply(res$pvalue, bins, function(p) mean(p < .01, na.rm=TRUE))
# plot these ratios
barplot(ratios, xlab="mean normalized count", ylab="ratio of small p values")
```
```{r}
options(repr.plot.width=15, repr.plot.height=15)
#padj<0.001
p1<-EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'padj',
title = paste0(condition1,' vs ',condition2,' in FC > 1'),
pCutoff = 0.01,
FCcutoff = 1,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1
)
print(p1)
p1_1<-EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'padj',
title = paste0(condition1,' vs ',condition2,' in FC > 0.75'),
pCutoff = 0.01,
FCcutoff = 0.75,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1
)
print(p1_1)
```
```{r}
res.list[[condition]][[paste0(condition1,"_",condition2)]] <- res
```
```{r}
write.csv(res, paste0("DGE_",condition,"_neuron_",condition1,"_",condition2,"_20250915.csv"), row.names=TRUE)
```
```{r}
res <- res[complete.cases(res$padj),]
res_filtered <- res[(res$padj < 0.05 & res$baseMean > 1),]
HVGs.list[[condition]][["upregulated"]][[paste0(condition1,"_",condition2)]] <- rownames(res_filtered[res_filtered$log2FoldChange>=0,])
HVGs.list[[condition]][["downregulated"]][[paste0(condition1,"_",condition2)]] <- rownames(res_filtered[res_filtered$log2FoldChange<0,])
```
## D14 vs Ctrl
```{r}
condition1 = "D14"
condition2 = "Ctrl"
```
```{r}
contrast <- c("group_id", condition1, condition2)
# resultsNames(dds)
res <- results(dds, contrast = contrast, alpha = 0.05)
```
```{r}
summary(res)
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
hist(res$pvalue, breaks=0:20/20, col="grey50", border="white", xlim=c(0,1), main="Histogram of p-values", xlab="p-value")
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
hist(res$padj, breaks=0:20/20, col="grey50", border="white", xlim=c(0,1), main="Histogram of adj p-values", xlab="adj p-value")
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
qs <- c(0, quantile(res$baseMean[res$baseMean > 0], 0:7/7))
# cut the genes into the bins
bins <- cut(res$baseMean, qs)
# rename the levels of the bins using the middle point
levels(bins) <- paste0("~",round(.5*qs[-1] + .5*qs[-length(qs)]))
# calculate the ratio of $p$ values less than .01 for each bin
ratios <- tapply(res$pvalue, bins, function(p) mean(p < .01, na.rm=TRUE))
# plot these ratios
barplot(ratios, xlab="mean normalized count", ylab="ratio of small p values")
```
```{r}
options(repr.plot.width=15, repr.plot.height=15)
#padj<0.001
p2 <- EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'padj',
title = paste0(condition1,' vs ',condition2,' in FC > 1'),
pCutoff = 0.01,
FCcutoff = 1,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1)
print(p2)
p2_1 <- EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'padj',
title = paste0(condition1,' vs ',condition2,' in FC > 0.75'),
pCutoff = 0.01,
FCcutoff = 0.75,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1)
print(p2_1)
```
```{r}
res.list[[condition]][[paste0(condition1,"_",condition2)]] <- res
```
```{r}
write.csv(res, paste0("DGE_",condition,"_neuron_",condition1,"_",condition2,"_20250915.csv"), row.names=TRUE)
```
```{r}
res <- res[complete.cases(res$padj),]
res_filtered <- res[(res$padj < 0.05 & res$baseMean > 1),]
HVGs.list[[condition]][["upregulated"]][[paste0(condition1,"_",condition2)]] <- rownames(res_filtered[res_filtered$log2FoldChange>=0,])
HVGs.list[[condition]][["downregulated"]][[paste0(condition1,"_",condition2)]] <- rownames(res_filtered[res_filtered$log2FoldChange<0,])
```
## D14 vs D7
```{r}
condition1 = "D14"
condition2 = "D7"
```
```{r}
contrast <- c("group_id", condition1, condition2)
# resultsNames(dds)
res <- results(dds, contrast = contrast, alpha = 0.05)
```
```{r}
summary(res)
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
hist(res$pvalue, breaks=0:20/20, col="grey50", border="white", xlim=c(0,1), main="Histogram of p-values", xlab="p-value")
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
hist(res$padj, breaks=0:20/20, col="grey50", border="white", xlim=c(0,1), main="Histogram of adj p-values", xlab="adj p-value")
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
qs <- c(0, quantile(res$baseMean[res$baseMean > 0], 0:7/7))
# cut the genes into the bins
bins <- cut(res$baseMean, qs)
# rename the levels of the bins using the middle point
levels(bins) <- paste0("~",round(.5*qs[-1] + .5*qs[-length(qs)]))
# calculate the ratio of $p$ values less than .01 for each bin
ratios <- tapply(res$pvalue, bins, function(p) mean(p < .01, na.rm=TRUE))
# plot these ratios
barplot(ratios, xlab="mean normalized count", ylab="ratio of small p values")
```
```{r}
options(repr.plot.width=15, repr.plot.height=15)
#padj<0.001
p3 <- EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'padj',
title = paste0(condition1,' vs ',condition2,' in FC > 1'),
pCutoff = 0.01,
FCcutoff = 1,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1)
print(p3)
p3_1 <- EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'padj',
title = paste0(condition1,' vs ',condition2,' in FC > 0.75'),
pCutoff = 0.01,
FCcutoff = 0.75,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1,
)
print(p3_1)
p3_2 <- EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'padj',
title = paste0(condition1,' vs ',condition2,' in FC > 0.5'),
pCutoff = 0.05,
FCcutoff = 0.5,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1,
)
print(p3_2)
```
```{r}
res.list[[condition]][[paste0(condition1,"_",condition2)]] <- res
```
```{r}
write.csv(res, paste0("DGE_",condition,"_neuron_",condition1,"_",condition2,"_20250915.csv"), row.names=TRUE)
```
```{r}
res <- res[complete.cases(res$padj),]
res_filtered <- res[(res$padj < 0.05 & res$baseMean > 1),]
HVGs.list[[condition]][["upregulated"]][[paste0(condition1,"_",condition2)]] <- rownames(res_filtered[res_filtered$log2FoldChange>=0,])
HVGs.list[[condition]][["downregulated"]][[paste0(condition1,"_",condition2)]] <- rownames(res_filtered[res_filtered$log2FoldChange<0,])
```
```{r}
# export the svg file
svglite::svglite(file = "volcano_plot_neuron_D7vsCtrl_FC_1.svg", width = 15, height = 15)
print(p1)
dev.off()
svglite::svglite(file = "volcano_plot_neuron_D7vsCtrl_FC_075.svg", width = 15, height = 15)
print(p1_1)
dev.off()
svglite::svglite(file = "volcano_plot_neuron_D14vsCtrl_FC_1.svg", width = 15, height = 15)
print(p2)
dev.off()
svglite::svglite(file = "volcano_plot_neuron_D14vsCtrl_FC_075.svg", width = 15, height = 15)
print(p2_1)
dev.off()
svglite::svglite(file = "volcano_plot_neuron_D14vsD7_FC_1.svg", width = 15, height = 15)
print(p3)
dev.off()
svglite::svglite(file = "volcano_plot_neuron_D14vsD7_FC_075.svg", width = 15, height = 15)
print(p3_1)
dev.off()
```
## GO on the res
```{r}
# load the different DEG list
res_D7_Ctrl <- read.csv(paste0("~/rds/rds-rds-karadottir-qqstKNzWvy8/Omar_Spatial/", "DGE_IO_neuron_","D7","_","Ctrl","_20250915.csv"), row.names = 1)
res_D14_Ctrl <- read.csv(paste0("~/rds/rds-rds-karadottir-qqstKNzWvy8/Omar_Spatial/", "DGE_IO_neuron_","D14","_","Ctrl","_20250915.csv"), row.names = 1)
res_D14_D7 <- read.csv(paste0("~/rds/rds-rds-karadottir-qqstKNzWvy8/Omar_Spatial/", "DGE_IO_neuron_","D14","_","D7","_20250915.csv"), row.names = 1)
# clean up the na
res_D7_Ctrl <- res_D7_Ctrl[!is.na(res_D7_Ctrl$padj),]
res_D14_Ctrl <- res_D14_Ctrl[!is.na(res_D14_Ctrl$padj),]
res_D14_D7 <- res_D14_D7[!is.na(res_D14_D7$padj),]
# Filter the res for giving different lists
res_D7_Ctrl_FC1_up <- rownames(res_D7_Ctrl[(res_D7_Ctrl$padj < 0.05) & (res_D7_Ctrl$log2FoldChange > 1),])
res_D7_Ctrl_FC1_down <- rownames(res_D7_Ctrl[(res_D7_Ctrl$padj < 0.05) & (res_D7_Ctrl$log2FoldChange < -1),])
res_D7_Ctrl_FC075_up <- rownames(res_D7_Ctrl[(res_D7_Ctrl$padj < 0.05) & (res_D7_Ctrl$log2FoldChange > 0.75),])
res_D7_Ctrl_FC075_down <- rownames(res_D7_Ctrl[(res_D7_Ctrl$padj < 0.05) & (res_D7_Ctrl$log2FoldChange < -0.75),])
res_D7_Ctrl_FC05_up <- rownames(res_D7_Ctrl[(res_D7_Ctrl$padj < 0.05) & (res_D7_Ctrl$log2FoldChange > 0.5),])
res_D7_Ctrl_FC05_down <- rownames(res_D7_Ctrl[(res_D7_Ctrl$padj < 0.05) & (res_D7_Ctrl$log2FoldChange < -0.5),])
res_D14_Ctrl_FC1_up <- rownames(res_D14_Ctrl[(res_D14_Ctrl$padj < 0.05) & (res_D14_Ctrl$log2FoldChange > 1),])
res_D14_Ctrl_FC1_down <- rownames(res_D14_Ctrl[(res_D14_Ctrl$padj < 0.05) & (res_D14_Ctrl$log2FoldChange < -1),])
res_D14_Ctrl_FC075_up <- rownames(res_D14_Ctrl[(res_D14_Ctrl$padj < 0.05) & (res_D14_Ctrl$log2FoldChange > 0.75),])
res_D14_Ctrl_FC075_down <- rownames(res_D14_Ctrl[(res_D14_Ctrl$padj < 0.05) & (res_D14_Ctrl$log2FoldChange < -0.75),])
res_D14_Ctrl_FC05_up <- rownames(res_D14_Ctrl[(res_D14_Ctrl$padj < 0.05) & (res_D14_Ctrl$log2FoldChange > 0.5),])
res_D14_Ctrl_FC05_down <- rownames(res_D14_Ctrl[(res_D14_Ctrl$padj < 0.05) & (res_D14_Ctrl$log2FoldChange < -0.5),])
res_D14_D7_FC1_up <- rownames(res_D14_D7[(res_D14_D7$padj < 0.05) & (res_D14_D7$log2FoldChange > 1),])
res_D14_D7_FC1_down <- rownames(res_D14_D7[(res_D14_D7$padj < 0.05) & (res_D14_D7$log2FoldChange < -1),])
res_D14_D7_FC075_up <- rownames(res_D14_D7[(res_D14_D7$padj < 0.05) & (res_D14_D7$log2FoldChange > 0.75),])
res_D14_D7_FC075_down <- rownames(res_D14_D7[(res_D14_D7$padj < 0.05) & (res_D14_D7$log2FoldChange < -0.75),])
res_D14_D7_FC05_up <- rownames(res_D14_D7[(res_D14_D7$padj < 0.05) & (res_D14_D7$log2FoldChange > 0.5),])
res_D14_D7_FC05_down <- rownames(res_D14_D7[(res_D14_D7$padj < 0.05) & (res_D14_D7$log2FoldChange < -0.5),])
```
```{r}
# Load the libraries
library(dplyr)
library(writexl)
library(tibble)
# Filter and rank the D7 vs Ctrl results, then add a column for the index
res_D7_Ctrl_ranked <- res_D7_Ctrl %>%
filter(abs(log2FoldChange) > 0.5) %>%
arrange(padj) %>%
rownames_to_column("GeneName")
# Filter and rank the D14 vs Ctrl results
res_D14_Ctrl_ranked <- res_D14_Ctrl %>%
filter(abs(log2FoldChange) > 0.5) %>%
arrange(padj) %>%
rownames_to_column("GeneName")
# Filter and rank the D14 vs D7 results
res_D14_D7_ranked <- res_D14_D7 %>%
filter(abs(log2FoldChange) > 0.5) %>%
arrange(padj) %>%
rownames_to_column("GeneName")
# Create a named list of the filtered data frames
filtered_results_list <- list(
D7_vs_Ctrl = res_D7_Ctrl_ranked,
D14_vs_Ctrl = res_D14_Ctrl_ranked,
D14_vs_D7 = res_D14_D7_ranked
)
# Export the list to a single Excel file with three sheets
write_xlsx(filtered_results_list, "DGE_IO_neuron_ConditionComparison.xlsx")
```
```{r}
dbs <- c("GO_Molecular_Function_2023","GO_Cellular_Component_2023","GO_Biological_Process_2023","KEGG_2019_Mouse") # ,"HDSigDB_Mouse_2021","Tabula_Muris","WikiPathways_2024_Mouse")
```
```{r}
enrichedR.list <- list()
enrichedR.list[["IO_D7_Ctrl_FC1_up"]] <- enrichr(res_D7_Ctrl_FC1_up, dbs)
enrichedR.list[["IO_D7_Ctrl_FC1_down"]] <- enrichr(res_D7_Ctrl_FC1_down, dbs)
enrichedR.list[["IO_D7_Ctrl_FC075_up"]] <- enrichr(res_D7_Ctrl_FC075_up, dbs)
enrichedR.list[["IO_D7_Ctrl_FC075_down"]] <- enrichr(res_D7_Ctrl_FC075_down, dbs)
enrichedR.list[["IO_D7_Ctrl_FC05_up"]] <- enrichr(res_D7_Ctrl_FC05_up, dbs)
enrichedR.list[["IO_D7_Ctrl_FC05_down"]] <- enrichr(res_D7_Ctrl_FC05_down, dbs)
enrichedR.list[["IO_D14_Ctrl_FC1_up"]] <- enrichr(res_D14_Ctrl_FC1_up, dbs)
enrichedR.list[["IO_D14_Ctrl_FC1_down"]] <- enrichr(res_D14_Ctrl_FC1_down, dbs)
enrichedR.list[["IO_D14_Ctrl_FC075_up"]] <- enrichr(res_D14_Ctrl_FC075_up, dbs)
enrichedR.list[["IO_D14_Ctrl_FC075_down"]] <- enrichr(res_D14_Ctrl_FC075_down, dbs)
enrichedR.list[["IO_D14_Ctrl_FC05_up"]] <- enrichr(res_D14_Ctrl_FC05_up, dbs)
enrichedR.list[["IO_D14_Ctrl_FC05_down"]] <- enrichr(res_D14_Ctrl_FC05_down, dbs)
enrichedR.list[["IO_D14_D7_FC1_up"]] <- enrichr(res_D14_D7_FC1_up, dbs)
enrichedR.list[["IO_D14_D7_FC1_down"]] <- enrichr(res_D14_D7_FC1_down, dbs)
enrichedR.list[["IO_D14_D7_FC075_up"]] <- enrichr(res_D14_D7_FC075_up, dbs)
enrichedR.list[["IO_D14_D7_FC075_down"]] <- enrichr(res_D14_D7_FC075_down, dbs)
enrichedR.list[["IO_D14_D7_FC05_up"]] <- enrichr(res_D14_D7_FC05_up, dbs)
enrichedR.list[["IO_D14_D7_FC05_down"]] <- enrichr(res_D14_D7_FC05_down, dbs)
```
```{r}
# Create an empty list to store the data frames for export
export_list <- list()
# Loop through your enrichedR.list and extract a specific database for each cluster
for (cluster_name in names(enrichedR.list)) {
# For example, extract the "GO_Biological_Process_2023" results
export_list[[cluster_name]] <- enrichedR.list[[cluster_name]]$GO_Biological_Process_2023
}
# Export the prepared list to a single Excel file
write.xlsx(export_list, file = "enrichment_IO_calb1_neuron_DayComparison_GO_KEGG.xlsx")
```
```{r}
pdf("enrichment_plot_IO_Neuron_DayComp.pdf", width=12, height=8)
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
p <- plotEnrich(
enr_df,
showTerms = 50,
numChar = 100,
y = "Count",
orderBy = "FDR",
title = paste0(dge," in ", go_db),
xlab = ""
)
print(
p + theme(
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 6),
plot.title = element_text(size = 10),
plot.margin = margin(10, 10, 10, 10)
)
)
}
}
dev.off()
```
```{r}
pdf("enrichment_plot_IO_Neuron_DayComp_filtered.pdf", width=8, height=4)
# Define your list of keywords
keywords <- c("Glycolytic Process", "Pyruvate Metabolic Process", "Carbohydrate Catabolic Process", "Synaptic Vesicle Recycling", "metabo")
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
# Filter for keywords and FDR < 0.1
filtered_df <- enr_df[
grepl(paste(keywords, collapse = "|"), enr_df$Term, ignore.case = TRUE) & enr_df$Adjusted.P.value < 0.1,
]
# Check if there are any results after filtering
if (nrow(filtered_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no terms found after filtering)")
next
}
# Use the filtered data for plotting
p <- plotEnrich(
filtered_df,
showTerms = 50,
numChar = 100,
y = "Count",
orderBy = "FDR",
title = paste0(dge, " in ", go_db, " (Filtered)"),
xlab = ""
)
print(
p + theme(
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 6),
plot.title = element_text(size = 10),
plot.margin = margin(10, 10, 10, 10)
)
)
}
}
dev.off()
```
----------------------------------------
## Combining day 7 and day 14 to compare against control (IO neuron)
```{r}
condition <- "IO"
HVGs.list <- list()
res.list <- list()
HVGs.list[[condition]] <- list()
res.list[[condition]] <- list()
```
```{r}
DefaultAssay(seurat_obj_IO_calb1) <- 'RNA'
table(seurat_obj_IO_calb1$orig.ident_merge)
```
```{r}
# Create a copy of the Seurat object to work with
seurat_obj_IO_calb1_copy <- seurat_obj_IO_calb1
# Access the orig.ident_merge metadata column and convert it to a character vector.
# This allows for the assignment of new values that are not already present in the factor levels.
seurat_obj_IO_calb1_copy$lesion_status <- as.character(seurat_obj_IO_calb1_copy$orig.ident_merge)
# Now, update the D7 and D14 values to "lesion"
seurat_obj_IO_calb1_copy$lesion_status[seurat_obj_IO_calb1_copy$lesion_status %in% c("D7", "D14")] <- "lesion"
# Finally, convert the new column back to a factor and set the levels.
# This step is crucial for downstream analysis to define your groups.
seurat_obj_IO_calb1_copy$lesion_status <- factor(seurat_obj_IO_calb1_copy$lesion_status, levels = c("Ctrl", "lesion"))
# Verify the changes
table(seurat_obj_IO_calb1_copy$lesion_status)
seurat_obj_sce_calb1 <- as.SingleCellExperiment(seurat_obj_IO_calb1_copy, assay="RNA")
seurat_obj_sce_calb1 <- prepSCE(seurat_obj_sce_calb1, kid = "InfOlive", gid = "lesion_status", sid = "orig.ident", drop = FALSE)
kids <- purrr::set_names(levels(seurat_obj_sce_calb1$cluster_id))
# Total number of clusters
nk <- length(kids)
# Named vector of sample names
sids <- purrr::set_names(levels(seurat_obj_sce_calb1$sample_id))
# Total number of samples
ns <- length(sids)
kids
sids
```
```{r}
## Determine the number of cells per sample
table(seurat_obj_sce_calb1$sample_id)
## Turn named vector into a numeric vector of number of cells per sample
n_cells <- as.numeric(table(seurat_obj_sce_calb1$sample_id))
## Determine how to reoder the samples (rows) of the metadata to match the order of sample names in sids vector
m <- match(sids, seurat_obj_sce_calb1$sample_id)
## Create the sample level metadata by combining the reordered metadata with the number of cells corresponding to each sample.
ei <- data.frame(colData(seurat_obj_sce_calb1)[m, ], n_cells, row.names = NULL) %>% dplyr::select(-"cluster_id")
ei
```
```{r}
# Aggregate the counts per sample_id and cluster_id
# Subset metadata to only include the cluster and sample IDs to aggregate across
groups <- colData(seurat_obj_sce_calb1)[, c("cluster_id", "sample_id")]
# Aggregate across cluster-sample groups
pb <- aggregate.Matrix(t(counts(seurat_obj_sce_calb1)), groupings = groups, fun = "sum")
splitf <- sapply(stringr::str_split(rownames(pb), pattern = "_", n = 2), `[`, 1)
pb <- split.data.frame(pb, factor(splitf)) %>% lapply(function(u) magrittr::set_colnames(t(u), stringr::str_extract(rownames(u), "(?<=_)[:alnum:]+")))
class(pb)
# Explore the different components of list
str(pb)
```
```{r}
options(width = 100)
table(seurat_obj_sce_calb1$cluster_id, seurat_obj_sce_calb1$group_id)
```
```{r}
options(width = 100)
table(seurat_obj_sce_calb1$cluster_id, seurat_obj_sce_calb1$seurat_clusters)
```
```{r}
# prep. data.frame for plotting
get_sample_ids <- function(x){pb[[x]] %>% colnames()}
de_samples <- purrr::map(1:length(kids), get_sample_ids) %>% unlist()
samples_list <- purrr::map(1:length(kids), get_sample_ids)
get_cluster_ids <- function(x){rep(names(pb)[x],each = length(samples_list[[x]]))}
de_cluster_ids <- purrr::map(1:length(kids), get_cluster_ids) %>% unlist()
gg_df <- data.frame(cluster_id = de_cluster_ids, sample_id = de_samples)
gg_df <- left_join(gg_df, ei[, c("sample_id", "group_id")])
metadata <- gg_df %>% dplyr::select(cluster_id, sample_id, group_id)
levels(metadata$cluster_id) <- unique(metadata$cluster_id)
# Generate vector of cluster IDs
clusters <- levels(metadata$cluster_id)
clusters
cluster_choice <- "core"
table_sample_low <- table(seurat_obj_sce_calb1$cluster_id, seurat_obj_sce_calb1$sample_id)[cluster_choice,][table(seurat_obj_sce_calb1$cluster_id, seurat_obj_sce_calb1$sample_id)[cluster_choice,] <= 5]
table_sample_low
sample_keep <- names(orig.ident_colors)[!names(orig.ident_colors) %in% names(table_sample_low)]
table_group_low <-table(seurat_obj_sce_calb1$cluster_id, seurat_obj_sce_calb1$group_id)[cluster_choice,][table(seurat_obj_sce_calb1$cluster_id, seurat_obj_sce_calb1$group_id)[cluster_choice,] <= 30]
table_group_low
group_keep <- names(orig.ident_merge_colors)[!names(orig.ident_merge_colors) %in% names(table_group_low)]
group_keep <- c("Ctrl", "lesion")
cluster_metadata <- metadata[which(metadata$cluster_id == cluster_choice), ]
# Assign the rownames of the metadata to be the sample IDs
rownames(cluster_metadata) <- cluster_metadata$sample_id
#Remove sample and timepoint with low counts (sample<=5, timepoint<=30)
cluster_metadata <- cluster_metadata[(cluster_metadata$sample_id %in% sample_keep & cluster_metadata$group_id %in% group_keep),]
new_group_levels <- levels(cluster_metadata$group_id)[levels(cluster_metadata$group_id) %in% group_keep]
cluster_metadata$group_id <- droplevels(cluster_metadata$group_id)
levels(cluster_metadata$group_id) <- new_group_levels
# Subset the counts to only this cluster
counts <- pb[[cluster_choice]]
cluster_counts <- data.frame(counts[, which(colnames(counts) %in% rownames(cluster_metadata))])
# Check that all of the row names of the metadata are the same and in the same order as the column names of the counts in order to use as input to DESeq2
all(rownames(cluster_metadata) == colnames(cluster_counts))
```
```{r}
dds <- DESeqDataSetFromMatrix(cluster_counts, colData = cluster_metadata, design = ~ group_id)
```
```{r}
# Transform counts for data visualization
rld <- rlog(dds, blind=TRUE)
```
```{r}
# Plot PCA
options(repr.plot.width=8, repr.plot.height=5)
DESeq2::plotPCA(rld, intgroup = "group_id")
```
```{r}
# Extract the rlog matrix from the object and compute pairwise correlation values
rld_mat <- assay(rld)
rld_cor <- cor(rld_mat)
```
```{r}
# Plot heatmap
options(repr.plot.width=10, repr.plot.height=8)
# pheatmap(rld_cor, annotation = cluster_metadata[, "group_id", drop=F], annotation_colors = list(group_id=orig.ident_merge_colors))
```
```{r}
dds <- DESeq(dds)
```
```{r}
libsize = data.frame(x=sizeFactors(dds), y=colSums(assay(dds)))
options(repr.plot.width=7, repr.plot.height=7)
ggplot(data=libsize, aes(x=x, y=y)) + geom_point() + geom_smooth(method="lm") + xlab("Estimated size factor") + ylab("Library size")
```
```{r}
data.frame(colData(dds))
```
```{r}
# Plot dispersion estimates
options(repr.plot.width=12, repr.plot.height=12)
plotDispEsts(dds)
```
```{r}
HVGs.list[[condition]][["upregulated"]] <- list()
HVGs.list[[condition]][["downregulated"]] <- list()
```
```{r}
condition1 = "lesion"
condition2 = "Ctrl"
contrast <- c("group_id", condition1, condition2)
# resultsNames(dds)
res <- results(dds, contrast = contrast, alpha = 0.05)
summary(res)
options(repr.plot.width=10, repr.plot.height=7)
hist(res$pvalue, breaks=0:20/20, col="grey50", border="white", xlim=c(0,1), main="Histogram of p-values", xlab="p-value")
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
hist(res$padj, breaks=0:20/20, col="grey50", border="white", xlim=c(0,1), main="Histogram of adj p-values", xlab="adj p-value")
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
qs <- c(0, quantile(res$baseMean[res$baseMean > 0], 0:7/7))
# cut the genes into the bins
bins <- cut(res$baseMean, qs)
# rename the levels of the bins using the middle point
levels(bins) <- paste0("~",round(.5*qs[-1] + .5*qs[-length(qs)]))
# calculate the ratio of $p$ values less than .01 for each bin
ratios <- tapply(res$pvalue, bins, function(p) mean(p < .01, na.rm=TRUE))
# plot these ratios
barplot(ratios, xlab="mean normalized count", ylab="ratio of small p values")
```
```{r}
options(repr.plot.width=15, repr.plot.height=15)
#padj<0.001
genes_to_highlight = c("Gapdh", "Cyc1", "ND6", "Spp1")
p1_2<-EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'padj',
title = paste0(condition1,' vs ',condition2,' in FC > 0.5'),
pCutoff = 0.05,
FCcutoff = 0.5,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1,
selectLab = genes_to_highlight,
drawConnectors = TRUE
)
print(p1_2)
```
```{r}
res.list[[condition]][[paste0(condition1,"_",condition2)]] <- res
```
```{r}
write.csv(res, paste0("DGE_",condition,"_neuron_",condition1,"_",condition2,"_20250915.csv"), row.names=TRUE)
```
```{r}
res <- res[complete.cases(res$padj),]
res_filtered <- res[(res$padj < 0.05 & res$baseMean > 1),]
HVGs.list[[condition]][["upregulated"]][[paste0(condition1,"_",condition2)]] <- rownames(res_filtered[res_filtered$log2FoldChange>=0,])
HVGs.list[[condition]][["downregulated"]][[paste0(condition1,"_",condition2)]] <- rownames(res_filtered[res_filtered$log2FoldChange<0,])
```
```{r}
# export the svg file
svglite::svglite(file = "volcano_plot_neuron_LesionvsCtrl_FC_05.svg", width = 8, height = 8)
print(p1_2)
dev.off()
```
### GO on the res
```{r}
# load the different DEG list
res_lesion_Ctrl <- read.csv(paste0("~/rds/rds-rds-karadottir-qqstKNzWvy8/Omar_Spatial/", "DGE_IO_neuron_lesion_Ctrl_20250915.csv"), row.names = 1)
# clean up the na
res_lesion_Ctrl <- res_lesion_Ctrl[!is.na(res_lesion_Ctrl$padj),]
# Filter the res for giving different lists
res_lesion_Ctrl_FC05_up <- rownames(res_lesion_Ctrl[(res_lesion_Ctrl$padj < 0.05) & (res_lesion_Ctrl$log2FoldChange > 0.5),])
res_lesion_Ctrl_FC05_down <- rownames(res_lesion_Ctrl[(res_lesion_Ctrl$padj < 0.05) & (res_lesion_Ctrl$log2FoldChange < -0.5),])
```
```{r}
# 1. Filter the data based on your criteria
filtered_data <- res_lesion_Ctrl[(res_lesion_Ctrl$padj < 0.05) & (res_lesion_Ctrl$log2FoldChange > 0.5),]
# 2. Rank the filtered data by 'padj' in ascending order
ranked_data <- filtered_data[order(filtered_data$padj),]
write.csv(ranked_data, "DGE_IO_neuron_lesion_ctrl_FC05_ranked.csv")
# 1. Filter the data based on your criteria
filtered_data <- res_lesion_Ctrl[(res_lesion_Ctrl$padj < 0.05) & (res_lesion_Ctrl$log2FoldChange < -0.5),]
# 2. Rank the filtered data by 'padj' in ascending order
ranked_data <- filtered_data[order(filtered_data$padj),]
write.csv(ranked_data, "DGE_IO_neuron_lesion_ctrl_FC05_ranked_down.csv")
```
```{r}
dbs <- c("GO_Molecular_Function_2023","GO_Cellular_Component_2023","GO_Biological_Process_2023","KEGG_2019_Mouse")
```
```{r}
enrichedR.list <- list()
enrichedR.list[["IO_Lesion_Ctrl_FC05_up"]] <- enrichr(res_lesion_Ctrl_FC05_up, dbs)
enrichedR.list[["IO_Lesion_Ctrl_FC05_down"]] <- enrichr(res_lesion_Ctrl_FC05_down, dbs)
```
```{r}
# Create an empty list to store the data frames for export
export_list <- list()
# Loop through your enrichedR.list and extract a specific database for each cluster
for (cluster_name in names(enrichedR.list)) {
# For example, extract the "KEGG" results
export_list[[cluster_name]] <- enrichedR.list[[cluster_name]]$KEGG_2019_Mouse
}
# Export the prepared list to a single Excel file
openxlsx::write.xlsx(export_list, file = "enrichment_IO_calb1_neuron_LesionVsCtrl_KEGG_20250924.xlsx")
```
```{r}
pdf("enrichment_plot_IO_Neuron_MergedLesionComp.pdf", width=12, height=8)
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
p <- plotEnrich(
enr_df,
showTerms = 50,
numChar = 100,
y = "Count",
orderBy = "FDR",
title = paste0(dge," in ", go_db),
xlab = ""
)
print(
p + theme(
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 6),
plot.title = element_text(size = 10),
plot.margin = margin(10, 10, 10, 10)
)
)
}
}
dev.off()
```
```{r}
pdf("enrichment_plot_IO_Neuron_MergedLesionComp_filtered.pdf", width=8, height=4)
# Define your list of keywords
keywords <- c("Oxidative phosphorylation", "Parkinson disease", "Retrograde endocannabinoid signaling", "Alzheimer disease")
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
# Filter for keywords and FDR < 0.1
filtered_df <- enr_df[
grepl(paste(keywords, collapse = "|"), enr_df$Term, ignore.case = TRUE) & enr_df$Adjusted.P.value < 0.1,
]
# Check if there are any results after filtering
if (nrow(filtered_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no terms found after filtering)")
next
}
# Use the filtered data for plotting
p <- plotEnrich(
filtered_df,
showTerms = 50,
numChar = 100,
y = "Count",
orderBy = "FDR",
title = paste0(dge, " in ", go_db, " (Filtered)"),
xlab = ""
)
print(
p + theme(
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 6),
plot.title = element_text(size = 10),
plot.margin = margin(10, 10, 10, 10)
)
)
}
}
dev.off()
```
## Comparison on Clusters
```{r}
library(Matrix)
library(dplyr)
# Grab metadata
meta <- colData(seurat_obj_sce_calb1)[, c("sample_id", "seurat_clusters")]
# counts = genes × cells
cts <- counts(seurat_obj_sce_calb1)
# build group labels: sample_cluster
groups <- paste(meta$sample_id, meta$seurat_clusters, sep = "_")
# aggregate counts across all cells belonging to each (sample, cluster) pair
pb <- aggregate.Matrix(t(cts), groupings = groups, fun = "sum")
pb <- t(pb) # back to genes × pseudobulks
dim(pb)
head(colnames(pb))
```
```{r}
# split sample_id and cluster from colnames
colinfo <- do.call(rbind, strsplit(colnames(pb), "_"))
coldata <- data.frame(
sample_id = colinfo[,1],
cluster = colinfo[,2],
row.names = colnames(pb)
)
head(coldata)
```
```{r}
dds <- DESeqDataSetFromMatrix(
countData = as.matrix(pb),
colData = coldata,
design = ~ sample_id + cluster
)
# filter low counts
dds <- dds[rowSums(counts(dds)) > 10, ]
# run DESeq2
dds <- DESeq(dds)
```
```{r}
# results: cluster0 (day7) vs cluster1 (ctrl)
res <- results(dds, contrast = c("cluster", "0", "1"))
res <- res[order(res$padj), ]
head(res)
```
```{r}
deg <- as.data.frame(res)
write.csv(deg, "DEG_neuron_cluster0_day7_vs_cluster1_day0_pseudobulk.csv")
```
```{r}
# Volcano plot
library(EnhancedVolcano)
p1 <- EnhancedVolcano(
deg,
lab = rownames(deg),
x = 'log2FoldChange',
y = 'pvalue',
title = "Cluster 0 (Day7) vs Cluster 1 (Day0)",
pCutoff = 0.01,
FCcutoff = 1,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1
)
print(p1)
```
```{r}
# load & filter
res_0_vs_1 <- read.csv("DEG_neuron_cluster0_day7_vs_cluster1_day0_pseudobulk.csv", row.names = 1)
res_0_vs_1 <- res_0_vs_1[!is.na(res_0_vs_1$padj),]
# FC > 0.5 (up and down)
res_0_vs_1_FC05_up <- rownames(res_0_vs_1[(res_0_vs_1$padj < 0.05) & (res_0_vs_1$log2FoldChange > 0.5),])
res_0_vs_1_FC05_down <- rownames(res_0_vs_1[(res_0_vs_1$padj < 0.05) & (res_0_vs_1$log2FoldChange < -0.5),])
```
```{r}
# Assuming res_lesion_Ctrl is your data frame
# 1. Filter the data based on your criteria
filtered_data <- res_0_vs_1[(res_0_vs_1$padj < 0.05) & (res_0_vs_1$log2FoldChange > 0.5),]
# 2. Rank the filtered data by 'padj' in ascending order
ranked_data <- filtered_data[order(filtered_data$padj),]
# 3. Export the ranked data to a CSV file
write.csv(ranked_data, "DGE_IO_neuron_Clus0vsClus1_FC05_ranked_up.csv")
# 1. Filter the data based on your criteria
filtered_data <- res_0_vs_1[(res_0_vs_1$padj < 0.05) & (res_0_vs_1$log2FoldChange < -0.5),]
# 2. Rank the filtered data by 'padj' in ascending order
ranked_data <- filtered_data[order(filtered_data$padj),]
# 3. Export the ranked data to a CSV file
write.csv(ranked_data, "DGE_IO_neuron_Clus0vsClus1_FC05_ranked_down.csv")
```
```{r}
dbs <- c("GO_Molecular_Function_2023",
"GO_Cellular_Component_2023",
"GO_Biological_Process_2023",
"KEGG_2019_Mouse")
enrichedR.list <- list()
enrichedR.list[["Cluster0_vs_1_FC05_up"]] <- enrichr(res_0_vs_1_FC05_up, dbs)
enrichedR.list[["Cluster0_vs_1_FC05_down"]] <- enrichr(res_0_vs_1_FC05_down, dbs)
```
```{r}
# Create an empty list to store the data frames for export
export_list <- list()
# Loop through your enrichedR.list and extract a specific database for each cluster
for (cluster_name in names(enrichedR.list)) {
# For example, extract the "GO_Biological_Process_2023" results
export_list[[cluster_name]] <- enrichedR.list[[cluster_name]]$GO_Biological_Process_2023
}
# Export the prepared list to a single Excel file
write.xlsx(export_list, file = "enrichment_IO_calb1_neuron_Cluster0vs1_DESeq_20250915.xlsx")
```
```{r}
pdf("enrichment_plot_IO_neuron_clus0_day7_vs_clus1_day0_filter.pdf", width=12, height=8)
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
p <- plotEnrich(
enr_df,
showTerms = 50,
numChar = 100,
y = "Count",
orderBy = "FDR",
title = paste0(dge," in ", go_db),
xlab = ""
)
print(
p + theme(
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 6),
plot.title = element_text(size = 10),
plot.margin = margin(10, 10, 10, 10)
)
)
}
}
dev.off()
```
```{r}
pdf("enrichment_plot_IO_neuron_clus0_day7_vs_clus1_day0.pdf", width=12, height=8)
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
p <- plotEnrich(
enr_df,
showTerms = 50,
numChar = 100,
y = "Count",
orderBy = "FDR",
title = paste0(dge," in ", go_db),
xlab = ""
)
print(
p + theme(
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 6),
plot.title = element_text(size = 10),
plot.margin = margin(10, 10, 10, 10)
)
)
}
}
dev.off()
```
```{r}
pdf("enrichment_plot_IO_neuron_clus0_day7_vs_clus1_day0_filter_20250905.pdf", width=6, height=4)
# Define your list of keywords
keywords <- c("matergic synapse", "Thermogenesis", "TCA", "oxidative phosphor")
# keywords <- c("Gluta")
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
# Filter for keywords and FDR < 0.1
filtered_df <- enr_df[
grepl(paste(keywords, collapse = "|"), enr_df$Term, ignore.case = TRUE) & enr_df$Adjusted.P.value < 0.1,
]
# Check if there are any results after filtering
if (nrow(filtered_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no terms found after filtering)")
next
}
# Use the filtered data for plotting
p <- plotEnrich(
filtered_df,
showTerms = 4,
numChar = 100,
y = "Count",
orderBy = "FDR",
title = paste0(dge, " in ", go_db, " (Filtered)"),
xlab = ""
)
print(
p + theme(
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 6),
plot.title = element_text(size = 10),
plot.margin = margin(10, 10, 10, 10)
)
)
}
}
dev.off()
```
### Day 14 vs Ctrl
```{r}
# results: cluster2 (day14) vs cluster1 (ctrl)
res <- results(dds, contrast = c("cluster", "2", "1"))
res <- res[order(res$padj), ]
head(res)
```
```{r}
deg <- as.data.frame(res)
write.csv(deg, "DEG_neuron_cluster2_day14_vs_cluster1_day0_pseudobulk.csv")
```
```{r}
# Volcano plot
library(EnhancedVolcano)
p2 <- EnhancedVolcano(
deg,
lab = rownames(deg),
x = 'log2FoldChange',
y = 'pvalue',
title = "Cluster 2 (Day14) vs Cluster 1 (Day0)",
pCutoff = 0.01,
FCcutoff = 1,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1
)
print(p2)
```
```{r}
# load & filter
res_2_vs_1 <- read.csv("DEG_neuron_cluster2_day14_vs_cluster1_day0_pseudobulk.csv", row.names = 1)
res_2_vs_1 <- res_2_vs_1[!is.na(res_2_vs_1$padj),]
# extract up/down
# FC > 0.5 (up and down)
res_2_vs_1_FC05_up <- rownames(res_2_vs_1[(res_2_vs_1$padj < 0.05) & (res_2_vs_1$log2FoldChange > 0.5),])
res_2_vs_1_FC05_down <- rownames(res_2_vs_1[(res_2_vs_1$padj < 0.05) & (res_2_vs_1$log2FoldChange < -0.5),])
```
```{r}
dbs <- c("GO_Molecular_Function_2023",
"GO_Cellular_Component_2023",
"GO_Biological_Process_2023",
"KEGG_2019_Mouse")
enrichedR.list <- list()
enrichedR.list[["Cluster2_vs_1_FC05_up"]] <- enrichr(res_2_vs_1_FC05_up, dbs)
enrichedR.list[["Cluster2_vs_1_FC05_down"]] <- enrichr(res_2_vs_1_FC05_down, dbs)
```
```{r}
# Create an empty list to store the data frames for export
export_list <- list()
# Loop through your enrichedR.list and extract a specific database for each cluster
for (cluster_name in names(enrichedR.list)) {
# For example, extract the "GO_Biological_Process_2023" results
export_list[[cluster_name]] <- enrichedR.list[[cluster_name]]$GO_Biological_Process_2023
}
# Export the prepared list to a single Excel file
write.xlsx(export_list, file = "enrichment_IO_calb1_neuron_Cluster2vs1_DESeq_20250915.xlsx")
```
```{r}
pdf("enrichment_plot_IO_neuron_clus2_day14_vs_clus1_day0.pdf", width=12, height=8)
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
p <- plotEnrich(
enr_df,
showTerms = 50,
numChar = 100,
y = "Count",
orderBy = "FDR",
title = paste0(dge," in ", go_db),
xlab = ""
)
print(
p + theme(
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 6),
plot.title = element_text(size = 10),
plot.margin = margin(10, 10, 10, 10)
)
)
}
}
dev.off()
```
### D14 vs D7
```{r}
# results: cluster 2 (day 14) cluster0 (day7)
res <- results(dds, contrast = c("cluster", "2", "0"))
res <- res[order(res$padj), ]
head(res)
```
```{r}
deg <- as.data.frame(res)
write.csv(deg, "DEG_neuron_cluster2_day14_vs_cluster0_day7_pseudobulk.csv")
```
```{r}
# Volcano plot
library(EnhancedVolcano)
p3 <- EnhancedVolcano(
deg,
lab = rownames(deg),
x = 'log2FoldChange',
y = 'pvalue',
title = "Cluster 2 (Day14) vs Cluster 0 (Day7)",
pCutoff = 0.01,
FCcutoff = 1,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1
)
print(p3)
```
```{r}
# load & filter
res_2_vs_0 <- read.csv("DEG_neuron_cluster2_day14_vs_cluster0_day7_pseudobulk.csv", row.names = 1)
res_2_vs_0 <- res_2_vs_0[!is.na(res_2_vs_0$padj),]
# FC > 0.5 (up and down)
res_2_vs_0_FC05_up <- rownames(res_2_vs_0[(res_2_vs_0$padj < 0.05) & (res_2_vs_0$log2FoldChange > 0.5),])
res_2_vs_0_FC05_down <- rownames(res_2_vs_0[(res_2_vs_0$padj < 0.05) & (res_2_vs_0$log2FoldChange < -0.5),])
```
```{r}
# Assuming res_lesion_Ctrl is your data frame
# 1. Filter the data based on your criteria
filtered_data <- res_2_vs_0[(res_2_vs_0$padj < 0.05) & (res_2_vs_0$log2FoldChange > 0.5),]
# 2. Rank the filtered data by 'padj' in ascending order
ranked_data <- filtered_data[order(filtered_data$padj),]
# 3. Export the ranked data to a CSV file
write.csv(ranked_data, "DGE_IO_neuron_Clus2vsClus0_FC05_ranked_up.csv")
# 1. Filter the data based on your criteria
filtered_data <- res_2_vs_0[(res_2_vs_0$padj < 0.05) & (res_2_vs_0$log2FoldChange < -0.5),]
# 2. Rank the filtered data by 'padj' in ascending order
ranked_data <- filtered_data[order(filtered_data$padj),]
# 3. Export the ranked data to a CSV file
write.csv(ranked_data, "DGE_IO_neuron_Clus2vsClus0_FC05_ranked_down.csv")
```
```{r}
dbs <- c("GO_Molecular_Function_2023",
"GO_Cellular_Component_2023",
"GO_Biological_Process_2023",
"KEGG_2019_Mouse")
enrichedR.list <- list()
enrichedR.list[["Cluster2_vs_0_FC05_up"]] <- enrichr(res_2_vs_0_FC05_up, dbs)
enrichedR.list[["Cluster2_vs_0_FC05_down"]] <- enrichr(res_2_vs_0_FC05_down, dbs)
```
```{r}
# Create an empty list to store the data frames for export
export_list <- list()
# Loop through your enrichedR.list and extract a specific database for each cluster
for (cluster_name in names(enrichedR.list)) {
# For example, extract the "GO_Biological_Process_2023" results
export_list[[cluster_name]] <- enrichedR.list[[cluster_name]]$KEGG_2019_Mouse
}
# Export the prepared list to a single Excel file
openxlsx::write.xlsx(export_list, file = "enrichment_IO_calb1_neuron_Cluster2vs0_DESeq2_KEGG_20250915.xlsx")
```
```{r}
pdf("enrichment_plot_IO_neuron_clus2_day14_vs_clus0_day7.pdf", width=12, height=8)
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
p <- plotEnrich(
enr_df,
showTerms = 50,
numChar = 100,
y = "Count",
orderBy = "FDR",
title = paste0(dge," in ", go_db),
xlab = ""
)
print(
p + theme(
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 6),
plot.title = element_text(size = 10),
plot.margin = margin(10, 10, 10, 10)
)
)
}
}
dev.off()
```
```{r}
pdf("enrichment_plot_IO_neuron_clus2_day14_vs_clus0_day7_filter.pdf", width=6, height=4)
# Define your list of keywords
keywords <- c("potentiation", "calcium signaling pathway", "depression", "synapse")
# keywords <- c("Gluta")
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
# Filter for keywords and FDR < 0.1
filtered_df <- enr_df[
grepl(paste(keywords, collapse = "|"), enr_df$Term, ignore.case = TRUE) & enr_df$Adjusted.P.value < 0.1,
]
# Check if there are any results after filtering
if (nrow(filtered_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no terms found after filtering)")
next
}
# Use the filtered data for plotting
p <- plotEnrich(
filtered_df,
showTerms = 4,
numChar = 100,
y = "Count",
orderBy = "FDR",
title = paste0(dge, " in ", go_db, " (Filtered)"),
xlab = ""
)
print(
p + theme(
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 6),
plot.title = element_text(size = 10),
plot.margin = margin(10, 10, 10, 10)
)
)
}
}
dev.off()
```
```{r}
NeuronClusterMarker <- c("COX1", "COX2", "COX3", "ATP6", "ND4", "ND5", "Caln1", "Calm1", "Camk2a", "Kcnab1", "Kcnc1", "Kcnip4", "Kcna1", "Kcnd3", "Cacna1c", "Cdh11", "Gria1", "Grin1", "Grin2a","Ajap1", "Wnt5a" )
p <- DotPlot(
seurat_obj_IO_calb1,
features = NeuronClusterMarker,
group.by = "seurat_clusters",
assay = "RNA_repaired"
) +
scale_color_gradient2(
low = "blue",
mid = "lightyellow",
high = "red",
midpoint = 0,
limits = c(-1,1),
oob = scales::squish
) +
ggtitle("IO Calb1 Neuron Cluster Markers") +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
axis.text.y = element_text(size = 10),
text = element_text(size = 12)
)
print(p)
ggsave("Dotplot_IO_Calb1_neuron_Cluster_markers.svg" , plot = p, width = 7, height = 4, units="in")
```
--------------------------------------------
# Microglia in the IO
*** obsolete codes due to loading .rds
```{r}
micro_Balazs_markers = c("Aif1", "Maf", "Spi1", "Csf1r")
seurat_obj_IO_calb1 <- AddModuleScore(
seurat_obj_IO_calb1,
features=list(micro_Balazs_markers),
nbin = 24,
ctrl = 100,
assay = "RNA_repaired",
slot = "data",
name = "AMS",
)
colnames([email hidden]) <- c(colnames([email hidden])[1:(length(colnames([email hidden]))-1)], "micro_Balazs_markers")
options(repr.plot.width=6, repr.plot.height=5)
FeaturePlot(seurat_obj_IO_calb1, reduction = "umap", features = "micro_Balazs_markers", order=TRUE, pt.size=1, raster=FALSE, min.cutoff=0)
seurat_obj_IO_calb1_micro <- subset(seurat_obj_IO_calb1, subset= micro_Balazs_markers > 0.05)
```
```{r}
# load the .rds file shared on 23rd Aug 2025
seurat_obj_IO_calb1_micro <- readRDS("~/rds/rds-rds-karadottir-qqstKNzWvy8/Omar_Spatial/seurat_obj_IO_calb1_micro_DIET.rds")
```
```{r}
options(repr.plot.width=11, repr.plot.height=9)
DimPlot(seurat_obj_IO_calb1, cells.highlight = colnames(seurat_obj_IO_calb1_micro), pt.size=0.5,cols.highlight = c("red", "grey50"), raster=FALSE)
```
```{r}
percent_micro <- (as.vector(table(seurat_obj_IO_calb1_micro$orig.ident))[1:12]*100)/as.vector(table(seurat_obj_IO$orig.ident))[1:12]
mat <- tibble(
group = c("Ctrl","Ctrl","Ctrl","Ctrl","D7","D7","D7","D7","D14","D14","D14","D14"),
measurement = percent_micro
)
mat$group <- factor(mat$group, levels = c("Ctrl", "D7", "D14"))
options(repr.plot.width=4, repr.plot.height=6)
barplot_with_error_and_dots(mat, x = group, y = measurement, error = "se", title="MG > 10% pixels (% of IO)")
```
```{r}
DefaultAssay(seurat_obj_IO_calb1_micro) <- "RNA"
seurat_obj_IO_calb1_micro <- NormalizeData(seurat_obj_IO_calb1_micro, normalization.method = "LogNormalize", scale.factor = 10000)
DefaultAssay(seurat_obj_IO_calb1_micro) <- "RNA_repaired"
seurat_obj_IO_calb1_micro <- NormalizeData(seurat_obj_IO_calb1_micro, normalization.method = "LogNormalize", scale.factor = 10000)
seurat_obj_IO_calb1_micro <- FindVariableFeatures(seurat_obj_IO_calb1_micro, selection.method = "vst", nfeatures = 1000)
options(repr.plot.width=10, repr.plot.height=7)
VariableFeaturePlot(object = seurat_obj_IO_calb1_micro, selection.method = "vst", log = TRUE)
```
```{r}
seurat_obj_IO_calb1_micro <- ScaleData(seurat_obj_IO_calb1_micro, features = rownames(seurat_obj_IO_calb1_micro))
```
```{r}
seurat_obj_IO_calb1_micro <- RunPCA(seurat_obj_IO_calb1_micro)
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
ElbowPlot(seurat_obj_IO_calb1_micro, reduction = "pca", ndims = 50)
```
```{r}
nb_pcs = 15
seurat_obj_IO_calb1_micro <- FindNeighbors(seurat_obj_IO_calb1_micro, reduction = "pca", dims = 1:nb_pcs, prune.SNN = 0)
cluster_res = 0.7
seurat_obj_IO_calb1_micro <- FindClusters(seurat_obj_IO_calb1_micro, resolution = cluster_res)
```
```{r}
table(seurat_obj_IO_calb1_micro$seurat_clusters)
```
```{r}
seurat_clusters_IO_calb1_micro_color <- list()
n_clusters <- 5
seurat_clusters_IO_calb1_micro_color <- color_palette_30[1:n_clusters]
names(seurat_clusters_IO_calb1_micro_color) <- seq(1:n_clusters)-1
seurat_obj_IO_calb1_micro <- RunUMAP(seurat_obj_IO_calb1_micro, reduction = "pca", dims = 1:nb_pcs)
```
```{r}
options(repr.plot.width=25, repr.plot.height=8)
DimPlot(seurat_obj_IO_calb1_micro, group.by="orig.ident", split.by="orig.ident_merge", cols=orig.ident_colors, shuffle=TRUE, seed=seed, ncol=4, pt.size=1, label=FALSE, raster=FALSE)
```
```{r}
options(repr.plot.width=11, repr.plot.height=9)
for(timepoint in names(orig.ident_merge_colors)){
highlight_cells <- rownames([email hidden][[email hidden]$orig.ident_merge==timepoint,])
print(DimPlot(seurat_obj_IO_calb1_micro, cells.highlight = highlight_cells, pt.size=0.5, cols.highlight = orig.ident_merge_colors[timepoint], raster=FALSE))
}
```
```{r}
options(repr.plot.width=11, repr.plot.height=9)
DimPlot(seurat_obj_IO_calb1_micro, group.by="seurat_clusters", cols=seurat_clusters_IO_calb1_micro_color, shuffle=TRUE, seed=seed, pt.size=1, label=FALSE, raster=FALSE)
```
```{r}
options(repr.plot.width=11, repr.plot.height=9)
p <- DimPlot(seurat_obj_IO_calb1_micro, group.by="seurat_clusters", cols=seurat_clusters_IO_calb1_micro_color, shuffle=TRUE, seed=seed, pt.size=1, label=FALSE, raster=FALSE) +NoAxes()
print(p)
ggsave("IO_calb1_microglia_Clusters_UMAP.svg", height = 3, width = 4, units = "in")
```
```{r}
options(repr.plot.width=25, repr.plot.height=8)
DimPlot(seurat_obj_IO_calb1_micro, group.by="seurat_clusters", split.by="orig.ident_merge", cols=seurat_clusters_IO_calb1_micro_color, shuffle=TRUE, seed=seed, ncol=4, pt.size=1, label=FALSE, raster=FALSE)
```
```{r}
options(repr.plot.width=15, repr.plot.height=15)
set.seed(115)
mat = as.matrix(table(seurat_obj_IO_calb1_micro$orig.ident_merge, seurat_obj_IO_calb1_micro$seurat_clusters))
for (row in 1:nrow(mat)) mat[row,] <- as.integer(mat[row,]/(sum(mat[row,])/min(rowSums(mat))))
circos.clear()
par(cex = 0.8)
chordDiagram(mat, annotationTrack = c("grid", "axis"), preAllocateTracks = list(track.height = max(strwidth(unlist(dimnames(mat))))/3), big.gap = 20, small.gap = 2, order = c(colnames(mat), names(orig.ident_merge_colors)), grid.col=c(rev(seurat_clusters_IO_calb1_micro_color), rev(orig.ident_merge_colors)))
circos.track(track.index = 1, panel.fun = function(x, y) {circos.text(CELL_META$xcenter, CELL_META$ylim[1], CELL_META$sector.index, facing = "clockwise", niceFacing = TRUE, adj = c(-0.3, 0.5))}, bg.border = NA)
```
```{r}
pdf("Chord_plot_IO_MgClus_20250906.pdf", width=8, height=8)
options(repr.plot.width=15, repr.plot.height=15)
set.seed(115)
mat = as.matrix(table(seurat_obj_IO_calb1_micro$orig.ident_merge, seurat_obj_IO_calb1_micro$seurat_clusters))
for (row in 1:nrow(mat)) mat[row,] <- as.integer(mat[row,]/(sum(mat[row,])/min(rowSums(mat))))
circos.clear()
par(cex = 0.8)
chordDiagram(mat, annotationTrack = c("grid", "axis"), preAllocateTracks = list(track.height = max(strwidth(unlist(dimnames(mat))))/3), big.gap = 20, small.gap = 2, order = c(colnames(mat), names(orig.ident_merge_colors)), grid.col=c(rev(seurat_clusters_IO_calb1_micro_color), rev(orig.ident_merge_colors)))
circos.track(track.index = 1, panel.fun = function(x, y) {circos.text(CELL_META$xcenter, CELL_META$ylim[1], CELL_META$sector.index, facing = "clockwise", niceFacing = TRUE, adj = c(-0.3, 0.5))}, bg.border = NA)
dev.off()
```
```{r}
t <- table(seurat_obj_IO_calb1_micro$orig.ident, seurat_obj_IO_calb1_micro$seurat_clusters)
t_prop <- t/as.vector(table(seurat_obj_IO$orig.ident))*100
rownames(t_prop) <- c("Ctrl","Ctrl","Ctrl","Ctrl","D7","D7","D7","D7","D14","D14","D14","D14")
df_table <- as.data.frame(t_prop)
colnames(df_table) <- c("group", "cluster", "count")
df_summary <- df_table %>%
group_by(group, cluster) %>%
summarise(
mean = mean(count),
se = sd(count) / sqrt(n()),
.groups = "drop"
)
```
```{r}
options(repr.plot.width=10, repr.plot.height=10)
ggplot(df_summary, aes(x = cluster, y = mean, fill = group)) +
# Bars
geom_bar(stat = "identity", position = position_dodge(width = 0.9), color = "black") +
# Error bars
geom_errorbar(aes(ymin = mean - se, ymax = mean + se),
position = position_dodge(width = 0.9),
width = 0.2) +
# Manual fill colors
scale_fill_manual(values = orig.ident_merge_colors) +
# Points dodged and jittered
geom_point(
data = df_table,
aes(x = cluster, y = count, fill = group), # <- fill included so dodge works
position = position_jitterdodge(jitter.width = 0.4, dodge.width = 0.9),
size = 2, shape = 21, stroke = 0.5,
color = "black", # border
show.legend = FALSE
) +
theme_minimal() +
labs(x = "", y = "% of pixels in the IO") +
ggtitle("Number of microglia pixels associated calbindin neurons divided by each IO for each sample") +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank())
```
```{r}
t <- table(Cluster=seurat_obj_IO_calb1_micro$seurat_clusters, Batch=[email hidden][["orig.ident_merge"]])
t <- t[,rev(names(orig.ident_merge_colors))]
options(repr.plot.width=10, repr.plot.height=10)
barplot(t(t), xlab = "Cluster", ylab="Number of cells", legend = TRUE, ylim = c(0, round_any(as.integer(max(rowSums(t)))+100, 20, f = ceiling)), col = rev(orig.ident_merge_colors), args.legend = list(bty = "n", x = "top", ncol = 3))
```
```{r}
t <- table(Cluster=seurat_obj_IO_calb1_micro$seurat_clusters, Batch=[email hidden][["orig.ident"]])
t <- t[,rev(names(orig.ident_colors))]
options(repr.plot.width=10, repr.plot.height=10)
barplot(t(t), xlab = "Cluster", ylab="Number of cells", legend = TRUE, ylim = c(0, round_any(as.integer(max(rowSums(t)))+200, 20, f = ceiling)), col = rev(orig.ident_colors), args.legend = list(bty = "n", x = "top", ncol = 3))
```
```{r}
DefaultAssay(seurat_obj_IO_calb1_micro) <- "RNA"
Idents(object = seurat_obj_IO_calb1_micro) <- "seurat_clusters"
seurat_obj_IO_calb1_micro.markers.seurat_clusters <- FindAllMarkers(seurat_obj_IO_calb1_micro, slot = "data", min.pct = 0.1, logfc.threshold = 0, only.pos = TRUE)
```
```{r}
write.csv(seurat_obj_IO_calb1_micro.markers.seurat_clusters, "IO_calb1_microglia_clusterMarkers_20250916.csv")
```
```{r}
seurat_obj_IO_calb1_micro.markers.seurat_clusters_pval <- seurat_obj_IO_calb1_micro.markers.seurat_clusters[seurat_obj_IO_calb1_micro.markers.seurat_clusters$p_val_adj < 0.05,]
seurat_obj_IO_calb1_micro.markers.seurat_clusters_best <- seurat_obj_IO_calb1_micro.markers.seurat_clusters_pval %>%
filter(avg_log2FC > 0.5) %>%
filter(pct.1 > 0.3) %>%
arrange(cluster, desc(avg_log2FC)) %>%
group_by(cluster)
markers_heatmap <- seurat_obj_IO_calb1_micro.markers.seurat_clusters_best %>% top_n(n = 10, wt = avg_log2FC)
topDiffCluster <- seurat_obj_IO_calb1_micro.markers.seurat_clusters_best %>% top_n(n = 200, wt = avg_log2FC)
sample_df <- as.data.frame(lapply(split(topDiffCluster, topDiffCluster$cluster), function(x) c(x$gene,rep("None", 200-length(x$gene)))))
colnames(sample_df) <- levels(topDiffCluster$cluster)
sample_df
```
```{r}
# keep all significant markers without top_n
seurat_obj_IO_calb1_micro.markers.seurat_clusters_pval <-
seurat_obj_IO_calb1_micro.markers.seurat_clusters %>%
filter(p_val_adj < 0.05)
seurat_obj_IO_calb1_micro.markers.seurat_clusters_best <-
seurat_obj_IO_calb1_micro.markers.seurat_clusters_pval %>%
filter(avg_log2FC > 0.5, pct.1 > 0.3) %>%
arrange(cluster, desc(avg_log2FC)) %>%
group_by(cluster)
# keep full list per cluster
topDiffCluster <- seurat_obj_IO_calb1_micro.markers.seurat_clusters_best
# split into clusters
gene_lists <- split(topDiffCluster$gene, topDiffCluster$cluster)
# find max list length
max_len <- max(lengths(gene_lists))
# pad shorter clusters with NA
sample_df <- as.data.frame(
lapply(gene_lists, function(x) c(x, rep(NA, max_len - length(x))))
)
# set cluster names as column names
colnames(sample_df) <- names(gene_lists)
# write to CSV
write.csv(sample_df, "~/rds/rds-rds-karadottir-qqstKNzWvy8/Omar_Spatial/Bastien20250827/ClusterMarker_IO_Calb1_MG_20250906.csv", row.names = FALSE)
```
```{r}
options(repr.plot.width=30, repr.plot.height=20)
DoMultiBarHeatmap(
subset(seurat_obj_IO_calb1_micro, downsample = 200),
features = unique(markers_heatmap$gene),
cells = NULL,
group.by = "seurat_clusters",
additional.group.by = c("orig.ident_merge"),
additional.group.sort.by = c("orig.ident_merge"),
cols.use = list(seurat_clusters=seurat_clusters_IO_calb1_micro_color, orig.ident_merge=orig.ident_merge_colors),
group.bar = TRUE,
disp.min = -2.5,
disp.max = NULL,
layer = "scale.data",
assay = "RNA_repaired",
label = TRUE,
size = 5.5,
hjust = 0,
angle = 45,
raster = TRUE,
draw.lines = TRUE,
lines.width = NULL,
group.bar.height = 0.02,
combine = TRUE
) + theme(text = element_text(size = 8))
```
```{r}
pdf("TissueMapping_IO_microglia_cluster_20250906.pdf", width = 5.3, height = 4.6)
for (sample_name in names(orig.ident_colors)){
annot_plot <- [email hidden][seurat_obj_IO_calb1_micro$orig.ident==sample_name,]
annot_plot_plus <- [email hidden][(seurat_obj_IO$orig.ident==sample_name & seurat_obj_IO$InfOlive=="core"),]
annot_plot_plus <- annot_plot_plus[!(rownames(annot_plot_plus) %in% rownames(annot_plot)),]
annot_plot_plus$seurat_clusters <- 99
annot_plot <- annot_plot[,c("x","y","seurat_clusters")]
annot_plot_plus <- annot_plot_plus[,c("x","y","seurat_clusters")]
annot_plot$seurat_clusters <- as.character(annot_plot$seurat_clusters)
annot_plot_plus$seurat_clusters <- as.character(annot_plot_plus$seurat_clusters)
annot_plot <- rbind(annot_plot_plus,annot_plot)
annot_plot$seurat_clusters <- factor(x = annot_plot$seurat_clusters, levels = c(names(seurat_clusters_IO_calb1_micro_color),99))
options(repr.plot.width=5.3, repr.plot.height=4.6)
df_tmp <- data.frame(
x = annot_plot$x,
y = annot_plot$y,
value = annot_plot$seurat_clusters
)
tmp_colors <- c(seurat_clusters_IO_calb1_micro_color,"lightgrey")
names(tmp_colors) <- c(names(seurat_clusters_IO_calb1_micro_color),99)
print(ggplot(df_tmp, aes(x, y)) +
geom_point(aes(colour = value),shape=15, size = 1)+
scale_color_manual(values=tmp_colors) +
scale_x_continuous(expand=c(-0.005,0), lim=c(0,101)) +
scale_y_continuous(expand=c(-0.006,0), lim=c(0,101)) +
theme_void() +
ggtitle(sample_name))
}
dev.off()
```
```{r}
# saveRDS(seurat_obj_IO_calb1_micro, "seurat_obj_IO_calb1_micro.rds")
```
```{r}
# load the 5 csv files
MC1_Volcano_df <- read.csv("goncaloGeneLists/MC1_filtered_genes.csv")
MC3_Volcano_df <- read.csv("goncaloGeneLists/MC3_filtered_genes.csv")
MC1_NonDE_df <- read.csv("goncaloGeneLists/MC1_Non_DE.csv")
MC2_NonDE_df <- read.csv("goncaloGeneLists/MC2_Non_DE.csv")
MC3_NonDE_df <- read.csv("goncaloGeneLists/MC3_Non_DE.csv")
```
```{r}
# Filter to make them all into lists
FDRCutOff <- 0.1
MC1_Volcano_df_fil <- MC1_Volcano_df[MC1_Volcano_df$FDR.x < FDRCutOff , ]
MC1_Volcano_list <- MC1_Volcano_df_fil$gene
MC3_Volcano_df_fil <- MC3_Volcano_df[MC3_Volcano_df$FDR.x < FDRCutOff , ]
MC3_Volcano_list <- MC3_Volcano_df_fil$gene
MC1_NonDE_df_fil <- MC1_NonDE_df[MC1_NonDE_df$FDR < FDRCutOff , ]
MC1_NonDE_df_list <- MC1_NonDE_df_fil$name
MC2_NonDE_df_fil <- MC2_NonDE_df[MC2_NonDE_df$FDR < FDRCutOff , ]
MC2_NonDE_df_list <- MC2_NonDE_df_fil$name
MC3_NonDE_df_fil <- MC3_NonDE_df[MC3_NonDE_df$FDR < FDRCutOff , ]
MC3_NonDE_df_list <- MC3_NonDE_df_fil$name
```
```{r}
# Perform gene scoring and plot the violin plots for different clusters
seurat_obj_IO_calb1_micro <- AddModuleScore(
object = seurat_obj_IO_calb1_micro,
features = list(MC1_Volcano_list, MC3_Volcano_list),
name = c("MC1_Score", "MC3_Score")
)
# violin plots
pMC1 <- VlnPlot(seurat_obj_IO_calb1_micro, features = c("MC1_Score1"), group.by = "seurat_clusters") +
ggtitle("IO microglia MC1 Gene Module Score")
print(pMC1)
pMC3 <- VlnPlot(seurat_obj_IO_calb1_micro, features = c("MC3_Score2"), group.by = "seurat_clusters") +
ggtitle("IO microglia MC3 Gene Module Score")+
scale_y_continuous(limits = c(-0.1, 0.7), breaks = c(0, 0.2, 0.4, 0.6))
print(pMC3)
FeaturePlot(seurat_obj_IO_calb1_micro, features = c("MC1_Score1", "MC3_Score2"))
# export the violin plots
ggsave("IO_calb1microglia_MC1score_20250909.svg", plot = pMC1)
ggsave("IO_calb1microglia_MC3score_20250909.svg", plot = pMC3)
```
# DAM score
```{r}
# load the 5 csv files
DAM_geneList_df <- read.csv("DAM_geneList.csv")
KS_DAM_geneList_df <- read.csv("Keren-Shaul DAM list.csv")
```
```{r}
# Filter to make them all into lists
DAM_geneList <- DAM_geneList_df$DAM_geneList
KS_DAM_geneList_DAM <- KS_DAM_geneList_df$Keren_shaul_DAM_up
KS_DAM_geneList_HM <- KS_DAM_geneList_df$Keren_shaul_DAM_down
KS_DAM_geneList_HM <- KS_DAM_geneList_HM[KS_DAM_geneList_HM!=""]
```
```{r}
# Perform gene scoring and plot the violin plots for different clusters
seurat_obj_IO_calb1_micro <- AddModuleScore(
object = seurat_obj_IO_calb1_micro,
features = list(DAM_geneList, KS_DAM_geneList_DAM, KS_DAM_geneList_HM),
name = c("DAM_geneList", "KS_DAM_geneList_DAM", "KS_DAM_geneList_HM")
)
# violin plots
VlnPlot(seurat_obj_IO_calb1_micro, features = c("DAM_geneList1"), group.by = "seurat_clusters") +
ggtitle("DAM_geneList VolcanoPlot Gene Module Score")
VlnPlot(seurat_obj_IO_calb1_micro, features = c("KS_DAM_geneList_DAM2"), group.by = "seurat_clusters") +
ggtitle("KS_DAM_geneList_DAM VolcanoPlot Gene Module Score")
VlnPlot(seurat_obj_IO_calb1_micro, features = c("KS_DAM_geneList_HM3"), group.by = "seurat_clusters") +
ggtitle("KS_DAM_geneList_DAM VolcanoPlot Gene Module Score")
VlnPlot(seurat_obj_IO_calb1_micro, features = c("DAM_geneList1"), group.by = "orig.ident_merge") +
ggtitle("DAM_geneList VolcanoPlot Gene Module Score")
VlnPlot(seurat_obj_IO_calb1_micro, features = c("KS_DAM_geneList_DAM2"), group.by = "orig.ident_merge") +
ggtitle("KS_DAM_geneList_DAM VolcanoPlot Gene Module Score")
VlnPlot(seurat_obj_IO_calb1_micro, features = c("KS_DAM_geneList_HM3"), group.by = "orig.ident_merge") +
ggtitle("KS_DAM_geneList_DAM VolcanoPlot Gene Module Score")
FeaturePlot(seurat_obj_IO_calb1_micro, features = c("MC1_Score1", "MC3_Score2"))
```
```{r}
pdf("DAM_IO_microglia_calb1.pdf", width=12, height=8)
# violin plots
VlnPlot(seurat_obj_IO_calb1_micro, features = c("DAM_geneList1"), group.by = "seurat_clusters") +
ggtitle("DAM geneList (from 3dpl) Violin Plot Gene Module Score")
VlnPlot(seurat_obj_IO_calb1_micro, features = c("KS_DAM_geneList_DAM2"), group.by = "seurat_clusters") +
ggtitle("KS DAM Violin Plot Gene Module Score")
VlnPlot(seurat_obj_IO_calb1_micro, features = c("KS_DAM_geneList_HM3"), group.by = "seurat_clusters") +
ggtitle("KS homeostasis Violin Plot Gene Module Score")
VlnPlot(seurat_obj_IO_calb1_micro, features = c("DAM_geneList1"), group.by = "orig.ident_merge") +
ggtitle("DAM geneList (from 3dpl) Violin Plot Gene Module Score")
VlnPlot(seurat_obj_IO_calb1_micro, features = c("KS_DAM_geneList_DAM2"), group.by = "orig.ident_merge") +
ggtitle("KS DAM Violin Plot Gene Module Score")
VlnPlot(seurat_obj_IO_calb1_micro, features = c("KS_DAM_geneList_HM3"), group.by = "orig.ident_merge") +
ggtitle("KS homeostasis Violin Plot Gene Module Score")
dev.off()
```
## NeuroProtective score
```{r}
# load the .csv files
Pro_GeneList_df <- read.csv("NeuProtectiveGeneList.csv")
```
```{r}
# get the different list
Protective_Thora_list <- Pro_GeneList_df$Thora_List
Protective_Full_list <- Pro_GeneList_df$Full_List
Protective_BV_list <- Pro_GeneList_df$BV_List
Protective_Thora_list <- Protective_Thora_list[Protective_Thora_list!=""]
Protective_Thora_list<- paste0(
toupper(substr(Protective_Thora_list, 1, 1)),
tolower(substr(Protective_Thora_list, 2, nchar(Protective_Thora_list)))
)
Protective_BV_list <- Protective_BV_list[Protective_BV_list!=""]
Protective_BV_list<- paste0(
toupper(substr(Protective_BV_list, 1, 1)),
tolower(substr(Protective_BV_list, 2, nchar(Protective_BV_list)))
)
Protective_Full_list <- Protective_Full_list[Protective_Full_list!=""]
Protective_Full_list<- paste0(
toupper(substr(Protective_Full_list, 1, 1)),
tolower(substr(Protective_Full_list, 2, nchar(Protective_Full_list)))
)
Trimmed_List <- c("Bdnf", "Igf1", "Fgf13")
# Trimmed_List_cl3 <-
```
```{r}
# Perform gene scoring and plot the violin plots for different clusters
seurat_obj_IO_calb1_micro <- AddModuleScore(
object = seurat_obj_IO_calb1_micro,
features = list(Trimmed_List),
name = c( "Trimmed_List")
)
# violin plots
VlnPlot(seurat_obj_IO_calb1_micro, features = c("Trimmed_List1"), group.by = "seurat_clusters") +
ggtitle("Bdnf and Igf1 Violin Plot Gene Module Score")
VlnPlot(seurat_obj_IO_calb1_micro, features = c("Trimmed_List1"), group.by = "orig.ident_merge") +
ggtitle("Bdnf, Igf1, Fgf13 Violin Plot Gene Module Score")
```
-----------------------
```{r}
# Extract metadata
df <- [email hidden] %>%
dplyr::select(seurat_clusters, orig.ident_merge)
# Count cells per cluster and per day
df_summary <- df %>%
group_by(orig.ident_merge, seurat_clusters) %>%
summarise(n = n(), .groups = "drop") %>%
group_by(orig.ident_merge) %>%
mutate(freq = n / sum(n)) # normalize to fractions
# Plot stacked barplot
ggplot(df_summary, aes(x = orig.ident_merge, y = freq, fill = seurat_clusters)) +
geom_bar(stat = "identity", position = "fill") +
ylab("Fraction of Cells") +
xlab("Day") +
scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggplot(df_summary, aes(x = orig.ident_merge, y = freq, fill = seurat_clusters)) +
geom_bar(stat = "identity", position = "fill") +
geom_text(aes(label = scales::percent(freq, accuracy = 1)),
position = position_stack(vjust = 0.5), size = 3) +
ylab("Fraction of Cells") +
xlab("Day") +
scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggplot(df_summary, aes(x = orig.ident_merge, y = freq, fill = seurat_clusters)) +
geom_bar(stat = "identity", position = "fill") +
geom_text(aes(label = scales::percent(freq, accuracy = 1)),
position = position_stack(vjust = 0.5), size = 3) +
ylab("Fraction of Cells") +
xlab("Day") +
scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))+
scale_fill_manual(values = seurat_clusters_IO_calb1_color)
```
```{r}
# export the bar plot
barPlot <- ggplot(df_summary, aes(x = orig.ident_merge, y = freq, fill = seurat_clusters)) +
geom_bar(stat = "identity", position = "fill") +
geom_text(aes(label = scales::percent(freq, accuracy = 1)),
position = position_stack(vjust = 0.5), size = 3) +
ylab("Fraction of Cells") +
xlab("Day") +
scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))+
scale_fill_manual(values = seurat_clusters_IO_calb1_color)
ggsave(filename = "IO_calb1Microglia_Cluster_barplot_20250907.svg", plot = barPlot, width = 10, height = 8, units = "in")
```
## Use the Marker genes for GO
```{r}
subset_list_len <- 200
# extract gene lists from sample_df
Clus0_markers <- sample_df[1:subset_list_len,1]
Clus1_markers <- sample_df[1:subset_list_len,2]
Clus2_markers <- sample_df[1:subset_list_len,3]
Clus3_markers <- sample_df[1:subset_list_len,4]
Clus4_markers <- sample_df[1:subset_list_len,5]
# dbs set up
dbs <- c("GO_Molecular_Function_2023","GO_Cellular_Component_2023","GO_Biological_Process_2023","KEGG_2019_Mouse")
#
enrichedR.list <- list()
enrichedR.list[["Clus0_markers"]] <- enrichr(Clus0_markers, dbs)
enrichedR.list[["Clus1_markers"]] <- enrichr(Clus1_markers, dbs)
enrichedR.list[["Clus2_markers"]] <- enrichr(Clus2_markers, dbs)
enrichedR.list[["Clus3_markers"]] <- enrichr(Clus3_markers, dbs)
enrichedR.list[["Clus4_markers"]] <- enrichr(Clus4_markers, dbs)
```
```{r}
# Create an empty list to store the data frames for export
export_list <- list()
# Loop through your enrichedR.list and extract a specific database for each cluster
for (cluster_name in names(enrichedR.list)) {
# For example, extract the "GO_Biological_Process_2023" results
export_list[[cluster_name]] <- enrichedR.list[[cluster_name]]$GO_Biological_Process_2023
}
# Export the prepared list to a single Excel file
openxlsx::write.xlsx(export_list, file = "enrichment_IO_calb1_microglia_ClusterMarker_GO_KEGG_20250915.xlsx")
```
```{r}
pdf("enrichment_plot_IO_Calb1MgClus_filSample_20250915.pdf", width=12, height=8)
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
p <- plotEnrich(
enr_df,
showTerms = 60,
numChar = 100,
y = "Count",
orderBy = "FDR",
title = paste0(dge," in ", go_db),
xlab = ""
)
print(
p + theme(
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 6),
plot.title = element_text(size = 10),
plot.margin = margin(10, 10, 10, 10)
)
)
}
}
dev.off()
```
```{r}
pdf("enrichment_plot_IO_Calb1MgClus_filterGO_20250915.pdf", width=8, height=3)
# Define your list of keywords
# for D14 clus 0
clus0_keywords <- c("Positive Regulation Of Receptor", "0030198",
"Cellular Response To Growth Factor Stimulus", "Positive Regulation Of Cell Migration")
# for D7
clus2_keywords <- c("Positive Regulation of Calcium Ion Transmembrane Transporter Activity", "Neuron Projection Morphogenesis", "Neurotransmitter Secretion", "Positive Regulation Of Cation Channel Activity")
keywords <- c(clus0_keywords, clus2_keywords)
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
# Filter for keywords and FDR < 0.1
filtered_df <- enr_df[
grepl(paste(keywords, collapse = "|"), enr_df$Term, ignore.case = TRUE) & enr_df$Adjusted.P.value < 0.1,
]
# Check if there are any results after filtering
if (nrow(filtered_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no terms found after filtering)")
next
}
# Use the filtered data for plotting
p <- plotEnrich(
filtered_df,
showTerms = 4,
numChar = 100,
y = "Count",
orderBy = "FDR",
title = paste0(dge, " in ", go_db, " (Filtered)"),
xlab = ""
)
print(
p + theme(
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 6),
plot.title = element_text(size = 10),
plot.margin = margin(10, 10, 10, 10)
)
)
}
}
dev.off()
```
```{r}
pdf("enrichment_dotplot_IO_Calb1MgClus_filterGO_20250916.pdf", width=8, height=3)
clus0_keywords <- c("Positive Regulation Of Receptor", "0030198",
"Cellular Response To Growth Factor Stimulus", "Positive Regulation Of Cell Migration")
# for D7
clus2_keywords <- c("Positive Regulation of Calcium Ion Transmembrane Transporter Activity", "Neuron Projection Morphogenesis", "Neurotransmitter Secretion", "Positive Regulation Of Cation Channel Activity")
keywords <- c(clus0_keywords, clus2_keywords)
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
# ---- Prepare data ----
# Adjusted P value column in enrichR is usually "Adjusted.P.value"
if (!("Adjusted.P.value" %in% colnames(enr_df))) {
stop("No Adjusted.P.value column found in enrichment results for ", dge, " / ", go_db)
}
# Significance (FDR)
enr_df$FDR <- as.numeric(enr_df$Adjusted.P.value)
# Overlap is in "k/n" format -> extract numerator and denominator
tmp <- do.call(rbind, strsplit(enr_df$Overlap, "/"))
enr_df$Count <- as.numeric(tmp[,1])
enr_df$GeneRatio <- as.numeric(tmp[,1]) / as.numeric(tmp[,2])
# Keep top 50 terms by FDR
enr_df <- enr_df[order(enr_df$FDR), ][1:min(50, nrow(enr_df)), ]
# Filter for keywords and FDR < 0.1
filtered_df <- enr_df[
grepl(paste(keywords, collapse = "|"), enr_df$Term, ignore.case = TRUE) & enr_df$Adjusted.P.value < 0.1,
]
# Check if there are any results after filtering
if (nrow(filtered_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no terms found after filtering)")
next
}
enr_df <- filtered_df
# ---- Dotplot ----
p <- ggplot(enr_df, aes(
x = GeneRatio,
y = reorder(Term, -FDR), # most significant (lowest FDR) at top
size = Count,
color = -log10(FDR)
)) +
geom_point() +
scale_color_gradient(low = "blue", high = "red") +
scale_size_continuous(range = c(2, 5)) +
labs(
title = paste0(dge," in ", go_db),
x = "Gene ratio",
y = "GO term",
color = "-log10(FDR)",
size = "Gene count"
) +
theme_bw() +
theme(
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 6),
plot.title = element_text(size = 10),
plot.margin = margin(10, 10, 10, 10)
)
print(p)
}
}
dev.off()
```
```{r}
library(GOplot)
pdf("enrichment_GOBubbleplot_IO_Calb1MgClus_filterGO_20250916.pdf", width=12, height=8)
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
# Prepare enrichment data
tmp <- do.call(rbind, strsplit(enr_df$Overlap, "/"))
Count <- as.numeric(tmp[,1])
ID <- grepl(" ", enr_df$Term, value = 1 )
GO_data <- data.frame(
category = go_db,
ID = paste0("GO:", seq_len(nrow(enr_df))), # fake GO IDs
term = enr_df$Term,
genes = enr_df$Genes, # enrichR column listing genes
adj_pval = enr_df$Adjusted.P.value,
zscore = NA,
Count = Count
)
# Split the genes into a vector, trim spaces
all_genes <- unique(trimws(unlist(strsplit(enr_df$Genes, ";"))))
# Dummy DEG table with matching names
DEG_data <- data.frame(ID = all_genes, logFC = 1)
# Merge for GOplot
circ <- circle_dat(GO_data, DEG_data)
# Plot GO bubble
p <- GOBubble(circ, display = "single", title = paste0(dge, " in ", go_db))
print(p)
}
}
dev.off()
```
# 9. DGE
```{r}
condition <- "IO"
HVGs.list <- list()
res.list <- list()
HVGs.list[[condition]] <- list()
res.list[[condition]] <- list()
```
```{r}
DefaultAssay(seurat_obj_IO_calb1_micro) <- 'RNA'
table(seurat_obj_IO_calb1_micro$orig.ident_merge)
```
```{r}
seurat_obj_sce_calb1_micro <- as.SingleCellExperiment(seurat_obj_IO_calb1_micro, assay="RNA")
seurat_obj_sce_calb1_micro <- prepSCE(seurat_obj_sce_calb1_micro, kid = "InfOlive", gid = "orig.ident_merge", sid = "orig.ident", drop = FALSE)
kids <- purrr::set_names(levels(seurat_obj_sce_calb1_micro$cluster_id))
# Total number of clusters
nk <- length(kids)
# Named vector of sample names
sids <- purrr::set_names(levels(seurat_obj_sce_calb1_micro$sample_id))
# Total number of samples
ns <- length(sids)
kids
sids
```
```{r}
## Determine the number of cells per sample
table(seurat_obj_sce_calb1_micro$sample_id)
## Turn named vector into a numeric vector of number of cells per sample
n_cells <- as.numeric(table(seurat_obj_sce_calb1_micro$sample_id))
## Determine how to reoder the samples (rows) of the metadata to match the order of sample names in sids vector
m <- match(sids, seurat_obj_sce_calb1_micro$sample_id)
## Create the sample level metadata by combining the reordered metadata with the number of cells corresponding to each sample.
ei <- data.frame(colData(seurat_obj_sce_calb1_micro)[m, ], n_cells, row.names = NULL) %>% dplyr::select(-"cluster_id")
ei
```
```{r}
# Aggregate the counts per sample_id and cluster_id
# Subset metadata to only include the cluster and sample IDs to aggregate across
groups <- colData(seurat_obj_sce_calb1_micro)[, c("cluster_id", "sample_id")]
# Aggregate across cluster-sample groups
pb <- aggregate.Matrix(t(counts(seurat_obj_sce_calb1_micro)), groupings = groups, fun = "sum")
splitf <- sapply(stringr::str_split(rownames(pb), pattern = "_", n = 2), `[`, 1)
pb <- split.data.frame(pb, factor(splitf)) %>% lapply(function(u) magrittr::set_colnames(t(u), stringr::str_extract(rownames(u), "(?<=_)[:alnum:]+")))
class(pb)
# Explore the different components of list
str(pb)
```
```{r}
options(width = 100)
table(seurat_obj_sce_calb1_micro$cluster_id, seurat_obj_sce_calb1_micro$group_id)
```
```{r}
options(width = 100)
table(seurat_obj_sce_calb1_micro$cluster_id, seurat_obj_sce_calb1_micro$seurat_clusters)
```
```{r}
# prep. data.frame for plotting
get_sample_ids <- function(x){pb[[x]] %>% colnames()}
de_samples <- purrr::map(1:length(kids), get_sample_ids) %>% unlist()
samples_list <- purrr::map(1:length(kids), get_sample_ids)
get_cluster_ids <- function(x){rep(names(pb)[x],each = length(samples_list[[x]]))}
de_cluster_ids <- purrr::map(1:length(kids), get_cluster_ids) %>% unlist()
gg_df <- data.frame(cluster_id = de_cluster_ids, sample_id = de_samples)
gg_df <- left_join(gg_df, ei[, c("sample_id", "group_id")])
metadata <- gg_df %>% dplyr::select(cluster_id, sample_id, group_id)
levels(metadata$cluster_id) <- unique(metadata$cluster_id)
# Generate vector of cluster IDs
clusters <- levels(metadata$cluster_id)
clusters
cluster_choice <- "core"
table_sample_low <- table(seurat_obj_sce_calb1_micro$cluster_id, seurat_obj_sce_calb1_micro$sample_id)[cluster_choice,][table(seurat_obj_sce_calb1_micro$cluster_id, seurat_obj_sce_calb1_micro$sample_id)[cluster_choice,] <= 5]
table_sample_low
sample_keep <- names(orig.ident_colors)[!names(orig.ident_colors) %in% names(table_sample_low)]
table_group_low <-table(seurat_obj_sce_calb1_micro$cluster_id, seurat_obj_sce_calb1_micro$group_id)[cluster_choice,][table(seurat_obj_sce_calb1_micro$cluster_id, seurat_obj_sce_calb1_micro$group_id)[cluster_choice,] <= 30]
table_group_low
group_keep <- names(orig.ident_merge_colors)[!names(orig.ident_merge_colors) %in% names(table_group_low)]
cluster_metadata <- metadata[which(metadata$cluster_id == cluster_choice), ]
# Assign the rownames of the metadata to be the sample IDs
rownames(cluster_metadata) <- cluster_metadata$sample_id
#Remove sample and timepoint with low counts (sample<=5, timepoint<=30)
cluster_metadata <- cluster_metadata[(cluster_metadata$sample_id %in% sample_keep & cluster_metadata$group_id %in% group_keep),]
new_group_levels <- levels(cluster_metadata$group_id)[levels(cluster_metadata$group_id) %in% group_keep]
cluster_metadata$group_id <- droplevels(cluster_metadata$group_id)
levels(cluster_metadata$group_id) <- new_group_levels
# Subset the counts to only this cluster
counts <- pb[[cluster_choice]]
cluster_counts <- data.frame(counts[, which(colnames(counts) %in% rownames(cluster_metadata))])
# Check that all of the row names of the metadata are the same and in the same order as the column names of the counts in order to use as input to DESeq2
all(rownames(cluster_metadata) == colnames(cluster_counts))
```
```{r}
dds <- DESeqDataSetFromMatrix(cluster_counts, colData = cluster_metadata, design = ~ group_id)
```
```{r}
# Transform counts for data visualization
rld <- rlog(dds, blind=TRUE)
```
```{r}
# Plot PCA
options(repr.plot.width=8, repr.plot.height=5)
DESeq2::plotPCA(rld, intgroup = "group_id")
```
```{r}
# Extract the rlog matrix from the object and compute pairwise correlation values
rld_mat <- assay(rld)
rld_cor <- cor(rld_mat)
```
```{r}
# Plot heatmap
options(repr.plot.width=10, repr.plot.height=8)
pheatmap(rld_cor, annotation = cluster_metadata[, "group_id", drop=F], annotation_colors = list(group_id=orig.ident_merge_colors))
```
```{r}
dds <- DESeq(dds)
```
```{r}
libsize = data.frame(x=sizeFactors(dds), y=colSums(assay(dds)))
options(repr.plot.width=7, repr.plot.height=7)
ggplot(data=libsize, aes(x=x, y=y)) + geom_point() + geom_smooth(method="lm") + xlab("Estimated size factor") + ylab("Library size")
```
```{r}
data.frame(colData(dds))
```
```{r}
# Plot dispersion estimates
options(repr.plot.width=12, repr.plot.height=12)
plotDispEsts(dds)
```
```{r}
HVGs.list[[condition]][["upregulated"]] <- list()
HVGs.list[[condition]][["downregulated"]] <- list()
```
## 9.1 Contrasts
```{r}
condition1 = "D7"
condition2 = "Ctrl"
contrast <- c("group_id", condition1, condition2)
# resultsNames(dds)
res <- results(dds, contrast = contrast, alpha = 0.05)
summary(res)
options(repr.plot.width=10, repr.plot.height=7)
hist(res$pvalue, breaks=0:20/20, col="grey50", border="white", xlim=c(0,1), main="Histogram of p-values", xlab="p-value")
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
hist(res$padj, breaks=0:20/20, col="grey50", border="white", xlim=c(0,1), main="Histogram of adj p-values", xlab="adj p-value")
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
qs <- c(0, quantile(res$baseMean[res$baseMean > 0], 0:7/7))
# cut the genes into the bins
bins <- cut(res$baseMean, qs)
# rename the levels of the bins using the middle point
levels(bins) <- paste0("~",round(.5*qs[-1] + .5*qs[-length(qs)]))
# calculate the ratio of $p$ values less than .01 for each bin
ratios <- tapply(res$pvalue, bins, function(p) mean(p < .01, na.rm=TRUE))
# plot these ratios
barplot(ratios, xlab="mean normalized count", ylab="ratio of small p values")
```
```{r}
options(repr.plot.width=15, repr.plot.height=15)
#padj<0.001
p1<-EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'padj',
title = paste0(condition1,' vs ',condition2,' in FC > 1'),
pCutoff = 0.01,
FCcutoff = 1,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1
)
print(p1)
p1_1<-EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'padj',
title = paste0(condition1,' vs ',condition2,' in FC > 0.75'),
pCutoff = 0.01,
FCcutoff = 0.75,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1
)
print(p1_1)
```
```{r}
res.list[[condition]][[paste0(condition1,"_",condition2)]] <- res
```
```{r}
write.csv(res, paste0("DGE_",condition,"_calb1_mg_",condition1,"_",condition2,"_20250915.csv"), row.names=TRUE)
```
```{r}
res <- res[complete.cases(res$padj),]
res_filtered <- res[(res$padj < 0.05 & res$baseMean > 1),]
HVGs.list[[condition]][["upregulated"]][[paste0(condition1,"_",condition2)]] <- rownames(res_filtered[res_filtered$log2FoldChange>=0,])
HVGs.list[[condition]][["downregulated"]][[paste0(condition1,"_",condition2)]] <- rownames(res_filtered[res_filtered$log2FoldChange<0,])
```
## D14 vs Ctrl
```{r}
condition1 = "D14"
condition2 = "Ctrl"
```
```{r}
contrast <- c("group_id", condition1, condition2)
# resultsNames(dds)
res <- results(dds, contrast = contrast, alpha = 0.05)
```
```{r}
summary(res)
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
hist(res$pvalue, breaks=0:20/20, col="grey50", border="white", xlim=c(0,1), main="Histogram of p-values", xlab="p-value")
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
hist(res$padj, breaks=0:20/20, col="grey50", border="white", xlim=c(0,1), main="Histogram of adj p-values", xlab="adj p-value")
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
qs <- c(0, quantile(res$baseMean[res$baseMean > 0], 0:7/7))
# cut the genes into the bins
bins <- cut(res$baseMean, qs)
# rename the levels of the bins using the middle point
levels(bins) <- paste0("~",round(.5*qs[-1] + .5*qs[-length(qs)]))
# calculate the ratio of $p$ values less than .01 for each bin
ratios <- tapply(res$pvalue, bins, function(p) mean(p < .01, na.rm=TRUE))
# plot these ratios
barplot(ratios, xlab="mean normalized count", ylab="ratio of small p values")
```
```{r}
options(repr.plot.width=15, repr.plot.height=15)
#padj<0.001
p2 <- EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'padj',
title = paste0(condition1,' vs ',condition2,' in FC > 1'),
pCutoff = 0.01,
FCcutoff = 1,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1)
print(p2)
p2_1 <- EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'padj',
title = paste0(condition1,' vs ',condition2,' in FC > 0.75'),
pCutoff = 0.01,
FCcutoff = 0.75,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1)
print(p2_1)
```
```{r}
res.list[[condition]][[paste0(condition1,"_",condition2)]] <- res
```
```{r}
write.csv(res, paste0("DGE_",condition,"_calb1_mg_",condition1,"_",condition2,"_20250915.csv"), row.names=TRUE)
```
```{r}
res <- res[complete.cases(res$padj),]
res_filtered <- res[(res$padj < 0.05 & res$baseMean > 1),]
HVGs.list[[condition]][["upregulated"]][[paste0(condition1,"_",condition2)]] <- rownames(res_filtered[res_filtered$log2FoldChange>=0,])
HVGs.list[[condition]][["downregulated"]][[paste0(condition1,"_",condition2)]] <- rownames(res_filtered[res_filtered$log2FoldChange<0,])
```
## D14 vs D7
```{r}
condition1 = "D14"
condition2 = "D7"
```
```{r}
contrast <- c("group_id", condition1, condition2)
# resultsNames(dds)
res <- results(dds, contrast = contrast, alpha = 0.05)
```
```{r}
summary(res)
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
hist(res$pvalue, breaks=0:20/20, col="grey50", border="white", xlim=c(0,1), main="Histogram of p-values", xlab="p-value")
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
hist(res$padj, breaks=0:20/20, col="grey50", border="white", xlim=c(0,1), main="Histogram of adj p-values", xlab="adj p-value")
```
```{r}
options(repr.plot.width=10, repr.plot.height=7)
qs <- c(0, quantile(res$baseMean[res$baseMean > 0], 0:7/7))
# cut the genes into the bins
bins <- cut(res$baseMean, qs)
# rename the levels of the bins using the middle point
levels(bins) <- paste0("~",round(.5*qs[-1] + .5*qs[-length(qs)]))
# calculate the ratio of $p$ values less than .01 for each bin
ratios <- tapply(res$pvalue, bins, function(p) mean(p < .01, na.rm=TRUE))
# plot these ratios
barplot(ratios, xlab="mean normalized count", ylab="ratio of small p values")
```
```{r}
options(repr.plot.width=15, repr.plot.height=15)
#padj<0.001
p3 <- EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'padj',
title = paste0(condition1,' vs ',condition2,' in FC > 1'),
pCutoff = 0.01,
FCcutoff = 1,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1)
print(p3)
p3_1 <- EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'padj',
title = paste0(condition1,' vs ',condition2,' in FC > 0.5'),
pCutoff = 0.01,
FCcutoff = 0.5,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1)
print(p3_1)
# Define the genes you want to highlight
genes_to_highlight <- c("Hexa", "Hexb", "Apoe", "Clu", "B2m", "Megf10", "Dock1")
p3_1 <- EnhancedVolcano(
res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'padj',
title = paste0(condition1,' vs ',condition2,' in FC > 0.5'),
pCutoff = 0.05,
FCcutoff = 0.5,
pointSize = 3.0,
labSize = 6.0,
col = c('black', 'grey', 'grey', 'blue'),
colAlpha = 1,
# Force-highlight these genes
selectLab = genes_to_highlight,
drawConnectors = TRUE
)
print(p3_1)
```
```{r}
res.list[[condition]][[paste0(condition1,"_",condition2)]] <- res
```
```{r}
write.csv(res, paste0("DGE_",condition,"_calb1_mg_",condition1,"_",condition2,"_20250915.csv"), row.names=TRUE)
```
```{r}
res <- res[complete.cases(res$padj),]
res_filtered <- res[(res$padj < 0.05 & res$baseMean > 1),]
HVGs.list[[condition]][["upregulated"]][[paste0(condition1,"_",condition2)]] <- rownames(res_filtered[res_filtered$log2FoldChange>=0,])
HVGs.list[[condition]][["downregulated"]][[paste0(condition1,"_",condition2)]] <- rownames(res_filtered[res_filtered$log2FoldChange<0,])
```
```{r}
# export the svg file
svglite::svglite(file = "volcano_plot_D7vsCtrl_FC_1.svg", width = 15, height = 15)
print(p1)
dev.off()
svglite::svglite(file = "volcano_plot_D7vsCtrl_FC_075.svg", width = 15, height = 15)
print(p1_1)
dev.off()
svglite::svglite(file = "volcano_plot_D14vsCtrl_FC_1.svg", width = 15, height = 15)
print(p2)
dev.off()
svglite::svglite(file = "volcano_plot_D14vsCtrl_FC_075.svg", width = 15, height = 15)
print(p2_1)
dev.off()
svglite::svglite(file = "volcano_plot_D14vsD7_FC_1.svg", width = 15, height = 15)
print(p3)
dev.off()
svglite::svglite(file = "volcano_plot_D14vsD7_FC_05_small.svg", width = 7, height = 7)
print(p3_1)
dev.off()
```
## GO on the res
```{r}
# load the different DEG list
res_D7_Ctrl <- read.csv(paste0("~/rds/rds-rds-karadottir-qqstKNzWvy8/Omar_Spatial/", "DGE_IO_calb1_mg_","D7","_","Ctrl","_20250915.csv"), row.names = 1)
res_D14_Ctrl <- read.csv(paste0("~/rds/rds-rds-karadottir-qqstKNzWvy8/Omar_Spatial/", "DGE_IO_calb1_mg_","D14","_","Ctrl","_20250915.csv"), row.names = 1)
res_D14_D7 <- read.csv(paste0("~/rds/rds-rds-karadottir-qqstKNzWvy8/Omar_Spatial/", "DGE_IO_calb1_mg_","D14","_","D7","_20250915.csv"), row.names = 1)
# clean up the na
res_D7_Ctrl <- res_D7_Ctrl[!is.na(res_D7_Ctrl$padj),]
res_D14_Ctrl <- res_D14_Ctrl[!is.na(res_D14_Ctrl$padj),]
res_D14_D7 <- res_D14_D7[!is.na(res_D14_D7$padj),]
# Filter the res for giving different lists
res_D7_Ctrl_FC05_up <- rownames(res_D7_Ctrl[(res_D7_Ctrl$padj < 0.05) & (res_D7_Ctrl$log2FoldChange > 0.5),])
res_D7_Ctrl_FC05_down <- rownames(res_D7_Ctrl[(res_D7_Ctrl$padj < 0.05) & (res_D7_Ctrl$log2FoldChange < -0.5),])
res_D14_Ctrl_FC05_up <- rownames(res_D14_Ctrl[(res_D14_Ctrl$padj < 0.05) & (res_D14_Ctrl$log2FoldChange > 0.5),])
res_D14_Ctrl_FC05_down <- rownames(res_D14_Ctrl[(res_D14_Ctrl$padj < 0.05) & (res_D14_Ctrl$log2FoldChange < -0.5),])
res_D14_D7_FC05_up <- rownames(res_D14_D7[(res_D14_D7$padj < 0.05) & (res_D14_D7$log2FoldChange > 0.5),])
res_D14_D7_FC05_down <- rownames(res_D14_D7[(res_D14_D7$padj < 0.05) & (res_D14_D7$log2FoldChange < -0.5),])
```
```{r}
dbs <- c("GO_Molecular_Function_2023","GO_Cellular_Component_2023","GO_Biological_Process_2023","KEGG_2019_Mouse")
```
```{r}
enrichedR.list <- list()
enrichedR.list[["IO_D7_Ctrl_FC05_up"]] <- enrichr(res_D7_Ctrl_FC05_up, dbs)
enrichedR.list[["IO_D7_Ctrl_FC05_down"]] <- enrichr(res_D7_Ctrl_FC05_down, dbs)
enrichedR.list[["IO_D14_Ctrl_FC05_up"]] <- enrichr(res_D14_Ctrl_FC05_up, dbs)
enrichedR.list[["IO_D14_Ctrl_FC05_down"]] <- enrichr(res_D14_Ctrl_FC05_down, dbs)
enrichedR.list[["IO_D14_D7_FC05_up"]] <- enrichr(res_D14_D7_FC05_up, dbs)
enrichedR.list[["IO_D14_D7_FC05_down"]] <- enrichr(res_D14_D7_FC05_down, dbs)
```
```{r}
# Create an empty list to store the data frames for export
export_list <- list()
# Loop through your enrichedR.list and extract a specific database for each cluster
for (cluster_name in names(enrichedR.list)) {
# For example, extract the "GO_Biological_Process_2023" results
export_list[[cluster_name]] <- enrichedR.list[[cluster_name]]$GO_Biological_Process_2023
}
openxlsx::write.xlsx(export_list, file = "enrichment_IO_calb1_microglia_DayComparison_GO_KEGG_20250915.xlsx")
```
```{r}
pdf("enrichment_plot_IO_calb1_mg_DayComparison_20250915.pdf", width=12, height=8)
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
p <- plotEnrich(
enr_df,
showTerms = 50,
numChar = 100,
y = "Count",
orderBy = "FDR",
title = paste0(dge," in ", go_db),
xlab = ""
)
print(
p + theme(
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 6),
plot.title = element_text(size = 10),
plot.margin = margin(10, 10, 10, 10)
)
)
}
}
dev.off()
```
```{r}
pdf("enrichment_plot_IO_calb1_mg_DayComparison_filtered_20250915.pdf", width=8, height=3)
# Define your list of keywords
keywords <- c("Positive Regulation Of Endocytosis", "Positive Regulation Of Peptidase Activity", "Phagocytosis, Engulfment", "Protein Processing", "Positive Regulation Of Receptor−Mediated Endocytosis")
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
# Filter for keywords and FDR < 0.1
filtered_df <- enr_df[
grepl(paste(keywords, collapse = "|"), enr_df$Term, ignore.case = TRUE) & enr_df$Adjusted.P.value < 0.1,
]
# Check if there are any results after filtering
if (nrow(filtered_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no terms found after filtering)")
next
}
# Use the filtered data for plotting
p <- plotEnrich(
filtered_df,
showTerms = 20,
numChar = 100,
y = "Count",
orderBy = "FDR",
title = paste0(dge, " in ", go_db, " (Filtered)"),
xlab = ""
)
print(
p + theme(
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 6),
plot.title = element_text(size = 10),
plot.margin = margin(10, 10, 10, 10)
)
)
}
}
dev.off()
```
# Comparison on Clusters
```{r}
library(Matrix)
library(dplyr)
# Grab metadata
meta <- colData(seurat_obj_sce_calb1_micro)[, c("sample_id", "seurat_clusters")]
# counts = genes × cells
cts <- counts(seurat_obj_sce_calb1_micro)
# build group labels: sample_cluster
groups <- paste(meta$sample_id, meta$seurat_clusters, sep = "_")
# aggregate counts across all cells belonging to each (sample, cluster) pair
pb <- aggregate.Matrix(t(cts), groupings = groups, fun = "sum")
pb <- t(pb) # back to genes × pseudobulks
dim(pb)
head(colnames(pb))
```
```{r}
# split sample_id and cluster from colnames
colinfo <- do.call(rbind, strsplit(colnames(pb), "_"))
coldata <- data.frame(
sample_id = colinfo[,1],
cluster = colinfo[,2],
row.names = colnames(pb)
)
head(coldata)
```
```{r}
dds <- DESeqDataSetFromMatrix(
countData = as.matrix(pb),
colData = coldata,
design = ~ sample_id + cluster
)
# filter low counts
dds <- dds[rowSums(counts(dds)) > 10, ]
# run DESeq2
dds <- DESeq(dds)
# results: cluster3 vs cluster0
res <- results(dds, contrast = c("cluster", "3", "0"))
res <- res[order(res$padj), ]
head(res)
```
```{r}
deg <- as.data.frame(res)
write.csv(deg, "DEG_calb1_mg_cluster3_vs_cluster0_pseudobulk_20250916.csv")
```
```{r}
# Volcano plot
library(EnhancedVolcano)
p1 <- EnhancedVolcano(
deg,
lab = rownames(deg),
x = 'log2FoldChange',
y = 'pvalue',
title = "Cluster 3 vs Cluster 0",
pCutoff = 0.01,
FCcutoff = 1,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1
)
print(p1)
```
```{r}
# load & filter
res_3_vs_0 <- read.csv("DEG_calb1_mg_cluster3_vs_cluster0_pseudobulk_20250916.csv", row.names = 1)
res_3_vs_0 <- res_3_vs_0[!is.na(res_3_vs_0$padj),]
# extract up/down
res_3_vs_0_FC1_up <- rownames(res_3_vs_0[(res_3_vs_0$padj < 0.05) & (res_3_vs_0$log2FoldChange > 1),])
res_3_vs_0_FC1_down <- rownames(res_3_vs_0[(res_3_vs_0$padj < 0.05) & (res_3_vs_0$log2FoldChange < -1),])
# FC > 0.75 (up and down)
res_3_vs_0_FC075_up <- rownames(res_3_vs_0[(res_3_vs_0$padj < 0.05) & (res_3_vs_0$log2FoldChange > 0.75),])
res_3_vs_0_FC075_down <- rownames(res_3_vs_0[(res_3_vs_0$padj < 0.05) & (res_3_vs_0$log2FoldChange < -0.75),])
```
```{r}
dbs <- c("GO_Molecular_Function_2023",
"GO_Cellular_Component_2023",
"GO_Biological_Process_2023",
"KEGG_2019_Mouse")
enrichedR.list <- list()
enrichedR.list[["Cluster3_vs_0_FC1_up"]] <- enrichr(res_3_vs_0_FC1_up, dbs)
enrichedR.list[["Cluster3_vs_0_FC1_down"]] <- enrichr(res_3_vs_0_FC1_down, dbs)
enrichedR.list[["Cluster3_vs_0_FC075_up"]] <- enrichr(res_3_vs_0_FC075_up, dbs)
enrichedR.list[["Cluster3_vs_0_FC075_down"]] <- enrichr(res_3_vs_0_FC075_down, dbs)
```
```{r}
pdf("enrichment_plot_calb1_mg_cluster3_vs_0_20250916.pdf", width=12, height=8)
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
p <- plotEnrich(
enr_df,
showTerms = 50,
numChar = 100,
y = "Count",
orderBy = "FDR",
title = paste0(dge," in ", go_db),
xlab = ""
)
print(
p + theme(
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 6),
plot.title = element_text(size = 10),
plot.margin = margin(10, 10, 10, 10)
)
)
}
}
dev.off()
```
```{r}
pdf("enrichment_dotplot_calb1_mg_cluster3_vs_0_20250916.pdf", width=12, height=8)
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
# ---- Prepare data ----
# Adjusted P value column in enrichR is usually "Adjusted.P.value"
if (!("Adjusted.P.value" %in% colnames(enr_df))) {
stop("No Adjusted.P.value column found in enrichment results for ", dge, " / ", go_db)
}
# Significance (FDR)
enr_df$FDR <- as.numeric(enr_df$Adjusted.P.value)
# Overlap is in "k/n" format -> extract numerator and denominator
tmp <- do.call(rbind, strsplit(enr_df$Overlap, "/"))
enr_df$Count <- as.numeric(tmp[,1])
enr_df$GeneRatio <- as.numeric(tmp[,1]) / as.numeric(tmp[,2])
# Keep top 50 terms by FDR
enr_df <- enr_df[order(enr_df$FDR), ][1:min(50, nrow(enr_df)), ]
# ---- Dotplot ----
p <- ggplot(enr_df, aes(
x = GeneRatio,
y = reorder(Term, -FDR), # most significant (lowest FDR) at top
size = Count,
color = -log10(FDR)
)) +
geom_point() +
scale_color_gradient(low = "blue", high = "red") +
labs(
title = paste0(dge," in ", go_db),
x = "Gene ratio",
y = "GO term",
color = "-log10(FDR)",
size = "Gene count"
) +
theme_bw() +
theme(
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 6),
plot.title = element_text(size = 10),
plot.margin = margin(10, 10, 10, 10)
)
print(p)
}
}
dev.off()
```
## Cluster 2 compared to cluster 0
```{r}
# results: cluster2 vs cluster0
res <- results(dds, contrast = c("cluster", "2", "0"))
res <- res[order(res$padj), ]
head(res)
```
```{r}
deg <- as.data.frame(res)
write.csv(deg, "DEG_calb1_mg_cluster2_vs_cluster0_pseudobulk.csv")
```
```{r}
# Volcano plot
library(EnhancedVolcano)
p1 <- EnhancedVolcano(
deg,
lab = rownames(deg),
x = 'log2FoldChange',
y = 'pvalue',
title = "Cluster 2 vs Cluster 0",
pCutoff = 0.05,
FCcutoff = 1,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1
)
print(p1)
```
```{r}
# load & filter
res_2_vs_0 <- read.csv("DEG_calb1_mg_cluster2_vs_cluster0_pseudobulk.csv", row.names = 1)
res_2_vs_0 <- res_2_vs_0[!is.na(res_2_vs_0$padj),]
# extract up/down
res_2_vs_0_FC1_up <- rownames(res_2_vs_0[(res_2_vs_0$padj < 0.05) & (res_2_vs_0$log2FoldChange > 1),])
res_2_vs_0_FC1_down <- rownames(res_2_vs_0[(res_2_vs_0$padj < 0.05) & (res_2_vs_0$log2FoldChange < -1),])
# FC > 0.75 (up and down)
res_2_vs_0_FC075_up <- rownames(res_2_vs_0[(res_2_vs_0$padj < 0.05) & (res_2_vs_0$log2FoldChange > 0.75),])
res_2_vs_0_FC075_down <- rownames(res_2_vs_0[(res_2_vs_0$padj < 0.05) & (res_2_vs_0$log2FoldChange < -0.75),])
# FC > 0.75 (up and down)
res_2_vs_0_FC05_up <- rownames(res_2_vs_0[(res_2_vs_0$padj < 0.05) & (res_2_vs_0$log2FoldChange > 0.5),])
res_2_vs_0_FC05_down <- rownames(res_2_vs_0[(res_2_vs_0$padj < 0.05) & (res_2_vs_0$log2FoldChange < -0.5),])
```
```{r}
dbs <- c("GO_Molecular_Function_2023",
"GO_Cellular_Component_2023",
"GO_Biological_Process_2023",
"KEGG_2019_Mouse")
enrichedR.list <- list()
enrichedR.list[["Cluster2_vs_0_FC1_up"]] <- enrichr(res_2_vs_0_FC1_up, dbs)
enrichedR.list[["Cluster2_vs_0_FC1_down"]] <- enrichr(res_2_vs_0_FC1_down, dbs)
enrichedR.list[["Cluster2_vs_0_FC075_up"]] <- enrichr(res_2_vs_0_FC075_up, dbs)
enrichedR.list[["Cluster2_vs_0_FC075_down"]] <- enrichr(res_2_vs_0_FC075_down, dbs)
enrichedR.list[["Cluster2_vs_0_FC05_up"]] <- enrichr(res_2_vs_0_FC05_up, dbs)
enrichedR.list[["Cluster2_vs_0_FC05_down"]] <- enrichr(res_2_vs_0_FC05_down, dbs)
```
```{r}
pdf("enrichment_plot_mg_cluster2_vs_0.pdf", width=12, height=8)
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
p <- plotEnrich(
enr_df,
showTerms = 50,
numChar = 100,
y = "Count",
orderBy = "FDR",
title = paste0(dge," in ", go_db),
xlab = ""
)
print(
p + theme(
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 6),
plot.title = element_text(size = 10),
plot.margin = margin(10, 10, 10, 10)
)
)
}
}
dev.off()
```
## Cluster 1 compared to cluster 0
```{r}
# results: cluster1 vs cluster0
res <- results(dds, contrast = c("cluster", "1", "0"))
res <- res[order(res$padj), ]
head(res)
```
```{r}
deg <- as.data.frame(res)
write.csv(deg, "DEG_mg_cluster1_vs_cluster0_pseudobulk.csv")
```
```{r}
# Volcano plot
library(EnhancedVolcano)
p1 <- EnhancedVolcano(
deg,
lab = rownames(deg),
x = 'log2FoldChange',
y = 'pvalue',
title = "Cluster 1 vs Cluster 0",
pCutoff = 0.01,
FCcutoff = 1,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1
)
print(p1)
```
```{r}
# load & filter
res_1_vs_0 <- read.csv("DEG_mg_cluster1_vs_cluster0_pseudobulk.csv", row.names = 1)
res_1_vs_0 <- res_1_vs_0[!is.na(res_1_vs_0$padj),]
# extract up/down
res_1_vs_0_FC1_up <- rownames(res_1_vs_0[(res_1_vs_0$padj < 0.05) & (res_1_vs_0$log2FoldChange > 1),])
res_1_vs_0_FC1_down <- rownames(res_1_vs_0[(res_1_vs_0$padj < 0.05) & (res_1_vs_0$log2FoldChange < -1),])
# FC > 0.75 (up and down)
res_1_vs_0_FC075_up <- rownames(res_1_vs_0[(res_1_vs_0$padj < 0.05) & (res_1_vs_0$log2FoldChange > 0.75),])
res_1_vs_0_FC075_down <- rownames(res_1_vs_0[(res_1_vs_0$padj < 0.05) & (res_1_vs_0$log2FoldChange < -0.75),])
```
```{r}
dbs <- c("GO_Molecular_Function_2023",
"GO_Cellular_Component_2023",
"GO_Biological_Process_2023",
"KEGG_2019_Mouse")
enrichedR.list <- list()
enrichedR.list[["Cluster1_vs_0_FC1_up"]] <- enrichr(res_1_vs_0_FC1_up, dbs)
enrichedR.list[["Cluster1_vs_0_FC1_down"]] <- enrichr(res_1_vs_0_FC1_down, dbs)
enrichedR.list[["Cluster1_vs_0_FC075_up"]] <- enrichr(res_1_vs_0_FC075_up, dbs)
enrichedR.list[["Cluster1_vs_0_FC075_down"]] <- enrichr(res_1_vs_0_FC075_down, dbs)
```
```{r}
pdf("enrichment_plot_mg_cluster1_vs_0.pdf", width=12, height=8)
for (dge in names(enrichedR.list)) {
for (go_db in dbs) {
enr_df <- enrichedR.list[[dge]][[go_db]]
if (is.null(enr_df) || nrow(enr_df) == 0) {
message("Skipping ", dge, " in ", go_db, " (no enrichment results)")
next
}
p <- plotEnrich(
enr_df,
showTerms = 50,
numChar = 100,
y = "Count",
orderBy = "FDR",
title = paste0(dge," in ", go_db),
xlab = ""
)
print(
p + theme(
axis.text.x = element_text(size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 6),
plot.title = element_text(size = 10),
plot.margin = margin(10, 10, 10, 10)
)
)
}
}
dev.off()
```
```{r}
# Save DEGs
deg <- as.data.frame(res)
write.csv(deg, "DEG_mg_cluster3_vs_cluster0.csv")
# Volcano plot
library(EnhancedVolcano)
p1 <- EnhancedVolcano(
deg,
lab = rownames(deg),
x = 'log2FoldChange',
y = 'pvalue',
title = "Cluster 3 vs Cluster 0",
pCutoff = 0.01,
FCcutoff = 1,
pointSize = 3.0,
labSize = 6.0,
col=c('black', 'grey', 'grey', 'blue'),
colAlpha = 1
)
print(p1)
```
```{r}
# Use the gene list for the GO
# load the 3 vs 0 DE results
res_3_vs_0 <- read.csv("DEG_mg_cluster3_vs_cluster0.csv", row.names = 1)
# clean up NA
res_3_vs_0 <- res_3_vs_0[!is.na(res_3_vs_0$padj),]
# FC > 1 (up and down)
res_3_vs_0_FC1_up <- rownames(res_3_vs_0[(res_3_vs_0$padj < 0.05) & (res_3_vs_0$log2FoldChange > 1),])
res_3_vs_0_FC1_down <- rownames(res_3_vs_0[(res_3_vs_0$padj < 0.05) & (res_3_vs_0$log2FoldChange < -1),])
# FC > 0.75 (up and down)
res_3_vs_0_FC075_up <- rownames(res_3_vs_0[(res_3_vs_0$padj < 0.05) & (res_3_vs_0$log2FoldChange > 0.75),])
res_3_vs_0_FC075_down <- rownames(res_3_vs_0[(res_3_vs_0$padj < 0.05) & (res_3_vs_0$log2FoldChange < -0.75),])
```
---------------------------------------------------------------
# plotting
## plotting all microglia pixels
```{r}
micro_Balazs_markers = c("Aif1", "Maf", "Spi1", "Csf1r")
seurat_obj_IO <- AddModuleScore(
seurat_obj_IO,
features=list(micro_Balazs_markers),
nbin = 24,
ctrl = 100,
assay = "RNA_repaired",
slot = "data",
name = "AMS",
)
colnames([email hidden]) <- c(colnames([email hidden])[1:(length(colnames([email hidden]))-1)], "micro_Balazs_markers")
options(repr.plot.width=6, repr.plot.height=5)
FeaturePlot(seurat_obj_IO, reduction = "umap", features = "micro_Balazs_markers", order=TRUE, pt.size=1, raster=FALSE, min.cutoff=0)
seurat_obj_IO_micro <- subset(seurat_obj_IO, subset= micro_Balazs_markers > -0.03)
```
```{r}
DefaultAssay(seurat_obj_IO_micro) <- "RNA"
seurat_obj_IO_micro <- NormalizeData(seurat_obj_IO_micro, normalization.method = "LogNormalize", scale.factor = 10000)
DefaultAssay(seurat_obj_IO_micro) <- "RNA_repaired"
seurat_obj_IO_micro <- NormalizeData(seurat_obj_IO_micro, normalization.method = "LogNormalize", scale.factor = 10000)
seurat_obj_IO_micro <- FindVariableFeatures(seurat_obj_IO_micro, selection.method = "vst", nfeatures = 1000)
seurat_obj_IO_micro <- ScaleData(seurat_obj_IO_micro, features = rownames(seurat_obj_IO_micro))
seurat_obj_IO_micro <- RunPCA(seurat_obj_IO_micro)
options(repr.plot.width=10, repr.plot.height=7)
ElbowPlot(seurat_obj_IO_micro, reduction = "pca", ndims = 50)
```
```{r}
nb_pcs = 20
seurat_obj_IO_micro <- FindNeighbors(seurat_obj_IO_micro, reduction = "pca", dims = 1:nb_pcs, prune.SNN = 0)
cluster_res = 0.8
seurat_obj_IO_micro <- FindClusters(seurat_obj_IO_micro, resolution = cluster_res)
table(seurat_obj_IO_micro$seurat_clusters)
```
```{r}
seurat_obj_IO_calb1_color <- list()
n_clusters <- 6
seurat_obj_IO_calb1_color <- color_palette_30[1:n_clusters]
names(seurat_obj_IO_calb1_color) <- seq(1:n_clusters)-1
seurat_obj_IO_micro <- RunUMAP(seurat_obj_IO_micro, reduction = "pca", dims = 1:nb_pcs)
```
```{r}
options(repr.plot.width=25, repr.plot.height=8)
DimPlot(seurat_obj_IO_micro, group.by="orig.ident", split.by="orig.ident_merge", cols=orig.ident_colors, shuffle=TRUE, seed=seed, ncol=4, pt.size=1, label=FALSE, raster=FALSE)
options(repr.plot.width=11, repr.plot.height=9)
for(timepoint in names(orig.ident_merge_colors)){
highlight_cells <- rownames([email hidden][[email hidden]$orig.ident_merge==timepoint,])
print(DimPlot(seurat_obj_IO_micro, cells.highlight = highlight_cells, pt.size=0.5, cols.highlight = orig.ident_merge_colors[timepoint], raster=FALSE))
}
```
```{r}
options(repr.plot.width=11, repr.plot.height=9)
DimPlot(seurat_obj_IO_micro, group.by="seurat_clusters", cols=seurat_obj_IO_calb1_color, shuffle=TRUE, seed=seed, pt.size=1, label=FALSE, raster=FALSE)
```
```{r}
options(repr.plot.width=25, repr.plot.height=8)
DimPlot(seurat_obj_IO_calb1, group.by="seurat_clusters", split.by="orig.ident_merge", cols=seurat_obj_IO_calb1_color, shuffle=TRUE, seed=seed, ncol=4, pt.size=1, label=FALSE, raster=FALSE)
```
```{r}
options(repr.plot.width=15, repr.plot.height=15)
set.seed(115)
mat = as.matrix(table(seurat_obj_IO_micro$orig.ident_merge, seurat_obj_IO_micro$seurat_clusters))
for (row in 1:nrow(mat)) mat[row,] <- as.integer(mat[row,]/(sum(mat[row,])/min(rowSums(mat))))
circos.clear()
par(cex = 0.8)
chordDiagram(mat, annotationTrack = c("grid", "axis"), preAllocateTracks = list(track.height = max(strwidth(unlist(dimnames(mat))))/3), big.gap = 20, small.gap = 2, order = c(colnames(mat), names(orig.ident_merge_colors)), grid.col=c(rev(seurat_obj_IO_calb1_color), rev(orig.ident_merge_colors)))
circos.track(track.index = 1, panel.fun = function(x, y) {circos.text(CELL_META$xcenter, CELL_META$ylim[1], CELL_META$sector.index, facing = "clockwise", niceFacing = TRUE, adj = c(-0.3, 0.5))}, bg.border = NA)
```
```{r}
df_table <- as.data.frame(table(seurat_obj_IO_micro$orig.ident, seurat_obj_IO_micro$seurat_clusters))
colnames(df_table) <- c("group", "cluster", "count")
t <- table(seurat_obj_IO_micro$orig.ident, seurat_obj_IO_micro$seurat_clusters)
t_prop <- t/as.vector(table(seurat_obj_IO$orig.ident))*100
rownames(t_prop) <- c("Ctrl","Ctrl","Ctrl","Ctrl","D7","D7","D7","D7","D14","D14","D14","D14")
df_table <- as.data.frame(t_prop)
colnames(df_table) <- c("group", "cluster", "count")
df_summary <- df_table %>%
group_by(group, cluster) %>%
summarise(
mean = mean(count),
se = sd(count) / sqrt(n()),
.groups = "drop"
)
ggplot(df_summary, aes(x = cluster, y = mean, fill = group)) +
# Bars
geom_bar(stat = "identity", position = position_dodge(width = 0.9), color = "black") +
# Error bars
geom_errorbar(aes(ymin = mean - se, ymax = mean + se),
position = position_dodge(width = 0.9),
width = 0.2) +
# Manual fill colors
scale_fill_manual(values = orig.ident_merge_colors) +
# Points dodged and jittered
geom_point(
data = df_table,
aes(x = cluster, y = count, fill = group), # <- fill included so dodge works
position = position_jitterdodge(jitter.width = 0.4, dodge.width = 0.9),
size = 2, shape = 21, stroke = 0.5,
color = "black", # border
show.legend = FALSE
) +
theme_minimal() +
labs(x = "", y = "% of pixels in the IO") +
ggtitle("Number of calbindin neurons pixels divided by each IO for each sample") +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank())
t <- table(Cluster=seurat_obj_IO_micro$seurat_clusters, Batch=[email hidden][["orig.ident_merge"]])
t <- t[,rev(names(orig.ident_merge_colors))]
options(repr.plot.width=10, repr.plot.height=10)
barplot(t(t), xlab = "Cluster", ylab="Number of cells", legend = TRUE, ylim = c(0, round_any(as.integer(max(rowSums(t)))+100, 20, f = ceiling)), col = rev(orig.ident_merge_colors), args.legend = list(bty = "n", x = "top", ncol = 3))
t <- table(Cluster=seurat_obj_IO_micro$seurat_clusters, Batch=[email hidden][["orig.ident"]])
t <- t[,rev(names(orig.ident_colors))]
options(repr.plot.width=10, repr.plot.height=10)
barplot(t(t), xlab = "Cluster", ylab="Number of cells", legend = TRUE, ylim = c(0, round_any(as.integer(max(rowSums(t)))+150, 20, f = ceiling)), col = rev(orig.ident_colors), args.legend = list(bty = "n", x = "top", ncol = 3))
```
```{r}
DefaultAssay(seurat_obj_IO_micro) <- "RNA"
Idents(object = seurat_obj_IO_micro) <- "seurat_clusters"
seurat_obj_IO_micro.markers.seurat_clusters <- FindAllMarkers(seurat_obj_IO_micro, slot = "data", min.pct = 0.1, logfc.threshold = 0, only.pos = TRUE)
```
```{r}
write.csv(seurat_obj_IO_micro.markers.seurat_clusters,"IO_microglia_clustermarkers_20250916.csv",row.names = FALSE)
```
```{r}
seurat_obj_IO_micro.markers.seurat_clusters_pval <- seurat_obj_IO_micro.markers.seurat_clusters[seurat_obj_IO_micro.markers.seurat_clusters$p_val_adj < 0.05,]
seurat_obj_IO_micro.markers.seurat_clusters_best <- seurat_obj_IO_micro.markers.seurat_clusters_pval %>%
filter(avg_log2FC > 0.5) %>%
filter(pct.1 > 0.3) %>%
arrange(cluster, desc(avg_log2FC)) %>%
group_by(cluster)
markers_heatmap <- seurat_obj_IO_micro.markers.seurat_clusters_best %>% top_n(n = 10, wt = avg_log2FC)
topDiffCluster <- seurat_obj_IO_micro.markers.seurat_clusters_best %>% top_n(n = 50, wt = avg_log2FC)
sample_df <- as.data.frame(lapply(split(topDiffCluster, topDiffCluster$cluster), function(x) c(x$gene,rep("None", 50-length(x$gene)))))
colnames(sample_df) <- levels(topDiffCluster$cluster)
sample_df
```
```{r}
options(repr.plot.width=30, repr.plot.height=20)
DoMultiBarHeatmap(
subset(seurat_obj_IO_micro, downsample = 200),
features = unique(markers_heatmap$gene),
cells = NULL,
group.by = "seurat_clusters",
additional.group.by = c("orig.ident_merge"),
additional.group.sort.by = c("orig.ident_merge"),
cols.use = list(seurat_clusters=seurat_obj_IO_calb1_color, orig.ident_merge=orig.ident_merge_colors),
group.bar = TRUE,
disp.min = -2.5,
disp.max = NULL,
layer = "scale.data",
assay = "RNA_repaired",
label = TRUE,
size = 5.5,
hjust = 0,
angle = 45,
raster = TRUE,
draw.lines = TRUE,
lines.width = NULL,
group.bar.height = 0.02,
combine = TRUE
) + theme(text = element_text(size = 8))
```
---------------------------------------------------------
```{r}
# load the 5 csv files
MC1_Volcano_df <- read.csv("goncaloGeneLists/MC1_filtered_genes.csv")
MC3_Volcano_df <- read.csv("goncaloGeneLists/MC3_filtered_genes.csv")
MC1_NonDE_df <- read.csv("goncaloGeneLists/MC1_Non_DE.csv")
MC2_NonDE_df <- read.csv("goncaloGeneLists/MC2_Non_DE.csv")
MC3_NonDE_df <- read.csv("goncaloGeneLists/MC3_Non_DE.csv")
```
```{r}
# Filter to make them all into lists
FDRCutOff <- 0.1
MC1_Volcano_df_fil <- MC1_Volcano_df[MC1_Volcano_df$FDR.x < FDRCutOff , ]
MC1_Volcano_list <- MC1_Volcano_df_fil$gene
MC3_Volcano_df_fil <- MC3_Volcano_df[MC3_Volcano_df$FDR.x < FDRCutOff , ]
MC3_Volcano_list <- MC3_Volcano_df_fil$gene
MC1_NonDE_df_fil <- MC1_NonDE_df[MC1_NonDE_df$FDR < FDRCutOff , ]
MC1_NonDE_df_list <- MC1_NonDE_df_fil$name
MC2_NonDE_df_fil <- MC2_NonDE_df[MC2_NonDE_df$FDR < FDRCutOff , ]
MC2_NonDE_df_list <- MC2_NonDE_df_fil$name
MC3_NonDE_df_fil <- MC3_NonDE_df[MC3_NonDE_df$FDR < FDRCutOff , ]
MC3_NonDE_df_list <- MC3_NonDE_df_fil$name
```
### load the Keren Shaul DAM file
```{r}
Keren_df <- read.csv("Keren-Shaul DAM list.csv")
Keren_up_list <- Keren_df$Keren_shaul_DAM_up
Keren_down_list <- Keren_df$Keren_shaul_DAM_down
Keren_down_list <- Keren_down_list[Keren_down_list!=""]
# create the neuroprotective list
NeuPro_list <- c("Igf1", "Bdnf", "Fgf13")
# Igf1
Igf1_list <- c("Igf1")
# create the additional neuroprotective list
NeuPro_Clus0_list <- c("Igf2", "Igfbp7", "Igf2r1", "Igf1r")
NeuPro_Clus2_list <- c("Bdnf", "Igf1", "Fbxw7", "Fgf13")
NeuPro_Clus02_list <- c(NeuPro_Clus0_list, NeuPro_Clus2_list)
```
```{r}
# calculate the gene score for the Keren up and down list
seurat_obj_IO_micro <- AddModuleScore(
object = seurat_obj_IO_micro,
features = list(Keren_up_list, Keren_down_list, NeuPro_list, MC1_Volcano_list, MC3_Volcano_list, NeuPro_Clus0_list, NeuPro_Clus2_list, NeuPro_Clus02_list),
name = c("KerenDAM", "KerenHomeo", "NeuPro", "MC1Score", "MC3Score", "NeuPro0Clus", "NeuPro2Clus", "NeuPro02Clus")
)
# Get normalized expression of Igf1 for all cells
seurat_obj_IO_micro$Igf1Ex4 <- FetchData(
seurat_obj_IO_micro,
vars = "Igf1"
)[,1]
```
```{r}
pdf("TissueMapping_IO_microglia_Loc_20250910.pdf", width = 5.3, height = 4.6)
# Make sure KerenDAM1 values exist in the full object
# Initialize with NA for all cells
seurat_obj_IO$Microglia <- NA
# Fill in values for microglia subset
seurat_obj_IO$Microglia[colnames(seurat_obj_IO_micro)] <- 1
# Now plot from the full seurat_obj_IO
for (sample_name in names(orig.ident_colors)) {
annot_plot <- [email hidden][seurat_obj_IO$orig.ident == sample_name, ]
options(repr.plot.width=5.3, repr.plot.height=4.6)
df_tmp <- data.frame(
x = annot_plot$x,
y = annot_plot$y,
value = annot_plot$Microglia
)
p <- ggplot(df_tmp, aes(x, y)) +
# Non-microglia (NA) cells plotted in lightgrey first
geom_point(data = subset(df_tmp, is.na(value)),
colour = "lightgrey", shape = 15, size = 1) +
# Microglia (with KerenDAM1 scores) plotted on top with gradient
geom_point(data = subset(df_tmp, !is.na(value)),
aes(colour = value), shape = 15, size = 1) +
scale_colour_gradientn(
colours = c("blue", "blue", "lightyellow", "red", "red"),
values = c(0, 0.25, 0.5, 0.75, 1),
limits = c(-0.1, 0.2),
na.value = "lightgrey",
oob = scales::squish
) +
scale_x_continuous(expand = c(-0.005, 0), lim = c(0, 101)) +
scale_y_continuous(expand = c(-0.006, 0), lim = c(0, 101)) +
theme_void() +
ggtitle(sample_name)
print(p)
}
dev.off()
```
```{r}
pdf("TissueMapping_IO_microglia_Loc_type2_20250910.pdf", width = 5.3, height = 4.6)
# Make sure KerenDAM1 values exist in the full object
# Initialize with NA for all cells
seurat_obj_IO$Microglia <- NA
# Fill in values for microglia subset
seurat_obj_IO$Microglia[colnames(seurat_obj_IO_mg_trial)] <- 1
# Now plot from the full seurat_obj_IO
for (sample_name in names(orig.ident_colors)) {
annot_plot <- [email hidden][seurat_obj_IO$orig.ident == sample_name, ]
options(repr.plot.width=5.3, repr.plot.height=4.6)
df_tmp <- data.frame(
x = annot_plot$x,
y = annot_plot$y,
value = annot_plot$Microglia
)
p <- ggplot(df_tmp, aes(x, y)) +
# Non-microglia (NA) cells plotted in lightgrey first
geom_point(data = subset(df_tmp, is.na(value)),
colour = "lightgrey", shape = 15, size = 1) +
# Microglia (with KerenDAM1 scores) plotted on top with gradient
geom_point(data = subset(df_tmp, !is.na(value)),
aes(colour = value), shape = 15, size = 1) +
scale_colour_gradientn(
colours = c("blue", "blue", "lightyellow", "red", "red"),
values = c(0, 0.25, 0.5, 0.75, 1),
limits = c(-0.1, 0.2),
na.value = "lightgrey",
oob = scales::squish
) +
scale_x_continuous(expand = c(-0.005, 0), lim = c(0, 101)) +
scale_y_continuous(expand = c(-0.006, 0), lim = c(0, 101)) +
theme_void() +
ggtitle(sample_name)
print(p)
}
dev.off()
```
```{r}
pdf("TissueMapping_IO_microglia_KerenDAM_20250910.pdf", width = 5.3, height = 4.6)
# Initialize with NA for all cells
seurat_obj_IO$KerenDAM1 <- NA
# Fill in values for microglia subset
seurat_obj_IO$KerenDAM1[colnames(seurat_obj_IO_micro)] <- seurat_obj_IO_micro$KerenDAM1
# Now plot from the full seurat_obj_IO
for (sample_name in names(orig.ident_colors)) {
annot_plot <- [email hidden][seurat_obj_IO$orig.ident == sample_name, ]
options(repr.plot.width=5.3, repr.plot.height=4.6)
df_tmp <- data.frame(
x = annot_plot$x,
y = annot_plot$y,
value = annot_plot$KerenDAM1
)
p <- ggplot(df_tmp, aes(x, y)) +
# Non-microglia (NA) cells plotted in lightgrey first
geom_point(data = subset(df_tmp, is.na(value)),
colour = "lightgrey", shape = 15, size = 1) +
# Microglia (with KerenDAM1 scores) plotted on top with gradient
geom_point(data = subset(df_tmp, !is.na(value)),
aes(colour = value), shape = 15, size = 1) +
scale_colour_gradientn(
colours = c("blue", "blue", "lightyellow", "red", "red"),
values = c(0, 0.25, 0.5, 0.75, 1),
limits = c(-0.1, 0.25),
na.value = "lightgrey",
oob = scales::squish
) +
scale_x_continuous(expand = c(-0.005, 0), lim = c(0, 101)) +
scale_y_continuous(expand = c(-0.006, 0), lim = c(0, 101)) +
theme_void() +
ggtitle(sample_name)
print(p)
}
dev.off()
```
```{r}
# pdf("TissueMapping_IO_microglia_KerenHomeo_20250910.pdf", width = 5.3, height = 4.6)
# Initialize with NA for all cells
seurat_obj_IO$KerenHomeo2 <- NA
# Fill in values for microglia subset
seurat_obj_IO$KerenHomeo2[colnames(seurat_obj_IO_micro)] <- seurat_obj_IO_micro$KerenHomeo2
# Now plot from the full seurat_obj_IO
for (sample_name in names(orig.ident_colors)) {
annot_plot <- [email hidden][seurat_obj_IO$orig.ident == sample_name, ]
options(repr.plot.width=5.3, repr.plot.height=4.6)
df_tmp <- data.frame(
x = annot_plot$x,
y = annot_plot$y,
value = annot_plot$KerenHomeo2
)
p <- ggplot(df_tmp, aes(x, y)) +
# Non-microglia (NA) cells plotted in lightgrey first
geom_point(data = subset(df_tmp, is.na(value)),
colour = "lightgrey", shape = 15, size = 1) +
# Microglia (with KerenDAM1 scores) plotted on top with gradient
geom_point(data = subset(df_tmp, !is.na(value)),
aes(colour = value), shape = 15, size = 1) +
scale_colour_gradientn(
colours = c("darkgreen", "yellow", "red"),
values = c(0,0.5, 1),
limits = c(-0.2, 0.3),
na.value = "lightgrey",
oob = scales::squish
) +
scale_x_continuous(expand = c(-0.005, 0), lim = c(0, 101)) +
scale_y_continuous(expand = c(-0.006, 0), lim = c(0, 101)) +
theme_void() +
ggtitle(sample_name)
print(p)
}
# dev.off()
```
```{r}
pdf("TissueMapping_IO_microglia_NeuProtective_20250910.pdf", width = 5.3, height = 4.6)
# Initialize with NA for all cells
seurat_obj_IO$NeuPro3 <- NA
# Fill in values for microglia subset
seurat_obj_IO$NeuPro3[colnames(seurat_obj_IO_micro)] <- seurat_obj_IO_micro$NeuPro3
# Now plot from the full seurat_obj_IO
for (sample_name in names(orig.ident_colors)) {
annot_plot <- [email hidden][seurat_obj_IO$orig.ident == sample_name, ]
options(repr.plot.width=5.3, repr.plot.height=4.6)
df_tmp <- data.frame(
x = annot_plot$x,
y = annot_plot$y,
value = annot_plot$NeuPro3
)
p <- ggplot(df_tmp, aes(x, y)) +
# Non-microglia (NA) cells plotted in lightgrey first
geom_point(data = subset(df_tmp, is.na(value)),
colour = "lightgrey", shape = 15, size = 1) +
# Microglia (with KerenDAM1 scores) plotted on top with gradient
geom_point(data = subset(df_tmp, !is.na(value)),
aes(colour = value), shape = 15, size = 1) +
scale_colour_gradientn(
colours = c("blue", "lightyellow", "red"),
values = c(0, 0.5, 1),
limits = c(-0.5, 0.5),
na.value = "lightgrey",
oob = scales::squish
) +
scale_x_continuous(expand = c(-0.005, 0), lim = c(0, 101)) +
scale_y_continuous(expand = c(-0.006, 0), lim = c(0, 101)) +
theme_void() +
ggtitle(sample_name)
print(p)
}
dev.off()
```
```{r}
pdf("TissueMapping_IO_microglia_MC1_20250910.pdf", width = 5.3, height = 4.6)
# Initialize with NA for all cells
seurat_obj_IO$MC1Score4 <- NA
# Fill in values for microglia subset
seurat_obj_IO$MC1Score4[colnames(seurat_obj_IO_micro)] <- seurat_obj_IO_micro$MC1Score4
# Now plot from the full seurat_obj_IO
for (sample_name in names(orig.ident_colors)) {
annot_plot <- [email hidden][seurat_obj_IO$orig.ident == sample_name, ]
options(repr.plot.width=5.3, repr.plot.height=4.6)
df_tmp <- data.frame(
x = annot_plot$x,
y = annot_plot$y,
value = annot_plot$MC1Score4
)
p <- ggplot(df_tmp, aes(x, y)) +
# Non-microglia (NA) cells plotted in lightgrey first
geom_point(data = subset(df_tmp, is.na(value)),
colour = "lightgrey", shape = 15, size = 1) +
# Microglia (with KerenDAM1 scores) plotted on top with gradient
geom_point(data = subset(df_tmp, !is.na(value)),
aes(colour = value), shape = 15, size = 1) +
scale_colour_gradientn(
colours = c("blue", "yellow", "red"),
values = c(0, 0.5, 1),
limits = c(0, 0.8),
na.value = "lightgrey",
oob = scales::squish
) +
scale_x_continuous(expand = c(-0.005, 0), lim = c(0, 101)) +
scale_y_continuous(expand = c(-0.006, 0), lim = c(0, 101)) +
theme_void() +
ggtitle(sample_name)
print(p)
}
dev.off()
```
```{r}
pdf("TissueMapping_IO_microglia_MC3_20250910.pdf", width = 5.3, height = 4.6)
# Initialize with NA for all cells
seurat_obj_IO$MC3Score5 <- NA
# Fill in values for microglia subset
seurat_obj_IO$MC3Score5[colnames(seurat_obj_IO_micro)] <- seurat_obj_IO_micro$MC3Score5
# Now plot from the full seurat_obj_IO
for (sample_name in names(orig.ident_colors)) {
annot_plot <- [email hidden][seurat_obj_IO$orig.ident == sample_name, ]
options(repr.plot.width=5.3, repr.plot.height=4.6)
df_tmp <- data.frame(
x = annot_plot$x,
y = annot_plot$y,
value = annot_plot$MC3Score5
)
p <- ggplot(df_tmp, aes(x, y)) +
# Non-microglia (NA) cells plotted in lightgrey first
geom_point(data = subset(df_tmp, is.na(value)),
colour = "lightgrey", shape = 15, size = 1) +
# Microglia (with KerenDAM1 scores) plotted on top with gradient
geom_point(data = subset(df_tmp, !is.na(value)),
aes(colour = value), shape = 15, size = 1) +
scale_colour_gradientn(
colours = c("blue", "yellow", "red"),
values = c(0, 0.5, 1),
limits = c(-0.05, 0.3),
na.value = "lightgrey",
oob = scales::squish
) +
scale_x_continuous(expand = c(-0.005, 0), lim = c(0, 101)) +
scale_y_continuous(expand = c(-0.006, 0), lim = c(0, 101)) +
theme_void() +
ggtitle(sample_name)
print(p)
}
dev.off()
```
```{r}
pdf("TissueMapping_IO_microglia_Igf1Exp_20250910.pdf", width = 5.3, height = 4.6)
# Initialize with NA for all cells
seurat_obj_IO$Igf1Ex4 <- NA
# Fill in values for microglia subset
seurat_obj_IO$Igf1Ex4[colnames(seurat_obj_IO_micro)] <- seurat_obj_IO_micro$Igf1Ex4
# Now plot from the full seurat_obj_IO
for (sample_name in names(orig.ident_colors)) {
annot_plot <- [email hidden][seurat_obj_IO$orig.ident == sample_name, ]
options(repr.plot.width=5.3, repr.plot.height=4.6)
df_tmp <- data.frame(
x = annot_plot$x,
y = annot_plot$y,
value = annot_plot$Igf1Ex4
)
p <- ggplot(df_tmp, aes(x, y)) +
# Non-microglia (NA) cells plotted in lightgrey first
geom_point(data = subset(df_tmp, is.na(value)),
colour = "lightgrey", shape = 15, size = 1) +
# Microglia (with KerenDAM1 scores) plotted on top with gradient
geom_point(data = subset(df_tmp, !is.na(value)),
aes(colour = value), shape = 15, size = 1) +
scale_colour_gradientn(
colours = c("lightyellow","red", "purple", "purple"),
values = c(0, 0.5, 0.95, 1),
# limits = c(0, 2),
na.value = "lightgrey",
oob = scales::squish
) +
scale_x_continuous(expand = c(-0.005, 0), lim = c(0, 101)) +
scale_y_continuous(expand = c(-0.006, 0), lim = c(0, 101)) +
theme_void() +
ggtitle(sample_name)
print(p)
}
# dev.off()
```
## cluster 0, popular in D14
```{r}
pdf("TissueMapping_IO_microglia_NeuProtectiveClus0_20250910.pdf", width = 5.3, height = 4.6)
# Initialize with NA for all cells
seurat_obj_IO$NeuPro0Clus6 <- NA
# Fill in values for microglia subset
seurat_obj_IO$NeuPro0Clu

ST_Analysis_IO_20250929_final.Rmd at commit 75f7289, no license · at the source

Overview

Authors: Omar de Faria Jr¹, Stavros Vagionitis¹, Andrea Lopez-Lopez^1,2, Michael Perry^1,3, Joseph Jo Yin Wong¹, Leslie Rodríguez-Kirby⁴, Bastien Hervé⁴, Balazs Viktor Varga¹, Eneritz Agirre⁴, Sabrina Ghosh¹, Sebastian Timmler¹, Mert Yucel¹, Andrew T. Setley⁵, Kimberley Anne Evans¹, Tanja Mist Birgisdóttir⁶, Sindri Gíslason⁶, Yan Ting Ng¹, Courtney Kremler¹, Helene O. B. Gautier¹, Yasmine Kamen¹

and 10 other authors

Helena Pivonkova¹, Katrin Volbracht¹, Felix Hildebrand¹, Christian A. Cepeda¹, Javier Rueda-Carrasco⁷, Soyon Hong⁷, George Malliaras⁵, Sabine Dietmann⁸, Gonçalo Castelo-Branco⁴, Ragnhildur Thóra Káradóttir^1,6

ORCID iDs: Omar de Faria Jr, Stavros Vagionitis, Andrea Lopez-Lopez, Bastien Hervé, Balazs Viktor Varga, Eneritz Agirre, Sebastian Timmler, Andrew T. Setley, Yan Ting Ng, Yasmine Kamen, Helena Pivonkova, Katrin Volbracht, Javier Rueda-Carrasco, Soyon Hong, George Malliaras, Gonçalo Castelo-Branco, Ragnhildur Thóra Káradóttir

Institutions: University of Cambridge (United Kingdom); Universidade de Santiago de Compostela (Spain); University College London (United Kingdom); Karolinska Institutet (Sweden); University of Iceland (Iceland); Washington University in St. Louis (United States)

Journal: Nature, volume 654, issue 8120, pages 1033-1043

Dates: received 16 August 2024; accepted 13 March 2026; published online 22 April 2026; in print 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1038/s41586-026-10414-w · PMID 42020752 · PMCID PMC13293868 · OpenAlex W7155174467

Open access: hybrid, a free copy (OpenAlex)

Status: code verified

Categories: structural MRI / diffusion (modality), mouse (organism), rat (organism), other condition (population), Alzheimer's / dementia (population), multiple sclerosis (population), cellular / molecular (subfield)

Methods: Connectivity, Spectral & time-frequency, Statistics, Smoothing, state filtering, decompositions, Preprocessing, Evoked potentials, Machine learning, fMRI & imaging, Single-unit activity, calcium imaging

Keywords: Multiple sclerosis, Alzheimer's disease, Neuroimmunology

MeSH: Gray Matter*, Inflammation*, Neuroinflammatory Diseases*, Synapses*, White Matter*, Animals, Female, Male, Mice, Microglia, Myelin Sheath, Neurodegenerative Diseases, Rats, Regeneration (* major topic)

Topic: Neuroinflammation and Neurodegeneration Mechanisms (Neurology, Neuroscience), according to OpenAlex

Funding: Multiple Sclerosis Society (207, 132); Wellcome Trust (218481/Z/19/Z, 320224/Z/24/Z)

Citations: cited by 2 papers (Europe PMC); 79 references in the paper

Abstract

Focal white matter lesions occur in most neurodegenerative disorders1–3. Despite occurring early in disease, white matter lesions are considered to be independent of, or secondary to, grey matter neuroinflammation, synapse loss and altered neuronal activity4–7. Notably, their functional effect on neuronal circuits remains understudied. To address this, we generated a focal white matter lesion in the rat brain within a clinically relevant, anatomically well-defined circuit, in which these lesions occur in many neurodegenerative disorders8–10. Here we show that focal white matter lesions evoke transient neuronal activity changes and microgliosis, with subsequent synapse loss and increased microglial engulfment in the grey matter, which is reversed if myelin regeneration completes. Grey matter microgliosis is often considered to be detrimental; however, we show that it is an integral part of regeneration and is conserved across three distinct mouse circuits and lesioning methods. Preventing these transient changes in the grey matter blocks myelin regeneration in the white matter. Conversely, inducing myelin regeneration failure leads to chronic grey matter neuroinflammation. This recapitulates the low-grade inflammation considered to be a dominant mechanism underlying neurodegeneration7,11,12. Our findings reveal a form of regenerative plasticity coupling white matter integrity to grey matter function, which may underlie multiple neurodegenerative conditions, and highlight the potential of targeting myelin regeneration to prevent chronic neuroinflammation.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 10 matches between paragraphs and lines of code.

Castelo-Branco-lab/Karadottir_DBiT_2025

License: none: the authors keep all their rights

State: the link answers, verified on 29 September 2026

Evidence: files inventoried

Commit: 75f72895e43d0bdfc73de65e3b25d1fe026e7404, 29 September 2025

Languages: R (5)

Size: 7 files, 5 scripts

Software Heritage: not archived

Found in: “Code availability”

Holds: README, 5 notebooks

Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: tidyverse (5 files), Seurat (4 files), circlize (3 files), DESeq2 (3 files), ggplot2 (3 files), igraph (3 files), patchwork (3 files), pheatmap (3 files), reticulate (3 files), ComplexHeatmap (2 files), cowplot (2 files), data.table (1 file), Harmony (1 file), SingleCellExperiment (1 file)

Availability: 1 check, the latest on 29 September 2026: the link answers

29 September 2026: the link answers

6 files

Bulk_deconvolution/Bulk_RNor2Mm_deconv_2025_EA_1.Rmd, R, 1,024 lines, 2 matches
Bulk_deconvolution/Bulk_RNor2Mm_deconv_markerlists_2025_EA_2.Rmd, R, 583 lines
DBiTseq_analysis/ST_Analysis_IO_20250929_final.Rmd, R, 5,962 lines, 4 matches
DBiTseq_analysis/ST_Analysis_Lesion_20250929.Rmd, R, 1,913 lines, 3 matches
DBiTseq_analysis/ST_analysis_Healthy.rmd, R, 1,531 lines, 1 match
README.md, Text, 1 line

At the source: github.com/Castelo-Branco-lab/Karadottir_DBiT_2025

Code availability

The codes used for RNA-seq deconvolution and DBit-seq analysis are available at GitHub (https://github.com/Castelo-Branco-lab/Karadottir_DBiT_2025).

Reproduced under the paper's license (CC BY), from the paper cited above.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
5 scripts, each with its path and the digest of its content;
10 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

geo:GSE311781, at NCBI GEO; found in “Data availability”

Data availability

All sequencing data have been deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO series accession number GSE274050 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE274050) (bulk RNA data) and series accession number GSE311781 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE311781) (DBiT-seq). The R. norvegicus reference genome used to align bulk RNA-seq and DBiT-seq reads (mRatBN7.2, RefSeq annotation, NCBI: GCF_015227675.2 (https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_015227675.2)) is available online. The scRNA-seq dataset used to deconvolute the bulk RNA-seq experiment is available online (https://cells.ucsc.edu/?ds=mouse-nervous-system). Source data are provided with this paper.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 29 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 30 authors, 3 keywords, 14 MeSH terms, 2 funders, 77 references.

Cite

This paper

de Faria, O., Vagionitis, S., Lopez-Lopez, A., Perry, M., Wong, J. J. Y., Rodríguez-Kirby, L., Hervé, B., Varga, B. V., Agirre, E., Ghosh, S., Timmler, S., Yucel, M., Setley, A. T., Evans, K. A., Birgisdóttir, T. M., Gíslason, S., Ng, Y. T., Kremler, C., Gautier, H. O. B., . . . Káradóttir, R. T. (2026). Focal white matter lesions drive grey matter inflammation and synapse loss. Nature, 654(8120), 1033-1043. https://doi.org/10.1038/s41586-026-10414-w

BibTeX

@article{defaria2026focal, author = {de Faria, Omar and Vagionitis, Stavros and Lopez-Lopez, Andrea and Perry, Michael and Wong, Joseph Jo Yin and Rodríguez-Kirby, Leslie and Hervé, Bastien and Varga, Balazs Viktor and Agirre, Eneritz and Ghosh, Sabrina and Timmler, Sebastian and Yucel, Mert and Setley, Andrew T. and Evans, Kimberley Anne and Birgisdóttir, Tanja Mist and Gíslason, Sindri and Ng, Yan Ting and Kremler, Courtney and Gautier, Helene O. B. and Kamen, Yasmine and Pivonkova, Helena and Volbracht, Katrin and Hildebrand, Felix and Cepeda, Christian A. and Rueda-Carrasco, Javier and Hong, Soyon and Malliaras, George and Dietmann, Sabine and Castelo-Branco, Gonçalo and Káradóttir, Ragnhildur Thóra}, title = {{Focal white matter lesions drive grey matter inflammation and synapse loss}}, journal = {Nature}, year = {2026}, month = apr, volume = {654}, number = {8120}, pages = {1033--1043}, publisher = {Nature Portfolio}, issn = {0028-0836}, doi = {10.1038/s41586-026-10414-w}, url = {https://doi.org/10.1038/s41586-026-10414-w}, pmid = {42020752}, pmcid = {PMC13293868} }

RIS

TY - JOUR AU - de Faria, Omar AU - Vagionitis, Stavros AU - Lopez-Lopez, Andrea AU - Perry, Michael AU - Wong, Joseph Jo Yin AU - Rodríguez-Kirby, Leslie AU - Hervé, Bastien AU - Varga, Balazs Viktor AU - Agirre, Eneritz AU - Ghosh, Sabrina AU - Timmler, Sebastian AU - Yucel, Mert AU - Setley, Andrew T. AU - Evans, Kimberley Anne AU - Birgisdóttir, Tanja Mist AU - Gíslason, Sindri AU - Ng, Yan Ting AU - Kremler, Courtney AU - Gautier, Helene O. B. AU - Kamen, Yasmine AU - Pivonkova, Helena AU - Volbracht, Katrin AU - Hildebrand, Felix AU - Cepeda, Christian A. AU - Rueda-Carrasco, Javier AU - Hong, Soyon AU - Malliaras, George AU - Dietmann, Sabine AU - Castelo-Branco, Gonçalo AU - Káradóttir, Ragnhildur Thóra TI - Focal white matter lesions drive grey matter inflammation and synapse loss T2 - Nature J2 - Nature PY - 2026 DA - 2026/04/22 VL - 654 IS - 8120 SP - 1033 EP - 1043 SN - 0028-0836 PB - Nature Portfolio DO - 10.1038/s41586-026-10414-w UR - https://doi.org/10.1038/s41586-026-10414-w LA - en ER -

CSL-JSON

{ "id": "10.1038/s41586-026-10414-w", "type": "article-journal", "title": "Focal white matter lesions drive grey matter inflammation and synapse loss", "container-title": "Nature", "author": [ { "family": "de Faria", "given": "Omar" }, { "family": "Vagionitis", "given": "Stavros" }, { "family": "Lopez-Lopez", "given": "Andrea" }, { "family": "Perry", "given": "Michael" }, { "family": "Wong", "given": "Joseph Jo Yin" }, { "family": "Rodríguez-Kirby", "given": "Leslie" }, { "family": "Hervé", "given": "Bastien" }, { "family": "Varga", "given": "Balazs Viktor" }, { "family": "Agirre", "given": "Eneritz" }, { "family": "Ghosh", "given": "Sabrina" }, { "family": "Timmler", "given": "Sebastian" }, { "family": "Yucel", "given": "Mert" }, { "family": "Setley", "given": "Andrew T." }, { "family": "Evans", "given": "Kimberley Anne" }, { "family": "Birgisdóttir", "given": "Tanja Mist" }, { "family": "Gíslason", "given": "Sindri" }, { "family": "Ng", "given": "Yan Ting" }, { "family": "Kremler", "given": "Courtney" }, { "family": "Gautier", "given": "Helene O. B." }, { "family": "Kamen", "given": "Yasmine" }, { "family": "Pivonkova", "given": "Helena" }, { "family": "Volbracht", "given": "Katrin" }, { "family": "Hildebrand", "given": "Felix" }, { "family": "Cepeda", "given": "Christian A." }, { "family": "Rueda-Carrasco", "given": "Javier" }, { "family": "Hong", "given": "Soyon" }, { "family": "Malliaras", "given": "George" }, { "family": "Dietmann", "given": "Sabine" }, { "family": "Castelo-Branco", "given": "Gonçalo" }, { "family": "Káradóttir", "given": "Ragnhildur Thóra" } ], "container-title-short": "Nature", "volume": "654", "issue": "8120", "page": "1033-1043", "DOI": "10.1038/s41586-026-10414-w", "PMID": "42020752", "PMCID": "PMC13293868", "ISSN": "0028-0836", "publisher": "Nature Portfolio", "URL": "https://doi.org/10.1038/s41586-026-10414-w", "language": "en", "issued": { "date-parts": [ [ 2026, 4, 22 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41593-026-02367-0 [code]: A reproducible three-dimensional model of human brain tissue to investigate physiological and disease-associated microglia phenotypes.
Journal: Nature neuroscience
In common: Harmony, SingleCellExperiment, reticulate, 11 other tools, Alzheimer's / dementia, cellular / molecular, 1 reference
[2] doi:10.1016/j.celrep.2026.117500 [code]: Spatio-molecular gene expression reflects dorsal anterior cingulate cortex structure and function in the human brain.
Journal: Cell reports
In common: Harmony, SingleCellExperiment, reticulate, 9 other tools, cellular / molecular, 3 references
[3] doi:10.1038/s41586-026-10214-2 [code]: Multidimensional profiling of heterogeneity in supratentorial ependymomas.
Journal: Nature
In common: Harmony, SingleCellExperiment, reticulate, 11 other tools, other condition, mouse
[4] doi:10.1002/imt2.70163 [code]: Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: Harmony, SingleCellExperiment, igraph, 10 other tools, mouse, cellular / molecular, 1 reference
[5] doi:10.1016/j.xcrm.2026.102766 [code]: A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: Harmony, SingleCellExperiment, igraph, 10 other tools, other condition, cellular / molecular, 1 reference
[6] doi:10.1016/j.xcrm.2026.102651 [code]: Integrative CSF profiling identifies disease-specific immune responses in leptomeningeal disease.
Journal: Cell reports. Medicine
In common: Harmony, SingleCellExperiment, reticulate, 10 other tools, other condition, cellular / molecular
[7] doi:10.1038/s41593-026-02384-z [code]: cGAS-mediated type I IFN signaling contributes to disease progression in drug-refractory epilepsy.
Journal: Nature neuroscience
In common: SingleCellExperiment, igraph, circlize, 8 other tools, mouse, cellular / molecular, 3 references
[8] doi:10.1016/j.isci.2026.115573 [code]: Female cortical cellular mosaicism underlies shared MeCP2 and PCB impacted gene pathways.
Journal: iScience
In common: SingleCellExperiment, igraph, circlize, 9 other tools, mouse, cellular / molecular, 1 reference
[9] doi:10.1038/s41586-026-10629-x [code]: Whole-genome duplication shaped cell-type evolution in the vertebrate brain.
Journal: Nature
In common: Harmony, reticulate, igraph, 8 other tools, mouse, cellular / molecular, 2 references
[10] doi:10.1038/s41514-026-00397-3 [code]: Nasal administration of Protollin enhances monocyte phagocytosis and decreases CD8<sup>+</sup> T cell cytotoxicity in subjects with early Alzheimer's disease: a Phase 1 clinical trial.
Journal: npj aging
In common: Harmony, reticulate, circlize, 8 other tools, Alzheimer's / dementia, 2 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.

Focal white matter lesions drive grey matter inflammation and synapse loss.

Code ↔ Paper

Paper

The authors' code

Overview

Abstract

Repository

Castelo-Branco-lab/Karadottir_DBiT_2025

Code availability

Tracing map

Data

Datasets cited

Data availability

Versions

Version 1, 29 September 2026: the first record

Cite

This paper

Similar papers

Contribute

Claim this paper

Correct its record

Validate its tracing map

Add the badge to its README

Request its removal

Discussion, reproductions, activity