Ischemic injury triggers a protective microglial phenotype in models of Aβ pathology.

Code ↔ Paper

3 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 3 matches

[1] § Results › Augmentation of Apoe expression is a defining feature of microglial stroke response in comorbidity ↔ analysis/scRNA_Analyses_Candlishetal.Rmd, lines 416–441 · score 0.93 · Cx3cr1, P2ry12, H2 Aa, H2 Ab1, Clec7a, response microglia
[2] § Results › Augmentation of Apoe expression is a defining feature of microglial stroke response in comorbidity ↔ analysis/scRNA_Analyses_Candlishetal.Rmd, lines 416–441 · score 0.80 · H2 Aa, H2 Ab1, Clec7a, response microglia, Oasl2, Ifit3
[3] § Results › Augmentation of Apoe expression is a defining feature of microglial stroke response in comorbidity ↔ analysis/scRNA_Analyses_Candlishetal.Rmd, lines 3766–3820 · score 0.58 · H2 D1, Siglech, Ccl6, B2m, Lipe, Cst7

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 · 3,920 lines · 127 KB · no license · 3 matches

---
title: "Analyses for Publication"
author: "AGC"
date: "2023-09-27"
output:
#pdf_document: default
html_document: default
---
```{r setup, include=FALSE}
#options(Seurat.object.assay.version = "v5")
outputdir="./output/Res_202602/"
if(!dir.exists(outputdir)) dir.create(outputdir, recursive = T)
## don't update matrixStats
knitr::opts_chunk$set(
echo = FALSE,
warning = FALSE,
message = FALSE)
#library(limma)
#library(DESeq2)
library(tidyverse)
library(scCustomize)
library(Seurat)
library(harmony)
library(SingleR)
library(data.table)
library(compareGroups)
library(openxlsx)
library(viridis)
library(DT)
library(slingshot)
library(SingleCellExperiment)
library(ggpubr)
library(clustree)
library(RCurl)
library(AnnotationHub)
library(WGCNA)
library(plotly)
library(org.Mm.eg.db)
library(knitr)
library(ggrepel)
set.seed(1324567)
getwd()
```
```{r load_annotations, cache=F, include=FALSE}
source ("./code/custom_functions.R")
cc_file <- getURL("https://raw.githubusercontent.com/hbc/tinyatlas/master/cell_cycle/Mus_musculus.csv")
cell_cycle_genes <- read.csv(text = cc_file)
ah <- AnnotationHub::AnnotationHub()
# Access the Ensembl database for organism
ahDb <- query(ah,
pattern = c("Mus musculus", "EnsDb"),
ignore.case = TRUE)
# Acquire the latest annotation files
id <- ahDb %>%
mcols() %>%
rownames() %>%
tail(n = 1)
# Download the appropriate Ensembldb database
edb <- ah[[id]]
# Extract gene-level information from database
annotations <- genes(edb,
return.type = "data.frame")
# Select annotations of interest
annotations <- annotations %>%
dplyr::select(gene_id, gene_name, seq_name, gene_biotype, description)
```
# Load batch 2026
## Metadata additional samples
```{r load_add_meta, cache=F}
# Read Sheet2
df_raw <- read.xlsx("/files/scSeq_Hefendehl/data/20260126/534_plates_wSexAge.xlsx", sheet = "Sheet2")
# Clean column names
colnames(df_raw) <- c("User_Sample_Name", "Plate", "Container_Type",
"Well_Location", "Protocol", "Species", "Mouse_ID",
"Tissue_or_Organ", "Cell_Type", "Genotype", "Treatment",
"Gender", "Age")
# Filter to actual sample entries (remove header rows / template rows)
df <- df_raw %>%
dplyr::filter(!is.na(User_Sample_Name),
!grepl("^<", User_Sample_Name), # remove <TABLE HEADER> etc.
Container_Type == "386 well plate") # keep only plate entries
# Function to expand a well range like "A1 - P12" into all individual wells
expand_wells <- function(well_range) {
# Parse "A1 - P12" format
parts <- trimws(strsplit(well_range, "-")[[1]])
start <- trimws(parts[1])
end <- trimws(parts[2])
# Extract row letters and column numbers
start_row <- substr(start, 1, 1)
start_col <- as.integer(substring(start, 2))
end_row <- substr(end, 1, 1)
end_col <- as.integer(substring(end, 2))
rows <- LETTERS[which(LETTERS == start_row):which(LETTERS == end_row)]
cols <- start_col:end_col
# Generate all combinations
expand.grid(Row = rows, Col = cols, stringsAsFactors = FALSE) %>%
mutate(Well = paste0(Row, Col)) %>%
pull(Well)
}
# Expand each sample entry to individual well positions
df_expanded <- df %>%
rowwise() %>%
mutate(Well = list(expand_wells(Well_Location))) %>%
unnest(Well) %>%
dplyr::select(Plate,
Well,
Species,
Mouse_ID,
Tissue_or_Organ,
Cell_Type,
Genotype,
Treatment,
Age, Gender)
df_expanded <- df_expanded %>%
mutate(
Plate = gsub("_011", "_11", Plate),
Plate = gsub("_012", "_12", Plate),
Cell_ID = paste0(Plate, "_", Well),
Gender = ifelse(Gender=="female", "f", "m"),
Genotype_Treatment = paste0(Genotype, "_", Treatment),
) %>%
dplyr::select(
Plate,
Cell_ID,
Mouse_ID,
Genotype,
Treatment,
Genotype_Treatment,
Age,
Gender) %>% as.data.frame()
rownames(df_expanded) <- df_expanded$Cell_ID
# View result
display_tab(df_expanded)
# Optional: save to CSV
write.csv(df_expanded,
"/files/scSeq_Hefendehl/data/20260126/534_Plate_positions_all.csv",
row.names = FALSE)
```
## load additional samples second GT
```{r load_add_samples 2nd, cache=F, dependson="load_add_meta"}
files <- list.files("/files/scSeq_Hefendehl/data/20260126/",
pattern = ".*Kallisto.*.csv",
full.names = T)
countlist = list()
for(f in files){
counts <- read.csv(f, row.names = 1, header = TRUE)
matched <- annotations$gene_name[match(
gsub("\\.[0-9]+", "", rownames(counts)),
annotations$gene_id)]
gene_names <- ifelse(is.na(matched), rownames(counts), matched)
dups <- duplicated(gene_names) | duplicated(gene_names, fromLast = TRUE)
dup_counts <- counts[dups, ] %>%
as.data.frame() %>%
mutate(gene = matched[dups]) %>%
group_by(gene) %>%
summarise(across(everything(), sum)) %>%
column_to_rownames("gene")
rownames(counts)[!dups] <- gene_names[!dups]
counts <- rbind(counts[!dups, ], dup_counts)
counts <- counts[rownames(counts) != "", ]
counts <- as((as.matrix(counts)), "sparseMatrix")
meta_counts <- df_expanded[colnames(counts), ]
counts <- CreateSeuratObject(counts = counts,
project = unique(meta_counts$Plate),
min.cells = 3,
min.features = 150,
meta.data =meta_counts)
Notannotgenes <- grepl("^ENSMUSG", rownames(counts))
counts <- counts[!Notannotgenes,]
print(unique(counts$Plate))
counts[["orig.ident"]] <- counts$Cell_ID
countlist[[unique(counts$Plate)]] <- counts
}
```
# Load batch 2025
```{r load_add_samples, cache=F, dependson="load_add_samples 2nd"}
Sample_files <- list.files("/files/scSeq_Hefendehl/data/20230203/Sample_Tables/",
pattern = "Sample.*.csv",
full.names = T)
counts_run1 <- read.csv("/files/scSeq_Hefendehl/data/20230203/Alignments/Run1/Kallisto/counts.csv", row.names = 1, header = TRUE) %>% rownames_to_column(var="GENE")
counts_run2 <- read.csv("/files/scSeq_Hefendehl/data/20230203/Alignments/Run2/Kallisto/counts.csv", row.names = 1, header = TRUE) %>% rownames_to_column(var="GENE")
counts <- counts_run1 %>% full_join(counts_run2, by = "GENE")
Platenames= c("025"="SS2_20_025_9",
"022"="SS2_20_022_7",
"024"="SS2_20_024_8",
"048"="SS2_21_048_6")
Samples_meta = data.frame()
for(s in Sample_files){
P = substr(basename(s), 14,16)
P = Platenames[P]
Samples <- read.csv(s, header = TRUE) %>%
mutate(
Plate=as.character(P),
Genotype = ifelse(Genotype=="WT", "WT", "APPPS1"),
Genotype_Treatment = paste0(Genotype, "_", Treatment),
Cell_ID=paste0(P,"_", Cell_ID),
Age = as.numeric(gsub(" weeks", "", Age))) %>%
dplyr::select(
"Plate",
"Cell_ID",
"Mouse_ID",
"Genotype",
"Treatment",
"Genotype_Treatment",
"Gender",
"Age")
Samples_meta <- rbind(Samples_meta, Samples)
}
counts_sel <- counts %>% dplyr::select(any_of(Samples_meta$Cell_ID), "GENE") %>% as.data.frame() %>% column_to_rownames("GENE")
counts_sel[is.na(counts_sel)]<-0
matched <- annotations$gene_name[match(
gsub("\\.[0-9]+", "", rownames(counts_sel)),
annotations$gene_id)]
gene_names <- ifelse(is.na(matched), rownames(counts_sel), matched)
dups <- duplicated(gene_names) | duplicated(gene_names, fromLast = TRUE)
dup_counts <- counts_sel[dups, ] %>%
as.data.frame() %>%
mutate(gene = matched[dups]) %>%
group_by(gene) %>%
summarise(across(everything(), sum)) %>%
column_to_rownames("gene")
rownames(counts_sel)[!dups] <- gene_names[!dups]
counts_sel <- rbind(counts_sel[!dups, ], dup_counts)
counts_sel <- counts_sel[rownames(counts_sel) != "", ]
counts_sel <- as((as.matrix(counts_sel)), "sparseMatrix")
Cells <- intersect(colnames(counts_sel), Samples_meta$Cell_ID)
counts_sel <- counts_sel[,Cells]
rownames(Samples_meta) <- Samples_meta$Cell_ID
Samples_meta <- Samples_meta[Cells,]
counts <- CreateSeuratObject(counts = counts_sel,
project = "Batch1",
min.cells = 3,
min.features = 150,
meta.data =Samples_meta)
Notannotgenes <- grepl("^ENSMUSG", rownames(counts))
counts <- counts[!Notannotgenes,]
counts_list_firstbatch <- SplitObject(counts, split.by = "Plate")
```
```{r}
seurat_list = c(counts_list_firstbatch, countlist)
seurat_obj <- merge(
seurat_list[[1]],
y = seurat_list[-1],
add.cell.ids = names(seurat_list))
# Exclusde not gene_name annotated genes
```
```{r}
rm(list=c("meta_counts",
"samples.integrated",
"seurat_list",
"tmp_counts",
"tmp_meta",
"tmp_seurat",
"ah",
"ahDb",
"countlist",
"counts",
"counts_list_firstbatch",
"edb",
"df", "df_expanded",
"df_raw",
"dup_counts"))
gc()
```
```{r, fig.width=18, fig.height=5}
Idents(seurat_obj) <- "Plate"
DefaultAssay(seurat_obj) <- "RNA"
seurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj, pattern = "^mt-")
seurat_obj <- JoinLayers(seurat_obj)
seurat_obj[["RNA"]] <- JoinLayers(seurat_obj[["RNA"]])
seurat_obj <- seurat_obj %>%
NormalizeData() %>%
FindVariableFeatures() %>%
ScaleData() %>%
RunPCA() %>%
FindNeighbors(dims = 1:30, reduction = "pca") %>%
FindClusters(cluster.name = "unintegrated_clusters") %>%
RunUMAP(dims = 1:30, reduction = "pca", reduction.name = "umap.unintegrated")
p <- DimPlot(seurat_obj, reduction = "umap.unintegrated", group.by = c("Plate", "Mouse_ID", "unintegrated_clusters"))
p
ggsave(paste0(outputdir, "DataLoad_umap.unintegrated.pdf"), p, width = 18, height=5)
```
```{r fig.width=12, fig.height=4}
VlnPlot(seurat_obj, features = c("nCount_RNA", "nFeature_RNA", "percent.mt"),
group.by = "Plate")
seurat.obj.filt <- subset(seurat_obj,
subset = nFeature_RNA > 100 &
nFeature_RNA < 7000 &
nCount_RNA <3e6 &
nCount_RNA > 400 &
percent.mt < 5)
counts <- JoinLayers(seurat.obj.filt, assay = "RNA") %>% GetAssayData(assay = "RNA", layer="data")
genes.percent.expression <- rowMeans(counts > 0)
keep.genes <- names(genes.percent.expression[genes.percent.expression >= 0.005])
# filter for expression in at least 0.5% of cells
seurat.obj.filt <- seurat.obj.filt[keep.genes, ]
VlnPlot(seurat.obj.filt, features = c("nCount_RNA", "nFeature_RNA", "percent.mt"),
group.by = "Plate")
# 7000 or less than 150 and cells with >7% mitochondrial counts.
samples.integrated <- seurat.obj.filt
```
```{r}
#samples.integrated = readRDS("/files/scSeq_Hefendehl/data/20230203/seu-SCT-harmony_02-02-23.rds")
samples.integrated[["Genotype_Treatment"]] = factor(samples.integrated$Genotype_Treatment,
levels = c("WT_Ctrl",
"WT_Stroke",
"APPPS1_Ctrl",
"APPPS1_Stroke"))
[email hidden] %>%
dplyr::select(Mouse_ID, Genotype_Treatment, Gender, Age) %>%
group_by(Mouse_ID) %>% mutate(
nCells = n()
) %>% ungroup() %>%
dplyr::distinct() -> Samplelist
openxlsx::write.xlsx(Samplelist, "/files/scSeq_Hefendehl/data/Samplelist.xlsx",rowNames = TRUE)
display_tab(Samplelist)
```
## Reclustering of datasets
### Define Gene lists
```{r}
DAM_up <- c("Itgax","Cst7","Ccl4","Clec7a","Lpl","Siglec1","Spp1","Axl","Apoe","Trem2","Csf1","Lyz2","H2-D1","Tyrobp","Cd74","Gpnmb","B2m","Cd9","Ctsb","Ctsd") # genes up regulated in DAMs versus homestatics
# additions
DAM1_up <- c("Tyrobp", "Ctsb", "Ctsd", "Apoe", "B2m", "Fth1", "Lyz2") # genes up regulated in DAM1 vs Homeostasis
DAM1_down <- c("Cx3cr1", "P2ry12", "Tmem119")# genes down regulated in DAM1 vs homeostasis
DAM2_up <- c("Trem2", "Axl", "Cst7", "Ctsl", "Lpl", "Cd9", "Csf1", "Ccl6", "Itgax", "Clec7a", "Lilrb4", "Timp2")
DAM1vsDAM2 = c("Trem2", "Cst7", "Spp1", "Lpl", "Itgax", "Ctsl", "Ctsz", "Cd68", "Axl", "Cd9", "Ccl6", "Csf1")
homeostat.genes=c("P2ry12","P2ry13","Tmem119","Cx3cr1","Selplg","Cd33")
INF_resp_microglia <- c("Ifit2", "Ifit3", "Ifitm3", "Irf7", "Oasl2")
Cycling_microglia <- c("Top2a", "Mcm2", "Tubb5", "Mki67", "Cdk1")
Act_resp_microglia <- c("Cd74", "H2-Ab1", "H2-Aa", "Ctsb", "Ctsd")
Axon_tract_microglia <- c("Spp1", "Gpnmb", "Igf1", "Lgals3", "Fabp5", "Lpl", "Lgals1", "Ctsl", "Anxa5")
cluster.genes=unique(c(DAM_up, DAM1vsDAM2, homeostat.genes))
```
```{r}
ngenes = 50
distal_dat = openxlsx::read.xlsx(
"/files/scSeq_Hefendehl/data/Copy of DEG_control_vs_stroke_vs_distal_original.xlsx",
sheet = 1)
distal_dat <- distal_dat[order(distal_dat$score, decreasing = T)[1:ngenes],]
distal_genes <- distal_dat$X1
control_dat = openxlsx::read.xlsx("/files/scSeq_Hefendehl/data/Copy of DEG_control_vs_stroke_vs_distal_original.xlsx",sheet = 2)
control_dat <- control_dat[order(control_dat$score, decreasing = T)[1:ngenes],]
control_genes <- control_dat$X1
Stroke_dat = openxlsx::read.xlsx("/files/scSeq_Hefendehl/data/Copy of DEG_control_vs_stroke_vs_distal_original.xlsx",sheet = 3)
Stroke_dat <- Stroke_dat[order(Stroke_dat$score, decreasing = T)[1:ngenes],]
Stroke_genes <- Stroke_dat$X1
```
## inspect data
### SCT CCA integration all cells
```{r, fig.width=8,fig.height=4}
DefaultAssay(samples.integrated) = "RNA"
samples.integrated <- UpdateSeuratObject(samples.integrated)
samples.integrated$Run <-substr(samples.integrated$Plate, 1, 6)
samples.integrated[["RNA"]] <- split(samples.integrated[["RNA"]], f= samples.integrated$Run)
samples.integrated <- SCTransform(samples.integrated, vst.flavor = "v2")
samples.integrated <- RunPCA(samples.integrated, npcs = 30, verbose = FALSE)
# one-liner to run Integration
DefaultAssay(samples.integrated) <- "RNA"
samples.integrated <- IntegrateLayers(object = samples.integrated, method = CCAIntegration,
orig.reduction = "pca", new.reduction = 'integrated.cca')
samples.integrated <- FindNeighbors(samples.integrated, reduction = "integrated.cca", dims = 1:30)
samples.integrated <- FindClusters(samples.integrated, resolution = 0.6)
samples.integrated <- RunUMAP(samples.integrated, reduction = "integrated.cca", dims = 1:30, reduction.name = "umap.integrated.cca")
a = DimPlot(samples.integrated, reduction = "umap.integrated.cca", group.by = "Plate")
b = FeaturePlot_scCustom(samples.integrated, reduction = "umap.integrated.cca", features = "Aif1", colors_use = viridis_dark_high)
b2 = FeaturePlot_scCustom(samples.integrated, reduction = "umap.integrated.cca", features = "Ptprc", colors_use = viridis_dark_high)
c = RidgePlot(samples.integrated, features = "Aif1", group.by = "Mouse_ID")
c2 = RidgePlot(samples.integrated, features = "Ptprc", group.by = "Mouse_ID")
d = DimPlot(samples.integrated, reduction = "umap.integrated.cca", split.by = "Genotype_Treatment", pt.size = 1)
e = FeaturePlot_scCustom(samples.integrated, reduction = "umap.integrated.cca", features = "Aif1", split.by = "Genotype_Treatment", colors_use = viridis_dark_high)
e2 = FeaturePlot_scCustom(samples.integrated, reduction = "umap.integrated.cca", features = "Ptprc", split.by = "Genotype_Treatment", colors_use = viridis_dark_high)
```
```{r, fig.width=12,fig.height=5}
a|b|b2
ggarrange(c,c2, common.legend=T, legend = "none")
ggsave(paste0(outputdir,"cluster_cca_integrated_nofilter.pdf"), a)
```
```{r, fig.width=12,fig.height=4}
d
e
e2
```
### Sex specific genes
```{r, fig.width=12,fig.height=4}
Idents(samples.integrated) <- "Mouse_ID"
Sexgenes = c("Ddx3y", #y linked and microglia expressed
"Xist") #X inactivation
# cluster all cells for visualization purpose prior to exclusion
a = DimPlot(samples.integrated, reduction = "umap.integrated.cca")
b = FeaturePlot_scCustom(samples.integrated, reduction = "umap.integrated.cca",
features = "Xist",
colors_use = viridis_dark_high)
c = RidgePlot(samples.integrated, features = Sexgenes)
d = DimPlot(samples.integrated,
reduction = "umap.integrated.cca",
group.by = "Plate",
split.by = "Mouse_ID",
pt.size = 1)
e = FeaturePlot_scCustom(samples.integrated,
reduction = "umap.integrated.cca",
features = Sexgenes,
split.by = "Mouse_ID",
colors_use = viridis_dark_high)
```
```{r, fig.width=10,fig.height=4}
a|b
c
ggsave(paste0(outputdir, "cluster_Xist_harmony_integrated_nofilter.pdf"), a|b)
ggsave(paste0(outputdir, "cluster_Xist_Ridge_integrated_nofilter.pdf"), c)
```
```{r, fig.width=30,fig.height=4}
d
```
```{r fig.width=6, fig.height=6}
DimPlot(samples.integrated, reduction="umap.integrated.cca", group.by = "Genotype_Treatment")
```
```{r fig.width=80, fig.height=15}
e
ggsave(paste0(outputdir, "cluster_Sexgenes_integrated_nofilter.pdf"), e, limitsize = F)
```
```{r, include=FALSE}
reads = JoinLayers(samples.integrated, assay = "RNA") %>% GetAssayData(assay = "RNA", layer="data")
reads = t(reads) %>% as.data.frame()
metadata <- [email hidden]
metadata$Aif1neg = NA
metadata$Aif1neg <- paste0("Aif1neg_",as.character(reads$Aif1 ==0))
metadata -> [email hidden]
```
```{r}
rm( list=c("reads", "seurat_obj", "seurat.obj.filt"))
```
### remove sex specific and not-annotated genes
```{r, fig.width=4,fig.height=4}
samples.integrated.filt <- samples.integrated
sex.genes <- rownames(samples.integrated.filt) %in% annotations$gene_name[grepl("X|Y", annotations$seq_name)]
Gmgenes <- grepl("^Gm", rownames(samples.integrated.filt))
Rikgenes <- grepl(".*Rik$", rownames(samples.integrated.filt))
samples.integrated.filt<- samples.integrated.filt[!(sex.genes|Gmgenes|Rikgenes),]
```
# Clean integration
## Pre cell filtering integration
```{r, fig.width=4,fig.height=4}
DefaultAssay(samples.integrated.filt)<- "RNA"
samples.integrated.filt <- samples.integrated.filt %>% RunHarmony("Run")
SDs=Stdev(object = samples.integrated.filt, reduction = "pca")
varexpl <- SDs/sum(SDs)
ndims <- which(cumsum(varexpl)>0.75)[1]
ElbowPlot(samples.integrated.filt, reduction = "pca", ndims = 30) + geom_vline(xintercept =ndims)
samples.integrated.filt <- samples.integrated.filt %>%
FindNeighbors(reduction = "integrated.cca", dims = 1:ndims) %>%
FindClusters(cluster.name = "integrated_clusters") %>%
RunUMAP(dims = 1:ndims, reduction = "integrated.cca")
```
```{r, fig.width=12,fig.height=6}
DimPlot(samples.integrated.filt, reduction = "umap",
group.by = c("Plate", "Mouse_ID", "integrated_clusters", "Genotype_Treatment"))
```
## map cell types
```{r}
ref.sce.mm <- celldex::MouseRNAseqData()
sce.filt<- as.SingleCellExperiment(JoinLayers(samples.integrated.filt, assay = "RNA"), assay = "RNA")
pred.mouseRNA <- SingleR(test =sce.filt, ref = ref.sce.mm, assay.type.test=1,
labels = ref.sce.mm$label.fine)
table(pred.mouseRNA$labels) #fine tuned cell types
#table(pred.mouseRNA$pruned.labels) #cell type after pruning
samples.integrated.filt$Immgen_sc_labels_agc <- NA
samples.integrated.filt$Immgen_sc_labels_agc <- pred.mouseRNA$labels
```
```{r}
pheat=plotScoreHeatmap(pred.mouseRNA, show.labels = TRUE,
annotation_col=data.frame(cluster=samples.integrated.filt$seurat_clusters,
row.names=rownames(pred.mouseRNA)))
Idents(samples.integrated.filt)<- "seurat.integrated"
a<- DimPlot(samples.integrated.filt, reduction="umap", group.by="Immgen_sc_labels_agc",
label = F)
b<- DimPlot(samples.integrated.filt, reduction="umap", group.by="Immgen_sc_labels_agc",
label = F, split.by = "Plate")
c<- DimPlot(samples.integrated.filt, reduction="umap", group.by="Immgen_sc_labels_agc", split.by = "Genotype_Treatment",
label = F)
```
```{r}
print(pheat)
ggsave(filename = paste0(outputdir, "CelltypesHeatmap_Opt_All_genes_all_cells_Clusters_integrated.pdf"), pheat)
```
```{r}
a
ggsave(filename = paste0(outputdir, "Celltypes_Opt_All_genes_all_cells_Clusters_integrated.pdf"), a)
```
```{r, fig.width=18, fig.height=4}
b
ggsave(filename = paste0(outputdir, "Celltypes_Opt_All_genes_all_cells_Clusters_byPlate_integrated.pdf"), b)
```
```{r, fig.width=12, fig.height=4}
c
ggsave(filename = paste0(outputdir, "Celltypes_Opt_All_genes_all_cells_Clusters_byCondition_integrated.pdf"), c)
```
## filter and recluster only Microglia
### Re-optimize clustering
```{r}
samples.integrated.filt <- samples.integrated.filt[,grepl("Microglia",samples.integrated.filt$Immgen_sc_labels_agc)]
DefaultAssay(samples.integrated.filt)<- "RNA"
samples.integrated.filt <- samples.integrated.filt %>%
NormalizeData() %>%
ScaleData() %>%
RunPCA() %>%
RunHarmony("Run")
```
#### Layer integration across Runs using CCA Integration
```{r}
reslist= seq(0.2,1.2,0.2)
samples.integrated.filt <- samples.integrated.filt %>%
IntegrateLayers(method = CCAIntegration,
orig.reduction = "pca",
new.reduction = 'integrated.cca') %>%
FindNeighbors(reduction = "integrated.cca", dims = 1:ndims) %>%
FindClusters(resolution = reslist) %>%
RunUMAP(reduction = "integrated.cca", min.dist = 0.1, spread = 1, dims=1:ndims)
clustreeplot<-clustree(samples.integrated.filt, prefix = "RNA_snn_res.")
samples.integrated.filt <- samples.integrated.filt %>% RunUMAP(reduction = "integrated.cca", min.dist = 0.1, spread = 1, dims=1:ndims, reduction.name = "umap")
```
### Define resoultion for clustering
```{r, fig.width=15, fig.height=12}
plotlist=list()
for( cl in (grep("RNA_snn", colnames([email hidden]), value=T))){
plotlist[[cl]] <- DimPlot(samples.integrated.filt, reduction="umap", group.by = cl)
}
plotlist <- plotlist[sort(names(plotlist))]
ggarrange(plotlist=plotlist)
```
```{r, fig.width=5, fig.height=6}
clustreeplot
ggsave(filename = paste0(outputdir, "All_genes_Clustertree_integrated.pdf"), clustreeplot)
```
```{r}
library(mclust)
samples.integrated.filt.meta <- [email hidden]
res_cols <- sort(grep("RNA_snn_res", colnames(samples.integrated.filt.meta), value = TRUE))
ari_results <- data.frame()
for (i in 1:(length(res_cols) - 1)) {
col1 <- res_cols[i]
col2 <- res_cols[i + 1]
ari <- adjustedRandIndex(samples.integrated.filt.meta[[col1]],
samples.integrated.filt.meta[[col2]])
ari_results <- rbind(ari_results, data.frame(res1 = col1, res2 = col2, ARI = ari))
}
best_res <- ari_results$res1[which.max(ari_results$ARI)]
best_res <- as.numeric(gsub("RNA_snn_res.", "", best_res))
ari_results
best_res <- 0.6
```
clustering resolultion was set to `r best_res`
```{r, warning=F, }
DefaultAssay(samples.integrated.filt) <- "RNA"
samples.integrated.filt <- samples.integrated.filt %>%
FindNeighbors(reduction = "integrated.cca", dims = 1:ndims) %>%
FindClusters(resolution = best_res, cluster.name = "seurat.integrated") %>%
RunUMAP(reduction = "integrated.cca", min.dist = 0.1, spread = 1, dims=1:ndims)
n = length(unique(samples.integrated.filt$seurat.integrated))
Idents(samples.integrated.filt) <- "seurat.integrated"
samples.integrated.filt$labels<- factor(Idents(samples.integrated.filt),
levels = c(0:(n-1)), labels=paste0("Microglia", 0:(n-1))
)
Idents(samples.integrated.filt) = "labels"
DefaultAssay(samples.integrated.filt) <- "RNA"
samples.integrated.filt[["RNA"]] <- samples.integrated.filt[["RNA"]] %>% JoinLayers()
samples.integrated.filt[["RNA"]]$scale.data <- NULL
samples.integrated.filt<- samples.integrated.filt %>% NormalizeData() %>% ScaleData()
a<-DimPlot(samples.integrated.filt, , reduction = "umap",
group.by = "seurat.integrated", pt.size =1, label = T) + NoLegend()
b<-DimPlot(samples.integrated.filt, reduction = "umap", split.by = "Plate", pt.size = 1)
c<-DimPlot(samples.integrated.filt, reduction = "umap", split.by = "Genotype_Treatment", pt.size =1)
d <- DimPlot(samples.integrated.filt, reduction = "umap", split.by = "Mouse_ID", pt.size =1)
e<-FeaturePlot_scCustom(samples.integrated.filt,reduction = "umap", layer = "scale.data",
features = homeostat.genes, colors_use = viridis_dark_high)
```
```{r, fig.width=5, fig.height=5}
a
ggsave(filename =
paste0(outputdir,
"Opt_All_genes_all_cells_Clusters_integrated.pdf"), a)
```
```{r, fig.width=30, fig.height=5}
b
ggsave(filename =
paste0(outputdir,
"Opt_All_genes_all_cells_Clusters_byPlate_integrated.pdf"),
b)
```
```{r, fig.width=12, fig.height=4}
c
ggsave(filename = paste0(
outputdir,
"Opt_All_genes_all_cells_Clusters_byCondition_integrated.pdf"),
c)
```
```{r, fig.width=30, fig.height=3}
d
ggsave(filename = paste0(
outputdir,
"Opt_All_genes_all_cells_Clusters_byMouseID_integrated.pdf"),
d)
```
```{r, fig.width=8, fig.height=12}
e
ggsave(filename = paste0(outputdir,
"All_genes_Clusters_byHomeostatgenes_integrated.pdf"), e)
```
```{r}
a <- DimPlot(samples.integrated.filt, reduction="umap", label = T)
b <- DimPlot(samples.integrated.filt, reduction="umap", split.by = "Plate", label = F)
c <- DimPlot(samples.integrated.filt, reduction="umap", split.by = "Genotype_Treatment", label = F)
d <- DimPlot(samples.integrated.filt, reduction = "umap", group.by = "Genotype_Treatment")
e <- DimPlot(samples.integrated.filt, reduction = "umap", group.by = "Gender")
g <- DimPlot(samples.integrated.filt, reduction = "umap", group.by = "Aif1neg")
h <- DimPlot(samples.integrated.filt, reduction = "umap", group.by = "Immgen_sc_labels_agc")
samples.integrated.filt[["Immgen_sc_label_short"]] <- ifelse(
grepl("Microglia", samples.integrated.filt$Immgen_sc_labels_agc),
"Microglia", NA)
i <- DimPlot(samples.integrated.filt, reduction = "umap", group.by = "Immgen_sc_label_short")
```
```{r, fig.width=8, fig.height=7}
a
ggsave(filename = paste0(outputdir, "Microglia_Opt_All_genes_Clusters_integrated.pdf"), a)
d
ggsave(filename = paste0(outputdir, "Microglia_Opt_All_genes_Condition_integrated.pdf"), d)
e
ggsave(filename = paste0(outputdir, "Microglia_Opt_All_genes_Gender_integrated.pdf"), e)
g
ggsave(filename = paste0(outputdir, "Microglia_Opt_All_genes_Aif1neg_integrated.pdf"), g)
h
ggsave(filename = paste0(outputdir, "Microglia_Opt_All_genes_Immgen_sc_labels_agc_integrated.pdf"), h)
i
ggsave(filename = paste0(outputdir, "Microglia_Opt_All_genes_Immgen_sc_labels_summed_integrated.pdf"), i)
```
```{r, fig.width=12, fig.height=4}
c
ggsave(filename = paste0(outputdir, "Microglia_Opt_All_genes_Clusters_byCondition.pdf"), c)
```
```{r, fig.width=18, fig.height=3}
b
ggsave(filename = paste0(outputdir, "Microglia_Opt_All_genes_Clusters_byPlate_integrated.pdf"), b)
```
```{r, fig.height=18, fig.width=10}
f <- FeaturePlot_scCustom(samples.integrated.filt, reduction ="umap",
features = sort(cluster.genes),order =T,
colors_use = viridis_dark_high)
f
ggsave(filename = paste0(outputdir, "Microglia_Opt_all_genes_MicrogliaGenesintegrated.pdf"), f)
```
```{r, fig.width=8, fig.height=8}
genelist <- c("Alpk1", "Mmp12", "Igf1", "S100a6", "Apoc4", "Dock10", "Nanog", "Olfr344")
f <- FeaturePlot_scCustom(samples.integrated.filt, reduction = "umap",
layer = "scale.data",
features = sort(genelist),order =T,
colors_use = viridis_dark_high)
f
ggsave(filename = paste0(outputdir, "Microglia_Opt_all_genes_AddGenelistintegrated.pdf"), f)
```
# Map initial clusters batch 1
```{r, fig.width=8, fig.height=7}
samples.integrated.filt_batch1 <- readRDS("./output/SeuratObjectafterProcessing.rds")
meta <- [email hidden]
meta_B1 <- [email hidden]
lkup <- c(SS2_20_22="SS2_20_022",
SS2_20_24="SS2_20_024",
SS2_20_25="SS2_20_025",
SS2_21_48="SS2_21_048")
meta_B1$plate_new <- lkup[meta_B1$Plate]
meta$label_B1<- meta_B1[match(meta$Cell_ID, paste0(meta_B1$plate_new, "_",meta_B1$Cell_ID)), "labels"]
[email hidden] <- meta
p<- DimPlot(samples.integrated.filt,
group.by = "label_B1",
reduction = "umap", order = T) + scale_color_discrete(na.value = "#F1F1F1")
ggsave(paste0(outputdir, "Overlay_Initital_clusters.pdf"), p)
p
table(samples.integrated.filt$labels, samples.integrated.filt$label_B1)
ggsave(filename = paste0(outputdir, "Microglia_Opt_all_genes_Batch1_Old_labels_integrated.pdf"), p)
rm("samples.integrated.filt_batch1")
```
```{r, fig.width=7, fig.height=7}
DefaultAssay(samples.integrated.filt) <- "RNA"
FeaturePlot_scCustom(samples.integrated.filt, reduction="umap", features = "P2ry1", colors_use = viridis_dark_high, layer = "scale.data")
```
# Markers for Clusters
```{r}
metadata<[email hidden]
res = compareGroups(Genotype_Treatment~labels, data = metadata)
compareGroups::createTable(res, show.p.mul = T)
compareGroups::createTable(res, show.p.mul = T) %>% export2xls(paste0(outputdir, "Cluster_distribution_genotype.xlsx"))
res = compareGroups(Genotype_Treatment~Mouse_ID, data = metadata, max.xlev = 70, max.ylev = 20)
compareGroups::createTable(res, show.p.mul = T)
compareGroups::createTable(res, show.p.mul = T) %>% export2xls(paste0(outputdir, "Cellcount_distribution_MouseBygenotype.xlsx"))
res = compareGroups(labels~Mouse_ID, data = metadata, max.xlev = 70, max.ylev = 20)
compareGroups::createTable(res, show.p.mul = T)
compareGroups::createTable(res, show.p.mul = T) %>% export2xls(paste0(outputdir, "Cellcount_distribution_MouseByCluster.xlsx"))
```
```{r, fig.height=18, fig.width=24}
DefaultAssay(samples.integrated.filt) <- "RNA"
f <- RidgePlot(samples.integrated.filt,
features = sort(cluster.genes),
assay = "RNA",
layer = "data", combine = F)
f <- ggarrange(plotlist = f, common.legend = T, legend = "right")
f
ggsave(filename = paste0(outputdir, "Microglia_Opt_all_genes_MicrogliaGenesRidgePlotintegrated.pdf"), f)
```
### all clusters
```{r}
Idents(samples.integrated.filt) <- "labels"
DefaultAssay(samples.integrated.filt)<- "RNA"
#samples.integrated.filt <- PrepSCTFindMarkers(samples.integrated.filt)
# remove not annotated genes
geneunivers<-rownames(samples.integrated.filt)
Markers=FindAllMarkers(samples.integrated.filt,
logfc.threshold = 1,
assay = "RNA", slot = "data",
features = geneunivers,
only.pos = T)
display_tab(Markers)
openxlsx::write.xlsx(Markers, paste0(outputdir, "Clusters_Markers_all.xlsx"),rowNames = TRUE)
genelist = list()
for ( i in unique(Markers$cluster)) genelist[[i]] <- Markers$gene[Markers$cluster==i]
GOterm=getGOresults(genelist, genereference = geneunivers, organism = "mmusculus")
display_tab(GOterm$result)
openxlsx::write.xlsx(GOterm$result, paste0(outputdir, "Clusters_GOterms_Markers_all.xlsx"),rowNames = TRUE)
```
```{r}
metadata=[email hidden]
restab = table(GT_Treat = metadata$Genotype_Treatment, DAMs=metadata$labels)
restab
prop.table(restab, margin = 1)
chisq.test(restab)
proplot <- ggplot(as.data.frame(restab),
aes(fill=DAMs, x=GT_Treat, y=Freq))+
geom_bar(position="fill", stat="identity", color="black", linewidth=0.5)+theme_classic()+ggtitle("Micorglia cluster proportions")
proplot
ggsave(filename = paste0(outputdir, "Microglia_Clusterproportions_per_condition_all_integrated.pdf"), proplot)
library(chisq.posthoc.test)
chisq.posthoc.test(restab)
#To Do Pairwise
```
```{r, fig.width=10, fig.height=12}
plotlist = list()
mat <- GetAssayData(samples.integrated.filt, assay = "RNA", slot = "data")
temp_obj <- samples.integrated.filt
for (m in names(genelist)){
#reads = rowSums(samples.integrated.filt@assays$RNA$counts)
targets <- Markers[Markers$cluster==m, ] %>% arrange(desc(pct.1))
targets <- targets$gene[1:10]
mat_sel <- mat[targets, , drop = FALSE]
mat_norm <- t(apply(mat_sel, 1, function(x) (x - min(x)) / (max(x) - min(x) + 1e-9))) %>% as.matrix()
temp_obj@assays$RNA@layers$data[match(targets, rownames(temp_obj)), ] <- mat_norm[targets, colnames(temp_obj)]
p<-DoHeatmap(temp_obj, features = targets, disp.max = 1,
group.by = "labels", assay = "RNA", slot = "data",
label = F, combine=F)
p <- p[[1]] + scale_fill_viridis()
ggsave(p, file=paste0(outputdir, "Top10GenesByExpr_",m,"_Heatmaps_byCluster_integrated.pdf"))
plotlist[[m]]<-p
}
p<- ggarrange(plotlist = plotlist, labels = names(genelist), common.legend = T)
p
ggsave(file=paste0(outputdir, "Top10Genes_Heatmaps_byCluster_Row_Scaled_integrated.pdf"), p)
```
```{r, fig.width=15, fig.height=12}
plotlist = list()
mat <- GetAssayData(samples.integrated.filt, assay = "RNA", slot = "data")
for (m in names(genelist)){
#reads = rowSums(samples.integrated.filt@assays$RNA$counts)
targets <- Markers[Markers$cluster==m, ] %>% arrange(desc(pct.1))
targets <- targets$gene[1:10]
mat_sel <- mat[targets, , drop = FALSE]
max_plot <- mean(mat_sel)+2*sd(mat_sel)
p<-DoHeatmap(samples.integrated.filt, features = targets, disp.max = ceiling(max_plot),
group.by = "labels", assay = "RNA", slot = "data",
label = F, combine=F)
p <- p[[1]] + scale_fill_viridis()
ggsave(p, file=paste0(outputdir, "Top10GenesByExpr_",m,"_Heatmaps_byCluster_integrated.pdf"))
plotlist[[m]]<-p
}
p<- ggarrange(plotlist = plotlist, labels = names(genelist), common.legend = F)
p
ggsave(file=paste0(outputdir, "Top10Genes_Heatmaps_byCluster_integrated.pdf"), p)
```
```{r, fig.width=14, fig.height=7}
FeaturePlot_scCustom(samples.integrated.filt, features = c("Lipe1", "Rfx7", "Olfr344", "S100a8"),
layer="data", colors_use = viridis_dark_high, , reduction = "umap")
```
```{r, fig.width=7, fig.height=7}
DimPlot(samples.integrated.filt, reduction = "umap")
```
## Comparison of DEGs in DAMs between treatments/genotypes
### define DAM
DAM profiles based on Cell. 2017 Jun 15;169(7):1276-1290.e17. doi: 10.1016/j.cell.2017.05.018)
##### DAM1
```{r, fig.width=12, fig.height=4}
#load("./data/DAM_genelists.RData")
DefaultAssay(samples.integrated.filt) = "RNA"
samples.integrated.filt<- AddModuleScore(object = samples.integrated.filt,
assay="RNA", slot="data",
features = list(DAM1_up), name = "DAM_score")
#FeaturePlot_scCustom(samples.integrated.filt, features = "DAM_score1")+scale_color_viridis_c()
samples.integrated.filt <- AddModuleScore(object = samples.integrated.filt,
assay="RNA", slot="data",
features = list(DAM_up), name = "DAM_score")
p4 <- FeaturePlot_scCustom(seurat_object = samples.integrated.filt, reduction = "umap",
features = "DAM_score1", colors_use = viridis_dark_high)
p4a <- DimPlot(object = samples.integrated.filt, reduction = "umap")
DAMcutoff = 1.5
phist = ggplot(data = [email hidden], aes(x=DAM_score1))+
geom_density()+geom_vline(xintercept = DAMcutoff)+
theme_classic()
samples.integrated.filt$DAM_binary = as.factor(samples.integrated.filt$DAM_score1>=DAMcutoff)
pbin <- DimPlot(object = samples.integrated.filt, group.by = "DAM_binary", reduction = "umap")
comb = phist | p4 | pbin| p4a
comb
prop.table(table(DAMs=samples.integrated.filt$DAM_binary, Clusters=samples.integrated.filt$labels), margin = 2)
dam_GT_p <-DimPlot(object = samples.integrated.filt, group.by = "DAM_binary", split.by = "Genotype_Treatment")
ggsave(filename = paste0(outputdir, "Microglia_DAM1_callingintegrated.pdf"), comb)
ggsave(filename = paste0(outputdir, "Microglia_DAM1_bygenotype_integrated.pdf"), dam_GT_p)
```
```{r}
phist
p4
pbin
ggsave(filename = paste0(outputdir, "Microglia_DAM1_thresholding_integrated.pdf"), phist)
ggsave(filename = paste0(outputdir, "Microglia_DAM1_scoring_integrated.pdf"), p4)
ggsave(filename = paste0(outputdir, "Microglia_DAM1_binary_integrated.pdf"), pbin)
```
##### DAM2
```{r, fig.width=12, fig.height=4}
samples.integrated.filt <- AddModuleScore(object = samples.integrated.filt,
features = list(DAM2_up), name = "DAM2_score")
p4 <- FeaturePlot_scCustom(seurat_object = samples.integrated.filt, reduction = "umap",
features = "DAM2_score1", colors_use = viridis_dark_high)
p4a <- DimPlot(object = samples.integrated.filt, reduction = "umap")
DAM2cutoff = 1.0
phist = ggplot(data = [email hidden], aes(x=DAM2_score1))+
geom_density()+geom_vline(xintercept = DAM2cutoff)+
theme_classic()
samples.integrated.filt$DAM2_binary = as.factor(samples.integrated.filt$DAM2_score1>=DAM2cutoff)
pbin <- DimPlot(object = samples.integrated.filt, reduction = "umap",
group.by = "DAM2_binary")
comb = phist | p4 | pbin| p4a
dam_GT_p <-DimPlot(object = samples.integrated.filt,reduction = "umap", group.by = "DAM2_binary", split.by = "Genotype_Treatment")
comb
dam_GT_p
ggsave(filename = paste0(outputdir, "Microglia_DAM2_calling_integrated.pdf"), comb)
ggsave(filename = paste0(outputdir, "Microglia_DAM2_bygenotype_integrated.pdf"), dam_GT_p)
prop.table(table(DAM2=samples.integrated.filt$DAM2_binary, Clusters=samples.integrated.filt$labels), margin = 2)
prop.table(table(DAM2=samples.integrated.filt$DAM2_binary, DAMs=samples.integrated.filt$DAM_binary), margin = 2)
```
```{r, fig.width=12}
a <- ggplot([email hidden], aes(x=DAM2_score1, y=DAM_score1))+
geom_density_2d_filled()+geom_point(col="#FFFFFF55")+theme_classic()+ggtitle("all cells")
b<- ggplot([email hidden][samples.integrated.filt$DAM_binary==TRUE,], aes(x=DAM2_score1, y=DAM_score1))+geom_density_2d_filled()+geom_point(col="#FFFFFF55")+theme_classic()+ggtitle("DAMs only")
combo <- a|b
combo
ggsave(filename = paste0(outputdir, "Microglia_DAMvsDAM2_contourplot_integrated.pdf"), combo)
```
```{r}
metadata=[email hidden]
restab = table(GT_Treat = metadata$Genotype_Treatment, DAMs=metadata$DAM_binary)
restab
prop.table(restab, margin = 1)
chisq.test(restab)
proplot <- ggplot(as.data.frame(restab),
aes(fill=DAMs, x=GT_Treat, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()+ggtitle("DAM proportions")
proplot
ggsave(filename = paste0(outputdir, "Microglia_DAMproportions_per_condition_all_integrated.pdf"), proplot)
library(chisq.posthoc.test)
chisq.posthoc.test(restab)
#To Do Pairwise
```
#### DAMs in APPPS1 vs. DAMs in APPPS1 + Stroke
```{r}
metadata=[email hidden][samples.integrated.filt$Genotype == "APPPS1", ]
restab = table(GT_Treat = as.character(metadata$Genotype_Treatment), DAMs=metadata$DAM_binary)
restab
prop.table(restab, margin = 1)
fisher.test(restab)
proplot<-ggplot(as.data.frame(restab),
aes(fill=DAMs, x=GT_Treat, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
proplot
ggsave(filename = paste0(outputdir, "Microglia_DAMproportions_per_condition_APPPS1_integrated.pdf"), proplot)
```
#### DAMs in WT vs. DAMs in WT + Stroke
```{r}
metadata=[email hidden][samples.integrated.filt$Genotype == "WT", ]
restab = table(GT_Treat = as.character(metadata$Genotype_Treatment), DAMs=metadata$DAM_binary)
restab
prop.table(restab, margin = 1)
fisher.test(restab)
proplot<-ggplot(as.data.frame(restab),
aes(fill=DAMs, x=GT_Treat, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
proplot
ggsave(filename = paste0(outputdir, "Microglia_DAMproportions_per_condition_WT_integrated.pdf"), proplot)
```
#### DAM1 vs DAM2
```{r}
metadata=[email hidden][samples.integrated.filt$DAM_binary==TRUE,]
restab = table(GT_Treat = metadata$Genotype_Treatment, DAM2=metadata$DAM2_binary)
restab
prop.table(restab, margin = 1)
chisq.test(restab)
proplot<-ggplot(as.data.frame(restab),
aes(fill=DAM2, x=GT_Treat, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
chisq.posthoc.test(restab)
proplot
ggsave(filename = paste0(outputdir, "Microglia_DAM2proportionsofDAM1_per_condition_integrated.pdf"), proplot)
```
```{r}
metadata=[email hidden][(samples.integrated.filt$DAM_binary==TRUE) & (samples.integrated.filt$Genotype=="APPPS1"),]
restab = table(GT_Treat = as.character(metadata$Genotype_Treatment), DAM2=metadata$DAM2_binary)
restab
prop.table(restab, margin = 1)
chisq.test(restab)
proplot<- ggplot(as.data.frame(restab),
aes(fill=DAM2, x=GT_Treat, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
proplot
ggsave(filename = paste0(outputdir, "Microglia_DAM2proportionsofDAM1_per_condition_APPPS1_integrated.pdf"), proplot)
```
```{r}
metadata=[email hidden][(samples.integrated.filt$DAM_binary==TRUE) & (samples.integrated.filt$Genotype=="WT"),]
restab = table(GT_Treat = as.character(metadata$Genotype_Treatment), DAM2=metadata$DAM2_binary)
restab
prop.table(restab, margin = 1)
chisq.test(restab)
proplot<-ggplot(as.data.frame(restab),
aes(fill=DAM2, x=GT_Treat, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
proplot
ggsave(filename = paste0(outputdir, "Microglia_DAM2proportionsofDAM1_per_condition_WT_integrated.pdf"), proplot)
```
#### Homeostatic Cells
```{r fig.width=9, fig.height=15}
# list from Jan Hofmann
a <- RidgePlot(samples.integrated.filt, features = homeostat.genes, ncol = 2)
a
b<- FeaturePlot_scCustom(samples.integrated.filt, reduction="umap",
features = homeostat.genes, colors_use = viridis_dark_high)
b
ggsave(filename = paste0(outputdir, "Microglia_HomeostatGenes_Cluster_Ridgeplot_integrated.pdf"), a)
ggsave(filename = paste0(outputdir, "Microglia_HomeostatGenes_Featureplot_integrated.pdf"), a)
```
```{r, fig.width=14, fig.height=4}
samples.integrated.filt <- AddModuleScore(object = samples.integrated.filt, layer ="scaled.data",
features = list(homeostat.genes), name = "Homeo_score")
p4 <- FeaturePlot_scCustom(seurat_object = samples.integrated.filt,
reduction = "umap",
features = "Homeo_score1",
colors_use = viridis_dark_high)
p4a <- DimPlot(object = samples.integrated.filt, reduction = "umap")
Homeocutff = 1.5
phist = ggplot(data = [email hidden], aes(x=Homeo_score1))+
geom_density()+geom_vline(xintercept = Homeocutff)+theme_classic()
samples.integrated.filt$Homeo_binary = as.factor(samples.integrated.filt$Homeo_score1>=Homeocutff)
pbin <- DimPlot(object = samples.integrated.filt, reduction = "umap", group.by = "Homeo_binary")
comb = phist | p4 | pbin| p4a
Hom_GT_p <-DimPlot(object = samples.integrated.filt, reduction = "umap",
group.by = "Homeo_binary", split.by = "Genotype_Treatment")
comb
Hom_GT_p
prop.table(table(Homeo=samples.integrated.filt$Homeo_binary,
Clusters=samples.integrated.filt$labels), margin = 2)
ggsave(filename = paste0(outputdir, "Microglia_Homeo_calling_integrated.pdf"), comb)
ggsave(filename = paste0(outputdir, "Microglia_Homeo_bygenotype_integrated.pdf"), Hom_GT_p)
```
```{r}
phist
p4
pbin
ggsave(filename = paste0(outputdir, "Microglia_Homeo_thresholding_integrated.pdf"), phist)
ggsave(filename = paste0(outputdir, "Microglia_Homeo_scoring_integrated.pdf"), p4)
ggsave(filename = paste0(outputdir, "Microglia_Homeo_binary_integrated.pdf"), pbin)
```
#### Homeo vs DAM
```{r}
metadata = [email hidden]
metadata$Homeo_DAM <- factor(1*as.logical(metadata$Homeo_binary) + 2*as.logical(metadata$DAM_binary), levels = c(0,1,2,3), labels = c("Neither", "Homeo", "DAM", "Unknown"))
[email hidden] <- metadata
a = DimPlot(samples.integrated.filt, reduction="umap", group.by = "Homeo_DAM")
a
ggsave(filename = paste0(outputdir, "Microglia_Homeo_DAM_integrated.pdf"), a)
```
```{r}
restab = table(GT_Treat = metadata$Genotype_Treatment, Homeo=metadata$Homeo_DAM)
restab
prop.table(restab, margin = 1)
chisq.test(restab)
proplot<-ggplot(as.data.frame(restab),
aes(fill=Homeo, x=GT_Treat, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
chisq.posthoc.test(restab)
proplot
ggsave(filename = paste0(outputdir, "Microglia_HomeoDAMproportions_per_condition_integrated.pdf"), proplot)
```
```{r}
restab = table(Clust = metadata$labels, Homeo=metadata$Homeo_DAM)
restab
prop.table(restab, margin = 1)
chisq.test(restab)
proplot<-ggplot(as.data.frame(restab),
aes(fill=Homeo, x=Clust, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
chisq.posthoc.test(restab)
proplot
ggsave(filename = paste0(outputdir, "Microglia_HomeoDAMproportions_per_ckuster_integrated.pdf"), proplot)
```
```{r}
restab = table(Clust = metadata$labels, Celltype=metadata$Mouse_ID)
restab
prop.table(restab, margin = 1)
chisq.test(restab)
proplot<-ggplot(as.data.frame(restab),
aes(fill=Clust, x=Celltype, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
chisq.posthoc.test(restab)
proplot <- proplot + theme(axis.text.x = element_text(angle=90))
proplot
ggsave(filename = paste0(outputdir, "Microglia_celltype_labels_integrated.pdf"), proplot)
```
```{r}
restab = table(gender=metadata$Gender, Clust = metadata$labels)
restab
prop.table(restab, margin = 1)
chisq.test(restab)
proplot<-ggplot(as.data.frame(restab),
aes(fill=gender, x=Clust, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
chisq.posthoc.test(restab)
proplot
ggsave(filename = paste0(outputdir, "Microglia_HomeoDAMproportions_per_gender_integrated.pdf"), proplot)
```
```{r, fig.width=12}
a <- ggplot([email hidden], aes(x=Homeo_score1, y=DAM_score1))+
geom_density_2d_filled()+geom_point(col="#FFFFFF55")+theme_classic()+ggtitle("all cells")
b<- ggplot([email hidden][samples.integrated.filt$DAM_binary==TRUE,], aes(x=Homeo_score1, y=DAM_score1))+geom_density_2d_filled()+geom_point(col="#FFFFFF55")+theme_classic()+ggtitle("DAMs only")
a|b
ggsave(filename = paste0(outputdir, "Microglia_DAMvsHomeo_integrated.pdf"), a|b)
```
#### Homeo in Genotype Treatment
proportion tables show percentages of Homeostatic cells per condition
```{r}
metadata=[email hidden]
restab = table(GT_Treat = metadata$Genotype_Treatment, Homeo=metadata$Homeo_binary)
restab
prop.table(restab, margin = 1)
chisq.test(restab)
proplot<-ggplot(as.data.frame(restab),
aes(fill=Homeo, x=GT_Treat, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
chisq.posthoc.test(restab)
proplot
ggsave(filename = paste0(outputdir, "Microglia_Homeoproportions_per_condition_integrated.pdf"), proplot)
```
```{r}
metadata=[email hidden][samples.integrated.filt$Genotype=="APPPS1",]
restab = table(GT_Treat = as.character(metadata$Genotype_Treatment), Homeo=metadata$Homeo_binary)
restab
prop.table(restab, margin = 1)
chisq.test(restab)
proplot<- ggplot(as.data.frame(restab),
aes(fill=Homeo, x=GT_Treat, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
proplot
ggsave(filename = paste0(outputdir, "Microglia_Homeoproportions_per_condition_APPPS1_integrated.pdf"), proplot)
```
```{r}
metadata=[email hidden][samples.integrated.filt$Genotype=="WT",]
restab = table(GT_Treat = as.character(metadata$Genotype_Treatment), Homeo=metadata$Homeo_binary)
restab
prop.table(restab, margin = 1)
chisq.test(restab)
proplot<-ggplot(as.data.frame(restab),
aes(fill=Homeo, x=GT_Treat, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
proplot
ggsave(filename = paste0(outputdir, "Microglia_Homeoproportions_per_condition_WT_integrated.pdf"), proplot)
```
### Spatial cell types
```{r, fig.height=12,fig.width=5}
Idents(samples.integrated.filt) = "Genotype_Treatment"
DoHeatmap(samples.integrated.filt, slot = "data", features = c(distal_genes, control_genes, Stroke_genes))+ggtitle("spatial genes all")
DoHeatmap(samples.integrated.filt, slot = "data",features = c(distal_genes))+ggtitle("spatial genes distal")
DoHeatmap(samples.integrated.filt, slot = "data",features = c(control_genes))+ggtitle("spatial genes control")
DoHeatmap(samples.integrated.filt, slot = "data",features = c(Stroke_genes))+ggtitle("spatial genes peri-ictus")
```
#### Distal Cell types
```{r, fig.width=12, fig.height=4}
#load("./data/DAM_genelists.RData")
DefaultAssay(samples.integrated.filt) = "RNA"
samples.integrated.filt <- AddModuleScore(object = samples.integrated.filt,
features = list(distal_genes), name = "Distal")
p4 <- FeaturePlot_scCustom(seurat_object = samples.integrated.filt, reduction = "umap",
features = "Distal1", colors_use = viridis_dark_high)
p4a <- DimPlot(object = samples.integrated.filt, reduction = "umap")
Distalcutoff = 0.125
phist = ggplot(data = [email hidden], aes(x=Distal1))+
geom_density()+geom_vline(xintercept = Distalcutoff)+
theme_classic()
samples.integrated.filt$Distal_binary = as.factor(samples.integrated.filt$Distal1>=Distalcutoff)
pbin <- DimPlot(object = samples.integrated.filt, reduction = "umap",group.by = "Distal_binary")
phist | p4 | pbin| p4a
prop.table(table(Distals=samples.integrated.filt$Distal_binary, Clusters=samples.integrated.filt$labels), margin = 2)
DimPlot(object = samples.integrated.filt, group.by = "Distal_binary", reduction = "umap", split.by = "Genotype_Treatment")
```
#### Control Cell types
```{r, fig.width=12, fig.height=4}
#load("./data/DAM_genelists.RData")
DefaultAssay(samples.integrated.filt) = "RNA"
samples.integrated.filt <- AddModuleScore(object = samples.integrated.filt,
features = list(control_genes), name = "Control")
p4 <- FeaturePlot_scCustom(seurat_object = samples.integrated.filt, reduction = "umap",
features = "Control1", colors_use = viridis_dark_high)
p4a <- DimPlot(object = samples.integrated.filt, reduction = "umap")
Controlcutoff = 0.1
phist = ggplot(data = [email hidden], aes(x=Control1))+
geom_density()+geom_vline(xintercept = Controlcutoff)+
theme_classic()
samples.integrated.filt$Control_binary = as.factor(samples.integrated.filt$Control1>=Controlcutoff)
pbin <- DimPlot(object = samples.integrated.filt, reduction = "umap",
group.by = "Control_binary")
phist | p4 | pbin| p4a
prop.table(table(Controls=samples.integrated.filt$Control_binary, Clusters=samples.integrated.filt$labels), margin = 2)
DimPlot(object = samples.integrated.filt, group.by = "Control_binary", reduction = "umap",
split.by = "Genotype_Treatment")
```
#### Peri_ictus Cell types
```{r, fig.width=12, fig.height=4}
#load("./data/DAM_genelists.RData")
DefaultAssay(samples.integrated.filt) = "RNA"
samples.integrated.filt <- AddModuleScore(object = samples.integrated.filt,
features = list(Stroke_genes), name = "Peri_ictus")
p4 <- FeaturePlot_scCustom(seurat_object = samples.integrated.filt, reduction = "umap",features = "Peri_ictus1", colors_use = viridis_dark_high)
p4a <- DimPlot(object = samples.integrated.filt, reduction = "umap")
Strokecutoff = 1.75
phist = ggplot(data = [email hidden], aes(x=Peri_ictus1))+
geom_density()+geom_vline(xintercept = Strokecutoff)+
theme_classic()+ggtitle("Peri_ictus")
samples.integrated.filt$Peri_ictus_binary = as.factor(samples.integrated.filt$Peri_ictus1>=Strokecutoff)
pbin <- DimPlot(object = samples.integrated.filt, reduction = "umap",group.by = "Peri_ictus_binary")
phist | p4 | pbin| p4a
prop.table(table(Strokes=samples.integrated.filt$Peri_ictus_binary, Clusters=samples.integrated.filt$labels), margin = 2)
DimPlot(object = samples.integrated.filt, reduction = "umap",group.by = "Peri_ictus_binary", split.by = "Genotype_Treatment")
```
```{r}
ggplot([email hidden], aes(x=Peri_ictus1, y=Control1))+
geom_density_2d_filled()+geom_point(col="#FFFFFF55")+theme_classic()+ggtitle("all cells")
ggplot([email hidden][samples.integrated.filt$Peri_ictus_binary==TRUE,], aes(x=Peri_ictus1, y=Control1))+geom_density_2d_filled()+geom_point(col="#FFFFFF55")+theme_classic()+ggtitle("Peri_ictus only")
```
proportion tables show percebtages of DAM cells per condition
#### Distals binary proportions
```{r}
metadata=[email hidden]
restab = table(GT_Treat = metadata$Genotype_Treatment, Distals=metadata$Distal_binary)
restab
prop.table(restab, margin = 1)
chisq.test(restab)
ggplot(as.data.frame(restab),
aes(fill=Distals, x=GT_Treat, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
library(chisq.posthoc.test)
chisq.posthoc.test(restab)
#To Do Pairwise
```
#### Distals in APPPS1
```{r}
metadata=[email hidden][samples.integrated.filt$Genotype == "APPPS1", ]
restab = table(GT_Treat = as.character(metadata$Genotype_Treatment), Distals=metadata$Distal_binary)
restab
prop.table(restab, margin = 1)
fisher.test(restab)
ggplot(as.data.frame(restab),
aes(fill=Distals, x=GT_Treat, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
```
#### Distals in WT
```{r}
metadata=[email hidden][samples.integrated.filt$Genotype == "WT", ]
restab = table(GT_Treat = as.character(metadata$Genotype_Treatment), Distals=metadata$Distal_binary)
restab
prop.table(restab, margin = 1)
fisher.test(restab)
ggplot(as.data.frame(restab),
aes(fill=Distals, x=GT_Treat, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
```
#### Controls binary proportions
```{r}
metadata=[email hidden]
restab = table(GT_Treat = metadata$Genotype_Treatment, Controls=metadata$Control_binary)
restab
prop.table(restab, margin = 1)
chisq.test(restab)
ggplot(as.data.frame(restab),
aes(fill=Controls, x=GT_Treat, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
library(chisq.posthoc.test)
chisq.posthoc.test(restab)
#To Do Pairwise
```
#### Controls in APPPS1
```{r}
metadata=[email hidden][samples.integrated.filt$Genotype == "APPPS1", ]
restab = table(GT_Treat = as.character(metadata$Genotype_Treatment), Controls=metadata$Control_binary)
restab
prop.table(restab, margin = 1)
fisher.test(restab)
ggplot(as.data.frame(restab),
aes(fill=Controls, x=GT_Treat, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
```
#### Controls in WT
```{r}
metadata=[email hidden][samples.integrated.filt$Genotype == "WT", ]
restab = table(GT_Treat = as.character(metadata$Genotype_Treatment), Controls=metadata$Control_binary)
restab
prop.table(restab, margin = 1)
fisher.test(restab)
ggplot(as.data.frame(restab),
aes(fill=Controls, x=GT_Treat, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
```
#### Peri_ictus binary proportions
```{r}
metadata=[email hidden]
restab = table(GT_Treat = metadata$Genotype_Treatment, Peri_Ictal=metadata$Peri_ictus_binary)
restab
prop.table(restab, margin = 1)
chisq.test(restab)
ggplot(as.data.frame(restab),
aes(fill=Peri_Ictal, x=GT_Treat, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
library(chisq.posthoc.test)
chisq.posthoc.test(restab)
#To Do Pairwise
```
#### Peri_ictus in APPPS1
```{r}
metadata=[email hidden][samples.integrated.filt$Genotype == "APPPS1", ]
restab = table(GT_Treat = as.character(metadata$Genotype_Treatment), Peri_Ictal=metadata$Peri_ictus_binary)
restab
prop.table(restab, margin = 1)
fisher.test(restab)
ggplot(as.data.frame(restab),
aes(fill=Peri_Ictal, x=GT_Treat, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
```
#### Peri_ictus WT
```{r}
metadata=[email hidden][samples.integrated.filt$Genotype == "WT", ]
restab = table(GT_Treat = as.character(metadata$Genotype_Treatment), Peri_Ictal=metadata$Peri_ictus_binary)
restab
prop.table(restab, margin = 1)
fisher.test(restab)
ggplot(as.data.frame(restab),
aes(fill=Peri_Ictal, x=GT_Treat, y=Freq))+
geom_bar(position="fill", color="black",stat="identity")+theme_classic()
```
# DEG comparing conditions overall
## WT Stroke vs no Stroke
avg_logFC: log fold-chage of the average expression between the two groups. Positive values indicate that the gene is more highly expressed in the first group
```{r, fig.width=6, fig.height=6}
genes_to_label <- openxlsx::read.xlsx("/files/scSeq_Hefendehl/data/Volcano pot labels on graph.xlsx",
sheet = "Wildtype control vs Wildtype st")
genes_to_label <- c(genes_to_label$Up, genes_to_label$Down)
Idents(samples.integrated.filt) <- "Genotype_Treatment"
DimPlot(samples.integrated.filt, reduction="umap")
DEG_Wt_CtrlvsStroke <-FindMarkers(samples.integrated.filt,
ident.1 = "WT_Stroke", ident.2 = "WT_Ctrl",
assay = "RNA", slot = "data",
test.use = "wilcox", logfc.threshold = 0)
genes_to_label[!genes_to_label %in% rownames(DEG_Wt_CtrlvsStroke)]
```
```{r, fig.width=6, fig.height=6}
p = EnhancedVolcano::EnhancedVolcano(DEG_Wt_CtrlvsStroke,
x = "avg_log2FC",
y = "p_val_adj",
pCutoff=0.05,
lab=rownames(DEG_Wt_CtrlvsStroke),
selectLab = genes_to_label,
title = "WT_Ctrl vs WT_Stroke")
display_tab(DEG_Wt_CtrlvsStroke)
p
write.xlsx(DEG_Wt_CtrlvsStroke, file = paste0(outputdir, "Microglia_DEX_WTCtrlvsWTStroke_integrated.xlsx"),rowNames = TRUE)
ggsave(filename = paste0(outputdir, "Microglia_Volcano_WTCtrlvsWTStroke_integrated.pdf"), p)
```
## APPPS1 Stroke vs APPPS1 no Stroke
avg_logFC: log fold-chage of the average expression between the two groups. Positive values indicate that the gene is more highly expressed in the first group
```{r fig.width=6, fig.height=6}
genes_to_label <- openxlsx::read.xlsx("/files/scSeq_Hefendehl/data/Volcano pot labels on graph.xlsx",
sheet = "APPPS1 control Vs APPPS1 stroke")
genes_to_label <- c(genes_to_label$Up, genes_to_label$Down)
DEG_APPPS1_CtrlvsStroke <-FindMarkers(samples.integrated.filt,
ident.1 = "APPPS1_Stroke", ident.2 = "APPPS1_Ctrl",
assay = "RNA", slot = "data",
test.use = "wilcox", logfc.threshold = 0)
display_tab(DEG_APPPS1_CtrlvsStroke)
genes_to_label[!genes_to_label %in% rownames(DEG_APPPS1_CtrlvsStroke)]
p<- EnhancedVolcano::EnhancedVolcano(DEG_APPPS1_CtrlvsStroke,
x = "avg_log2FC",
y = "p_val_adj",
lab=rownames(DEG_APPPS1_CtrlvsStroke),
selectLab = genes_to_label,
pCutoff = 0.05,
title = "APPPS1_Ctrl vs APPPS1_Stroke")
p
write.xlsx(DEG_APPPS1_CtrlvsStroke, file = paste0(outputdir, "Microglia_DEX_APPPS1CtrlvsAPPPS1Stroke_integrated.xlsx"),rowNames = TRUE)
ggsave(filename = paste0(outputdir, "Microglia_Volcano_APPPS1CtrlvsAPPPS1Stroke_integrated.pdf"), p)
```
## WT Stroke vs APPPS1Stroke
```{r fig.width=6, fig.height=6}
# plot here only DAM markers
genes_to_label <- unique(c(DAM_up, DAM1_down, DAM1_up, DAM1vsDAM2, DAM2_up))
DEG_Stroke_APPPS1vsWT <-FindMarkers(samples.integrated.filt,
ident.1 = "APPPS1_Stroke", ident.2 = "WT_Stroke",
assay = "RNA", slot = "data",
test.use = "wilcox", logfc.threshold = 0)
display_tab(DEG_Stroke_APPPS1vsWT)
genes_to_label[!genes_to_label %in% rownames(DEG_Stroke_APPPS1vsWT)]
p<- EnhancedVolcano::EnhancedVolcano(DEG_Stroke_APPPS1vsWT,
x = "avg_log2FC",
y = "p_val_adj",
lab=rownames(DEG_Stroke_APPPS1vsWT),
selectLab = genes_to_label,
pCutoff = 0.05,
title = "WT_stroke vs APPPS1_Stroke")
p
write.xlsx(DEG_Stroke_APPPS1vsWT, file = paste0(outputdir, "Microglia_DEX_Stroke_APPPS1vsWT_integrated.xlsx"),rowNames = TRUE)
ggsave(filename = paste0(outputdir, "Microglia_Volcano_WTstrokevsAPPPS1StrokeDAMgenes_integrated.pdf"), p)
```
## WT ctr; vs APPPS1 ctrl
avg_logFC: log fold-chage of the average expression between the two groups. Positive values indicate that the gene is more highly expressed in the first group
```{r fig.width=6, fig.height=6}
genes_to_label <- openxlsx::read.xlsx("/files/scSeq_Hefendehl/data/Volcano pot labels on graph.xlsx",
sheet = "Wild type control Vs APPPS1 con")
genes_to_label <- c(genes_to_label$Up, genes_to_label$Down)
DEG_APPPS1_CtrlvsWT_Ctrl <-FindMarkers(samples.integrated.filt,
ident.1 = "APPPS1_Ctrl", ident.2 = "WT_Ctrl",
assay = "RNA", slot = "data",
test.use = "wilcox", logfc.threshold = 0)
display_tab(DEG_APPPS1_CtrlvsWT_Ctrl)
genes_to_label[!genes_to_label %in% rownames(DEG_APPPS1_CtrlvsWT_Ctrl)]
p<- EnhancedVolcano::EnhancedVolcano(DEG_APPPS1_CtrlvsWT_Ctrl,
x = "avg_log2FC",
y = "p_val_adj",
lab=rownames(DEG_APPPS1_CtrlvsWT_Ctrl),
selectLab = genes_to_label,
pCutoff = 0.05,
title = "WT_Ctrl vs APPPS1_Ctrl")
p
write.xlsx(DEG_APPPS1_CtrlvsWT_Ctrl, file = paste0(outputdir, "Microglia_DEX_APPPS1_CtrlvsWT_Ctrl_integrated.xlsx"),rowNames = TRUE)
ggsave(filename = paste0(outputdir, "Microglia_Volcano_WTCtrlvsAPPPS1ctrl_integrated.pdf"), p)
```
## Scatterplot of DEGs of all four conditions
```{r}
genelist_both= intersect(rownames(DEG_Wt_CtrlvsStroke), rownames(DEG_APPPS1_CtrlvsStroke))
plotdata=data.frame(log2FC_WTCtrl_vs_WTStroke=DEG_Wt_CtrlvsStroke[genelist_both,]$avg_log2FC,
padj_WTCtrl_vs_WTStroke=DEG_Wt_CtrlvsStroke[genelist_both,]$p_val,
log2FC_APPPS1Ctrl_vs_APPPS1Stroke = DEG_APPPS1_CtrlvsStroke[genelist_both,]$avg_log2FC,
padj_APPPS1Ctrl_vs_APPPS1Stroke = DEG_APPPS1_CtrlvsStroke[genelist_both,]$p_val) %>%
mutate(
collab=ifelse((padj_WTCtrl_vs_WTStroke < 0.05) & (padj_APPPS1Ctrl_vs_APPPS1Stroke < 0.05), "Both",
ifelse((padj_WTCtrl_vs_WTStroke > 0.05) & (padj_APPPS1Ctrl_vs_APPPS1Stroke < 0.05), "APPPS1Stroke",
ifelse((padj_WTCtrl_vs_WTStroke < 0.05) & (padj_APPPS1Ctrl_vs_APPPS1Stroke > 0.05), "WTStroke", NA)))
)
rownames(plotdata) = genelist_both
plotdata$label = rownames(plotdata)
# the next lines plot the top 20 genes of either analyses
n=20
plotdata$label[! c(1:nrow(plotdata)) %in% c(
order(abs(plotdata$log2FC_WTCtrl_vs_WTStroke), decreasing = T)[1:n],
order(abs(plotdata$log2FC_APPPS1Ctrl_vs_APPPS1Stroke), decreasing = T)[1:n])]=""
#if you execute the next two lines the labels you choose will be plotted
# to execute set FALSE to TRUE and change the list of gene names
if(TRUE){
genestoplot = c("S100a6", "Apoc4", "C4b", "Atp1a2", "Apoc1", "Fn1","Ntrk2")
print(paste0(genestoplot[! genestoplot %in% rownames(plotdata)], " cannot be found"))
plotdata$label[! rownames(plotdata) %in% genestoplot] = ""
}
```
```{r}
plotdata <- plotdata %>%
arrange(desc(is.na(collab)))
p <- ggplot(plotdata, aes(x=log2FC_WTCtrl_vs_WTStroke, y=log2FC_APPPS1Ctrl_vs_APPPS1Stroke, label=label))+
geom_point(aes(colour = collab), alpha=0.7)+
geom_text_repel(box.padding = 0.5, max.overlaps = Inf, color="black")+
scale_color_discrete(palette=c( "#56B4E9", "#009E73","#F0E442"))+
#geom_abline(slope=-1, intercept=0, lty=2)+geom_abline(slope=1, intercept=0, lty=2)+
theme_classic()+ggtitle("log_FC")
p
ggsave(filename = paste0(outputdir, "Microglia_logFCCtrlvsStroke_APPPS1vsWT_all_integrated.pdf"), p)
```
```{r, fig.width=9, fig.height=7}
# sig only
genelist_both=
unique(c(rownames(DEG_Wt_CtrlvsStroke %>% dplyr::filter(p_val<0.05 & abs(avg_log2FC)>0.5)),
rownames(DEG_APPPS1_CtrlvsStroke%>% dplyr::filter(p_val<0.05 & abs(avg_log2FC)>0.5))))
plotdata=data.frame(log2FC_WTCtrl_vs_WTStroke=DEG_Wt_CtrlvsStroke[genelist_both,]$avg_log2FC,
padj_WTCtrl_vs_WTStroke=DEG_Wt_CtrlvsStroke[genelist_both,]$p_val,
log2FC_APPPS1Ctrl_vs_APPPS1Stroke = DEG_APPPS1_CtrlvsStroke[genelist_both,]$avg_log2FC,
padj_APPPS1Ctrl_vs_APPPS1Stroke = DEG_APPPS1_CtrlvsStroke[genelist_both,]$p_val) %>%
mutate(
collab=ifelse((padj_WTCtrl_vs_WTStroke < 0.05) & (padj_APPPS1Ctrl_vs_APPPS1Stroke < 0.05), "Both",
ifelse((padj_WTCtrl_vs_WTStroke > 0.05) & (padj_APPPS1Ctrl_vs_APPPS1Stroke < 0.05), "APPPS1Stroke",
ifelse((padj_WTCtrl_vs_WTStroke < 0.05) & (padj_APPPS1Ctrl_vs_APPPS1Stroke > 0.05), "WTStroke", NA)))
)
rownames(plotdata) = genelist_both
plotdata$label = rownames(plotdata)
#if you execute the next lines the labels you choose will be plotted
# if set to TRUE, change the list of gene names and they will be plotted
# if set to FALSE the top n=15 of both will be plotted
mode = "both" # alternative "selected" or "topN" or "both"
genestoplot = c("S100a6", "Apoc4", "C4b", "Atp1a2", "Apoc1", "Fn1","Ntrk2") # name here your gene of interest
topn = 15 # name the top number of genes to be plotted
### don't change anything below the line here
plotdata$label = rownames(plotdata)
if(mode=="selected"){
print(paste0(genestoplot[! genestoplot %in% rownames(plotdata)], " cannot be found"))
plotdata$label[which(! rownames(plotdata) %in% genestoplot)] = ""
} else if (mode=="topN"){
n=topn
genestoplotn15 <-c(1:nrow(plotdata)) %in%
c(
order(abs(plotdata$log2FC_WTCtrl_vs_WTStroke), decreasing = T)[1:n],
order(abs(plotdata$log2FC_APPPS1Ctrl_vs_APPPS1Stroke), decreasing = T)[1:n])
plotdata$label[!genestoplotn15]=""
} else if (mode =="both") {
print(paste0(genestoplot[! genestoplot %in% rownames(plotdata)], " cannot be found"))
genestoplotcmbn <-(c(1:nrow(plotdata)) %in%
c(
order(abs(plotdata$log2FC_WTCtrl_vs_WTStroke), decreasing = T)[1:n],
order(abs(plotdata$log2FC_APPPS1Ctrl_vs_APPPS1Stroke), decreasing = T)[1:n])) | (rownames(plotdata) %in% genestoplot)
plotdata$label[which(! genestoplotcmbn)] = ""
}
```
```{r, fig.width=9, fig.height=7}
p <- ggplot(plotdata, aes(x=log2FC_WTCtrl_vs_WTStroke, y=log2FC_APPPS1Ctrl_vs_APPPS1Stroke, label=label))+
geom_point(aes(colour = collab), alpha=0.7)+
geom_text_repel(box.padding = 0.5, max.overlaps = Inf,
color="black")+
scale_color_discrete(palette=c( "#56B4E9", "#009E73","#F0E442"))+
#geom_abline(slope=-1, intercept=0, lty=2)+geom_abline(slope=1, intercept=0, lty=2)+
theme_classic()+ggtitle("log_FC")
p
ggsave(filename = paste0(outputdir, "Microglia_logFCCtrlvsStroke_APPPS1vsWT_sig_only_integrated.pdf"), p)
```
```{r fig.width=8, fig.height=8}
reads_data = JoinLayers(samples.integrated.filt, assay = "RNA") %>%
GetAssayData(assay = "RNA", layer="data") # loaded data as used above
#reads_data = samples.integrated.filt@assays$RNA@counts # Raw counts/reads
#reads_data = samples.integrated.filt@assays$[email hidden] #scaled/reads
barplot_genes = function(gene, split.by =T){
require(ggbeeswarm)
plotdata = data.frame(log2Reads = log2(reads_data[gene,]+1),
Genotype_Treatment = [email hidden]$Genotype_Treatment,
Cluster = [email hidden]$labels)
if(split.by=="none" | split.by=="cluster") {
p=ggplot(plotdata, aes(x=Genotype_Treatment, y=log2Reads, fill=Genotype_Treatment))+
geom_violin()+
geom_beeswarm(alpha=0.4, method = "hex")+ggtitle(gene)+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
if(split.by=="cluster"){p=p+facet_wrap(~Cluster)}
} else if (split.by=="Microglia") {
p=ggplot(plotdata, aes(x=Cluster, y=log2Reads, fill=Cluster))+
geom_violin()+
geom_beeswarm(alpha=0.4, method = "hex")+
ggtitle(gene)+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
}
p
}
barplot_genes(gene="Ptprc", split.by ="cluster") # alternative split.by = "none" oder " Micorglia
barplot_genes(gene="Ptprc", split.by ="none") # just across the genotype_treatments
barplot_genes(gene="Ptprc", split.by ="Microglia") # add +coord_cartesian(ylim=c(0, 7)) if you want to change the ylims
barplot_genes(gene="Ptprc", split.by ="Microglia")+coord_cartesian(ylim=c(0, 2))
barplot_genes(gene="Ptprc", split.by ="Microglia") + facet_wrap(~Genotype_Treatment)
FeaturePlot_scCustom(samples.integrated.filt, features = sort(genestoplot), reduction="umap", colors_use = viridis_dark_high)
```
## DEG Stroke Cells only
```{r, eval=FALSE}
DEG_Strokebins_Wt_CtrlvsStroke <-FindMarkers(samples.integrated.filt, subset.ident = "Stroke_bin", subset="TRUE",
assay = "RNA", slot = "data",
ident.1 = "WT_Stroke", ident.2 = "WT_Ctrl",
test.use = "wilcox", logfc.threshold = 0)
display_tab(DEG_Strokebins_Wt_CtrlvsStroke)
DEG_Strokebins_Appps1_CtrlvsStroke <-FindMarkers(samples.integrated.filt, subset.ident = "Stroke_bin", subset="TRUE",
assay = "RNA", slot = "data",
ident.1 = "APPPS1_Stroke", ident.2 = "APPPS1_Ctrl",
test.use = "wilcox", logfc.threshold = 0)
display_tab(DEG_Strokebins_Appps1_CtrlvsStroke)
genelist_both= intersect(rownames(DEG_Strokebins_Appps1_CtrlvsStroke), rownames(DEG_Strokebins_Wt_CtrlvsStroke))
plotdata=data.frame(WT=DEG_Strokebins_Wt_CtrlvsStroke[genelist_both,]$avg_log2FC, APPPS1 = DEG_Strokebins_Appps1_CtrlvsStroke[genelist_both,]$avg_log2FC)
rownames(plotdata) = genelist_both
plotdata$label = rownames(plotdata)
n=5
plotdata$label[! c(1:nrow(plotdata)) %in%
c(
order(abs(plotdata$WT), decreasing = T)[1:n],
order(abs(plotdata$APPPS1), decreasing = T)[1:n])]=""
plotdata$col="WT"
apps1=abs(plotdata$APPPS1)>abs(plotdata$WT)
plotdata$col[ apps1]="APPPS1"
plotdata$col[(abs(plotdata$APPPS1)<0.25) & (abs(plotdata$WT)<0.25)] <- NA
ggplot(plotdata, aes(x=WT, y=APPPS1, color=col, label=label))+geom_point(alpha=0.7)+
geom_text(vjust=-0.7, color="black")+geom_abline(slope=-1, intercept=0, lty=2)+geom_abline(slope=1, intercept=0, lty=2)+theme_classic()+ggtitle("log_FC_Stroke_bins")
```
## DEG comparing Cluster
All Clusters versus All Genes versus all Conditions
### Define Main analysis pipeline
```{r}
compareCluster <- function(
ClusterA = "Microglia6",
ClusterB = "Microglia0",
seuratObject = samples.integrated.filt,
output = "./",
genelist = NULL,
FCcutoff = 1,
comorbid = FALSE,
evgenecode=F
) {
# ensure output exists
if (!dir.exists(output)) dir.create(outputdir, recursive = TRUE)
# name + logfile ALWAYS defined
resname <- if (isTRUE(comorbid)) {
paste0(ClusterA, "vs", ClusterB, "_comorbid")
} else {
paste0(ClusterA, "vs", ClusterB)
}
logfile <- file.path(output, paste0("log_", resname, ".txt"))
if (file.exists(logfile)) file.remove(logfile)
# helper to append to log
log_line <- function(level, msg) {
cat(
paste0(
format(Sys.time(), "%Y-%m-%d %H:%M:%S"),
" | ", level,
" | ", ClusterA, " vs ", ClusterB,
" | ", msg, "\n"
),
file = logfile,
append = TRUE
)
}
# FindMarkers call (different args if comorbid)
DEG_tmp <- tryCatch(
withCallingHandlers(
{
if (isTRUE(comorbid)) {
FindMarkers(
seuratObject,
assay = "RNA",
subset.ident = "Genotype_Treatment",
subset = "APPPS1_Stroke",
layer = "data",
ident.1 = ClusterA, ident.2 = ClusterB,
test.use = "wilcox", logfc.threshold = 0
)
} else {
FindMarkers(
seuratObject,
assay = "RNA",
layer = "data",
ident.1 = ClusterA, ident.2 = ClusterB,
test.use = "wilcox", logfc.threshold = 0
)
}
},
warning = function(w) {
log_line("WARNING", conditionMessage(w))
invokeRestart("muffleWarning")
}
),
error = function(e) {
log_line("ERROR", conditionMessage(e))
return(NULL)
}
)
# If FindMarkers failed, still return log + empty slots
if (is.null(DEG_tmp)) {
return(list(
Markers = NULL,
VulcanoPlot = NULL,
qqplot = NULL,
Goresults = NULL,
GOplots = NULL,
log = logfile
))
}
openxlsx::write.xlsx(
DEG_tmp %>% dplyr::filter(p_val_adj<0.1),
file = file.path(output, paste0("DEG_Gene_list_", resname, ".xlsx")),
rowNames = TRUE
)
# label genes that actually exist in DEG table
genes_to_label2 <- intersect(genelist, rownames(DEG_tmp))
q <- gg_qqplot(DEG_tmp$p_val)
p <- EnhancedVolcano::EnhancedVolcano(
DEG_tmp,
x = "avg_log2FC",
y = "p_val_adj",
lab = rownames(DEG_tmp),
selectLab = genes_to_label2,
pCutoff = 0.05,
title = resname
)
ggsave(
filename = file.path(output, paste0("DEG_Vulcano_", resname, "_integrated.pdf")),
plot = (p | q),
width = 10, height = 5
)
gene_univers <- rownames(seuratObject)
gogenes <- list(
genesigall = rownames(DEG_tmp)[DEG_tmp$p_val_adj < 0.05 & abs(DEG_tmp$avg_log2FC) > FCcutoff],
genesigup = rownames(DEG_tmp)[DEG_tmp$p_val_adj < 0.05 & DEG_tmp$avg_log2FC > FCcutoff],
genesigdown = rownames(DEG_tmp)[DEG_tmp$p_val_adj < 0.05 & DEG_tmp$avg_log2FC < -FCcutoff]
)
Gores <- getGOresults(
geneset = gogenes,
domain_scope = "known",
genereference = gene_univers,
return_genelist = evgenecode,
organism = "mmusculus"
)
openxlsx::write.xlsx(
Gores$result %>% dplyr::select(-parents),
file = file.path(output, paste0("DEG_GOlist_", resname, ".xlsx")),
rowNames = TRUE
)
plotlist <- list()
for (m in unique(Gores$result$query)) {
idx <- Gores$result$query == m & grepl("GO", Gores$result$source)
if (sum(idx) >= 1) {
plotlist[[m]] <- GOplot(Gores$result[idx, ], N = 10, Title = m)
} else {
plotlist[[m]] <- ggplot() +
annotate(
"text", x = 10, y = 10, size = 6,
label = "no significant GO term\ncheck excel sheet for other enrichments"
) +
theme_void() +
ggtitle(m)
}
}
ggsave(
filename = file.path(output, paste0("DEG_GO_plots_", resname, ".pdf")),
plot = ggarrange(plotlist = plotlist),
width = 10, height = 5
)
return(list(
Markers = DEG_tmp,
VulcanoPlot = p,
qqplot = q,
Goresults = Gores$result,
GOplots = plotlist,
log = logfile
))
}
```
### DEG run for each combination
```{r, warning=FALSE, message=FALSE}
if(T){
Idents(samples.integrated.filt) <- "labels"
labels_unique = sort(unique(samples.integrated.filt$labels))
label_pairs <- combn(labels_unique, 2, simplify = FALSE)
results_list <- list()
for (pair in label_pairs) {
print(paste(pair, collapse = "_vs_"))
print("all")
res <- compareCluster(
ClusterA = pair[1],
ClusterB = pair[2],
FCcutoff = 0.5,
output = outputdir,
genelist = genes_to_label
)
print("comorbid")
res_comorbid <- compareCluster(
ClusterA = pair[1],
ClusterB = pair[2],
FCcutoff = 0.5,
output = outputdir,
genelist = genes_to_label,
comorbid=T
)
results_list[[paste(pair, collapse = "_vs_")]] <- res
results_list[[paste0(paste(pair, collapse = "_vs_"),
"_comorbid")]] <- res_comorbid
}
}
```
```{r, results='asis', message=FALSE, warning=FALSE, fig.width=10, fig.height=5}
# results_list: named list where each element is res from compareCluster()
for (nm in names(results_list)) {
res <- results_list[[nm]]
cat("\n\n## ", nm, "\n\n", sep = "")
# --- 1) Volcano + QQ (patchwork-style)
# If you're using patchwork, this should work:
if (!is.null(res$VulcanoPlot) && !is.null(res$qqplot)) {
p1 <- res$VulcanoPlot | res$qqplot
print(p1)
} else {
cat("*No VulcanoPlot and/or qqplot found for this comparison.*\n\n")
}
# --- 2) GO plots (arranged)
if (!is.null(res$GOplots) && length(res$GOplots) > 0) {
p2 <- ggarrange(plotlist = res$GOplots, ncol = 2)
print(p2)
} else {
cat("*No GOplots found for this comparison.*\n\n")
}
# --- 3) Markers table
cat("\n\n### Markers\n\n")
if (!is.null(res$Markers) && nrow(res$Markers) > 0) {
print(kable(res$Markers %>% arrange(p_val) %>% slice_head(n = 100), format = "markdown"))
} else {
cat("*No markers found for this comparison.*\n\n")
}
}
```
### within cluster Wt vs WT Stroke
```{r, results='asis', message=FALSE, warning=FALSE, fig.width=10, fig.height=5}
Idents(samples.integrated.filt) <- "Genotype_Treatment"
for(clust in unique(samples.integrated.filt$labels)){
print(clust)
FCcutoff = 0.5
resname=paste0(clust, "_WT_Ctrl_vs_WT_Stroke_")
DEG_tmp<- FindMarkers(samples.integrated.filt[,samples.integrated.filt$labels==clust],
assay = "RNA",
layer = "data",
ident.1 = "WT_Ctrl",
ident.2 = "WT_Stroke",
test.use = "wilcox", logfc.threshold = 0
)
res_sig <- DEG_tmp %>% dplyr::filter(p_val<0.05)
openxlsx::write.xlsx(res_sig,
file = file.path(outputdir, paste0("DEG_Gene_list_", resname, ".xlsx")),
rowNames = TRUE
)
q <- gg_qqplot(DEG_tmp$p_val)
p <- EnhancedVolcano::EnhancedVolcano(
DEG_tmp,
x = "avg_log2FC",
y = "p_val_adj",
lab = rownames(DEG_tmp),
selectLab = genes_to_label,
pCutoff = 0.05,
title = resname
)
ggsave(
filename = file.path(outputdir, paste0("DEG_Vulcano_", resname, "integrated.pdf")),
plot = (p | q),
width = 10, height = 5
)
gene_univers <- rownames(samples.integrated.filt)
gogenes <- list(
genesigall = rownames(DEG_tmp)[DEG_tmp$p_val_adj < 0.05 & abs(DEG_tmp$avg_log2FC) > FCcutoff],
genesigup = rownames(DEG_tmp)[DEG_tmp$p_val_adj < 0.05 & DEG_tmp$avg_log2FC > FCcutoff],
genesigdown = rownames(DEG_tmp)[DEG_tmp$p_val_adj < 0.05 & DEG_tmp$avg_log2FC < -FCcutoff]
)
if(any(lapply(gogenes, length)>0)){
Gores <- getGOresults(
geneset = gogenes,
domain_scope = "known",
genereference = gene_univers,
return_genelist = F,
organism = "mmusculus"
)
} else {
print("no significant genes identified")
Gores <- NULL
}
if(length(Gores)!=0){
openxlsx::write.xlsx(
Gores$result %>% dplyr::select(-parents),
file = file.path(outputdir, paste0("DEG_GOlist_", resname, ".xlsx")),
rowNames = TRUE
)
plotlist <- list()
for (m in unique(Gores$result$query)) {
idx <- Gores$result$query == m & grepl("GO", Gores$result$source)
if (sum(idx) >= 1) {
plotlist[[m]] <- GOplot(Gores$result[idx, ], N = 10, Title = m)
} else {
plotlist[[m]] <- ggplot() +
annotate(
"text", x = 10, y = 10, size = 6,
label = "no significant GO term\ncheck excel sheet for other enrichments"
) +
theme_void() +
ggtitle(m)
}
}
ggsave(
filename = file.path(outputdir, paste0("DEG_GO_plots_", resname, ".pdf")),
plot = ggarrange(plotlist = plotlist),
width = 10, height = 5
) } else {
openxlsx::write.xlsx(data.frame(resname="No results identified"),
file = file.path(outputdir, paste0("DEG_GOlist_", resname, ".xlsx")))
p <- ggplot() +
annotate(
"text", x = 10, y = 10, size = 6,
label = "no significant enrichments"
) +
theme_void()
ggsave(
filename = file.path(outputdir, paste0("DEG_GO_plots_", resname, ".pdf")),
plot = p,
width = 10, height = 5
)
}
}
Idents(samples.integrated.filt) <- "labels"
```
### within cluster APPPS1_Stroke vs WT Stroke
```{r, results='asis', message=FALSE, warning=FALSE, fig.width=10, fig.height=5}
Idents(samples.integrated.filt) <- "Genotype_Treatment"
for(clust in unique(samples.integrated.filt$labels)){
print(clust)
FCcutoff = 0.5
resname=paste0(clust, "_WT_stroke_vs_APPPS1_Stroke_")
DEG_tmp<- FindMarkers(samples.integrated.filt[,samples.integrated.filt$labels==clust],
assay = "RNA",
layer = "data",
ident.1 = "APPPS1_Stroke",
ident.2 = "WT_Stroke",
test.use = "wilcox", logfc.threshold = 0
)
res_sig <- DEG_tmp %>% dplyr::filter(p_val<0.05)
openxlsx::write.xlsx(res_sig,
file = file.path(outputdir, paste0("DEG_Gene_list_", resname, ".xlsx")),
rowNames = TRUE
)
q <- gg_qqplot(DEG_tmp$p_val)
p <- EnhancedVolcano::EnhancedVolcano(
DEG_tmp,
x = "avg_log2FC",
y = "p_val_adj",
lab = rownames(DEG_tmp),
selectLab = genes_to_label,
pCutoff = 0.05,
title = resname
)
ggsave(
filename = file.path(outputdir, paste0("DEG_Vulcano_", resname, "integrated.pdf")),
plot = (p | q),
width = 10, height = 5
)
gene_univers <- rownames(samples.integrated.filt)
gogenes <- list(
genesigall = rownames(DEG_tmp)[DEG_tmp$p_val_adj < 0.05 & abs(DEG_tmp$avg_log2FC) > FCcutoff],
genesigup = rownames(DEG_tmp)[DEG_tmp$p_val_adj < 0.05 & DEG_tmp$avg_log2FC > FCcutoff],
genesigdown = rownames(DEG_tmp)[DEG_tmp$p_val_adj < 0.05 & DEG_tmp$avg_log2FC < -FCcutoff]
)
if(any(lapply(gogenes, length)>0)){
Gores <- getGOresults(
geneset = gogenes,
domain_scope = "known",
genereference = gene_univers,
return_genelist = F,
organism = "mmusculus"
)
} else {
print("no significant genes identified")
Gores <- NULL
}
if(length(Gores)!=0){
openxlsx::write.xlsx(
Gores$result %>% dplyr::select(-parents),
file = file.path(outputdir, paste0("DEG_GOlist_", resname, ".xlsx")),
rowNames = TRUE
)
plotlist <- list()
for (m in unique(Gores$result$query)) {
idx <- Gores$result$query == m & grepl("GO", Gores$result$source)
if (sum(idx) >= 1) {
plotlist[[m]] <- GOplot(Gores$result[idx, ], N = 10, Title = m)
} else {
plotlist[[m]] <- ggplot() +
annotate(
"text", x = 10, y = 10, size = 6,
label = "no significant GO term\ncheck excel sheet for other enrichments"
) +
theme_void() +
ggtitle(m)
}
}
ggsave(
filename = file.path(outputdir, paste0("DEG_GO_plots_", resname, ".pdf")),
plot = ggarrange(plotlist = plotlist),
width = 10, height = 5
) } else {
openxlsx::write.xlsx(data.frame(resname="No results identified"),
file = file.path(outputdir, paste0("DEG_GOlist_", resname, ".xlsx")))
p <- ggplot() +
annotate(
"text", x = 10, y = 10, size = 6,
label = "no significant enrichments"
) +
theme_void()
ggsave(
filename = file.path(outputdir, paste0("DEG_GO_plots_", resname, ".pdf")),
plot = p,
width = 10, height = 5
)
}
}
Idents(samples.integrated.filt) <- "labels"
```
### within cluster APPPS1 vs APPPS1 Stroke
```{r, results='asis', message=FALSE, warning=FALSE, fig.width=10, fig.height=5}
Idents(samples.integrated.filt) <- "Genotype_Treatment"
for(clust in unique(samples.integrated.filt$labels)){
print(clust)
FCcutoff = 0.5
resname=paste0(clust, "_APPPS1_Stroke_vs_APPPS1_Ctrl_")
DEG_tmp<- FindMarkers(samples.integrated.filt[,samples.integrated.filt$labels==clust],
assay = "RNA",
subset.ident = "labels",
subset = clust,
layer = "data",
ident.1 = "APPPS1_Stroke", ident.2 = "APPPS1_Ctrl",
test.use = "wilcox", logfc.threshold = 0
)
res_sig <- DEG_tmp %>% dplyr::filter(p_val<0.05)
openxlsx::write.xlsx(res_sig,
file = file.path(outputdir, paste0("DEG_Gene_list_", resname, ".xlsx")),
rowNames = TRUE
)
q <- gg_qqplot(DEG_tmp$p_val)
p <- EnhancedVolcano::EnhancedVolcano(
DEG_tmp,
x = "avg_log2FC",
y = "p_val_adj",
lab = rownames(DEG_tmp),
selectLab = genes_to_label,
pCutoff = 0.05,
title = resname
)
ggsave(
filename = file.path(outputdir, paste0("DEG_Vulcano_", resname, "integrated.pdf")),
plot = (p | q),
width = 10, height = 5
)
gene_univers <- rownames(samples.integrated.filt)
gogenes <- list(
genesigall = rownames(DEG_tmp)[DEG_tmp$p_val_adj < 0.05 & abs(DEG_tmp$avg_log2FC) > FCcutoff],
genesigup = rownames(DEG_tmp)[DEG_tmp$p_val_adj < 0.05 & DEG_tmp$avg_log2FC > FCcutoff],
genesigdown = rownames(DEG_tmp)[DEG_tmp$p_val_adj < 0.05 & DEG_tmp$avg_log2FC < -FCcutoff]
)
if(any(lapply(gogenes, length)>0)){
Gores <- getGOresults(
geneset = gogenes,
domain_scope = "known",
genereference = gene_univers,
return_genelist = F,
organism = "mmusculus"
)
} else {
print("no significant genes identified")
Gores <- NULL
}
if(length(Gores)!=0){
openxlsx::write.xlsx(
Gores$result %>% dplyr::select(-parents),
file = file.path(outputdir, paste0("DEG_GOlist_", resname, ".xlsx")),
rowNames = TRUE
)
plotlist <- list()
for (m in unique(Gores$result$query)) {
idx <- Gores$result$query == m & grepl("GO", Gores$result$source)
if (sum(idx) >= 1) {
plotlist[[m]] <- GOplot(Gores$result[idx, ], N = 10, Title = m)
} else {
plotlist[[m]] <- ggplot() +
annotate(
"text", x = 10, y = 10, size = 6,
label = "no significant GO term\ncheck excel sheet for other enrichments"
) +
theme_void() +
ggtitle(m)
}
}
ggsave(
filename = file.path(outputdir, paste0("DEG_GO_plots_", resname, ".pdf")),
plot = ggarrange(plotlist = plotlist),
width = 10, height = 5
) } else {
openxlsx::write.xlsx(data.frame(resname="No results identified"),
file = file.path(outputdir, paste0("DEG_GOlist_", resname, ".xlsx")))
p <- ggplot() +
annotate(
"text", x = 10, y = 10, size = 6,
label = "no significant enrichments"
) +
theme_void()
ggsave(
filename = file.path(outputdir, paste0("DEG_GO_plots_", resname, ".pdf")),
plot = p,
width = 10, height = 5
)
}
}
Idents(samples.integrated.filt) <- "labels"
```
```{r}
rm(list=c("seuratObject", "sce.filt", "reads", "reads_data", "comb",
"ref.sce.mm",
"e2",
"res",
"res_comorbid",
"samples.integrated"))
gc()
```
# Pseudotime analysis
## slinghsot starting Microglia 0
```{r, fig.width==7, fig.height=7}
samples.integrated.filt.meta = [email hidden] %>% dplyr::select(- contains("Pseudo"))
[email hidden] <- samples.integrated.filt.meta
sce.microglia <- as.SingleCellExperiment(samples.integrated.filt, assay = "RNA")
sce.microglia <- slingshot(sce.microglia,
clusterLabels = 'labels',
reducedDim = 'UMAP',
start.clus="Microglia0",
end.clus=c("Microglia3", "Microglia4", "Microglia5"),
omega=T,
omega_scale=4,
extend="n")
pt_col <- grep("slingPseudotime", names(colData(sce.microglia)), value = TRUE)
metasce <- data.frame(colData(sce.microglia)[,pt_col],
row.names = rownames(colData(sce.microglia)))
colnames(metasce) <- pt_col
samples.integrated.filt.meta <- cbind(samples.integrated.filt.meta,metasce)
samples.integrated.filt.meta -> [email hidden]
dataset = data.frame(reducedDim(sce.microglia, "UMAP"),
cluster = as.character(sce.microglia$labels),
DAMs = as.character(sce.microglia$DAM_binary))
curves <- slingCurves(sce.microglia, as.df=T)
p = ggplot(dataset, aes(x=umap_1, y=umap_2))+
geom_point(aes(col = cluster), size=2)+
geom_point(data=subset(dataset, dataset$DAMs==T), aes(fill=DAMs), pch=21, size=0.7)+
scale_fill_manual(values="black")
curves = curves %>% arrange(Lineage, Order)
p = p + geom_path(data = curves, aes(group = Lineage), lwd=1, lineend = "butt",)+theme_classic()
curves_arrows = curves %>% group_by(Lineage) %>% mutate(xend = lead(umap_1), yend = lead(umap_2))
curves_arrows = curves_arrows[curves_arrows$Order %in% seq(10,150, by=floor((150-10)/10)),]
p = p + geom_segment(data = curves_arrows,
aes(x = umap_1, y = umap_2, xend = xend, yend = yend),
linewidth=0.5,
arrow = arrow(length=unit(0.3, "cm"), type="closed"))
p
ggsave(filename = paste0(outputdir, "Pseudotime_linages_startMicroglia0_integrated.pdf"), p)
write_rds(p, paste0(outputdir, "Pseudotime_linages_startMicroglia0_integrated.rds"))
```
```{r}
library(TSCAN)
n.top.gene = 10
tscan_res <- list()
for (pt in pt_col){
print(pt)
tscan_res [[pt]]<- testPseudotime(sce.microglia, pseudotime = sce.microglia[[pt]])
}
```
Look for genes that are up/down-regulated along pseudotime, for the re-clustered cells (Same paper Fig 2G/H)
## Trajectory genes
By default, estimates of the spline coefficients are not returned as they are difficult to interpret. Rather, a log-fold change of expression along each path is estimated to provide some indication of the overall magnitude and direction of any change.
```{r}
reslist <- list()
reslist_feat <- list()
for (pt in names(tscan_res)){
scePT <- tscan_res[[pt]]
scePT<- scePT[order(scePT$FDR, decreasing = F),]
display_tab(as.data.frame(scePT)[1:20,])
reslist_feat[[pt]] <- FeaturePlot_scCustom(samples.integrated.filt, features = pt, reduction = "umap", colors_use = viridis_dark_high)
write.xlsx(as.data.frame(scePT), rowNames = TRUE,
file = paste0(outputdir, pt, "_genes.xslx"))
PT1_up = scePT[scePT$logFC>0,]
PT1_low = scePT[scePT$logFC<0,]
p <- scater::plotExpression(sce.microglia, features=rownames(PT1_up)[1:n.top.gene], x=pt, colour_by="labels")+scale_color_hue() + ggtitle(paste(pt, "up reg"))
p2 <- scater::plotExpression(sce.microglia, features=rownames(PT1_low)[1:n.top.gene], x=pt, colour_by="labels")+scale_color_hue()+ ggtitle(paste(pt, "down reg"))
reslist[[pt]]=p|p2
}
```
```{r, fig.height=15, fig.width=36}
a <- ggarrange(plotlist = reslist)
a
ggsave(filename = paste0(outputdir, "Pseudotime_Top10up_down_reg_Genes-correct_integrated.pdf"), a)
```
```{r, fig.height=8, fig.width=10}
a <- ggarrange(plotlist = reslist_feat)
a
ggsave(filename = paste0(outputdir, "Pseudotime_Umap_Genes-correct_integrated.pdf"), a)
```
# WGCNA
In the paper <https://www.nature.com/articles/s41586-018-0023-4> this analysis was used to generate different modules of gene lists that can then be mapped onto different treatment/genotypes within one experiment (see Figure 4). This is an analysis that can be done in Seurat <https://smorabit.github.io/tutorials/9_scWGCNA_tutorial/> Which could be very cool for our data. In the paper they also have 3 different conditions in 2 genotypes.
avoid running this section as it takes very long
```{r}
nclusters <- length(unique(samples.integrated.filt$labels))
genes.use <- rownames(samples.integrated.filt)
targets <- [email hidden]
group <- as.factor(samples.integrated.filt$Genotype_Treatment)
```
```{r fig.width=8, fig.height=6}
force.recalc = T
powers = c(seq(1,10,by=1), seq(12,20, by=2));
# skips recalculation of net object as it takes long...
reads <- JoinLayers(samples.integrated.filt, assay = "RNA") %>% GetAssayData(assay = "RNA", layer="data")
datExpr <- as.data.frame(reads[genes.use,])
rm("reads")
datExpr <- as.data.frame(t(datExpr))
datExpr <- datExpr[,goodGenes(datExpr)]
# expressed in more than 10% of cells
exprgenes = colSums(datExpr>0)>nrow(datExpr)*0.1
datExpr <- datExpr[,exprgenes]
#annotated only
mmacc=as.list(org.Mm.egALIAS2EG)
annotatedgenes = colnames(datExpr) %in% names(mmacc)
datExpr <- datExpr[,annotatedgenes]
powers = c(seq(1,10,by=1), seq(12,20, by=2));
enableWGCNAThreads(nThreads = 7)
if(force.recalc | (! file.exists(file=paste0(outputdir, "powertable.rds")))){
# Call the network topology analysis function for each set in turn
powerTable = list(
data = pickSoftThreshold(
datExpr,
powerVector=powers,
verbose = 100,
networkType="unsigned",
corFnc="cor"
)[[2]]
);
saveRDS(powerTable, file=paste0(outputdir, "powertable.rds"))
} else {
powerTable = readRDS(file=paste0(outputdir, "powertable.rds"))
}
colors = c("blue", "red","black")
# Will plot these columns of the returned scale free analysis tables
plotCols = c(2,5,6,7)
colNames = c("Scale Free Topology Model Fit", "Mean connectivity", "mean connectivity",
"Max connectivity");
# Get the minima and maxima of the plotted points
ylim = matrix(NA, nrow = 2, ncol = 4);
for (col in 1:length(plotCols)){
ylim[1, col] = min(ylim[1, col], powerTable$data[, plotCols[col]], na.rm = TRUE);
ylim[2, col] = max(ylim[2, col], powerTable$data[, plotCols[col]], na.rm = TRUE);
}
# Plot the quantities in the chosen columns vs. the soft thresholding power
par(mfcol = c(2,2));
par(mar = c(4.2, 4.2 , 2.2, 0.5))
cex1 = 0.7;
for (col in 1:length(plotCols)){
plot(powerTable$data[,1], -sign(powerTable$data[,3])*powerTable$data[,2],
xlab="Soft Threshold (power)",ylab=colNames[col],type="n", ylim = ylim[, col],
main = colNames[col]);
addGrid();
if (col==1){
text(powerTable$data[,1], -sign(powerTable$data[,3])*powerTable$data[,2],
labels=powers,cex=cex1,col=colors[1]);
} else
text(powerTable$data[,1], powerTable$data[,plotCols[col]],
labels=powers,cex=cex1,col=colors[1]);
if (col==1){
legend("bottomright", legend = 'Metacells', col = colors, pch = 20) ;
} else
legend("topright", legend = 'Metacells', col = colors, pch = 20) ;
}
softPower=3
nSets = 1
setLabels = 'OCS'
shortLabels = setLabels
multiExpr <- list()
multiExpr[['ODC']] <- list(data=datExpr)
checkSets(multiExpr) # check data size
# construct network
if(force.recalc | (! file.exists(file=paste0(outputdir, "net.rds")))){
net=blockwiseConsensusModules(multiExpr, blocks = NULL,
maxBlockSize = 30000, ## This should be set to a smaller size if the user has limited RAM
randomSeed = 12345,
corType = "pearson",
power = softPower,
consensusQuantile = 0.3,
networkType = "unsigned",
TOMType = "signed",
TOMDenom = "min",
scaleTOMs = TRUE, scaleQuantile = 0.8,
sampleForScaling = TRUE, sampleForScalingFactor = 1000,
useDiskCache = TRUE, chunkSize = NULL,
deepSplit = 4,
pamStage=FALSE,
detectCutHeight = 0.995, minModuleSize = 50,
mergeCutHeight = 0.2,
saveConsensusTOMs = TRUE,
consensusTOMFilePattern = paste0(outputdir,"ConsensusTOM-block.%b.rda"))
saveRDS(net, file=paste0(outputdir, "net.rds"))
} else {
net = readRDS(file=paste0(outputdir, "net.rds"))
}
consMEs = net$multiMEs;
moduleLabels = net$colors;
# Convert the numeric labels to color labels
moduleColors = as.character(moduleLabels)
consTree = net$dendrograms[[1]];
# module eigengenes
MEs=moduleEigengenes(multiExpr[[1]]$data, colors = moduleColors, nPC=1)$eigengenes
MEs=orderMEs(MEs)
meInfo<-data.frame(rownames(datExpr), MEs)
colnames(meInfo)[1]= "SampleID"
# intramodular connectivity
KMEs<-signedKME(datExpr, MEs,outputColumnName = "kME",corFnc = "bicor")
# compile into a module metadata table
geneInfo=as.data.frame(cbind(colnames(datExpr),moduleColors, KMEs))
# how many modules did we get?
nmodules <- length(unique(moduleColors))
# merged gene symbol column
colnames(geneInfo)[1]= "GeneSymbol"
colnames(geneInfo)[2]= "Initially.Assigned.Module.Color"
# save info
write.csv(geneInfo,file=paste0(outputdir,'/geneInfoSigned.csv'))
PCvalues=MEs
```
```{r}
p = WGCNA::plotDendroAndColors(consTree, moduleColors[net$goodGenes], "Module colors", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05,
main = paste0("Microglia gene dendrogram and module colors"))
ggsave(file=paste0(outputdir, "Microglia_WGCNA_Dendrocolplot_integrated.pdf"))
```
## WGCNA by Genotype_treatment
```{r, fig.width=7, fig.height=4}
plot_df <- cbind(dplyr::select(targets, c(Genotype, Treatment, labels)), PCvalues)
plot_df <- reshape2::melt(plot_df, id.vars = c('Genotype', 'Treatment', "labels"))
plot_df$Genotype <- factor(plot_df$Genotype, levels=c('WT','APPPS1'))
plot_df$Genotype_Treatment <- paste0(plot_df$Genotype, "_", plot_df$Treatment)
plot_df$Genotype_Treatment <- factor(plot_df$Genotype_Treatment, levels=c('WT_Ctrl','WT_Stroke', "APPPS1_Ctrl", "APPPS1_Stroke"))
colors <- sub('ME', '', as.character(levels(plot_df$variable)))
p <- ggplot(plot_df, aes(x=variable, y=value, fill=Treatment)) +
geom_boxplot(notch=FALSE) +
RotatedAxis() + ylab('Module Eigengene') + xlab('') +
theme(
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
)
p <- p + facet_wrap(~variable+Genotype, scales='free_x', ncol=nmodules)+theme_classic()
p
ggsave(file=paste0(outputdir, "Microgila_WGCNA_Eigenscore_byGenotype_integrated.pdf"), p)
```
## WGCNA violin by condition
```{r, fig.height=2, fig.width=8}
p2 <- ggplot(plot_df, aes(x=variable, y=value, fill=Genotype_Treatment)) +
geom_violin() +
# RotatedAxis() +
ylab('Module Eigengene') + xlab('') +
theme(
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
)
p2 <- p2 + facet_wrap(~variable, scales='free', ncol=nmodules)+theme_classic()
p2
ggsave(file=paste0(outputdir, "Microgila_WGCNA_Eigenscore_byGenotypeTreatment_integrated.pdf"), p2)
```
## WGCNA violin by cluster
```{r, fig.height=2, fig.width=10}
p2 <- ggplot(plot_df, aes(x=variable, y=value, fill=labels)) +
geom_violin() +
# RotatedAxis() +
ylab('Module Eigengene') + xlab('') +
theme(
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
)
p2 <- p2 + facet_wrap(~variable, scales='free', ncol=nmodules)+theme_classic()
p2
ggsave(file=paste0(outputdir, "Microgila_WGCNA_Eigenscore_byCluster_integrated.pdf"), p2)
```
## WGCNA Heatmap by condition
```{r, fig.width=6, fig.height=4}
reslist = list()
for(M in unique(plot_df$variable)){
Mdata = plot_df[plot_df$variable==M,]
fit= lm(value~Genotype*Treatment, data = Mdata)
resfit = summary(fit)
reslist[[M]] <- c(Genotype_Treatment = resfit$coefficients["GenotypeAPPPS1:TreatmentStroke", "Pr(>|t|)"])
}
ME_avg <- plot_df %>% dplyr::group_by(variable, Treatment, Genotype) %>% summarise(avg_Eigenmodule = mean(value),
sd_Eigenvalue = sd(value),
med_Eigenvalue=median(value))
ME_avg$pval=unlist(reslist[ME_avg$variable])
ME_avg$pval_ast=convertpvaltostars(ME_avg$pval)
ME_avg$Genotype_Treatment = paste0(ME_avg$Genotype, "_",ME_avg$Treatment)
ME_avg$Genotype_Treatment <- factor(ME_avg$Genotype_Treatment, levels=c('WT_Ctrl','WT_Stroke', "APPPS1_Ctrl", "APPPS1_Stroke"))
ME_avg$Module = paste0(ME_avg$pval_ast, ME_avg$variable)
ME_avg
openxlsx::write.xlsx(ME_avg, file=paste0(outputdir, "WGCNA_LinModStat_GTxTreatm.xlsx"),rowNames = TRUE)
```
```{r, fig.width=8, fig.height=5}
Meandata = MEs %>% group_by(samples.integrated.filt$Genotype_Treatment) %>% summarise_all(mean) %>% column_to_rownames("samples.integrated.filt$Genotype_Treatment") %>% t()
SDdata = MEs %>% group_by(samples.integrated.filt$Genotype_Treatment) %>% summarise_all(sd)%>% column_to_rownames("samples.integrated.filt$Genotype_Treatment") %>% t()
colLabels=colnames(Meandata)
rowLabels=rownames(Meandata)
rowLabels=paste0(convertpvaltostars(unlist(reslist[rowLabels])), " ",rowLabels=rownames(Meandata))
textlab = paste0(signif(Meandata,2), "\n", "(", signif(SDdata, 2), ")")
pdf(paste0(outputdir, "WGCNA_labeled_Heatmap_integrated.pdf"))
par(mar=c(5,8,5,5))
labeledHeatmap(Meandata, xLabels = colLabels,
yLabels = rowLabels,
yColorLabels = T,
colors = viridis(50),
setStdMargins = FALSE, textMatrix = textlab,
main = "Mean(SD) Eigengene-expression");
dev.off()
par(mar=c(5,8,5,5))
labeledHeatmap(Meandata, xLabels = colLabels,
yLabels = rowLabels,
yColorLabels = T,
colors = viridis(50),
setStdMargins = FALSE, textMatrix = textlab,
main = "Mean(SD) Eigengene-expression");
```
## WGCNA Heatmap by cluster
```{r, fig.width=6, fig.height=4}
reslist = list()
for(M in unique(plot_df$variable)){
Mdata = plot_df[plot_df$variable==M,]
fit= aov(value~labels, data = Mdata)
resfit = summary(fit)
reslist[[M]] <- c(Clusters = resfit[[1]]["labels", "Pr(>F)"])
}
ME_avg <- plot_df %>% dplyr::group_by(variable, labels) %>% summarise(avg_Eigenmodule = mean(value),
sd_Eigenvalue = sd(value),
med_Eigenvalue=median(value))
ME_avg$pval=unlist(reslist[ME_avg$variable])
ME_avg$pval_ast=convertpvaltostars(ME_avg$pval)
ME_avg$Module = paste0(ME_avg$pval_ast, ME_avg$variable)
ME_avg
openxlsx::write.xlsx(ME_avg, file=paste0(outputdir, "WGCNA_LinModStat_ByCluster.xlsx"),rowNames = TRUE)
```
```{r, fig.width=10, fig.height=5}
Meandata = MEs %>% group_by(samples.integrated.filt$labels) %>% summarise_all(mean) %>% column_to_rownames("samples.integrated.filt$labels") %>% t()
SDdata = MEs %>% group_by(samples.integrated.filt$labels) %>% summarise_all(sd)%>% column_to_rownames("samples.integrated.filt$labels") %>% t()
colLabels=colnames(Meandata)
rowLabels=rownames(Meandata)
rowLabels=rownames(Meandata)
rowLabels=paste0(convertpvaltostars(unlist(reslist[rowLabels])), " ",rowLabels=rownames(Meandata))
textlab = paste0(signif(Meandata,2), "\n", "(", signif(SDdata, 2), ")")
pdf(paste0(outputdir, "WGCNA_labeled_Heatmap_byCluster_integrated.pdf"))
par(mar=c(5,8,5,5))
labeledHeatmap(Meandata, xLabels = colLabels,
yLabels = rowLabels,
colorLabels = TRUE,
colors = viridis(50),
setStdMargins = FALSE, textMatrix = textlab,
main = "Mean(SD) Eigengene-expression");
dev.off()
par(mar=c(5,8,5,5))
labeledHeatmap(Meandata, xLabels = colLabels,
yLabels = rowLabels,
colorLabels = TRUE,
colors = viridis(50),
setStdMargins = FALSE, textMatrix = textlab,
main = "Mean(SD) Eigengene-expression");
```
## Heatmaps Genewise per module
```{r}
Modules = unique(moduleColors)
Modules = Modules[! Modules %in% "grey"] # drop grey one
for(m in Modules){
p = DoHeatmap(samples.integrated.filt, features = rownames(samples.integrated.filt)[moduleColors==m], slot="data",
group.by = "Genotype_Treatment")+scale_fill_viridis()
ggsave(paste0(paste0(outputdir, "WGCNAGeneWise_ExpressionModule", m, "_integrated.pdf")), p)
}
```
## Heatmpas top 50 genes per module Conditions wise
```{r, fig.width=14, fig.height=14}
genelist = split(names(net$colors), net$colors)
gene_univers = names(net$colors)
plotlist = list()
reads = rowSums(JoinLayers(samples.integrated.filt, assay = "RNA") %>% GetAssayData(assay = "RNA", layer="data") )
for (m in names(genelist)){
targets = genelist[[m]]
targets <- reads[targets] %>% sort(., decreasing = T)
targets <- names(targets)[1:50]
p<-DoHeatmap(samples.integrated.filt,features = targets,
group.by = "Genotype_Treatment", slot = "data")
ggsave(p, file=paste0(outputdir, "WGCNA_Module_",m,"_allgenes_Heatmaps_integrated.pdf"))
plotlist[[m]]<-p
}
p<- ggarrange(plotlist = plotlist, labels = names(genelist))
p
ggsave(file=paste0(outputdir, "WGCNA_Module_Heatmaps_integrated.pdf"), p)
names(plotlist)
```
```{r, fig.width=5, fig.height=5}
module = "green"
featureslist=c("Apoe", "Axl", "Cst7", "Tyrobp","Cd63", "Clec7a", "H2-D1", "Lyz2", "Npc2")
pdf(paste0(outputdir, "WGCNA_", module, "selectedgenes_heatmap_integrated.pdf"))
DoHeatmap(samples.integrated.filt,features = featureslist,
group.by = "Genotype_Treatment", slot = "data")
dev.off()
```
```{r, fig.width=5, fig.height=5}
module = "yellow"
featureslist=c("Trem2", "Csf1r", "C1qa", "C1qb", "C1qc", "Ctss",
"Ctsl", "Ctsh", "Cyba", "Lamp1")
pdf(paste0(outputdir, "WGCNA_", module, "selectedgenes_heatmap_integrated.pdf"))
DoHeatmap(samples.integrated.filt,features = featureslist,
group.by = "Genotype_Treatment", slot = "data")
dev.off()
```
```{r, fig.width=5, fig.height=5}
module = "turqouise"
featureslist=c("Olfr1307", "Olfr635", "Olfr1336", "Olfr46",
"Olfr344", "Fgd4", "Pkib", "Shisa5", "Senp7", "Snhg14")
pdf(paste0(outputdir, "WGCNA_", module, "selectedgenes_heatmap_integrated.pdf"))
DoHeatmap(samples.integrated.filt,features = featureslist,
group.by = "Genotype_Treatment", slot = "data")
dev.off()
```
```{r, fig.width=5, fig.height=5}
module = "grey"
featureslist=c("Adam10", "Abca1", "CD74","Arpc1b", "Eef2",
"Itgb2", "Pdia3", "Aldoa", "Cfh", "Ssh2")
pdf(paste0(outputdir, "WGCNA_", module, "selectedgenes_heatmap_integrated.pdf"))
DoHeatmap(samples.integrated.filt,features = featureslist,
group.by = "Genotype_Treatment", slot = "data")
dev.off()
```
## Heatmpas top 50 genes per module Cluster wise
```{r, fig.width=14, fig.height=14}
plotlist = list()
reads = rowSums(JoinLayers(samples.integrated.filt, assay = "RNA") %>% GetAssayData(assay = "RNA", layer="data") )
for (m in names(genelist)){
targets = genelist[[m]]
targets <- reads[targets] %>% sort(., decreasing = T)
targets <- names(targets)[1:50]
p<-DoHeatmap(samples.integrated.filt, features = targets,
group.by = "labels", slot = "data")
ggsave(p, file=paste0(outputdir, "WGCNA_Module_",m,"_Heatmaps_byCluster_integrated.pdf"))
plotlist[[m]]<-p
}
p<- ggarrange(plotlist = plotlist, labels = names(genelist))
p
ggsave(file=paste0(outputdir,"WGCNA_Module_Heatmaps_byCluster_integrated.pdf"), p)
names(plotlist)
```
```{r, fig.width=5, fig.height=5}
module = "green"
featureslist=c("Apoe", "Axl", "Cst7", "Tyrobp","Cd63", "Clec7a", "H2-D1", "Lyz2", "Npc2")
pdf(paste0(outputdir, "WGCNA_", module, "selectedgenes_heatmap_byCluster_integrated.pdf"))
DoHeatmap(samples.integrated.filt,features = featureslist,
group.by = "labels", slot = "data")
dev.off()
```
```{r, fig.width=5, fig.height=5}
module = "yellow"
featureslist=c("Trem2", "Csf1r", "C1qa", "C1qb", "C1qc", "Ctss",
"Ctsl", "Ctsh", "Cyba", "Lamp1")
pdf(paste0(outputdir, "WGCNA_", module, "selectedgenes_heatmap_byCluster_integrated.pdf"))
DoHeatmap(samples.integrated.filt,features = featureslist,
group.by = "labels", slot = "data")
dev.off()
```
```{r, fig.width=5, fig.height=5}
module = "turqouise"
featureslist=c("Olfr1307", "Olfr635", "Olfr1336", "Olfr46",
"Olfr344", "Fgd4", "Pkib", "Shisa5", "Senp7", "Snhg14")
pdf(paste0(outputdir, "WGCNA_", module, "selectedgenes_heatmap_byCluster_integrated.pdf"))
DoHeatmap(samples.integrated.filt,features = featureslist,
group.by = "labels", slot = "data")
dev.off()
```
```{r, fig.width=5, fig.height=5}
module = "grey"
featureslist=c("Adam10", "Abca1", "CD74","Arpc1b", "Eef2",
"Itgb2", "Pdia3", "Aldoa", "Cfh", "Ssh2")
pdf(paste0(outputdir, "WGCNA_", module, "selectedgenes_heatmap_byCluster_integrated.pdf"))
DoHeatmap(samples.integrated.filt,features = featureslist,
group.by = "labels", slot = "data")
dev.off()
```
## WGCNA GO terms
```{r, out.width='130%',out.height='100%'}
pvallist = unlist(reslist)
MOI = gsub("ME", "", names(pvallist)) %>% gsub(".Clusters", "", .)
# MOI = MOI[pvallist<0.05]
MOI <- MOI[MOI != "grey"]
# if you want to have the affected genes set return_genelist = T; this takes roughly 4-8 hours
Gores = getGOresults(geneset = genelist[MOI], domain_scope = "known",
genereference = gene_univers,
return_genelist = T,
organism = "mmusculus")
display_tab(Gores$result)
openxlsx::write.xlsx(Gores$result%>% dplyr::select(-parents), file=paste0(outputdir, "WGCNA_GO_terms.xlsx"),rowNames = TRUE)
plotlist <- list()
for(m in unique(Gores$result$query)){
p <- GOplot(Gores$result[Gores$result$query ==m,], N = 10, Title = m)
plotlist[[m]] <- p
}
```
```{r, fig.width=15, fig.height=6}
p<- ggarrange(plotlist = plotlist)
p
ggsave(file=paste0(outputdir,"WGCNA_GO_terms.pdf"), p)
```
# Exdendet Plots and DAta inspection
## Dot plot showing top markers for condition clusters, same paper Fig 2f
### MG related genes
```{r, fig.width=18, fig.height=5}
reads = JoinLayers(samples.integrated.filt, assay = "RNA") %>%
GetAssayData(assay = "RNA", layer="data")
reads.filt= reads[cluster.genes,]
datatoplot_mean = aggregate(t(reads.filt),
list(paste0(samples.integrated.filt$labels, "/",samples.integrated.filt$Genotype_Treatment)),
mean, na.rm=T) %>%
as.data.frame() %>%
column_to_rownames("Group.1") %>% scale()
datatoplot_mean <- datatoplot_mean %>% as.data.frame() %>%
rownames_to_column("Var1") %>%
pivot_longer(cols = -Var1,
names_to = "Genes",
values_to = "value")
datatoplot_mean$Microglia = factor(gsub("/.*", "", datatoplot_mean$Var1), levels=rev(levels(samples.integrated.filt$labels)))
datatoplot_mean$Genotype_Treatment = factor(gsub(".*/", "", datatoplot_mean$Var1), levels=c("WT_Ctrl", "WT_Stroke", "APPPS1_Ctrl","APPPS1_Stroke") )
colnames(datatoplot_mean) = c("Cluster_Condition", "Gene", "Scaled_mean", "Microglia", "Genotype_Treatment")
datatoplot_perc = aggregate(t(reads.filt), list(paste0(samples.integrated.filt$labels, "/",samples.integrated.filt$Genotype_Treatment)),
function(x){length(which(x>0))/length(x)})
datatoplot_perc <- datatoplot_perc %>% as.data.frame() %>%
pivot_longer(cols = -Group.1,
names_to = "Genes",
values_to = "perc")
datatoplot = cbind(datatoplot_mean, perc=datatoplot_perc$perc)
p <- ggplot(datatoplot, aes(x=Gene, y=Microglia, col=Scaled_mean, size=perc))+geom_point()+scale_size_continuous(range = c(0,3))+
theme_classic()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+facet_grid(~Genotype_Treatment)
p
ggsave(filename = paste0(outputdir, "Microglia_Dotplot_ClsuterCondition_Genelists_integrated.pdf"), p)
```
### Spatial_mutual genes
```{r}
Spatial_mutual <- c("S100a6", "Apoc4", "C4b", "Atp1a2", "Apoc1", "Fn1", "Ntrk2") # genes up in both spatial and APPPS1_stroke vs APPPS1
### Dot plot showing top markers for condition clusters, same paper Fig 2f
reads = JoinLayers(samples.integrated.filt, assay = "RNA") %>%
GetAssayData(assay = "RNA", layer="data")
reads.filt= reads[Spatial_mutual,]
datatoplot_mean = aggregate(t(reads.filt),
list(paste0(samples.integrated.filt$labels, "/",samples.integrated.filt$Genotype_Treatment)),
mean, na.rm=T) %>%
as.data.frame() %>%
column_to_rownames("Group.1") %>% scale() %>% as.data.frame() %>%
rownames_to_column("Var1")
datatoplot_mean=datatoplot_mean %>% pivot_longer(cols = -Var1, names_to = "Gene", values_to = "value")
datatoplot_mean$Microglia = factor(gsub("/.*", "", datatoplot_mean$Var1), levels=rev(levels(samples.integrated.filt$labels)))
datatoplot_mean$Genotype_Treatment = factor(gsub(".*/", "", datatoplot_mean$Var1), levels=c("WT_Ctrl", "WT_Stroke", "APPPS1_Ctrl","APPPS1_Stroke") )
colnames(datatoplot_mean) = c("Cluster_Condition", "Gene", "Scaled_mean", "Microglia", "Genotype_Treatment")
datatoplot_perc = aggregate(t(reads.filt), list(paste0(samples.integrated.filt$labels, "/",samples.integrated.filt$Genotype_Treatment)),
function(x){length(which(x>0))/length(x)})
datatoplot_perc=datatoplot_perc %>% pivot_longer(cols=-Group.1,
values_to = "perc",
names_to = "Gene")
datatoplot = cbind(datatoplot_mean, perc=datatoplot_perc$perc)
p <- ggplot(datatoplot, aes(x=Gene, y=Microglia, col=Scaled_mean, size=perc))+geom_point()+scale_size_continuous(range = c(0,3))+
theme_classic()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+facet_grid(~Genotype_Treatment)
p
ggsave(filename = paste0(outputdir, "Microglia_Dotplot_Up_spatial_and_APPPS1vsAPPPS1stroke_integrated.pdf"), p)
```
### Custom Genes
```{r}
Test_genes <- c("S100a6", "S100a9", "S100a8", "S100a10", "S100a11", "S100b", "S100a7a") # genes up in both spatial and APPPS1_stroke vs APPPS1
### Dot plot showing top markers for condition clusters, same paper Fig 2f
reads = JoinLayers(samples.integrated.filt, assay = "RNA") %>%
GetAssayData(assay = "RNA", layer="data")
reads.filt= reads[Test_genes,]
datatoplot_mean = aggregate(t(reads.filt),
list(paste0(samples.integrated.filt$labels, "/",samples.integrated.filt$Genotype_Treatment)),
mean, na.rm=T) %>%
as.data.frame() %>%
column_to_rownames("Group.1") %>% scale() %>% as.data.frame() %>%
rownames_to_column("Var1")
datatoplot_mean=datatoplot_mean %>% pivot_longer(cols=-Var1,
names_to = "Gene",
values_to = "value")
datatoplot_mean$Microglia = factor(gsub("/.*", "", datatoplot_mean$Var1), levels=rev(levels(samples.integrated.filt$labels)))
datatoplot_mean$Genotype_Treatment = factor(gsub(".*/", "", datatoplot_mean$Var1), levels=c("WT_Ctrl", "WT_Stroke", "APPPS1_Ctrl","APPPS1_Stroke") )
colnames(datatoplot_mean) = c("Cluster_Condition", "Gene", "Scaled_mean", "Microglia", "Genotype_Treatment")
datatoplot_perc = aggregate(t(reads.filt), list(paste0(samples.integrated.filt$labels, "/",samples.integrated.filt$Genotype_Treatment)),
function(x){length(which(x>0))/length(x)})
datatoplot_perc=datatoplot_perc %>% pivot_longer(
cols=-Group.1, values_to = "perc", names_to = "Gene")
datatoplot = cbind(datatoplot_mean, perc=datatoplot_perc$perc)
p <- ggplot(datatoplot, aes(x=Gene, y=Microglia, col=Scaled_mean, size=perc))+geom_point()+scale_size_continuous(range = c(0,3))+
theme_classic()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+facet_grid(~Genotype_Treatment)
p
ggsave(filename = paste0(outputdir, "Microglia_Dotplot_test_integrated.pdf"), p)
```
## cell cylce scoring plots
```{r}
cell_cycle_markers <- dplyr::left_join(cell_cycle_genes, annotations, by = c("geneID" = "gene_id"))
# Acquire the S phase genes
s_genes <- cell_cycle_markers %>%
dplyr::filter(phase == "S") %>%
pull("gene_name")
# Acquire the G2M phase genes
g2m_genes <- cell_cycle_markers %>%
dplyr::filter(phase == "G2/M") %>%
pull("gene_name")
samples.integrated.filt<-
CellCycleScoring(samples.integrated.filt, layer="scaled.data",
s.features = s_genes,
g2m.features = g2m_genes)
a<-DimPlot(samples.integrated.filt, reduction = "umap", group.by = "Phase", pt.size =1)
b<-DimPlot(samples.integrated.filt, reduction = "umap", group.by = "Phase", split.by = "Plate", pt.size =1)
c<-DimPlot(samples.integrated.filt, reduction = "umap", group.by = "Phase", split.by = "Genotype_Treatment", pt.size =1)
a
ggsave(filename = paste0(outputdir, "G2M_Opt_All_genes_all_cells_Clusters_integrated.pdf"), a)
```
```{r, fig.width=18, fig.height=4}
b
ggsave(filename = paste0(outputdir, "G2M_Opt_All_genes_all_cells_Clusters_byPlate_integrated.pdf"), b)
```
```{r fig.width=10, fig.height=4}
c
ggsave(filename = paste0(outputdir, "G2M_Opt_All_genes_all_cells_Clusters_byCondition_integrated.pdf"), c)
```
```{r, fig.width=7, fig.height=7}
p <- DimPlot(samples.integrated.filt, reduction = "umap", group.by = "Phase")
p
ggsave(file=paste0(outputdir,"CClustering_CellPhase.pdf_integrated.pdf"), p)
p <- FeaturePlot_scCustom(samples.integrated.filt,
reduction = "umap",
features = "S.Score",
colors_use = viridis_dark_high)
p
ggsave(file=paste0(outputdir,"CClustering_CellPhaseScore.pdf"), p)
p <- FeaturePlot_scCustom(samples.integrated.filt, reduction = "umap",features = "G2M.Score", colors_use = viridis_dark_high)
p
ggsave(file=paste0(outputdir,"CClustering_CellPhaseG2MScore.pdf"), p)
metadata <- [email hidden]
p<- ggplot(metadata, aes(x=labels, y= S.Score))+geom_violin(aes(fill=labels))+geom_boxplot(aes(alpha=1))+theme_classic()
p
ggsave(file=paste0(outputdir,"CClustering_CellPhaseScoreViolin.pdf"), p)
p<- ggplot(metadata, aes(x=labels, y= G2M.Score))+geom_violin(aes(fill=labels))+geom_boxplot(aes(alpha=1))+theme_classic()
p
ggsave(file=paste0(outputdir,"CClustering_CellPhaseG2MScoreViolin.pdf"), p)
pairwise.t.test(metadata$S.Score, metadata$labels)
pairwise.t.test(metadata$G2M.Score, metadata$labels)
```
## Additonal Microglia subtype scores
```{r, fig.width=6, fig.height=6}
DefaultAssay(samples.integrated.filt) <- "RNA"
BAM = c("Cd36", "Cd38","Lyve1","Cd206","Cd163", "Cd169")
CAM = c("Emilin2", "Pf4", "Msa4a", "Hp","F5","Mki67")
Lipid = c("Lipe", "Acnat1", "Lratd1", "Or8b44", "Vmn1r13")
samples.integrated.filt<- AddModuleScore(object = samples.integrated.filt,
assay="RNA", slot="data",
features = list(INF_resp_microglia),
name = "MG_score_INF_resp_microglia")
samples.integrated.filt<- AddModuleScore(object = samples.integrated.filt,
assay="RNA", slot="data",
features = list(Cycling_microglia),
name = "MG_score_Cycling_microglia")
samples.integrated.filt<- AddModuleScore(object = samples.integrated.filt,
assay="RNA", slot="data",
features = list(Act_resp_microglia),
name = "MG_score_Act_resp_microglia")
samples.integrated.filt<- AddModuleScore(object = samples.integrated.filt,
assay="RNA", slot="data",
features = list(Axon_tract_microglia),
name = "MG_score_Axon_tract_microglia")
samples.integrated.filt <- AddModuleScore(object = samples.integrated.filt,
assay = "RNA", layer="data",
features = list(BAM), name = "BAM_score")
samples.integrated.filt <- AddModuleScore(object = samples.integrated.filt,
assay = "RNA", layer="data",
features = list(CAM), name = "CAM_score")
samples.integrated.filt <- AddModuleScore(object = samples.integrated.filt,
assay = "RNA", layer="data",
features = list(Lipid), name = "Lipid_score")
featlist=c("MG_score_Axon_tract_microglia1",
"MG_score_Act_resp_microglia1",
"MG_score_Cycling_microglia1",
"MG_score_INF_resp_microglia1",
"BAM_score1",
"CAM_score1",
"Lipid_score1"
)
p <- FeaturePlot_scCustom(samples.integrated.filt,
na_cutoff = 0, reduction="umap",
features = featlist,
colors_use = viridis_dark_high)
```
```{r, fig.width=18, fig.height=14}
p
ggsave(paste0(outputdir, "DimPlots_MGSubtypes_Genes_integrated.pdf"), p)
```
```{r, fig.width=12, fig.height=16}
p <- RidgePlot(samples.integrated.filt,
features = c("MG_score_Axon_tract_microglia1",
"MG_score_Act_resp_microglia1",
"MG_score_Cycling_microglia1",
"MG_score_INF_resp_microglia1",
"BAM_score1",
"CAM_score1",
"Lipid_score1") ,
group.by = "labels", ncol = 2)
p
ggsave(paste0(outputdir, "Ridge_MGSubtypes_Genes_integrated.pdf"), p)
```
## Additional QC plots final dataset
```{r fig.width=12, fig.height=4}
DefaultAssay(samples.integrated.filt) <- "RNA"
metadata <- [email hidden]
metadata$nCount_RNA_finset <- colSums(samples.integrated.filt)
metadata$nFeature_RNA_finset <- colSums((JoinLayers(samples.integrated.filt, assay = "RNA") %>% GetAssayData(assay = "RNA", layer="data")) >0)
metadata -> [email hidden]
p <- VlnPlot(samples.integrated.filt, features = c("nCount_RNA_finset") , group.by = "labels") +geom_hline(yintercept = 200)
p
p2 <- VlnPlot(samples.integrated.filt, features = c("nFeature_RNA_finset") , group.by = "labels") +geom_hline(yintercept = 150)
p2
p3 <- VlnPlot(samples.integrated.filt, features = "percent.mt" , group.by = "labels") + ylim(c(0,10)) +geom_hline(yintercept = 5)
p3
pcomb <- ggarrange(p,p2,p3, common.legend = T, ncol = 3, legend = "right")
pcomb
ggsave(paste0(outputdir, "Qualityplots_integrated.pdf"), pcomb)
```
## Additional Plots DEX
```{r, fig.width=4, fig.height=6}
genes_to_label_dwn <- c("Siglech", "Ccl6" , "Cst7", "Tyrobp")
genes_to_label_up <- c("H2-D1", "B2m" , "C4b", "Lipe", "Apoe")
all_genes <- c(genes_to_label_up, genes_to_label_dwn)
all_genes <- c(genes_to_label_up, genes_to_label_dwn)
# Get remaining genes not in all_genes
rest_genes <- rownames(DEG_APPPS1_CtrlvsStroke)[!rownames(DEG_APPPS1_CtrlvsStroke) %in% all_genes]
# Combine: labeled genes first, then the rest
df_res <- DEG_APPPS1_CtrlvsStroke[rev(c(all_genes, rest_genes)), ]
df_res <- df_res %>%
mutate(sign = ifelse(df_res$avg_log2FC < 0, "down",
ifelse(df_res$avg_log2FC > 0, "up", "no.reg")),
sign = ifelse(df_res$p_val_adj > 0.05, "no.reg", sign),
sign = ifelse(rownames(df_res) %in% genes_to_label_up, "up", sign),
sign = ifelse(rownames(df_res) %in% genes_to_label_dwn, "down", sign),
goi = ifelse(rownames(df_res) %in% all_genes, 1, 0))
library(ggrepel)
panelA <- ggplot(df_res, aes(x=avg_log2FC, y=-log10(p_val_adj),
col=sign, alpha=goi, size=as.factor(goi))) +
geom_point() +
theme_minimal() +
geom_hline(yintercept = -log10(0.05), lwd=0.3, linetype="dashed") +
geom_vline(xintercept = 0, lwd=0.3, linetype="dashed") +
scale_color_manual(name = "",
values = c("up" = "red3", "no.reg" = "gray60", "down" = "royalblue")) +
guides(alpha = "none", size = "none") +
geom_label_repel(data = df_res[all_genes, ],
aes(label = rownames(df_res[all_genes, ])),
size = 3,
max.overlaps = Inf,
box.padding = 0.5,
point.padding = 0.3,
segment.color = "gray40",
segment.size = 0.3,
min.segment.length = 0,
force = 2,
show.legend = FALSE) +
theme(legend.position = "top",
legend.justification = "left") +
ggtitle("Bulk differential expression \nall micorglia integrated")
panelA
```
```{r}
reads = JoinLayers(samples.integrated.filt, assay = "RNA") %>%
GetAssayData(assay = "RNA", layer="data")
subset = samples.integrated.filt$Genotype =="APPPS1"
reads.filt= reads[all_genes, subset]
datatoplot_mean = aggregate(t(reads.filt),
list(paste0(samples.integrated.filt$labels[subset], "/",samples.integrated.filt$Genotype_Treatment[subset])),
mean, na.rm=T) %>%
as.data.frame() %>%
column_to_rownames("Group.1")
datatoplot_mean <- datatoplot_mean %>% as.data.frame() %>%
rownames_to_column("Var1") %>%
pivot_longer(cols = -Var1,
names_to = "Genes",
values_to = "value")
datatoplot_mean$Microglia = factor(gsub("/.*", "", datatoplot_mean$Var1), levels=rev(levels(samples.integrated.filt$labels)))
datatoplot_mean$Genotype_Treatment = factor(gsub(".*/", "", datatoplot_mean$Var1), levels=c("WT_Ctrl", "WT_Stroke", "APPPS1_Ctrl","APPPS1_Stroke") )
colnames(datatoplot_mean) = c("Cluster_Condition", "Gene", "Scaled_mean", "Microglia", "Genotype_Treatment")
datatoplot_perc = aggregate(t(reads.filt), list(paste0(samples.integrated.filt$labels[subset], "/",samples.integrated.filt$Genotype_Treatment[subset])),
function(x){length(which(x>0))/length(x)})
datatoplot_perc <- datatoplot_perc %>% as.data.frame() %>%
pivot_longer(cols = -Group.1,
names_to = "Genes",
values_to = "perc")
datatoplot = cbind(datatoplot_mean, perc=datatoplot_perc$perc)
datatoplot$Microglia <- as.character(datatoplot$Microglia)
datatoplot$Gene <- factor(datatoplot$Gene, levels=all_genes, labels=all_genes)
datatoplot <- datatoplot %>%
group_by(Microglia, Gene) %>%
mutate(Scaled_Diff = Scaled_mean - mean(Scaled_mean)) %>%
ungroup()
library(ggh4x)
Microglia_color_palette <- scales::hue_pal()(nlevels(factor(datatoplot$Microglia)))
# Define colors for each Microglia level
strip_colors <- setNames(Microglia_color_palette, levels(factor(datatoplot$Microglia)))
strip <- strip_themed(background_x = elem_list_rect(fill = strip_colors))
panelB <- ggplot(datatoplot, aes(x=Genotype_Treatment, y=Gene, col=Scaled_Diff, size=perc)) +
geom_point() +
scale_size_continuous(range = c(1,7)) +
theme_classic() +
geom_hline(yintercept = 5.5, lwd=0.5, linetype="dashed", color="gray40")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
strip.text = element_text(face = "bold", color = "white")) +
facet_grid2(~Microglia, strip = strip) +
scale_color_continuous(palette = colorRampPalette(rev(RColorBrewer::brewer.pal(5,"RdBu")))(100))+
ggtitle("Per cluster difference in expression \n normlaized to within cluster mean")
```
```{r, fig.width=15, fig.height=6}
p<- ggarrange(panelA,panelB, widths = c(1.2,4), labels = c("A", "B"))
ggsave(filename = paste0(outputdir, "APPPS1_Bulk_vs_Microglia_Dotplot_ClusterCondition_Genelists_integrated.pdf"), p, width = 15, height = 6)
```
```{r}
rm(list=c("reads", "counts", "counts_run2", "counts_run1",
"counts_sel", "results_list", "sce.microglia"))
gc()
```
```{r}
saveRDS(samples.integrated.filt, paste0(outputdir, "SeuratObjectafterProcessing.rds"))
# samples.integrated.filt <- readRDS( paste0(outputdir, "SeuratObjectafterProcessing.rds"))
```

scRNA_Analyses_Candlishetal.Rmd at commit e83d217, no license · at the source

Overview

Authors: Michael Candlish¹, Jan Hofmann¹, Desirée Brösamle^2,3,4,5, Annika Haessler⁶, Murphy DeMeglio¹, Angelos Skodras², Georgi Tushev⁷, Eloah S De Biasi¹, Stefan Günther⁸, René Wiegandt⁸, Heidi Theis⁹, Elena De Domenico⁹, Nina Sofia Hermann^2,3,4,5, Peter Breunig¹, Christina Sauerland¹, K Peter R Nilsson¹⁰, Marc D Beyer^9,11, Mario Looso⁸, Maike Windbergs⁶, Sigrun Roeber¹², Jochen Herms¹³, Jonas J Neher^2,3,4,5, Andreas G Chiocchetti¹⁴, Jasmin K Hefendehl¹

ORCID iDs: Andreas G Chiocchetti, Jasmin K Hefendehl

14 affiliations

Institutions: Goethe University Frankfurt (Germany); German Center for Neurodegenerative Diseases (Germany); Munich Cluster for Systems Neurology (Germany); Hertie Institute for Clinical Brain Research (Germany); University of Tübingen (Germany); Ludwig-Maximilians-Universität München (Germany); Max Planck Institute for Brain Research (Germany); Max Planck Institute for Heart and Lung Research (Germany); German Center for Lung Research (Germany); Cardio-Pulmonary Institute (Germany); Linköping University (Sweden); Institut für Neuropathologie (Germany); University Hospital Frankfurt (Germany)

Journal: Journal of neuroinflammation, volume 23, issue 1, article 213

Dates: received 29 August 2025; accepted 29 May 2026; published online 9 June 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1186/s12974-026-03897-x · PMID 42265753 · PMCID PMC13292430 · OpenAlex W4417255812

Open access: gold, a free copy (OpenAlex)

Preprint: doi:10.64898/2025.12.09.692939

Status: code verified

Categories: human (organism), mouse (organism), Alzheimer's / dementia (population), stroke (population), cellular / molecular (subfield)

Methods: Statistics, Smoothing, state filtering, decompositions

Keywords: Alzheimer’s disease, Stroke, Microglia, Co-morbidity

MeSH: Amyloid beta-Peptides*, Brain Ischemia*, Microglia*, Alzheimer Disease, Amyloid beta-Protein Precursor, Animals, Apolipoproteins E, Disease Models, Animal, Humans, Mice, Mice, Inbred C57BL, Mice, Transgenic, Phenotype (* major topic)

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

Funding: Johann Wolfgang Goethe-Universität, Frankfurt am Main

Citations: cited by 1 paper (Europe PMC); 39 references in the paper

Abstract

Microglia are highly plastic cells that are capable of integrating subsequent insults. As the majority of Alzheimer’s Disease (AD) patients also show cerebrovascular pathology, we here aimed to dissect the interactions between AD and ischemic brain injury on the microglial response to amyloid beta (Aβ) pathology. Unexpectedly, ischemic stroke in the context of cerebral β-amyloidosis drives the emergence of a neuroprotective microglial phenotype characterized by an ApoE-enriched transcriptional state and enhanced lipid handling. These microglia promote the rapid formation of highly compact Aβ plaques that are relatively inert and strikingly reminiscent of those observed in cognitively resilient AD patients. Our findings thus reveal that the microglial response to Aβ pathology is not a fixed trajectory toward dysfunction, but retains a capacity for beneficial reprogramming when engaged by the appropriate stimulus. Beyond characterizing this comorbid state, our data identify specific molecular pathways, centered on ApoE, complement activation, and lysosomal processing, that may be amenable to therapeutic targeting to promote protective microglial function in AD.

Supplementary Information: The online version contains supplementary material available at 10.1186/s12974-026-03897-x.

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

Repositories

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

KJPMolgenLab/scSeq_Hefendehl

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: e83d21755cf398796bb6dc512bcd7bc3548688da, 28 May 2026

Languages: JavaScript (27), R (7), Shell (1)

Size: 345 files, 35 scripts

Software Heritage: not archived

Found in: “Data availability”

Holds: README, documentation, 4 notebooks

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

Tools: tidyverse (2 files), data.table (1 file), DESeq2 (1 file), ggpubr (1 file), Harmony (1 file), limma (1 file), pheatmap (1 file), Plotly (1 file), reshape2 (1 file), Seurat (1 file), SingleCellExperiment (1 file), WGCNA (1 file)

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

27 September 2026: the link answers

36 files

analysis/about.Rmd, R, 20 lines
analysis/index.Rmd, R, 274 lines
analysis/license.Rmd, R, 30 lines
analysis/scRNA_Analyses_Candlishetal.Rmd, R, 3,920 lines, 3 matches
code/Installer_Pak.R, R, 71 lines
code/custom_functions.R, R, 325 lines
code/installer.sh, Shell, 72 lines
docs/site_libs/bootstrap-3.3.5/js/bootstrap.js, JavaScript, 2,363 lines
docs/site_libs/bootstrap-3.3.5/js/bootstrap.min.js, JavaScript, 7 lines
docs/site_libs/bootstrap-3.3.5/js/npm.js, JavaScript, 13 lines
docs/site_libs/bootstrap-3.3.5/shim/html5shiv.min.js, JavaScript, 7 lines
docs/site_libs/bootstrap-3.3.5/shim/respond.min.js, JavaScript, 8 lines
docs/site_libs/crosstalk-1.2.2/js/crosstalk.js, JavaScript, 1,474 lines
docs/site_libs/crosstalk-1.2.2/js/crosstalk.min.js, JavaScript, 2 lines
docs/site_libs/datatables-binding-0.34.0/datatables.js, JavaScript, 1,539 lines
docs/site_libs/dt-core-1.13.6/js/jquery.dataTables.min.js, JavaScript, 4 lines
docs/site_libs/dt-ext-buttons-1.13.6/js/buttons.colVis.min.js, JavaScript, 5 lines
docs/site_libs/dt-ext-buttons-1.13.6/js/buttons.html5.min.js, JavaScript, 8 lines
docs/site_libs/dt-ext-buttons-1.13.6/js/buttons.print.min.js, JavaScript, 5 lines
docs/site_libs/dt-ext-buttons-1.13.6/js/dataTables.buttons.min.js, JavaScript, 4 lines
docs/site_libs/header-attrs-2.30/header-attrs.js, JavaScript, 12 lines
docs/site_libs/highlightjs-9.12.0/highlight.js, JavaScript, 2 lines
docs/site_libs/htmlwidgets-1.6.4/htmlwidgets.js, JavaScript, 901 lines
docs/site_libs/jquery-3.6.0/jquery-3.6.0.js, JavaScript, 7,407 lines
docs/site_libs/jquery-3.6.0/jquery-3.6.0.min.js, JavaScript, 2 lines
docs/site_libs/jqueryui-1.13.2/jquery-ui.js, JavaScript, 7,423 lines
docs/site_libs/jqueryui-1.13.2/jquery-ui.min.js, JavaScript, 6 lines
docs/site_libs/jszip-1.13.6/jszip.min.js, JavaScript, 13 lines
docs/site_libs/navigation-1.1/codefolding.js, JavaScript, 76 lines
docs/site_libs/navigation-1.1/sourceembed.js, JavaScript, 12 lines
docs/site_libs/navigation-1.1/tabsets.js, JavaScript, 141 lines
docs/site_libs/nouislider-7.0.10/jquery.nouislider.min.js, JavaScript, 3 lines
docs/site_libs/selectize-0.12.0/selectize.min.js, JavaScript, 3 lines
docs/site_libs/tocify-1.9.1/jquery.tocify.js, JavaScript, 1,002 lines
wflowhelper.R, R, 12 lines
README.md, Text, 257 lines

At the source: github.com/KJPMolgenLab/scSeq_Hefendehl/

gitlab.mpcdf.mpg.de/mpibr/scic/roioverstack

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: c86198d87dcd2f6dea1896e12a607be5cf2be609, 6 May 2024

Languages: MATLAB (1)

Size: 4 files, 1 script

Software Heritage: not archived

Found in: “Code availability”

Holds: README, tests

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

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

27 September 2026: the link answers

2 files

test/findClosestPathToEdge.m, MATLAB, 119 lines
README.md, Text, 22 lines

At the source: gitlab.mpcdf.mpg.de/mpibr/scic/roioverstack

Code availability

The code for segmentation of qFTAA- and hFTAA- labelled Aβ plaques can be found here: https://gitlab.mpcdf.mpg.de/mpibr/scic/roioverstack. Any additional requests for code should be addressed to Prof. Jasmin Hefendehl.

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:

2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
36 scripts, each with its path and the digest of its content;
3 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

figshare:32787540, at figshare; found in DataCite
figshare:32787543, at figshare; found in DataCite
figshare:32787546, at figshare; found in DataCite
figshare:32787549, at figshare; found in DataCite
kjpmolgenlab.github.io/scseq_hefendehl, at kjpmolgenlab.github.io; found in “Data availability”

Data Availability Statement

Any additional requests for data should be addressed to Prof. Jasmin Hefendehl. Details regarding scRNAseq analysis and data can be found at the following websites respectively. https://github.com/KJPMolgenLab/scSeq_Hefendehl/ and https://kjpmolgenlab.github.io/scSeq_Hefendehl/.

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, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 24 authors, 4 keywords, 13 MeSH terms, 1 funder, 39 references.

Cite

This paper

Candlish, M., Hofmann, J., Brösamle, D., Haessler, A., DeMeglio, M., Skodras, A., Tushev, G., De Biasi, E. S., Günther, S., Wiegandt, R., Theis, H., De Domenico, E., Hermann, N. S., Breunig, P., Sauerland, C., Nilsson, K. P. R., Beyer, M. D., Looso, M., Windbergs, M., . . . Hefendehl, J. K. (2026). Ischemic injury triggers a protective microglial phenotype in models of Aβ pathology. Journal of neuroinflammation, 23(1), 213. https://doi.org/10.1186/s12974-026-03897-x

BibTeX

@article{candlish2026ischemic, author = {Candlish, Michael and Hofmann, Jan and Brösamle, Desirée and Haessler, Annika and DeMeglio, Murphy and Skodras, Angelos and Tushev, Georgi and De Biasi, Eloah S and Günther, Stefan and Wiegandt, René and Theis, Heidi and De Domenico, Elena and Hermann, Nina Sofia and Breunig, Peter and Sauerland, Christina and Nilsson, K Peter R and Beyer, Marc D and Looso, Mario and Windbergs, Maike and Roeber, Sigrun and Herms, Jochen and Neher, Jonas J and Chiocchetti, Andreas G and Hefendehl, Jasmin K}, title = {{Ischemic injury triggers a protective microglial phenotype in models of Aβ pathology}}, journal = {Journal of neuroinflammation}, year = {2026}, month = jun, volume = {23}, number = {1}, pages = {213}, publisher = {BMC}, issn = {1742-2094}, doi = {10.1186/s12974-026-03897-x}, url = {https://doi.org/10.1186/s12974-026-03897-x}, pmid = {42265753}, pmcid = {PMC13292430} }

RIS

TY - JOUR AU - Candlish, Michael AU - Hofmann, Jan AU - Brösamle, Desirée AU - Haessler, Annika AU - DeMeglio, Murphy AU - Skodras, Angelos AU - Tushev, Georgi AU - De Biasi, Eloah S AU - Günther, Stefan AU - Wiegandt, René AU - Theis, Heidi AU - De Domenico, Elena AU - Hermann, Nina Sofia AU - Breunig, Peter AU - Sauerland, Christina AU - Nilsson, K Peter R AU - Beyer, Marc D AU - Looso, Mario AU - Windbergs, Maike AU - Roeber, Sigrun AU - Herms, Jochen AU - Neher, Jonas J AU - Chiocchetti, Andreas G AU - Hefendehl, Jasmin K TI - Ischemic injury triggers a protective microglial phenotype in models of Aβ pathology T2 - Journal of neuroinflammation J2 - J Neuroinflammation PY - 2026 DA - 2026/06/09 VL - 23 IS - 1 SP - 213 SN - 1742-2094 PB - BMC DO - 10.1186/s12974-026-03897-x UR - https://doi.org/10.1186/s12974-026-03897-x LA - en ER -

CSL-JSON

{ "id": "10.1186/s12974-026-03897-x", "type": "article-journal", "title": "Ischemic injury triggers a protective microglial phenotype in models of Aβ pathology", "container-title": "Journal of neuroinflammation", "author": [ { "family": "Candlish", "given": "Michael" }, { "family": "Hofmann", "given": "Jan" }, { "family": "Brösamle", "given": "Desirée" }, { "family": "Haessler", "given": "Annika" }, { "family": "DeMeglio", "given": "Murphy" }, { "family": "Skodras", "given": "Angelos" }, { "family": "Tushev", "given": "Georgi" }, { "family": "De Biasi", "given": "Eloah S" }, { "family": "Günther", "given": "Stefan" }, { "family": "Wiegandt", "given": "René" }, { "family": "Theis", "given": "Heidi" }, { "family": "De Domenico", "given": "Elena" }, { "family": "Hermann", "given": "Nina Sofia" }, { "family": "Breunig", "given": "Peter" }, { "family": "Sauerland", "given": "Christina" }, { "family": "Nilsson", "given": "K Peter R" }, { "family": "Beyer", "given": "Marc D" }, { "family": "Looso", "given": "Mario" }, { "family": "Windbergs", "given": "Maike" }, { "family": "Roeber", "given": "Sigrun" }, { "family": "Herms", "given": "Jochen" }, { "family": "Neher", "given": "Jonas J" }, { "family": "Chiocchetti", "given": "Andreas G" }, { "family": "Hefendehl", "given": "Jasmin K" } ], "container-title-short": "J Neuroinflammation", "volume": "23", "issue": "1", "page": "213", "DOI": "10.1186/s12974-026-03897-x", "PMID": "42265753", "PMCID": "PMC13292430", "ISSN": "1742-2094", "publisher": "BMC", "URL": "https://doi.org/10.1186/s12974-026-03897-x", "language": "en", "issued": { "date-parts": [ [ 2026, 6, 9 ] ] } }

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.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: WGCNA, Harmony, SingleCellExperiment, 9 other tools, cellular / molecular
[2] doi:10.1002/imt2.70163 [code]: Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: Harmony, SingleCellExperiment, limma, 8 other tools, mouse, cellular / molecular, 1 reference
[3] 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: WGCNA, Harmony, SingleCellExperiment, 8 other tools, Alzheimer's / dementia, cellular / molecular
[4] doi:10.1038/s41467-026-73305-8 [code]: Comparative analysis of the cellular landscape in mammalian striatum.
Journal: Nature communications
In common: WGCNA, Harmony, SingleCellExperiment, 7 other tools, mouse, cellular / molecular
[5] doi:10.1038/s41467-026-70232-6 [code]: Gene expression dynamics of human and mouse craniofacial development at the single-cell level.
Journal: Nature communications
In common: Harmony, SingleCellExperiment, limma, 7 other tools, mouse, cellular / molecular
[6] doi:10.1038/s41586-026-10414-w [code]: Focal white matter lesions drive grey matter inflammation and synapse loss.
Journal: Nature
In common: Harmony, SingleCellExperiment, DESeq2, 4 other tools, Alzheimer's / dementia, mouse, cellular / molecular, 3 references
[7] doi:10.1038/s42003-026-10957-8 [code]: Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: WGCNA, limma, DESeq2, 7 other tools, mouse, cellular / molecular
[8] 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, limma, DESeq2, 7 other tools, Alzheimer's / dementia
[9] doi:10.1016/j.isci.2026.115657 [code]: Integration of machine learning to develop a disulfidptosis model for predicting glioma prognosis, immunotherapy response, and drug.
Journal: iScience
In common: WGCNA, Harmony, limma, 7 other tools
[10] doi:10.1016/j.isci.2026.115196 [code]: Transcriptional and cellular maturation of the chick spinal cord in the context of distinct neuromuscular circuits.
Journal: iScience
In common: WGCNA, SingleCellExperiment, limma, 7 other tools

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.

Ischemic injury triggers a protective microglial phenotype in models of Aβ pathology.

Code ↔ Paper

Paper

The authors' code

Overview

Abstract

Repositories

KJPMolgenLab/scSeq_Hefendehl

gitlab.mpcdf.mpg.de/mpibr/scic/roioverstack

Code availability

Tracing map

Data

Datasets cited

Data Availability Statement

Versions

Version 1, 27 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