OSCR

A reproducible three-dimensional model of human brain tissue to investigate physiological and disease-associated microglia phenotypes.

Code ↔ Paper

18 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 18 matches · 1 of them tie a paragraph to a whole file, not to given lines: a weak match, whose lines are not tinted
  1. [1] § Results › 3BTMs enable studying AD-relevant perturbations and drug responses in human brain cells ↔ analysis_3BTM_rhapsody/analysis_3BTM_downstream.Rmd, lines 6167–6187 · score 0.85 · apoptotic signaling pathways, MHC class, protein folding, negative regulation, neuron projection, stress
  2. [2] § Methods › scRNA-seq › Downstream analysis of the BD Rhapsody datasets ↔ install_bioc.R, the whole file · a weak match · score 0.76 · ComplexHeatmap, UCell, clusterProfiler, fgsea, MAST, seq
  3. [3] § Methods › scRNA-seq › Downstream analysis of the 10X dataset ↔ analysis_3BTM_10x/Figure5_scRNA_microglia_KI.Rmd, lines 787–822 · score 0.76 · downregulated DEGs, min.pct, FindMarkers, WT Resting, MAST, DEAs
  4. [4] § Methods › scRNA-seq › Downstream analysis of the BD Rhapsody datasets ↔ analysis_3BTM_10x/Figure5_scRNA_microglia_KI.Rmd, lines 82–219 · score 0.72 · FindMarkers, VlnPlot, ribosomal, Seurat, downstream, UMIs
  5. [5] § Methods › scRNA-seq › Downstream analysis of the BD Rhapsody datasets ↔ analysis_3BTM_rhapsody/analysis_3BTM_downstream.Rmd, lines 5642–5690 · score 0.72 · latent.vars, FindMarkers, Sample_univoque, MAST, Rhapsody, threshold
  6. [6] § Methods › scRNA-seq › Downstream analysis of the 10X dataset ↔ analysis_3BTM_10x/Figure3_scRNA_microglia_maturation_clean.Rmd, lines 246–299 · score 0.66 · k.weight, IntegrateData, CCA, transform, Seurat, Downstream
  7. [7] § Results › iMGs in AD-modeling 3BTMs recapitulate disease-related changes ↔ analysis_3BTM_10x/Figure5_scRNA_microglia_KI.Rmd, lines 1238–1261 · score 0.64 · AD risk genes, KI cluster, WT resting, BLNK, SORL1, TREM2
  8. [8] § Methods › scRNA-seq › Downstream analysis of the 10X dataset ↔ analysis_3BTM_10x/Figure5_scRNA_microglia_KI.Rmd, lines 82–219 · score 0.64 · standard Seurat, workflow, ribosomal, downstream, UMIs, mitochondrial
  9. [9] § Results › 3BTMs enable studying AD-relevant perturbations and drug responses in human brain cells ↔ analysis_3BTM_rhapsody/analysis_3BTM_downstream.Rmd, lines 6267–6333 · score 0.63 · MHC protein complex, innate immune response, GO term, phagocytosis, aducanumab, signal
  10. [10] § Results › iMGs, iNEs and iAS show reproducible, AD-related changes across cell lines ↔ analysis_3BTM_rhapsody/analysis_3BTM_downstream.Rmd, lines 6023–6049 · score 0.59 · unfolded protein response, endoplasmic reticulum, neurons, GO
  11. [11] § Results › iMGs mature in 3BTMs and adopt transcriptomic signatures similar to the human brain ↔ analysis_3BTM_10x/Figure3_scRNA_microglia_maturation_clean.Rmd, lines 2023–2060 · score 0.56 · adult human, signatures derived, scRNA, fetal, maturation, seq
  12. [12] § Results › Brain cells in 3BTMs show increased maturity, homeostasis and functionality compared to 2D co-cultures ↔ analysis_3BTM_3Dvs2D/analysis_3BTM_3Dvs2D_minimal.Rmd, lines 3259–3282 · score 0.56 · ATP1A2, SLC2A1, MT3
  13. [13] § Methods › scRNA-seq › Downstream analysis of the 10X dataset ↔ analysis_3BTM_10x/Figure3_scRNA_microglia_maturation_clean.Rmd, lines 595–616 · score 0.56 · ALDH1A1, SOX11, TYROBP, AQP4, canonical, downstream
  14. [14] § Results › CRISPR–Cas9-mediated insertion of synergistic APP mutations to model AD ↔ analysis_3BTM_10x/Figure5_scRNA_microglia_KI.Rmd, lines 241–274 · score 0.53 · causing mutations, stem cells, iPSC, APP, knocked, cultured
  15. [15] § Results › iMGs in AD-modeling 3BTMs recapitulate disease-related changes ↔ analysis_3BTM_10x/Figure5_scRNA_microglia_KI.Rmd, lines 1238–1261 · score 0.52 · AD risk genes, WT resting, Cell state, Violin, Figure 5, subset
  16. [16] § Results › Brain cells in 3BTMs show increased maturity, homeostasis and functionality compared to 2D co-cultures ↔ analysis_3BTM_rhapsody/analysis_3BTM_downstream.Rmd, lines 6167–6187 · score 0.52 · synaptic signaling, synapse organization, pathways, regulation, protein
  17. [17] § Results › iMGs, iNEs and iAS show reproducible, AD-related changes across cell lines ↔ analysis_3BTM_3Dvs2D/analysis_3BTM_3Dvs2D_minimal.Rmd, lines 2372–2397 · score 0.51 · BCL11B, mid, GABAergic, layer, SATB2, cortical
  18. [18] § Results › Brain cells in 3BTMs show increased maturity, homeostasis and functionality compared to 2D co-cultures ↔ analysis_3BTM_3Dvs2D/analysis_3BTM_3Dvs2D_minimal.Rmd, lines 2886–2912 · score 0.51 · synaptic signaling, synapse organization, regulation, protein

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 authors' code

R Markdown · 1,404 lines · 50 KB · no license · 6 matches

  1. ---
  3. title: "Characterization of the KI 3months model"
  4. author: ""
  5. date: "`r format(Sys.time(), '%d %B, %Y')`"
  6. output:
  7. html_document:
  8. code_download: yes
  9. df_print: kable
  10. theme: united
  11. toc: yes
  12. toc_depth: 6
  13. toc_float: yes
  14. output: html_document
  15. editor_options:
  16. chunk_output_type: inline
  17. ---
  20. # Load packages
  22. ```{r}
  23. library(Seurat)
  24. library(SeuratDisk)
  25. library(ggplot2)
  26. library(SingleR)
  27. library(biomaRt)
  28. library(org.Hs.eg.db)
  29. library(clusterProfiler)
  30. library(RColorBrewer)
  31. library(ggpubr)
  32. library(reshape2)
  33. library(org.Hs.eg.db)
  34. library(tidyverse)
  35. library(gridExtra)
  36. library(ComplexHeatmap)
  37. library(readxl)
  38. library(scCustomize)
  39. library(data.table)
  40. library(ggVennDiagram)
  41. library(slingshot)
  42. library(ggplotify)
  43. library(patchwork)
  44. library(AUCell)
  45. library(ggVennDiagram)
  46. library(nichenetr)
  47. library(viridis)
  48. library(dittoSeq)
  49. library(ggvolc)
  50. library(cowplot)
  51. library(reticulate)
  52. library(sceasy)
  53. library(MAST)
  54. library(UpSetR)
  55. library(sceasy)
  56. library(UCell)
  57. library(PCAtools)
  58. library(scIntegrationMetrics)
  59. library(topGO)
  60. library(harmony)
  61. library(genesorteR)
  62. library(ggVolcano)
  63. library(speckle)
  64. library(monocle3)
  65. library(monocle)
  66. library(SeuratWrappers)
  67. library(SingleCellExperiment)
  68. library(Nebulosa)
  69. library(org.Mm.eg.db)
  71. BiocManager::install("GeneOverlap")
  72. library("GeneOverlap")
  74. devtools::install_github("mw201608/SuperExactTest")
  75. library(SuperExactTest)
  77. sess.i <- sessionInfo()
  78. ```
  82. # Define custom functions
  84. ```{r}
  87. # barplot of DEGs from filtered FindMarkers output
  88. degs_barplot <- function(degs_df){
  89. h <- degs_df[degs_df$Regulation != "none",]
  90. f_u <- table(h$Regulation)
  91. f_uh <- data.frame(status = names(f_u),
  92. total = as.vector(f_u))
  93. degs_bp <- ggplot(f_uh, aes(x = status, y = total, fill = status))+
  94. geom_bar(stat= "identity", width = 0.9)+
  95. geom_text(aes(label = total), vjust= -0.2, size = 6)+
  96. theme_classic()+
  97. theme(aspect.ratio = 1.9,
  98. legend.position = 'none',
  99. axis.text.x = element_text(size=20, angle=0),
  100. axis.text.y = element_text(size=20, angle=0),
  101. axis.title = element_text(size = 20))+
  102. scale_fill_manual(values = c("#00798c","#cc0000"))+
  103. ylab("number of DEGs")+
  104. xlab(label= NULL)
  105. return(degs_bp)
  106. }
  108. # Filtering FindMarkers output to get up and down regulated genes
  110. up_down_degs <- function(degs_df){
  111. degs.list <- degs_df
  112. degs.list$genes <- rownames(degs_df)
  113. degs.list$Significance <- "none"
  114. degs.list$Significance <- ifelse(degs.list$p_val_adj < 0.01, "yes", degs.list$Significance)
  115. degs.list$Significance <- ifelse(degs.list$p_val_adj > 0.01, "no", degs.list$Significance)
  116. degs.list$Regulation <- "none"
  117. degs.list$Regulation <- ifelse(degs.list$avg_log2FC > 0.25 & degs.list$p_val_adj < 0.01, "up", degs.list$Regulation)
  118. degs.list$Regulation <- ifelse(degs.list$avg_log2FC < -0.25 & degs.list$p_val_adj < 0.01, "down", degs.list$Regulation)
  119. degs.list <- subset(degs.list, degs.list$Regulation != "none")
  120. return(degs.list)
  121. }
  125. # seuratv3 integration wrapper
  127. integr_wrap <- function(sobject.list, nfeatures, k_anchor, k_filter, k_score ){
  128. # select features that are repeatedly variable across datasets for integration
  129. features <- SelectIntegrationFeatures(object.list = sobject.list, nfeatures = 2000)
  130. features.to.integrate <- intersect(rownames(sobject.list[[2]]), rownames(sobject.list[[1]]))
  131. uglia.anchors <- FindIntegrationAnchors(object.list = sobject.list, anchor.features = features, k.anchor = k_anchor, k.filter = k_filter, k.score = k_score)
  132. # this command creates an 'integrated' data assay
  133. uglia.combined <- IntegrateData(anchorset = uglia.anchors, features.to.integrate = features.to.integrate)
  134. # specify that we will perform downstream analysis on the corrected data note that the
  135. # original unmodified data still resides in the 'RNA' assay
  136. DefaultAssay(uglia.combined) <- "integrated"
  137. # Run the standard workflow for visualization and clustering
  138. uglia.combined <- ScaleData(uglia.combined, verbose = FALSE)
  139. uglia.combined <- RunPCA(uglia.combined, npcs = 30, verbose = FALSE)
  140. uglia.combined <- RunUMAP(uglia.combined, reduction = "pca", dims = 1:30)
  141. uglia.combined <- FindNeighbors(uglia.combined, reduction = "pca", dims = 1:30)
  142. uglia.combined <- FindClusters(uglia.combined, resolution = 0.5)
  143. return(uglia.combined)
  144. }
  148. # dotplot to display GO enrichment results
  150. ego_dp <- function(df){
  151. plot <- df %>%
  152. arrange(Description) %>%
  153. mutate(status = factor(status, levels = c("up", "down"))) %>%
  154. ggplot(aes(x = status, y = Description, color = `qvalue`, size = Count)) +
  155. geom_point() +
  156. scale_color_gradient(low = "#006f62", high = "#b6d7a8") +
  157. ylab("") +
  158. xlab("")
  159. return(plot)
  160. }
  164. # QC plots
  166. plot_QC <- function(seurat){
  168. #Features (genes)
  169. n_feat_cl <- VlnPlot(seurat,
  170. features = c("nFeature_RNA"),
  171. pt.size = 0) +
  172. geom_hline(yintercept = 250, linetype = "dashed", colour = "red")+
  173. ylim(500, 10000)+
  174. NoLegend()
  175. #Counts (UMIs)
  176. n_count_cl <- VlnPlot(seurat,
  177. features = c("nCount_RNA"),
  178. pt.size = 0)+
  179. geom_hline(yintercept = 15000, linetype = "dashed", colour = "red")+
  180. ylim(500, 80000)+
  181. NoLegend()
  183. #Percent of mitochondrial reads
  185. # mito cutoff 10%
  186. mito_pt_cl <- VlnPlot(seurat,
  187. features = c("percent.mt"),
  188. pt.size = 0) +
  189. geom_hline(yintercept = 7, linetype = "dashed", colour = "red")+
  190. ylim(0, 40)+
  191. NoLegend()
  193. #Percent of ribosomal reads
  194. seurat[["percent.rb"]] <- PercentageFeatureSet(seurat, pattern = "^RP[SL]")
  195. ribo_pt_cl <- VlnPlot(seurat,
  196. features = c("percent.rb"),
  197. pt.size = 0) +
  198. geom_hline(yintercept = 25, linetype = "dashed", colour = "red")+
  199. ylim(0, 40)+
  200. NoLegend()
  202. plots <- list(n_feat_cl, n_count_cl, mito_pt_cl, ribo_pt_cl)
  203. print(plots)
  204. }
  207. # standard seurat pre-processing wrapper
  208. standard_preprocessing <- function(sobj){
  209. sobj <- NormalizeData(sobj, normalization.method = "LogNormalize", scale.factor = 10000)
  210. sobj <- FindVariableFeatures(sobj, selection.method = "vst", nfeatures = 2000)
  211. all.genes <- rownames(sobj)
  212. sobj <- ScaleData(sobj, features = all.genes_mancuso)
  213. sobj <- RunPCA(sobj, features = VariableFeatures(object = sobj))
  214. sobj <- FindNeighbors(sobj, reduction = "pca",dims = 1:50)
  215. sobj <- FindClusters(sobj, resolution = 0.5)
  216. sobj <- RunUMAP(sobj, reduction = "pca",dims = 1:50, seed.use = 42)
  217. return(sobj)
  218. }
  219. ```
  221. # Define colors
  223. ```{r}
  224. pal <- c("#E1BE6A","#004949","#009292","#ff6db6","#ffb6db",
  225. "#490092","#006ddb","#b66dff","#6db6ff","#b6dbff",
  226. "#920000","#924900","#db6d00","#24ff24","#ffff6d")
  230. pal_viridis <- viridis(n = 10, option = "D")
  232. sobj_3mo_pal <- c("WT Resting" = "#228b22",
  233. "KI Cluster 2" = "#ffb90f",
  234. "KI Cluster 1" = "#eae4b7",
  235. "Chemokine 3mo" = "#74d600",
  236. "Glycolytic 3mo" = "#b0bf1a",
  237. "Proliferative 3mo" = "#b2ffff")
  239. ```
  241. #_____________________________________________________________________________________________________________________________________________________
  242. #DATASET description
  243. #_____________________________________________________________________________________________________________________________________________________
  245. All 10X samples derive from human stem cell (iPSC-) derived 3D cultures made of cortical neurons (NE), astrocytes (AS) and microglia (MG). The cell types were differentiated separately, then aggregated into a spheroid-like culture, and grown for 1 month (1mo) or 3 months (3mo). For the experiment, wildtype (WT) and knock-in (KI) cultures carrying three familial Alzheimer’s disease-causing mutations in the APP gene are analyzed. The experiment was performed as indicated below:
  247. o Sample 3: WT 1mo (pooled 2 biological replicates with hashtags)
  248. o Sample 4: WT 3mo (pooled 2 biological replicates with hashtags)
  249. o Sample 5: KI 3mo (pooled 2 biological replicates with hashtags)
  252. #Load Data
  254. The code chunks relative to the loading of external objects(.csv, .xlsx) containing signature gene lists have been silenced. When loading the minimal environment it won't be necessary to load seurat objects and running pre-processing and differential expression aanlysis code chunks.
  256. ```{r}
  257. #load main seurat objects
  258. seurat_objects_figure6 <- readRDS("/path/seurat_figure6.rds")
  260. #Load environment:
  261. #In this environment are included:
  262. # - main seurat objects produced
  263. # - all signatures used for scoring
  264. # - differential analysis code
  266. #load("/path/figure6_env.RData")
  268. #extract each object from list (if environment is not loaded)
  269. #sobj_sub <- seurat_objects_figure6[[1]] # The code for pre-processing and subsetting this object are reported in Figure 3 script
  270. #sobj_mg <- seurat_objects_figure6[[2]] #seurat object filtered to have only microglia cells
  271. #sobj_3mo <- seurat_objects_figure6[[3]] # subset seurat object for 3 months old WT and KI neurocultures
  274. ```
  276. # Pre-processing
  278. ```{r}
  279. #Normalization
  280. sobj_mg <- NormalizeData(sobj_mg, normalization.method = "LogNormalize", scale.factor = 10000)
  281. #Variable features exploration
  282. sobj_mg <- FindVariableFeatures(sobj_mg, selection.method = "vst", nfeatures = 2000)
  283. # Identify the 10 most highly variable genes
  284. top10 <- head(VariableFeatures(sobj_mg), 10)
  285. #Scaling
  286. all.genes <- rownames(sobj_mg)
  287. sobj_mg <- ScaleData(sobj_mg, features = all.genes)
  288. #Dimensional reduction
  289. sobj_mg <- RunPCA(sobj_mg, features = VariableFeatures(object = sobj_mg))
  290. print(sobj_mg[["pca"]], dims = 1:5, nfeatures = 5)
  291. DimPlot(sobj_mg, reduction = "pca")
  292. #determine informative PCs
  293. ElbowPlot(sobj_mg)
  294. ```
  296. ```{r}
  297. set.seed(42)
  298. # Find clusters
  299. sobj_mg <- FindNeighbors(object=sobj_mg, dims=1:10)
  300. # Find clusters for several resolutions and choose appropriate one for downstream analysis
  301. sobj_mg_cl_tr <- FindClusters(object = sobj_mg, resolution = c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0))
  302. # Plot and save several resolutions
  303. sobj_mg_cl_tr <- RunUMAP(object=sobj_mg_cl_tr, dims=1:10, seed.use = 42)
  304. DimPlot(object=sobj_mg_cl_tr, reduction="umap", group.by="RNA_snn_res.0.1", label=T, seed = 42) + ggtitle("Res 0.1")
  305. DimPlot(object=sobj_mg_cl_tr, reduction="umap", group.by="RNA_snn_res.0.2", label=T, seed = 42) + ggtitle("Res 0.2")
  306. DimPlot(object=sobj_mg_cl_tr, reduction="umap", group.by="RNA_snn_res.0.3", label=T, seed = 42) + ggtitle("Res 0.3")
  307. DimPlot(object=sobj_mg_cl_tr, reduction="umap", group.by="RNA_snn_res.0.4", label=T, seed = 42) + ggtitle("Res 0.4")
  308. DimPlot(object=sobj_mg_cl_tr, reduction="umap", group.by="RNA_snn_res.0.5", label=T, seed = 42) + ggtitle("Res 0.5")
  309. DimPlot(object=sobj_mg_cl_tr, reduction="umap", group.by="RNA_snn_res.0.6", label=T, seed = 42) + ggtitle("Res 0.6")
  310. DimPlot(object=sobj_mg_cl_tr, reduction="umap", group.by="RNA_snn_res.0.8", label=T, seed = 42) + ggtitle("Res 0.8")
  311. DimPlot(object=sobj_mg_cl_tr, reduction="umap", group.by="RNA_snn_res.1", label=T, seed = 42) + ggtitle("Res 1.0")
  312. DimPlot(object=sobj_mg_cl_tr, reduction="umap", group.by="RNA_snn_res.1.2", label=T, seed = 42) + ggtitle("Res 1.2")
  313. DimPlot(object=sobj_mg_cl_tr, reduction="umap", group.by="RNA_snn_res.1.4", label=T, seed = 42) + ggtitle("Res 1.4")
  314. DimPlot(object=sobj_mg_cl_tr, reduction="umap", group.by="RNA_snn_res.0.7", label=T, seed = 42) + ggtitle("Res 0.7")
  316. #Here the resolution to be used was chosen empirically to avoid over-clustering. In case of re-analysis of published dataset the resolution was chosen accordingly to the information on the original publication if mentioned with the aim of replicating the published data
  317. ```
  319. ```{r}
  320. #Find clusters with chosen resolution
  321. sobj_mg <- FindClusters(object = sobj_mg, resolution=0.5)
  322. remove(sobj_mg_cl_tr) #removing object to limit memory usage
  323. ```
  325. ```{r}
  326. set.seed(42)
  327. #UMAP
  328. sobj_mg <- RunUMAP(sobj_mg, dims = 1:10, seed.use = 42)
  329. #plots
  330. DimPlot(sobj_mg, reduction = "umap", cols = pal, seed = 1) + theme(aspect.ratio=1)
  331. DimPlot(sobj_mg, split.by = "condition", cols = pal, seed = 1) + theme(aspect.ratio=1)
  332. #total number of cells and number of cells per cluster
  333. table(Idents(sobj_mg))
  334. ncol(sobj_mg)
  335. ```
  338. ```{r, fig.height=7, fig.width=20}
  339. #UMAP plot splitted by replicate
  340. DimPlot(sobj_mg, split.by = "orig.ident", cols = pal, seed = 1) + theme(aspect.ratio=1)
  341. ```
  343. # QC per cluster (only microglia)
  345. ```{r, fig.height=7, fig.width=20}
  346. #QC plots per condition
  347. Idents(sobj_mg) <- "condition"
  348. sobj_mg_QC_c <-plot_QC(sobj_mg)
  350. #QC plots per cluster
  351. Idents(sobj_mg) <- "seurat_clusters"
  352. sobj_mg_QC_cl <-plot_QC(sobj_mg)
  353. ```
  355. ## Subset to exclude sample enriched clusters
  357. some clusters may be sample specific and so not really informative in the downstream analysis. Here we excluded one cluster with such characteristics to avoid its inclusion in the downstream process.
  359. The seurat object was subject to pre-processing following the subset
  361. ```{r}
  362. set.seed(42)
  363. #subset
  364. clusters_to_keep <- c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9")
  365. sobj_mg_filt <- subset(sobj_mg, subset = seurat_clusters %in% clusters_to_keep)
  366. #Normalization
  367. sobj_mg_filt <- NormalizeData(sobj_mg_filt, normalization.method = "LogNormalize", scale.factor = 10000)
  368. #Variable features exploration
  369. sobj_mg_filt <- FindVariableFeatures(sobj_mg_filt, selection.method = "vst", nfeatures = 2000)
  370. # Identify the 10 most highly variable genes
  371. top10 <- head(VariableFeatures(sobj_mg_filt), 10)
  372. #Scaling
  373. all.genes <- rownames(sobj_mg_filt)
  374. sobj_mg_filt <- ScaleData(sobj_mg_filt, features = all.genes)
  375. #Dimensional reduction
  376. sobj_mg_filt <- RunPCA(sobj_mg_filt, features = VariableFeatures(object = sobj_mg_filt))
  377. print(sobj_mg_filt[["pca"]], dims = 1:5, nfeatures = 5)
  378. DimPlot(sobj_mg_filt, reduction = "pca")
  379. #determine informative PCs
  380. ElbowPlot(sobj_mg_filt)
  381. # Find clusters
  382. sobj_mg_filt <- FindNeighbors(object=sobj_mg_filt, dims=1:10)
  383. # Explore clustering resolutions
  384. sobj_mg_filt_t <- FindClusters(object = sobj_mg_filt, resolution = c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0))
  385. # Plot
  386. sobj_mg_filt_t <- RunUMAP(object=sobj_mg_filt_t, dims=1:10, seed.use =42)
  387. DimPlot(object=sobj_mg_filt_t, reduction="umap", group.by="RNA_snn_res.0.2", label=T,seed = 42) + ggtitle("Res 0.1")
  388. DimPlot(object=sobj_mg_filt_t, reduction="umap", group.by="RNA_snn_res.0.2", label=T, seed = 42) + ggtitle("Res 0.2")
  389. DimPlot(object=sobj_mg_filt_t, reduction="umap", group.by="RNA_snn_res.0.2", label=T, seed = 42) + ggtitle("Res 0.3")
  390. DimPlot(object=sobj_mg_filt_t, reduction="umap", group.by="RNA_snn_res.0.4", label=T, seed = 42) + ggtitle("Res 0.4")
  391. DimPlot(object=sobj_mg_filt_t, reduction="umap", group.by="RNA_snn_res.0.5", label=T, seed = 42) + ggtitle("Res 0.5")
  392. DimPlot(object=sobj_mg_filt_t, reduction="umap", group.by="RNA_snn_res.0.6", label=T, seed = 42) + ggtitle("Res 0.6")
  393. DimPlot(object=sobj_mg_filt_t, reduction="umap", group.by="RNA_snn_res.0.8", label=T, seed = 42) + ggtitle("Res 0.8")
  394. DimPlot(object=sobj_mg_filt_t, reduction="umap", group.by="RNA_snn_res.1", label=T, seed = 42) + ggtitle("Res 1.0")
  395. DimPlot(object=sobj_mg_filt_t, reduction="umap", group.by="RNA_snn_res.1.2", label=T, seed = 42) + ggtitle("Res 1.2")
  396. DimPlot(object=sobj_mg_filt_t, reduction="umap", group.by="RNA_snn_res.1.4", label=T, seed = 42) + ggtitle("Res 1.4")
  399. #Here the resolution to be used was chosen empirically to avoid over-clustering. In case of re-analysis of published dataset the resolution was chosen accordingly to the information on the original publication if mentioned with the aim of replicating the published data
  400. ```
  402. ```{r}
  403. set.seed(42)
  404. # Find clusters this time using the good resolution. Old resolution data is automatically overwritten.
  405. sobj_mg_filt <- FindClusters(object = sobj_mg_filt, resolution=0.5)
  406. #UMAP
  407. sobj_mg_filt <- RunUMAP(sobj_mg_filt, dims = 1:10, seed.use = 42)
  409. DimPlot(sobj_mg_filt, reduction = "umap", cols = pal, seed = 1) + theme(aspect.ratio=1)
  410. DimPlot(sobj_mg_filt, split.by = "condition", cols = pal, seed = 1) + theme(aspect.ratio=1)
  412. remove(sobj_mg_filt_t) #remove objects to limit memory usage
  413. ```
  415. ```{r, fig.height=7, fig.width=20}
  416. DimPlot(sobj_mg_filt, split.by = "orig.ident", cols = pal, seed = 1) + theme(aspect.ratio=1)
  417. ```
  419. # Cell type specific markers (feature plots)
  421. Check for the presence of only microglia in the subset.
  424. ```{r, fig.height=8, fig.width=8}
  425. tyrobp <- FeaturePlot(sobj_mg_filt, features = c("TYROBP")) + theme_void()
  426. sox <- FeaturePlot(sobj_mg_filt, features = c("SOX11")) + theme(aspect.ratio=1)
  427. astro <- FeaturePlot(sobj_mg_filt, features = c("ALDH1A1")) + theme(aspect.ratio=1)
  429. grid.arrange(sox, astro)
  430. grid.arrange(tyrobp)
  431. ```
  434. # Identify cluster biomarkers
  436. To identify different cell states in the microglia data the top 10 differentially expressed genes per cluster were calculated and their averaged expression shown in a dotplot
  438. ```{r, fig.height= 6, fig.width= 18}
  439. #set idents to clusters
  440. Idents(sobj_mg_filt) <- "seurat_clusters"
  441. #Find cluster enriched markers
  442. s.markers <- FindAllMarkers(sobj_mg_filt, assay = "RNA", only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
  443. #select top 10 markers per cluster to plot
  444. s.markers %>%
  445. group_by(cluster) %>%
  446. top_n(n = 10, wt = avg_log2FC) -> top10_s
  447. ```
  450. #___________________________________________________________________________________________________________________________
  451. # SUBSET FOR 3mo (KI, WT)
  452. #_____________________________________________________________________________________________________________________________________________________________________
  454. # Subset for 3mo cells
  456. If loading directly the sobj_3mo object at the beginning of the script running the subsetting and pre-processing code chunks is not necessary
  458. ```{r}
  459. #subset by time
  460. Idents(sobj_mg_filt) <- "time"
  461. sobj_3mo <- subset(sobj_mg_filt, idents= c("3mo"))
  462. #set and check subset idents
  463. Idents(sobj_3mo) <- "genotype"
  464. levels(Idents(sobj_3mo))
  465. ```
  467. ## Pre-processing
  469. ```{r}
  470. #Normalization
  471. sobj_3mo <- NormalizeData(sobj_3mo, normalization.method = "LogNormalize", scale.factor = 10000)
  472. #Find Variable features
  473. sobj_3mo <- FindVariableFeatures(sobj_3mo, selection.method = "vst", nfeatures = 2000)
  474. # Identify the 10 most highly variable genes
  475. top10 <- head(VariableFeatures(sobj_3mo), 10)
  476. plot1 <- VariableFeaturePlot(sobj_3mo)
  477. LabelPoints(plot = plot1, points = top10, repel = TRUE)
  478. ## Scaling
  479. all.genes <- rownames(sobj_3mo)
  480. sobj_3mo <- ScaleData(sobj_3mo, features = all.genes)
  481. ## PCA
  482. sobj_3mo <- RunPCA(sobj_3mo, features = VariableFeatures(object = sobj_3mo))
  483. print(sobj_3mo[["pca"]], dims = 1:5, nfeatures = 5)
  484. DimPlot(sobj_3mo, reduction = "pca")
  485. #determine more informative principal components
  486. ElbowPlot(sobj_3mo)
  488. ```
  492. ```{r}
  493. ## Clustering
  494. set.seed(42)
  495. # Find clusters
  496. sobj_3mo <- FindNeighbors(object = sobj_3mo, dims=1:10)
  497. #Plot different clustering resolutions
  498. sobj_3mo_t <- FindClusters(object = sobj_3mo, resolution = c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0), verbose = F)
  499. # Plots
  500. sobj_3mo_t <- RunUMAP(object= sobj_3mo_t, dims=1:10, seed.use = 42)
  501. DimPlot(object=sobj_3mo_t, reduction="umap", group.by="RNA_snn_res.0.1", label=T, seed = 42) + ggtitle("Res 0.1")
  502. DimPlot(object=sobj_3mo_t, reduction="umap", group.by="RNA_snn_res.0.2", label=T, seed = 42) + ggtitle("Res 0.2")
  503. DimPlot(object=sobj_3mo_t, reduction="umap", group.by="RNA_snn_res.0.3", label=T, seed = 42) + ggtitle("Res 0.3")
  504. DimPlot(object=sobj_3mo_t, reduction="umap", group.by="RNA_snn_res.0.4", label=T, seed = 42) + ggtitle("Res 0.4")
  505. DimPlot(object=sobj_3mo_t, reduction="umap", group.by="RNA_snn_res.0.5", label=T, seed = 42) + ggtitle("Res 0.5")
  506. DimPlot(object=sobj_3mo_t, reduction="umap", group.by="RNA_snn_res.0.6", label=T, seed = 42) + ggtitle("Res 0.6")
  507. DimPlot(object=sobj_3mo_t, reduction="umap", group.by="RNA_snn_res.0.8", label=T, seed = 42) + ggtitle("Res 0.8")
  508. DimPlot(object=sobj_3mo_t, reduction="umap", group.by="RNA_snn_res.1", label=T, seed = 42) + ggtitle("Res 1.0")
  509. DimPlot(object=sobj_3mo_t, reduction="umap", group.by="RNA_snn_res.1.2", label=T, seed = 42) + ggtitle("Res 1.2")
  510. DimPlot(object=sobj_3mo_t, reduction="umap", group.by="RNA_snn_res.1.4", label=T, seed = 42) + ggtitle("Res 1.4")
  511. ```
  515. ```{r}
  516. ## UMAP
  517. set.seed(42)
  519. #remove object to reduce memory usage
  520. rm(sobj_3mo_t)
  521. #choose dimensional reduction to use
  522. sobj_3mo <- FindNeighbors(object = sobj_3mo, dims=1:10)
  523. sobj_3mo <- FindClusters(object = sobj_3mo, resolution = c(0.4), verbose = F)
  524. sobj_3mo <- RunUMAP(object= sobj_3mo, dims=1:10)
  525. #reset idents
  526. Idents(sobj_3mo) <- "seurat_clusters"
  527. ```
  529. ```{r}
  530. #plot
  531. DimPlot(sobj_3mo, split.by = "condition", cols = pal, seed = 1) + theme(aspect.ratio=1)
  532. ```
  534. ```{r, fig.height=5, fig.width=10}
  535. #plot UMAP splitted by sample
  536. DimPlot(sobj_3mo, split.by = "orig.ident", cols = pal, seed = 1) + theme(aspect.ratio=1)
  537. ```
  539. # Characterization of the clusters
  541. ## Top markers
  543. ```{r}
  544. #set idents
  545. Idents(sobj_3mo) <- "seurat_clusters"
  546. #Find cluster enriched markers
  547. s3mo.markers <- FindAllMarkers(sobj_3mo, assay = "RNA", only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
  548. #select top 10 markers per cluster to plot
  549. s3mo.markers %>%
  550. group_by(cluster) %>%
  551. top_n(n = 10, wt = avg_log2FC) -> top10_s3mo
  552. ```
  554. ```{r, fig.height= 6, fig.width= 18}
  555. #extract genes to plot
  556. to_plot_s3mo <- unique(top10_s3mo$gene)
  557. #dotplot
  558. plot_mpc_s3mo <- DotPlot(sobj_3mo, features = to_plot_s3mo, group.by = "seurat_clusters", cols = c("RdBu"))+
  559. theme(axis.text.x = element_text(angle = 90))
  561. #cell number per cluster
  562. table(Idents(sobj_3mo))
  564. ```
  568. ## Functional enrichments in specific clusters
  571. ## Cell states annotation
  573. The annotation was given looking at the markers per cluster and in accordance with the annotation of the Figure 3 UMAP
  575. ```{r}
  576. Idents(sobj_3mo) <- "seurat_clusters"
  578. cell_states <- c( "0"="WT Resting",
  579. "1"="KI Cluster 1",
  580. "2"="KI Cluster 2",
  581. "3"="Chemokine 3mo",
  582. "4"="Glycolytic 3mo",
  583. "5"="Proliferative 3mo")
  585. names(cell_states) <- levels(sobj_3mo)
  586. sobj_3mo <- RenameIdents(sobj_3mo, cell_states)
  587. sobj_3mo$cell_states <- Idents(sobj_3mo)
  588. ```
  590. # QC per cluster
  592. ```{r, fig.height=7, fig.width=10}
  593. #QC plots per cluster
  594. #set idents
  595. Idents(sobj_3mo) <- "seurat_clusters"
  596. #plots
  597. sobj_3mo_QC_cl <- plot_QC(sobj_3mo)
  599. a <- sobj_3mo_QC_cl[[1]]
  600. b <- sobj_3mo_QC_cl[[2]]
  601. c <- sobj_3mo_QC_cl[[3]]
  602. d <- sobj_3mo_QC_cl[[4]]
  604. plot_grid(a, b, c, d)
  605. ```
  607. ```{r}
  608. #here I give names to metadata columns to store percentage of mitochondrial and ribosomal genes in the object before subsetting. The information will be retained in all subsets of the object
  610. #sobj_sub[["percent.mt"]] <- PercentageFeatureSet(sobj_sub, pattern = "^MT-")
  611. #sobj_sub[["percent.rb"]] <- PercentageFeatureSet(sobj_sub, pattern = "^RP[SL]")
  613. #here I visualize such parameters in all the clusters of the 3mo subset
  614. mito <- VlnPlot(object = sobj_3mo, features = "percent.mt", group.by = "seurat_clusters", pt.size = 0)+ theme(legend.position="none") + ggtitle("Percent of mitochondrial counts") + ylab("% of all counts") + xlab("")+ylim(0, 40)+geom_hline(yintercept = 10)
  616. ribo <- VlnPlot(object = sobj_3mo, features = "percent.rb", group.by = "seurat_clusters", pt.size = 0)+ theme(legend.position="none") + ggtitle("Percent of ribosomal counts") + ylab("% of all counts") + xlab("")+ylim(0, 40)+geom_hline(yintercept = 20)
  618. mito+ribo
  619. ```
  622. #______________________________________________________________________________
  623. # Figure
  625. Umap plots for the 3 mo subset
  627. ```{r}
  628. Idents(sobj_3mo) <- "cell_states"
  630. umap_3mo <-DimPlot(sobj_3mo, reduction = "umap", group.by = "cell_states", pt.size = 0.5,
  631. cols = c(sobj_3mo_pal))+
  632. theme(aspect.ratio = 1,
  633. legend.text = element_text(size = 15),
  634. axis.line = element_blank(),
  635. axis.title = element_blank(),
  636. axis.text = element_blank(),
  637. axis.ticks = element_blank(),
  638. legend.key.size = unit(0.5, "cm"),
  639. legend.title = element_text(color = "black", size = 18, face = 2),)+
  640. guides(color = guide_legend(override.aes = list(size = 6)))+
  641. ggtitle(element_blank())
  643. umap_3mo
  644. ```
  647. # Figure
  649. UMAP plot colored by genotype
  651. ```{r}
  652. #Highlight plots
  653. levels(Idents(sobj_3mo))
  654. Idents(sobj_3mo) <- "genotype"
  656. WT <- WhichCells(sobj_3mo, idents = c("WT"))
  658. cells_wt <- list("WT" = WT)
  660. highlight_wt <- Cell_Highlight_Plot(seurat_object = sobj_3mo, cells_highlight = cells_wt, highlight_color = c("#004D00"))+
  661. theme(aspect.ratio = 1,
  662. legend.text = element_text(size = 15),
  663. axis.line = element_blank(),
  664. axis.title = element_blank(),
  665. axis.text = element_blank(),
  666. axis.ticks = element_blank(),
  667. legend.key.size = unit(0.5, "cm"),
  668. legend.title = element_text(color = "black", size = 18, face = 2))+
  669. guides(color = guide_legend(override.aes = list(size = 6)))
  671. highlight_wt
  673. KI<- WhichCells(sobj_3mo, idents = c("KI"))
  675. cells_ki <- list("KI" = KI)
  677. highlight_ki <- Cell_Highlight_Plot(seurat_object = sobj_3mo, cells_highlight = cells_ki, highlight_color = c("#c99005"))+
  678. theme(aspect.ratio = 1,
  679. legend.text = element_text(size = 15),
  680. axis.line = element_blank(),
  681. axis.title = element_blank(),
  682. axis.text = element_blank(),
  683. axis.ticks = element_blank(),
  684. legend.key.size = unit(0.5, "cm"),
  685. legend.title = element_text(color = "black", size = 18, face = 2))+
  686. guides(color = guide_legend(override.aes = list(size = 6)))
  688. highlight_ki
  690. split_3mo_umap<- plot_grid(highlight_wt, highlight_ki, ncol= 2)
  691. split_3mo_umap
  692. ```
  693. # Figure
  694. Cell proportions
  696. ```{r}
  697. #adding replicate information to the metadata
  698. levels(as.factor([email hidden]$orig.ident))
  700. sobj_3mo$replicate <- ifelse([email hidden]$orig.ident == "JK211112_J2_WT_T1", "J1WT1", ifelse(
  701. [email hidden]$orig.ident == "JK211112_J2_WT_T2", "J1WT2", ifelse(
  702. [email hidden]$orig.ident == "JK211112_J3_KI_T1", "J3KI1", ifelse(
  703. [email hidden]$orig.ident == "JK211112_J3_KI_T2", "J3KI2", ""
  704. )
  705. )
  706. ))
  708. levels(as.factor([email hidden]$replicate))
  709. #DA
  710. propeller_3mo <- propeller(clusters = sobj_3mo$cell_states, sample = sobj_3mo$replicate, group = sobj_3mo$condition)
  712. # Plot cell type proportions
  713. propeller_prop_3mo <- plotCellTypeProps(clusters=sobj_3mo$cell_states, sample=sobj_3mo$genotype)
  715. ```
  717. ```{r}
  718. #barplot
  719. KI_cell_props <- propeller_prop_3mo+
  720. scale_fill_manual(values = c("WT Resting" = "#228b22",
  721. "KI Cluster 2" = "#ffb90f",
  722. "KI Cluster 1" = "#eae4b7",
  723. "Chemokine 3mo" = "#74d600",
  724. "Glycolytic 3mo" = "#b0bf1a",
  725. "Proliferative 3mo" = "#b2ffff"))+
  726. xlab("")+
  727. ylab("Frequency")+
  728. scale_x_discrete(limits =c("WT", "KI"))+
  729. theme(aspect.ratio = 1.5,
  730. panel.background = element_blank(),
  731. axis.line = element_line(colour = "black", linewidth = 0.5),
  732. axis.text.x = element_text(face = "bold"),
  733. axis.title.y = element_text(face = "bold"))
  736. KI_cell_props
  737. ```
  739. # Figure
  741. QC per cluster
  743. ```{r, fig.height=7, fig.width=10}
  744. #QC plots per cluster
  745. Idents(sobj_3mo) <- "cell_states"
  746. plot_grid(a, b, c, d)
  747. ```
  750. #Figure
  752. ```{r}
  753. Idents(sobj_3mo) <- "cell_states"
  755. #order sample
  756. sobj_3mo$sample <- ifelse([email hidden]$orig.ident == "JK211112_J2_WT_T1", "WT rep1", ifelse(
  757. [email hidden]$orig.ident == "JK211112_J2_WT_T2", "WT rep2", ifelse(
  758. [email hidden]$orig.ident == "JK211112_J3_KI_T1", "KI rep1", ifelse(
  759. [email hidden]$orig.ident == "JK211112_J3_KI_T2", "KI rep2", ""
  760. )
  761. )
  762. ))
  764. umap.3mo.sample <- DimPlot(sobj_3mo, reduction = "umap", group.by = "cell_states", split.by = "sample", pt.size = 0.5, cols = c(sobj_3mo_pal), ncol = 2)+
  765. theme(aspect.ratio = 1,
  766. legend.text = element_text(size = 15),
  767. axis.line = element_blank(),
  768. axis.title = element_blank(),
  769. axis.text = element_blank(),
  770. axis.ticks = element_blank(),
  771. legend.key.size = unit(0.5, "cm"),
  772. legend.title = element_text(color = "black", size = 18, face = 2),)+
  773. guides(color = guide_legend(override.aes = list(size = 6)))+
  774. ggtitle(element_blank())
  776. umap.3mo.sample
  777. ```
  778. #Figure
  780. ```{r, fig.height= 6, fig.width= 18}
  781. plot_mpc_s3mo
  782. ```
  787. #______________________________________________________________________________
  788. # DEA
  790. ## Cluster-wise comparison
  792. ## KI Cluster 2 vs WT resting
  794. ```{r}
  795. set.seed(42)
  796. #set and check idents
  797. Idents(sobj_3mo) <- "cell_states"
  798. levels(Idents(sobj_3mo))
  799. #subset the object to have only the clusters interested in the comparison
  800. sobj_3mo_sub1 <- subset(sobj_3mo, idents = c("KI Cluster 2", "WT Resting"))
  801. #assess differentially expressed genes with MAST
  802. ki.c2.wt.dea <- FindMarkers(sobj_3mo_sub1, min.pct = 0.25, logfc.threshold = 0, ident.1 = "KI Cluster 2", ident.2 = "WT Resting", test.use= "MAST")
  803. #add gene names to the output dataframe
  804. ki.c2.wt.dea$genes <- rownames(ki.c2.wt.dea)
  805. #select statistically significant genes with p.adjust < 0.01 and log2FC >< 0.25
  806. ki.c2.wt.dea_filt <- up_down_degs(ki.c2.wt.dea)
  807. #determination of markers to plot
  808. #select upregulated degs
  809. up_degs_ki.c2.wt.dea <- dplyr::filter(ki.c2.wt.dea_filt, Regulation == "up")
  810. #order in function of log2FC
  811. up_degs_ki.c2.wt.dea_ord <- dplyr::arrange(up_degs_ki.c2.wt.dea, desc(avg_log2FC))
  812. #subset for top 20 genes to display
  813. top20_up_degs_fc_ord_ki.c2.wt.dea <- up_degs_ki.c2.wt.dea_ord[1:20, ]
  814. #select downregulated degs
  815. dn_degs_ki.c2.wt.dea <- dplyr::filter(ki.c2.wt.dea_filt, Regulation == "down")
  816. #order in function of log2FC
  817. dn_degs_ki.c2.wt.dea_ord <- dplyr::arrange(dn_degs_ki.c2.wt.dea, avg_log2FC)
  818. #subset for top 20 genes to display
  819. top20_dn_degs_fc_ord_ki.c2.wt.dea <- dn_degs_ki.c2.wt.dea_ord[1:20, ]
  822. ```
  824. ## KI Cluster 1 vs WT resting
  826. ```{r}
  827. set.seed(42)
  828. #set and check idents
  829. Idents(sobj_3mo) <- "cell_states"
  830. levels(Idents(sobj_3mo))
  831. #subset the object to have only the clusters interested in the comparison
  832. sobj_3mo_sub2<- subset(sobj_3mo, idents = c("KI Cluster 1", "WT Resting"))
  833. #assess differentially expressed genes with MAST
  834. ki.c1.wt.dea <- FindMarkers(sobj_3mo_sub2, min.pct = 0.25, logfc.threshold = 0, ident.1 = "KI Cluster 1", ident.2 = "WT Resting", test.use= "MAST")
  835. ki.c1.wt.dea$genes <- rownames(ki.c1.wt.dea)
  836. #filter for statistically significant genes with p.adjust < 0.01 and log2FC >< 0.25
  837. ki.c1.wt.dea_filt <- up_down_degs(ki.c1.wt.dea)
  838. #determination of markers to plot
  839. up_degs_ki.c1.wt.dea <- dplyr::filter(ki.c1.wt.dea_filt, Regulation == "up")
  840. up_degs_ki.c1.wt.dea_ord <- dplyr::arrange(up_degs_ki.c1.wt.dea, desc(avg_log2FC))
  841. top20_up_degs_fc_ord_ki.c1.wt.dea <- up_degs_ki.c1.wt.dea_ord[1:20, ]
  843. dn_degs_ki.c1.wt.dea <- dplyr::filter(ki.c1.wt.dea_filt, Regulation == "down")
  844. dn_degs_ki.c1.wt.dea_ord <- dplyr::arrange(dn_degs_ki.c1.wt.dea, avg_log2FC)
  845. top20_dn_degs_fc_ord_ki.c1.wt.dea <- dn_degs_ki.c1.wt.dea_ord[1:20, ]
  848. ```
  850. ## Union of upregulated genes in KI clusters when compared with WT Resting
  852. ```{r}
  853. #filter for upregulated genes in KI Cluster 2
  854. up_ki.c2vswt_all <- subset(ki.c2.wt.dea_filt, Regulation == "up")
  855. #fiilter for upregulated genes in KI Cluster 1
  856. up_ki.c1vswt_all <- subset(ki.c1.wt.dea_filt, Regulation == "up")
  857. #union of upregulated genes in KI Cluster 1 and 2
  858. up_kivswt_union <- unique(c(up_ki.c2vswt_all$genes, up_ki.c1vswt_all$genes))
  859. #union of all DEGs in KI Cluster 1 and 2
  860. all_kivswt_union <- unique(c(ki.c2.wt.dea_filt$genes, ki.c1.wt.dea_filt$genes))
  861. ```
  863. # GO enrichment
  865. ## KI Cluster 2 vs WT resting
  867. ```{r}
  869. ## Define universe for GO enrichment
  870. univers <- unique(rownames(sobj_3mo@assays[["RNA"]]@counts))
  872. # Perform enrichment analysis for up-regulated genes
  873. ego_up <- enrichGO(gene = up_ki.c2vswt_all$genes,
  874. OrgDb = org.Hs.eg.db,
  875. keyType = "SYMBOL",
  876. ont = "BP",
  877. pvalueCutoff = 0.05,
  878. pAdjustMethod = "BH",
  879. universe = univers,
  880. qvalueCutoff = 0.2,
  881. minGSSize = 10,
  882. maxGSSize = 500,
  883. readable = FALSE)
  884. #add information to enrich GO result to later merge results for enriched terms relative to up- and downregulated genes
  885. ego_up_res_cmplt <- ego_up@result
  886. ego_up_res_cmplt$status <- "up"
  887. # Subset down-regulated genes
  888. dn_genes <- subset(ki.c2.wt.dea_filt, Regulation == "down")
  889. # Perform enrichment analysis for down-regulated genes
  890. ego_down <- enrichGO(gene = dn_genes$genes,
  891. OrgDb = org.Hs.eg.db,
  892. keyType = "SYMBOL",
  893. ont = "BP",
  894. pvalueCutoff = 0.05,
  895. pAdjustMethod = "BH",
  896. universe = univers,
  897. qvalueCutoff = 0.2,
  898. minGSSize = 10,
  899. maxGSSize = 500,
  900. readable = FALSE)
  901. #add information to enrich GO result to later merge results for enriched terms relative to up- and downregulated genes
  902. ego_down_res_cmplt <- ego_down@result
  903. ego_down_res_cmplt$status <- "down"
  904. # Combine the results and filter for the most relevant ones
  905. ego_combined_cmplt_ki.c2 <- dplyr::bind_rows(ego_down_res_cmplt, ego_up_res_cmplt)
  906. # filter for statistically significant enriched terms
  907. ego_combined_ki.c2 <- dplyr::filter(ego_combined_cmplt_ki.c2, qvalue < 0.2)
  909. ```
  911. ## KI Cluster 1 vs WT resting
  913. ```{r}
  914. # Define a function to perform the go enrichment analysis and generate the dot plot
  916. # Subset up-regulated genes
  917. up_genes_ki.c1 <- subset(ki.c1.wt.dea_filt, Regulation == "up")
  919. # Perform enrichment analysis for up-regulated genes
  920. ego_up_ki.c1 <- enrichGO(gene = up_genes_ki.c1$genes,
  921. OrgDb = org.Hs.eg.db,
  922. keyType = "SYMBOL",
  923. ont = "BP",
  924. pvalueCutoff = 0.05,
  925. pAdjustMethod = "BH",
  926. universe = univers,
  927. qvalueCutoff = 0.2,
  928. minGSSize = 10,
  929. maxGSSize = 500,
  930. readable = FALSE)
  932. ego_up_res_cmplt_ki.c1 <- ego_up_ki.c1@result
  933. ego_up_res_cmplt_ki.c1$status <- "up"
  935. # Subset up-regulated genes
  936. dn_genes_ki.c1 <- subset(ki.c1.wt.dea_filt, Regulation == "down")
  937. # Perform enrichment analysis for down-regulated genes
  938. ego_down_ki.c1 <- enrichGO(gene = dn_genes_ki.c1$genes,
  939. OrgDb = org.Hs.eg.db,
  940. keyType = "SYMBOL",
  941. ont = "BP",
  942. pvalueCutoff = 0.05,
  943. pAdjustMethod = "BH",
  944. universe = univers,
  945. qvalueCutoff = 0.2,
  946. minGSSize = 10,
  947. maxGSSize = 500,
  948. readable = FALSE)
  950. ego_down_res_cmplt_ki.c1 <- ego_down_ki.c1@result
  951. ego_down_res_cmplt_ki.c1$status <- "down"
  953. # Combine the results
  954. ego_combined_cmplt_ki.c1 <- dplyr::bind_rows(ego_down_res_cmplt_ki.c1, ego_up_res_cmplt_ki.c1)
  955. ego_combined_ki.c1 <- dplyr::filter(ego_combined_cmplt_ki.c1, qvalue < 0.2)
  956. ```
  959. #______________________________________________________________________________
  960. # Figure
  962. ## DEGs between KI Cluster 2 and WT resting
  964. ## Highlight plots
  966. ```{r}
  967. #Highlight plots
  968. levels(Idents(sobj_3mo))
  969. Idents(sobj_3mo) <- "cell_states"
  971. WT <- WhichCells(sobj_3mo, idents = c("WT Resting"))
  972. ki.c2 <- WhichCells(sobj_3mo, idents = c("KI Cluster 2"))
  974. cells_wtki.c2 <- list("WT" = WT,
  975. "KI Cluster 2" = ki.c2)
  977. highlight_ki.c2 <- Cell_Highlight_Plot(seurat_object = sobj_3mo, cells_highlight = cells_wtki.c2, highlight_color = c("#228b22", "#ffb90f"))+
  978. theme(aspect.ratio = 1,
  979. legend.text = element_text(size = 15),
  980. axis.line = element_blank(),
  981. axis.title = element_blank(),
  982. axis.text = element_blank(),
  983. axis.ticks = element_blank(),
  984. legend.key.size = unit(0.5, "cm"),
  985. legend.title = element_text(color = "black", size = 18, face = 2),
  986. legend.position.inside = c(0.65 ,0.25))+
  987. guides(color = guide_legend(override.aes = list(size = 6)))
  989. highlight_ki.c2
  990. ```
  992. ```{r, fig.height=4, fig.width=8}
  993. #input dataframe for barplot of DEGs
  994. ki.c2.wt.degs.count <- data.frame(status = names(table(ki.c2.wt.dea_filt$Regulation)),
  995. total = as.vector(table(ki.c2.wt.dea_filt$Regulation)))
  996. #barplot
  997. ki.c2.wt.dea_bp <- ggplot(ki.c2.wt.degs.count, aes(x = status, y = total ,fill = status))+
  998. geom_bar(stat= "identity", width = 0.9)+
  999. scale_fill_manual(values = c("#00798c","#cc0000"))+
  1000. geom_text(aes(label = total), vjust= -0.2, size = 6, hjust = 1)+
  1001. coord_flip()+
  1002. theme_classic()+
  1003. theme(aspect.ratio = 0.5,
  1004. legend.position = 'none',
  1005. axis.text.x = element_text(size=20, angle=0),
  1006. axis.text.y = element_text(size=20, angle=0),
  1007. axis.title = element_text(size = 20),
  1008. axis.line = element_line(linewidth = 1.5))+
  1009. scale_x_discrete(labels = c("up" = "KI Cluster 2 vs WT Resting up",
  1010. "down" = "KI Cluster 2 vs WT Resting down"))+
  1011. ylab("number of DEGs")+
  1012. xlab(label= NULL)
  1013. #number of cells in the respective clusters
  1014. table(Idents(sobj_3mo_sub1))
  1015. ki.c2.wt.dea_bp
  1016. ```
  1018. #Figure
  1019. ## DEGs heatmap
  1021. ```{r, fig.height=7, fig.width=15}
  1022. # create vector with the top 20 degs
  1023. ki2wt_markers_to_plot <- dplyr::bind_rows(top20_up_degs_fc_ord_ki.c2.wt.dea, top20_dn_degs_fc_ord_ki.c2.wt.dea)
  1024. ki2wt_markers_to_plot_vec <- ki2wt_markers_to_plot$gene
  1025. ```
  1027. ```{r, fig.height= 10, fig.width=5}
  1028. Idents(sobj_3mo_sub1) <- factor(Idents(sobj_3mo_sub1),
  1029. levels = c("WT Resting", "KI Cluster 2"))
  1031. #heatmap of top 20 DEGs expressed between two time points (assessed with MAST)
  1032. [email hidden]$HM <- Idents(sobj_3mo_sub1)
  1033. wtki2_heat <- DoHeatmap(sobj_3mo_sub1, features = ki2wt_markers_to_plot_vec, size = 4,
  1034. angle = 45, group.by = "HM", draw.lines = T, raster = FALSE)+
  1035. scale_fill_gradientn(colors = c("blue", "white", "red"))
  1037. wtki2_heat
  1038. ```
  1042. #Figure
  1043. ## GO dotplot
  1045. ```{r, fig.width=9, fig.height=3.5}
  1046. #split in up and down to select non redundant terms more easily
  1047. up_ki.ki2 <- subset(ego_combined_ki.c2, status== "up")
  1048. dn_ki.ki2 <- subset(ego_combined_ki.c2, status== "down")
  1051. #selection of non redundant terms
  1053. selected_got_ki.c2vswt0 <- subset(ego_combined_ki.c2, Description %in% c("humoral immune response", "immune response-regulating cell surface receptor signaling pathway", "regulation of myeloid leukocyte differentiation", "regulation of chemotaxis", "acute inflammatory response", "regulation of neuron death", "antigen processing and presentation of peptide antigen", "MHC class II protein complex assembly", "peptide antigen assembly with MHC protein complex"))
  1055. #modify aesthetics
  1056. selected_got_ki.c2vswt0$Description <- factor(selected_got_ki.c2vswt0$Description, levels = unique(c("humoral immune response", "immune response-regulating cell surface receptor signaling pathway", "regulation of myeloid leukocyte differentiation", "regulation of chemotaxis", "acute inflammatory response", "regulation of neuron death", "antigen processing and presentation of peptide antigen", "MHC class II protein complex assembly", "peptide antigen assembly with MHC protein complex")))
  1057. #dotplot
  1058. ki.c2.wt_goe_sel <- ego_dp(selected_got_ki.c2vswt0)
  1059. #modify aesthetics
  1060. ki.c2.wt_goe_sel+
  1061. scale_color_gradient(low ="#6a329f", high = "#cd99e7")+
  1062. theme(axis.text.x = element_text(size = 14, face = "bold"),
  1063. axis.text.y = element_text(size = 14, face = "bold"))
  1064. ```
  1066. #Figure
  1068. ## DEGs between KI Cluster 1 vs WT resting
  1071. ```{r}
  1072. #Highlight plots
  1073. levels(Idents(sobj_3mo))
  1074. Idents(sobj_3mo) <- "cell_states"
  1076. WT <- WhichCells(sobj_3mo, idents = c("WT Resting"))
  1077. ki.c1 <- WhichCells(sobj_3mo, idents = c("KI Cluster 1"))
  1079. cells_wtki.c1 <- list("WT" = WT,
  1080. "KI Cluster 1" = ki.c1)
  1082. highlight_ki.c1 <- Cell_Highlight_Plot(seurat_object = sobj_3mo, cells_highlight = cells_wtki.c1, highlight_color = c("#228b22", "#eae4b7"))+
  1083. theme(aspect.ratio = 1,
  1084. legend.text = element_text(size = 15),
  1085. axis.line = element_blank(),
  1086. axis.title = element_blank(),
  1087. axis.text = element_blank(),
  1088. axis.ticks = element_blank(),
  1089. legend.key.size = unit(0.5, "cm"),
  1090. legend.title = element_text(color = "black", size = 18, face = 2),
  1091. legend.position.inside = c(0.65 ,0.25))+
  1092. guides(color = guide_legend(override.aes = list(size = 6)))
  1094. highlight_ki.c1
  1095. ```
  1098. barplot
  1100. ```{r}
  1101. #input dataframe for barplot of DEGs
  1102. ki.c1.wt.degs.count <- data.frame(status = names(table(ki.c1.wt.dea_filt$Regulation)),
  1103. total = as.vector(table(ki.c1.wt.dea_filt$Regulation)))
  1105. #barplot
  1106. ki.c1.wt.dea_bp <- ggplot(ki.c1.wt.degs.count, aes(x = status, y = total ,fill = status))+
  1107. geom_bar(stat= "identity", width = 0.9)+
  1108. scale_fill_manual(values = c("#00798c","#cc0000"))+
  1109. geom_text(aes(label = total), vjust= -0.2, size = 6, hjust = 1)+
  1110. coord_flip()+
  1111. theme_classic()+
  1112. theme(aspect.ratio = 0.5,
  1113. legend.position = 'none',
  1114. axis.text.x = element_text(size=20, angle=0),
  1115. axis.text.y = element_text(size=20, angle=0),
  1116. axis.title = element_text(size = 20),
  1117. axis.line = element_line(linewidth = 1.5))+
  1118. scale_x_discrete(labels = c("up" = "KI Cluster 1 vs WT Resting up",
  1119. "down" = "KI Cluster 1 vs WT Resting down"))+
  1120. ylab("number of DEGs")+
  1121. xlab(label= NULL)
  1123. #number of cells in the respective clusters
  1124. table(Idents(sobj_3mo_sub2))
  1125. ki.c1.wt.dea_bp
  1126. ```
  1128. #Figure
  1129. ## DEGs heatmap
  1132. ```{r, fig.height= 10, fig.width=5}
  1133. Idents(sobj_3mo_sub2) <- factor(Idents(sobj_3mo_sub2),
  1134. levels = c("WT Resting", "KI Cluster 1"))
  1136. # create vector with the top 20 degs
  1137. ki1wt_markers_to_plot <- dplyr::bind_rows(top20_up_degs_fc_ord_ki.c1.wt.dea, top20_dn_degs_fc_ord_ki.c1.wt.dea)
  1138. ki1wt_markers_to_plot_vec <- ki1wt_markers_to_plot$gene
  1140. #heatmap of top 20 DEGs expressed between two time points (assessed with MAST)
  1141. [email hidden]$HM <- Idents(sobj_3mo_sub2)
  1142. wtki1_heat <- DoHeatmap(sobj_3mo_sub2, features = ki1wt_markers_to_plot_vec, size = 4,
  1143. angle = 45, group.by = "HM", draw.lines = T, raster = FALSE)+
  1144. scale_fill_gradientn(colors = c("blue", "white", "red"))
  1146. wtki1_heat
  1147. ```
  1149. #Figure
  1150. ## GO dotplot
  1152. ```{r, fig.width=6, fig.height=3.5}
  1153. #split in up and down to select non redundant terms more easily
  1154. up_ki.c1 <- subset(ego_combined_ki.c1, status== "up")
  1155. dn_ki.c1 <- subset(ego_combined_ki.c1, status== "down")
  1157. #selection of non redundant terms
  1159. selected_got_ki.c1vswt0 <- subset(ego_combined_ki.c1, Description %in% c("cell junction disassembly", "positive regulation of endocytosis", "inflammatory response to antigenic stimulus", "mononuclear cell migration", "regulation of lipid metabolic process", "response to nutrient levels"))
  1161. #modify aesthetics
  1162. selected_got_ki.c1vswt0$Description <- factor(selected_got_ki.c1vswt0$Description, levels = unique(c("cell junction disassembly", "positive regulation of endocytosis", "inflammatory response to antigenic stimulus", "mononuclear cell migration", "regulation of lipid metabolic process", "response to nutrient levels")))
  1163. #dotplot
  1164. ki.c1.wt_goe_sel <- ego_dp(selected_got_ki.c1vswt0)
  1165. #modify aestethics
  1166. ki.c1.wt_goe_sel+
  1167. scale_color_gradient(low ="#6a329f", high = "#cd99e7")+
  1168. theme(axis.text.x = element_text(size = 14, face = "bold"),
  1169. axis.text.y = element_text(size = 14, face = "bold"))
  1170. ```
  1176. #_____________________________________________________________________________________________________
  1177. ## Expression of AD-risk genes
  1179. Expression levels of genes curated from GWAS studies aimed at identifying risk factors in AD. The list of genes was taken from Mancuso et al., 2024 supplementary table 4
  1181. Publication:
  1182. Mancuso, R., Fattorelli, N., Martinez-Muriana, A. et al. Xenografted human microglia display diverse transcriptomic states in response to Alzheimer’s disease-related amyloid-β pathology. Nat Neurosci 27, 886–900 (2024). https://doi.org/10.1038/s41593-024-01600-y
  1183. https://www.nature.com/articles/s41593-024-01600-y?s=31
  1185. ```{r}
  1186. #load gene list extrapolated from Mancuso et al., supplementary
  1187. #Mancuso_2024_AD_risk_gene_list (import .xlsx file)
  1189. #get list of genes
  1190. ad_risk_genes <- mancuso_ad_risk$genes
  1192. length(ad_risk_genes)
  1193. ```
  1195. ```{r, fig.width= 30, fig.height= 8}
  1196. #are any of these genes DEG in the KI clusters of this dataset?
  1197. ad_risk_deg_paquet <- intersect(ad_risk_genes, all_kivswt_union)
  1199. ad_risk_deg_paquet
  1201. #select genes
  1202. ad_risk_genes_deg <- c("ANKH", "BLNK", "CLU", "SORL1", "TREM2", "HLA-DQA1", "HLA-DRB1")
  1203. ```
  1205. ## Signature scoring
  1207. Zhou et al.,2020 signature
  1209. Publication:
  1210. Zhou, Y., Song, W.M., Andhey, P.S. et al. Human and mouse single-nucleus transcriptomics reveal TREM2-dependent and TREM2-independent cellular responses in Alzheimer’s disease. Nat Med 26, 131–142 (2020). https://doi.org/10.1038/s41591-019-0695-9
  1211. https://www.nature.com/articles/s41591-019-0695-9
  1213. Dataset description:
  1214. Human and mouse single-nucleus transcriptomics reveal TREM2-dependent and -independent cellular responses in Alzheimer’s disease
  1216. source of signature: supplementary table 4-Micro 0 signature
  1218. ```{r}
  1219. #import zhou et al supplementary from excel
  1220. #zhou_supp4 <-read_excel("/path/zhou_tablesup4_micro0.xlsx")
  1222. #extract gene list
  1223. zhou_sig <- list("gene.symbol" = zhou_supp4$gene)
  1224. #dimension of signature
  1225. length(zhou_sig[[1]])
  1226. ```
  1230. ```{r}
  1231. #set idents
  1232. Idents(sobj_3mo) <- "genotype"
  1233. #scoring of signature
  1234. sobj_3mo <- AddModuleScore_UCell(sobj_3mo, features = zhou_sig, name = "zhou1")
  1235. ```
  1238. #________________________________________________________________________________________________
  1240. # Figure
  1242. ## AD risk genes in KI cultures
  1244. ```{r, fig.height=8, fig.width=3}
  1245. #stacked violin
  1246. Idents(sobj_3mo) <- "cell_states"
  1248. sobj_3mo_sub4 <- subset(sobj_3mo, idents = c("WT Resting", "KI Cluster 1", "KI Cluster 2"))
  1250. [email hidden]$id_factor <- factor([email hidden]$cell_states, levels = c("WT Resting", "KI Cluster 1", "KI Cluster 2"))
  1252. Idents(sobj_3mo_sub4) <- "id_factor"
  1255. #violin plot
  1256. Stacked_VlnPlot(seurat_object = sobj_3mo_sub4, features = c("ANKH", "BLNK", "SORL1", "TREM2", "HLA-DQA1", "HLA-DRB1", "CLU"),
  1257. x_lab_rotate = TRUE, idents = c("WT Resting", "KI Cluster 1", "KI Cluster 2"), colors_use = c("#228b22","#eae4b7", "#ffb90f"))+
  1258. theme(aspect.ratio = 1.5)
  1261. ```
  1263. # Fig
  1265. ## Enrichment for Zhou et al, signature
  1267. ```{r}
  1268. #WT resting vs KI Inflammatory vs KI 2
  1270. #set idents
  1271. Idents(sobj_3mo_sub4) <- "cell_states"
  1272. #subset for clusters of interest
  1273. v.cl <- VlnPlot(sobj_3mo_sub4, assay="RNA",
  1274. features = "gene.symbolzhou1" ,
  1275. pt.size = 0, cols = c("#228b22", "#eae4b7", "#ffb90f")) +
  1276. geom_boxplot(width = 0.4)+
  1277. ylim(0.1, 0.45)+
  1278. scale_x_discrete(limits = c("WT Resting", "KI Cluster 1", "KI Cluster 2"))+
  1279. NoLegend() +
  1280. geom_signif(test = "wilcox.test",
  1281. map_signif_level = F,
  1282. comparisons = list(c("WT Resting", "KI Cluster 1"),
  1283. c("WT Resting", "KI Cluster 2")),
  1284. y_position = c(0.37, 0.4))+
  1285. theme(aspect.ratio = 0.9,
  1286. title = element_blank(),
  1287. axis.text.x = element_text(size=14, angle=45, hjust = 1),
  1288. axis.text.y = element_text(size=14, angle=0),
  1289. axis.title.x = element_text(face = "bold", size = 14),
  1290. axis.title.y = element_text(face = "bold", size = 14))+
  1291. xlab("")+
  1292. ylab("ES")+
  1293. ggtitle("AD microglia signature", subtitle = "Zhou et al., 2020")
  1297. v.cl
  1299. #stat
  1300. es_df_ad <- [email hidden][, c("cell_states", "gene.symbolzhou1")]
  1302. #stat
  1303. es_df_ad <- [email hidden][, c("cell_states", "gene.symbolzhou1")]
  1305. #ordering
  1306. es_df_ad$cell_states <- ordered(es_df_ad$cell_states,
  1307. levels = c("WT Resting", "KI Cluster 1", "KI Cluster 2"))
  1309. kruskal.test(gene.symbolzhou1 ~ cell_states, data = es_df_ad)
  1312. ```
  1316. ```{r}
  1318. #levels(Idents(sobj_3mo_sub4))
  1319. # sobj_3mo_sub4$cell_states_subset <- Idents(sobj_3mo_sub4)
  1320. #Idents(sobj_3mo_sub4) <- "cell_states_subset"
  1321. #scoring of signature
  1322. # sobj_3mo_sub4 <- AddModuleScore_UCell(sobj_3mo_sub4, features = zhou_sig, name = "zhou_sub")
  1324. es_df <- FetchData(sobj_3mo_sub4, vars = c("cell_states_subset", "gene.symbolzhou_sub"))
  1326. summary(es_df)
  1327. table(es_df$cell_states_subset, useNA = "ifany")
  1328. table(es_df$gene.symbolzhou_sub, useNA = "ifany")
  1329. kruskal.test(gene.symbolzhou_sub ~ cell_states_subset, data = es_df)
  1331. p <- ggplot(es_df, aes(x= cell_states_subset, y = gene.symbolzhou_sub))+
  1332. geom_violin(trim = FALSE)
  1334. p <- p + geom_signif(comparisons = list(c("WT Resting", "KI Cluster 1"),
  1335. c("WT Resting", "KI Cluster 2")),
  1336. map_signif_level = FALSE,
  1337. test = "wilcox.test")
  1339. p
  1343. rm(p)
  1345. p <- ggplot(es_df, aes(x= cell_states_subset, y = gene.symbolzhou_sub))+
  1346. geom_violin(trim = FALSE)
  1348. p <- p + geom_signif(map_signif_level = FALSE,
  1349. test = "kruskal.test")
  1350. p
  1351. ```
  1353. #______________________________________________________________________________________________________________________________
  1355. #Figure
  1356. #Overlap between updegs in KI Cluster 2 and Zhou signature
  1358. ```{r}
  1360. #paquet: up_kivswt_union
  1361. #zhou:
  1362. zhou_chr <- zhou_sig[[1]]
  1363. #background: total of tested genes in differential expression analysis
  1364. total.degs <- union(x = ki.c2.wt.dea$genes, y = ki.c1.wt.dea$genes)
  1365. total.degs <- unique(total.degs)
  1366. total.degs.num <- length(total.degs)
  1367. #create gene overlap object
  1368. go.objB <- newGeneOverlap(listA = up_kivswt_union, listB = zhou_chr, genome.size = total.degs.num)
  1369. #stat
  1370. go.objB <- testGeneOverlap(go.objB)
  1371. #show odds ratio
  1372. print(go.objB)
  1374. ```
  1377. ```{r}
  1378. gene.set <- list("KI 3mo" = up_kivswt_union, "Zhou et al., 2020 signature" = zhou_chr)
  1380. resPZ <- supertest(x= gene.set, n= total.degs.num)
  1382. plot.pz <- as.ggplot(~plot(resPZ, Layout="landscape", sort.by="size", keep=FALSE,
  1383. show.elements=TRUE,
  1384. elements.cex=0.7,
  1385. elements.list=subset(summary(resPZ)$Table,Observed.Overlap <= 20),
  1386. show.expected.overlap=F,
  1387. color.expected.overlap='red',
  1388. color.on = NULL))
  1390. plot.pz
  1392. summary(resPZ)$Table
  1394. ```

Figure5_scRNA_microglia_KI.Rmd at commit 31c0793, no license · at the source

Overview

Authors: Julien Klimmt1,2, Carolina Cardoso Gonçalves1,2, Jessica Valentina Montgomery3, Stephan A. Müller4,5, Merle Bublitz1,2, Severin Filser1, Lars Paeger4,6, Brigitte Nuscher7, Angelika Dannert1,2, Sigrun Roeber6, Veronica Pravata8, Martina Schifferer4,9, Joshua J. Shrouder1, Nathalie Schulz1, Judit González-Gallego1,2, Silvia Cappello8,9,10, Thomas Misgeld4,9,11, Nikolaus Plesnila1,9, Eduardo Beltrán9,12,13, Jochen Herms4,6,9
and 7 other authorsElena De Domenico3,14, Marc D. Beyer3,14,15, Joachim L. Schultze3,14,16, Christian Haass4,7,9, Stefan F. Lichtenthaler4,5,9, Caterina Carraro3, Dominik Paquet1,9
ORCID iDs: Julien Klimmt, Carolina Cardoso Gonçalves, Stephan A. Müller, Angelika Dannert, Veronica Pravata, Martina Schifferer, Joshua J. Shrouder, Nathalie Schulz, Nikolaus Plesnila, Eduardo Beltrán, Jochen Herms, Elena De Domenico, Marc D. Beyer, Joachim L. Schultze, Christian Haass, Stefan F. Lichtenthaler, Caterina Carraro, Dominik Paquet
16 affiliations
  1. Institute for Stroke and Dementia Research (ISD), LMU University Hospital, LMU Medizin, LMU Munich,Munich, Germany
  2. Graduate School of Systemic Neuroscience (GSN), LMU Munich,Munich, Germany
  3. Systems Medicine, Deutsches Zentrum für Neurodegenerative Erkrankungen (DZNE) e.V.,Bonn, Germany
  4. German Center for Neurodegenerative Diseases (DZNE) Munich,Munich, Germany
  5. Neuroproteomics, School of Medicine and Health, TUM University Hospital, Technical University of Munich,Munich, Germany
  6. Institute of Neuropathology, LMU Medizin, LMU Munich,Munich, Germany
  7. Metabolic Biochemistry, Biomedical Center (BMC), Faculty of Medicine, LMU Munich,Munich, Germany
  8. Chair of Physiological Genomics, Biomedical Center (BMC), LMU Medizin, LMU Munich,Munich, Germany
  9. Munich Cluster for Systems Neurology (SyNergy),Munich, Germany
  10. Max-Planck-Institute of Psychiatry,Munich, Germany
  11. Institute of Neuronal Cell Biology, Technical University of Munich,Munich, Germany
  12. Institute of Clinical Neuroimmunology, LMU University Hospital, LMU Medizin, LMU Munich,Munich, Germany
  13. Biomedical Center (BMC), LMU Medizin, LMU Munich,Munich, Germany
  14. PRECISE Platform for Genomics and Epigenomics, Deutsches Zentrum für Neurodegenerative Erkrankungen (DZNE) e.V. and University of Bonn and West German Genome Center,Bonn, Germany
  15. Immunogenomics & Neurodegeneration, Deutsches Zentrum für Neurodegenerative Erkrankungen (DZNE) e.V.,Bonn, Germany
  16. Genomics and Immunoregulation, Life & Medical Sciences (LIMES) Institute, University of Bonn,Bonn, Germany
Institutions: Institut für Schlaganfall- und Demenzforschung (Germany); Ludwig-Maximilians-Universität München (Germany); German Center for Neurodegenerative Diseases (Germany); Technical University of Munich (Germany); Institut für Neuropathologie (Germany); Munich Cluster for Systems Neurology (Germany); Max Planck Institute of Psychiatry (Germany); Institut für Klinische Neuroimmunologie (Germany); University of Bonn (Germany)
Journal: Nature neuroscience, volume 29, issue 9, pages 2324-2340
Dates: received 17 October 2024; accepted 10 June 2026; published online 4 August 2026; in print 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41593-026-02367-0 · PMID 42552384 · PMCID PMC13533838 · OpenAlex W7172420001
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), Alzheimer's / dementia (population), cellular / molecular (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning, Connectivity, Single-unit activity, calcium imaging
Keywords: Neuroimmunology, Induced pluripotent stem cells, Alzheimer's disease, Neurodegeneration, Experimental models of disease
MeSH: Alzheimer Disease*, Brain*, Induced Pluripotent Stem Cells*, Microglia*, Amyloid beta-Peptides, Astrocytes, Humans, Neurons, Phenotype (* major topic)
Journal subjects: Technical Report
Topic: Neuroinflammation and Neurodegeneration Mechanisms (Neurology, Neuroscience), according to OpenAlex
Funding: BrightFocus Foundation (BrightFocus) (ADR AD2019604S); Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) (16LW047)
Citations: not cited yet (Europe PMC); 106 references in the paper
Research resources: male iPSC line GM05399 RRID:CVCL_7425, 3) KOLF2.1J cell line RRID:CVCL_B5P3, 409B2 line RRID:CVCL_K092, female iPSC line A18944 RRID:CVCL_RM92

Abstract

Stem-cell-based in vitro models offer promising potential to elucidate human brain cell functions and interactions, but limitations in reproducibility, maturation and cell-type diversity persist. Especially, prolonged incorporation of mature microglia and studies of neuroinflammation have proven challenging. Here, we developed a human induced pluripotent stem cell-based three-dimensional cortical brain tissue model (3BTM) containing neurons, astrocytes and microglia with high reproducibility, maturity and viability. 3BTMs show morphological, functional and proteomic maturation of all cell types, leading to high similarity to their in vivo counterparts. Incorporated microglia survive for over 6 months and display mature morphology, functions and gene expression. Importantly, when engineered to model Alzheimer’s disease pathology, 3BTMs recapitulate key disease hallmarks, including amyloid deposition, increased phospho-tau levels and neuroinflammation, with microglia shifting their transcriptional landscape to disease-relevant signatures. Treatment of Alzheimer’s disease 3BTMs with anti-Aβ immunotherapy cleared deposits and largely reversed disease signatures in glia. Together, our microglia-containing model provides a platform for studying physiological and pathological states of human brain tissue.

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 18 matches between paragraphs and lines of code.

jsschrepping/r_docker

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 187c49cdfa925ebd14226ede042ddc0122cf6215, 29 April 2024
Languages: R (10), Shell (1)
Size: 20 files, 11 scripts
Software Heritage: archived
Found in: the text, “Sample preparation for BD Rhapsody datasets”
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Monocle 3 (2 files), circlize (1 file), clusterProfiler (1 file), ComplexHeatmap (1 file), data.table (1 file), DESeq2 (1 file), edgeR (1 file), ggpubr (1 file), Harmony (1 file), igraph (1 file), limma (1 file), multcomp (1 file), patchwork (1 file), pheatmap (1 file), reshape2 (1 file), reticulate (1 file), Seurat (1 file), SingleCellExperiment (1 file), Subread (featureCounts) (1 file), tidyverse (1 file), WGCNA (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
13 files
At the source: github.com/jsschrepping/r_docker

gitlab.dzne.de/ag-beyer/singlecell_analysis_3btm

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 31c079315b8210b043d0b6d72b92ad5a18e0d168, 8 May 2026
Languages: R (6)
Size: 10 files, 6 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, 6 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: clusterProfiler (6 files), ComplexHeatmap (6 files), cowplot (6 files), ggpubr (6 files), Harmony (6 files), patchwork (6 files), reshape2 (6 files), Seurat (6 files), tidyverse (6 files), data.table (5 files), ggplot2 (5 files), reticulate (5 files), DESeq2 (4 files), edgeR (4 files), circlize (3 files), Monocle 3 (2 files), igraph (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
7 files
At the source: gitlab.dzne.de/ag-beyer/singlecell_analysis_3btm

Code availability

All original code to analyze the scRNA-seq data, and the tool and package information to fully reproduce the analysis is publicly available via GitLab at https://gitlab.dzne.de/ag-beyer/singlecell_analysis_3btm. Any additional information required to reanalyze the data reported in this article is available from the lead contact upon reasonable request.

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;
  • 17 scripts, each with its path and the digest of its content;
  • 18 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

Data availability

The scRNA-seq profiling datasets have been deposited with the European Genome-phenome Archive (EGA) (10X dataset: EGAS50000000469 (https://ega-archive.org/studies/EGAS50000000469); Rhapsody dataset 2D versus 3D: EGAS50000001392 (https://ega-archive.org/studies/EGAS50000001392); Rhapsody dataset SA2 + A18 WT versus knock-in and sAβ42 ± aducanumab treatment: EGAS50000001397 (https://ega-archive.org/studies/EGAS50000001397)). Access to these datasets will be granted upon request to the Data Access Committee to protect the rights of human donors. The MS proteomic datasets have been deposited with the ProteomeXchange Consortium via the PRIDE partner repository with dataset ID PXD055254 (https://www.ebi.ac.uk/pride/archive/projects/PXD055254) (dataset 1–3BTM maturation in SA2 line in comparison to NOs and human brain samples), PXD069984 (https://www.ebi.ac.uk/pride/archive/projects/PXD069984) (dataset 2–reproducibility of 3BTMs across iPSC lines and comparison to 2D co-cultures) and PXD071849 (https://www.ebi.ac.uk/pride/archive/projects/PXD071849) (dataset 3–comparison of WT and knock-in 3BTMs). The original microscopy data reported in this article will be shared by the lead contact upon request because of large file sizes. Source data are provided with this paper.

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

Versions

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

Version 3, 28 September 2026

  • Publisher: n/a → Nature Portfolio

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 27 authors, 5 keywords, 9 MeSH terms, 2 funders, 104 references, 4 RRIDs.

Cite

This paper

Klimmt, J., Cardoso Gonçalves, C., Montgomery, J. V., Müller, S. A., Bublitz, M., Filser, S., Paeger, L., Nuscher, B., Dannert, A., Roeber, S., Pravata, V., Schifferer, M., Shrouder, J. J., Schulz, N., González-Gallego, J., Cappello, S., Misgeld, T., Plesnila, N., Beltrán, E., . . . Paquet, D. (2026). A reproducible three-dimensional model of human brain tissue to investigate physiological and disease-associated microglia phenotypes. Nature neuroscience, 29(9), 2324-2340. https://doi.org/10.1038/s41593-026-02367-0

BibTeX

@article{klimmt2026reproducible,
author = {Klimmt, Julien and Cardoso Gonçalves, Carolina and Montgomery, Jessica Valentina and Müller, Stephan A. and Bublitz, Merle and Filser, Severin and Paeger, Lars and Nuscher, Brigitte and Dannert, Angelika and Roeber, Sigrun and Pravata, Veronica and Schifferer, Martina and Shrouder, Joshua J. and Schulz, Nathalie and González-Gallego, Judit and Cappello, Silvia and Misgeld, Thomas and Plesnila, Nikolaus and Beltrán, Eduardo and Herms, Jochen and De Domenico, Elena and Beyer, Marc D. and Schultze, Joachim L. and Haass, Christian and Lichtenthaler, Stefan F. and Carraro, Caterina and Paquet, Dominik},
title = {{A reproducible three-dimensional model of human brain tissue to investigate physiological and disease-associated microglia phenotypes}},
journal = {Nature neuroscience},
year = {2026},
month = aug,
volume = {29},
number = {9},
pages = {2324--2340},
publisher = {Nature Portfolio},
issn = {1097-6256},
doi = {10.1038/s41593-026-02367-0},
url = {https://doi.org/10.1038/s41593-026-02367-0},
pmid = {42552384},
pmcid = {PMC13533838}
}

RIS

TY - JOUR
AU - Klimmt, Julien
AU - Cardoso Gonçalves, Carolina
AU - Montgomery, Jessica Valentina
AU - Müller, Stephan A.
AU - Bublitz, Merle
AU - Filser, Severin
AU - Paeger, Lars
AU - Nuscher, Brigitte
AU - Dannert, Angelika
AU - Roeber, Sigrun
AU - Pravata, Veronica
AU - Schifferer, Martina
AU - Shrouder, Joshua J.
AU - Schulz, Nathalie
AU - González-Gallego, Judit
AU - Cappello, Silvia
AU - Misgeld, Thomas
AU - Plesnila, Nikolaus
AU - Beltrán, Eduardo
AU - Herms, Jochen
AU - De Domenico, Elena
AU - Beyer, Marc D.
AU - Schultze, Joachim L.
AU - Haass, Christian
AU - Lichtenthaler, Stefan F.
AU - Carraro, Caterina
AU - Paquet, Dominik
TI - A reproducible three-dimensional model of human brain tissue to investigate physiological and disease-associated microglia phenotypes
T2 - Nature neuroscience
J2 - Nat Neurosci
PY - 2026
DA - 2026/08/04
VL - 29
IS - 9
SP - 2324
EP - 2340
SN - 1097-6256
PB - Nature Portfolio
DO - 10.1038/s41593-026-02367-0
UR - https://doi.org/10.1038/s41593-026-02367-0
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41593-026-02367-0",
"type": "article-journal",
"title": "A reproducible three-dimensional model of human brain tissue to investigate physiological and disease-associated microglia phenotypes",
"container-title": "Nature neuroscience",
"author": [
{
"family": "Klimmt",
"given": "Julien"
},
{
"family": "Cardoso Gonçalves",
"given": "Carolina"
},
{
"family": "Montgomery",
"given": "Jessica Valentina"
},
{
"family": "Müller",
"given": "Stephan A."
},
{
"family": "Bublitz",
"given": "Merle"
},
{
"family": "Filser",
"given": "Severin"
},
{
"family": "Paeger",
"given": "Lars"
},
{
"family": "Nuscher",
"given": "Brigitte"
},
{
"family": "Dannert",
"given": "Angelika"
},
{
"family": "Roeber",
"given": "Sigrun"
},
{
"family": "Pravata",
"given": "Veronica"
},
{
"family": "Schifferer",
"given": "Martina"
},
{
"family": "Shrouder",
"given": "Joshua J."
},
{
"family": "Schulz",
"given": "Nathalie"
},
{
"family": "González-Gallego",
"given": "Judit"
},
{
"family": "Cappello",
"given": "Silvia"
},
{
"family": "Misgeld",
"given": "Thomas"
},
{
"family": "Plesnila",
"given": "Nikolaus"
},
{
"family": "Beltrán",
"given": "Eduardo"
},
{
"family": "Herms",
"given": "Jochen"
},
{
"family": "De Domenico",
"given": "Elena"
},
{
"family": "Beyer",
"given": "Marc D."
},
{
"family": "Schultze",
"given": "Joachim L."
},
{
"family": "Haass",
"given": "Christian"
},
{
"family": "Lichtenthaler",
"given": "Stefan F."
},
{
"family": "Carraro",
"given": "Caterina"
},
{
"family": "Paquet",
"given": "Dominik"
}
],
"container-title-short": "Nat Neurosci",
"volume": "29",
"issue": "9",
"page": "2324-2340",
"DOI": "10.1038/s41593-026-02367-0",
"PMID": "42552384",
"PMCID": "PMC13533838",
"ISSN": "1097-6256",
"publisher": "Nature Portfolio",
"URL": "https://doi.org/10.1038/s41593-026-02367-0",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
4
]
]
}
}

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: Monocle 3, WGCNA, Harmony, 17 other tools, cellular / molecular, 6 references
[2] doi:10.1038/s41586-026-10214-2 [code]
Multidimensional profiling of heterogeneity in supratentorial ependymomas.
Journal: Nature
In common: Monocle 3, Harmony, SingleCellExperiment, 16 other tools, 2 references
[3] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: Subread (featureCounts), Monocle 3, WGCNA, 16 other tools, cellular / molecular
[4] doi:10.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: multcomp, Subread (featureCounts), WGCNA, 16 other tools
[5] doi:10.1002/imt2.70163 [code]
Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: Monocle 3, Harmony, SingleCellExperiment, 16 other tools, cellular / molecular
[6] 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, edgeR, 13 other tools, cellular / molecular, 3 references
[7] doi:10.1016/j.isci.2026.115573 [code]
Female cortical cellular mosaicism underlies shared MeCP2 and PCB impacted gene pathways.
Journal: iScience
In common: WGCNA, SingleCellExperiment, edgeR, 13 other tools, cellular / molecular, 2 references
[8] 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, edgeR, 15 other tools
[9] doi:10.1016/j.xcrm.2026.102682 [code]
TET CpG sequence-context-specific DNA demethylation shapes progression of IDH-mutant gliomas.
Journal: Cell reports. Medicine
In common: multcomp, edgeR, limma, 13 other tools, 2 references
[10] doi:10.1016/j.celrep.2026.117500 [code]
Spatio-molecular gene expression reflects dorsal anterior cingulate cortex structure and function in the human brain.
Journal: Cell reports
In common: Harmony, SingleCellExperiment, edgeR, 12 other tools, cellular / molecular, 4 references

Contribute

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

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

Request its removal

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

Discussion, reproductions, activity

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

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

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