OSCR

Modelling synaptic dysfunction in childhood dementia using human iPSC-derived cortical networks.

Code ↔ Paper

13 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 13 matches
  1. [1] § Methods › Electrophysiological feature analysis › Gene set enrichment analysis (GSEA) ↔ Rscripts/Differential_Expression_Analysis_120D.Rmd, lines 694–780 · score 0.96 · simplifyGO, GO IDs, Cellular metabolism, avglog2FC, common words, downregulated genes
  2. [2] § Methods › Electrophysiological feature analysis › Single nuclei data processing ↔ Rscripts/QualityControl_120Day.Rmd, lines 255–289 · score 0.95 · generate cell clusters, Principal component, FindClusters, FindNeighbors, RunUMAP, SCTransform
  3. [3] § Methods › Electrophysiological feature analysis › Differential expression analysis ↔ Rscripts/Differential_Expression_Analysis_120D.Rmd, lines 782–901 · score 0.92 · FindMarkers, min.diff.pct, logfc.threshold, fold change, log normalised, excitatory neurons
  4. [4] § Methods › Electrophysiological feature analysis › Gene set enrichment analysis (GSEA) ↔ Rscripts/Differential_Expression_Analysis_120D.Rmd, lines 607–652 · score 0.90 · minSize, nPermSimple, gseaParam, maxSize, synaptic genes, inhibitory neurons
  5. [5] § Methods › Electrophysiological feature analysis › scMayoMap cell annotations ↔ Rscripts/Neuronal_Glial_Types_Analysis_120D.Rmd, lines 325–401 · score 0.85 · FindConservedMarkers, cluster_type2, scMayoMap, cluster_type1, cell annotations, database
  6. [6] § Methods › Electrophysiological feature analysis › Integration with neurotypical, human cortical biopsies ↔ Rscripts/QualityControl_120Day.Rmd, lines 255–289 · score 0.80 · SCT assay, variable features, RunPCA, Principal component, SCTransform, days
  7. [7] § Results › MPS IIIA neurons reveal dysregulated gene expression implicated in synaptic homeostasis ↔ Rscripts/Differential_Expression_Analysis_120D.Rmd, lines 694–780 · score 0.71 · Cellular metabolism, Ions, transport, inhibitory neurons, organising, upregulated
  8. [8] § Methods › Electrophysiological feature analysis › UMAP and Seurat cluster-based cell classification ↔ Rscripts/Neuronal_Glial_Types_Analysis_120D.Rmd, lines 572–609 · score 0.71 · FindClusters, FindNeighbors, RunUMAP, Seurat cluster, harmony, cells
  9. [9] § Methods › Electrophysiological feature analysis › Classification of cycling cells ↔ Rscripts/Neuronal_Glial_Types_Analysis_120D.Rmd, lines 211–278 · score 0.71 · cell cycle scores, NonCycling, AddModuleScore, G2M, G1S, gene
  10. [10] § Results › MPS IIIA neurons reveal dysregulated gene expression implicated in synaptic homeostasis ↔ Rscripts/Differential_Expression_Analysis_120D.Rmd, lines 506–563 · score 0.61 · expressed genes, NEUROD6, SATB2, excitatory neurons, inhibitory neurons, overlap
  11. [11] § Methods › Electrophysiological feature analysis › Calculation of cell type scores ↔ Rscripts/Neuronal_Glial_Types_Analysis_120D.Rmd, lines 150–209 · score 0.54 · AddModuleScore, Immature Neurons, Nowakowski, scores, OPCs, Astrocytes
  12. [12] § Methods › Electrophysiological feature analysis › Final cell type assignments ↔ Rscripts/Neuronal_Glial_Types_Analysis_120D.Rmd, lines 325–401 · score 0.51 · scMayoMap, cluster_type1, Seurat cluster, Inhibitory, Excitatory, cell
  13. [13] § Methods › Electrophysiological feature analysis › Unsupervised clustering analysis ↔ Rscripts/QualityControl_120Day.Rmd, lines 291–323 · score 0.50 · principal components, elbow, variance, PCA, clustering, filtered

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 · 929 lines · 54 KB · no license · 5 matches

  1. ---
  2. title: "{[Sanfilippo] Differential expression analysis - 120 Day"
  3. author: "Inushi De Silva & Paris Mazzachi"
  4. date: "2023-11-16"
  5. output: html_document
  6. editor_options:
  7. chunk_output_type: console
  8. chunk_output_type: console
  9. ---
  11. ```{r setup, include=FALSE}
  12. knitr::opts_chunk$set(echo = TRUE)
  13. ```
  15. 1. ***First we will load the required libraries and set workingdirectories, datadirectories and plotdirectories.***
  16. ```{r Load Libraries and Set Directory Paths}
  17. options(future.globals.maxSize=1000* 1024^3)
  19. script_path <- rstudioapi::getActiveDocumentContext()$path
  20. script_dir <- dirname(script_path)
  21. workingdirectory <- dirname(script_dir)
  22. setwd(workingdirectory)
  23. datadirectory <- paste0(workingdirectory,"/data")
  24. resultsdirectory <- paste0(workingdirectory, "/results")
  25. RData_directory <- paste0(datadirectory,"/RData")
  26. ##dir.create(resultsdirectory)
  28. ## Load required libraries
  29. if(!require(dplyr)){
  30. install.packages("dplyr")
  31. library(dplyr)
  32. }
  33. if(!require(Seurat)){
  34. install.packages("Seurat")
  35. library(Seurat)
  36. }
  37. if(!require(ggplot2)){
  38. install.packages("ggplot2")
  39. library(ggplot2)
  40. }
  41. if(!require(sctransform)){
  42. install.packages("sctransform")
  43. library(sctransform)
  44. }
  45. if(!require(RColorBrewer)){
  46. install.packages("RColorBrewer")
  47. library(RColorBrewer)
  48. }
  49. if(!require(ggpointdensity)){
  50. install.packages("ggpointdensity")
  51. library(ggpointdensity)
  52. }
  53. if(!require(devtools)){
  54. install.packages("devtools")
  55. library(devtools)
  56. }
  57. if(!require(BiocManager)){
  58. install.packages("BiocManager")
  59. library(BiocManager)
  60. }
  61. if(!require(patchwork)){
  62. install.packages("patchwork")
  63. library(patchwork)
  64. }
  65. ```
  67. ```{r Other functions required in this markdown}
  68. ## A function to calculate differentially and significantly expressed geens
  69. significantly_expressed_genes <- function(differential_genes_df){
  70. differential_genes_df <- differential_genes_df[order(-differential_genes_df$avg_log2FC, differential_genes_df$p_val_adj),]
  71. differential_genes_df$p_val_adj[differential_genes_df$p_val_adj == 0] <- 1e-304
  72. differential_genes_df$log10p_adj <- -log10(differential_genes_df$p_val_adj)
  73. differential_genes_df$rank <- (differential_genes_df$avg_log2FC)*(differential_genes_df$log10p_adj)
  74. differential_genes_df <- differential_genes_df[order(differential_genes_df$rank, decreasing = T),]
  75. top <- differential_genes_df %>% top_n(15, rank)
  76. bottom <- differential_genes_df %>% top_n(-15, rank)
  77. label_genes <- c(rownames(top), rownames(bottom))
  79. differential_genes_df$diffexpressed <- ifelse(differential_genes_df$p_val_adj < 0.05 & differential_genes_df$avg_log2FC > 0, "TOP UP", ifelse(differential_genes_df$p_val_ad < 0.05 & differential_genes_df$avg_log2FC < 0,"TOP DOWN", "NO"))
  81. differential_genes_df$delabel <- NA
  82. differential_genes_df$delabel[(rownames(differential_genes_df) %in% label_genes)] <- rownames(differential_genes_df)[(rownames(differential_genes_df) %in% label_genes)]
  84. differential_genes_df$diffexpressed[is.na(differential_genes_df$delabel)] <- "NO"
  86. df_name <- deparse(substitute(differential_genes_df))
  88. return(differential_genes_df)
  90. }
  92. ## Adjusted EnrichmentPlot
  93. EnrichmentPlot <- function (pathway, stats, gseaParam = 1, ticksSize = 0.2)
  94. {
  95. pd <- plotEnrichmentData(pathway = pathway, stats = stats,
  96. gseaParam = gseaParam)
  98. with(pd, ggplot(data = curve) + geom_line(aes(x = rank, y = ES),
  99. color = "deepskyblue") + geom_segment(data = ticks, mapping = aes(x = rank,
  100. y = -spreadES/16, xend = rank, yend = spreadES/16), linewidth = ticksSize) +
  101. #geom_hline(yintercept = posES, colour = "red", linetype = "dashed") +
  102. #geom_hline(yintercept = negES, colour = "red", linetype = "dashed") +
  103. geom_hline(yintercept = 0, colour = "black") + theme(panel.background = element_blank(),
  104. panel.grid.major = element_line(color = "grey92"), axis.text.x = element_text(size = 18, angle = 90, vjust = 0.5), axis.text.y = element_text(size = 18, hjust = 0.5), axis.title.x = element_text(size = 18), axis.title.y = element_text(size = 18)) +
  105. labs(x = "Gene rank", y = "Enrichment score")) # + coord_cartesian(ylim = c(-0.025,0.45))
  106. }
  108. ## Manually calculate percentage of cells
  110. PrctCellExpringGene <- function(object, genes, group.by = "all"){
  111. if(group.by == "all"){
  112. prct = unlist(lapply(genes,calc_helper, object=object))
  113. result = data.frame(Markers = genes, Cell_proportion = prct)
  114. return(result)
  115. }
  117. else{
  118. list = SplitObject(object, group.by)
  119. factors = names(list)
  121. results = lapply(list, PrctCellExpringGene, genes=genes)
  122. for(i in 1:length(factors)){
  123. results[[i]]$Feature = factors[i]
  124. }
  125. combined = do.call("rbind", results)
  126. return(combined)
  127. }
  128. }
  130. calc_helper <- function(object,genes){
  131. counts = object[['RNA']]@data
  132. ncells = ncol(counts)
  133. if(genes %in% row.names(counts)){
  134. sum(counts[genes,]>0)/ncells
  135. }else{return(NA)}
  136. }
  138. ## function detect the most common words per cluster
  139. most_common_word=function(x){
  141. #Split every word into single words for counting
  142. splitTest=unlist(strsplit(x," "))
  144. #Counting words
  145. count=table(splitTest)
  147. #Sorting to select only the highest value, which is the first one
  148. count=count[order(count, decreasing=TRUE)][1:50]
  150. #Return the desired character.
  151. #By changing this you can choose whether it show the number of times a word repeats
  152. return(names(count))
  153. }
  155. # from https://hbctraining.github.io/scRNA-seq/lessons/elbow_plot_metric.html
  156. significant_pcs <- function(x){
  157. stdv <- x@reductions$pca@stdev # stdev per PC
  158. sum.stdv <- sum(x@reductions$pca@stdev) # Sum of stdev
  159. percent.stdv <- (stdv/sum.stdv)*100 # percentage of variation associated with each PC
  160. cumulative <- cumsum(percent.stdv) # cumulative percent of each PC
  161. # Determine which PC exhibits cumulative percent greater than 90% and % variation associated with the PC as less than 5
  162. co1 <- which(cumulative > 90 & percent.stdv < 5)[1]
  163. # Determine the difference between variation of PC and subsequent PC
  164. co2 <- sort(which((percent.stdv[1:length(percent.stdv) - 1] - percent.stdv[2:length(percent.stdv)]) > 0.1),
  165. decreasing = T)[1] + 1
  166. # Minimum of the two calculations used for clustering
  167. min.pc <- min(co1, co2)
  169. return(min.pc)
  170. }
  172. ```
  174. 2. ***Load the required RData files that have previously been QC'd***
  175. ```{r Load RData files that have previously been QC'd}
  176. ## Load required datasets
  177. set.seed(1)
  178. file_names <- list.files(paste0(RData_directory,"/120D/revisions/"), pattern = "*RData",full.names = T)
  179. lapply(file_names[c(6, 8)], load, .GlobalEnv)
  180. rm(file_names)
  181. ```
  183. 3. ***We now load in genes from the Nowakowski and Liddelow datasets***
  184. ```{r Load in neuronal gene lists and create colour scheme}
  185. ##load updated neuronal gene lists
  186. ##The gene names have been updated previously.
  187. load(file = paste0(datadirectory, "/reference_sheets/updated_neuronal_genes_ls.RData"))
  189. library(readxl)
  190. synaptic_genes <- read_excel(paste0(datadirectory, "/reference_sheets/syngo_genes.xlsx")) # SynGO Portal v1.2 (2023-12-01)
  191. stress_genes <- read_excel(paste0(datadirectory, "/reference_sheets/response_to_stress_genes.xlsx"))
  192. neuronal_genes_ls <- updated_neuronal_genes_ls
  193. rm(updated_neuronal_genes_ls)
  195. ## We also create a colour palette for each of the cell types
  196. library(colorspace)
  197. colour_palette <- data.frame(matrix(nrow = 44, ncol = 2))
  198. colour_palette$X1 <- c(names(neuronal_genes_ls)[1:43], "Other")
  199. colour_palette$X2[grepl("Astrocyte", colour_palette$X1, ignore.case = T)] <- colorRampPalette(c("salmon", "red"))(3)
  200. colour_palette$X2[grepl("eN", colour_palette$X1, ignore.case = T)] <- colorRampPalette(c("chocolate4", "tan"))(18)
  201. colour_palette$X2[grepl("RG", colour_palette$X1, ignore.case = T)] <- colorRampPalette(c("limegreen", "palegreen"))(11)
  202. colour_palette$X2[grepl("MGE_newborn", colour_palette$X1, ignore.case = T)] <- colorRampPalette(c("plum", "plum4"))(5)
  203. colour_palette$X2[grepl("MGE_Ctx", colour_palette$X1, ignore.case = T)] <- colorRampPalette(c("purple", "purple3"))(2)
  204. colour_palette$X2[grepl("CGE_LGE", colour_palette$X1, ignore.case = T)] <- colorRampPalette(c("slateblue", "slateblue4"))(2)
  205. colour_palette$X2[grepl("OPC", colour_palette$X1, ignore.case = T)] <- colorRampPalette(c("royalblue", "royalblue"))(1)
  206. colour_palette$X2[grepl("Striatal", colour_palette$X1, ignore.case = T)] <- colorRampPalette(c("sienna", "sienna"))(1)
  207. colour_palette$X2[grepl("Microglia", colour_palette$X1, ignore.case = T)] <- colorRampPalette(c("aquamarine3", "aquamarine3"))(1)
  208. colour_palette$X2[grepl("Mural", colour_palette$X1, ignore.case = T)] <- colorRampPalette(c("khaki", "khaki"))(1)
  209. colour_palette$X2[grepl("Other", colour_palette$X1, ignore.case = T)] <- colorRampPalette(c("black", "black"))(1)
  211. ##convert to named vector
  212. colours <- as.vector(colour_palette$X2)
  213. names(colours) <- colour_palette$X1
  214. ```
  216. 4. ***We now perform a differential expression analysis between control and MPSIIIA***
  217. ```{r Differential Expression Analysis and Volcano Plots - ALL}
  218. ## To perform differential expression analysis we first merge the two Seurat objects together.
  219. ## We will use the datasets which do not contain any "Unassigned" values for the donor id
  220. Sanfilippo_120D <- merge(Control_filtered_ambientRNA_120D, MPSIIIA_filtered_ambientRNA_120D) # remove 2 unassigned from MPSIIIA seurat
  221. table(Sanfilippo_120D$orig.ident)
  222. tbl <- [email hidden] %>% group_by(orig.ident, BIAD_id, final_cell_type) %>% summarise(n = n())
  224. cells_to_remove <- c("unassigned", "doublet") # remove doublets and unassigned cells
  226. Sanfilippo_120D <- subset(Sanfilippo_120D, cells = which(Sanfilippo_120D$BIAD_id %in% cells_to_remove), invert = T) # need to removed 2 unassigned from sanfil + striatal neurons
  227. tbl <- [email hidden] %>% group_by(orig.ident, BIAD_id, final_cell_type) %>% summarise(n = n())
  228. [email hidden] %>% group_by(orig.ident) %>% summarise(n = n())
  229. Idents(Sanfilippo_120D) <- 'orig.ident'
  230. Idents(Sanfilippo_120D) <- 'final_cell_type'
  232. Sanfilippo_120D <- subset(Sanfilippo_120D, downsample = 6619) # Downsample to have same pool of cells in each genotype
  233. [email hidden] %>% group_by(orig.ident) %>% summarise(n = n())
  235. ## We now SCTransform, RunPCA and UMAP to normalise gene expressions.
  236. Sanfilippo_120D <- JoinLayers(Sanfilippo_120D, assay = 'RNA')
  237. Sanfilippo_120D <- SCTransform(Sanfilippo_120D, vars.to.regress = c("percent.mt", "BIAD_id"))
  239. ## Save the seurat object to save time in future without re-running SCTransform each time
  240. save(Sanfilippo_120D, file = paste0(datadirectory, "/RData/120D/revisions/Combined_Control_MPSIIIA_120D_NOambientRNA.RData"))
  241. load(file = paste0(datadirectory, "/RData/120D/revisions/Combined_Control_MPSIIIA_120D_NOambientRNA.RData")) # Once run SCTransform once, can load the save RData file
  243. ## Run PCA and clustering for visualisation of cell types
  244. Sanfilippo_120D <- RunPCA(Sanfilippo_120D, npcs = 100, assay = "SCT")
  245. Sanfilippo_120D <- NormalizeData(Sanfilippo_120D, assay = "RNA", normalization.method = "LogNormalize")
  247. min.pc <- significant_pcs(Sanfilippo_120D)
  248. Sanfilippo_120D <- FindNeighbors(Sanfilippo_120D, dims = 1:min.pc, reduction = "pca")
  249. Sanfilippo_120D <- FindClusters(Sanfilippo_120D, resolution = 0.1)
  250. Sanfilippo_120D <- RunUMAP(Sanfilippo_120D, reduction = "pca", dims= 1:min.pc, n.components = 2L)
  252. tbl <- [email hidden] %>%
  253. group_by(orig.ident, BIAD_id, final_cell_type) %>%
  254. summarise(n = n()) %>%
  255. mutate(freq = (n / sum(n))*100)
  256. tbl
  258. write.csv(tbl, file = paste0(resultsdirectory, "/sheets/revisions_120D/20250717_MPSIIIAvHealthy_CellTypes_PatientProportions.csv"))
  260. Sanfilippo_120D$final_cell_type <- factor(Sanfilippo_120D$final_cell_type, levels = c( "Astrocyte","OPC","RG","Immature Neurons","Inhibitory","Excitatory")) # Re-order cell types for visualisation
  262. ## Visualisation of different cell types in culture - Fig 1B ##
  263. umap <- DimPlot(Sanfilippo_120D, reduction = "umap", group.by = "final_cell_type", split.by = "orig.ident", cols = c("firebrick4", "darkorchid2", "orange", "forestgreen", "mediumblue","deeppink2"), pt.size = 0.5, alpha = 1.0, order = T, label = FALSE, repel = FALSE) + ##you can replace the group.by variable to any column on the metadata
  264. theme(legend.position = "bottom", legend.text = element_text(size = 10), legend.justification = "center",
  265. axis.text.x = element_text(size = 22), axis.text.y = element_text(size = 22),
  266. axis.title.x = element_text(size = 22), axis.title.y = element_text(size = 22),
  267. panel.border = element_rect(fill=NA, colour = "black", size=0.5), plot.title= element_blank()) +
  268. coord_cartesian(ylim = c(-10, 10))
  269. umap <- umap + NoAxes() #+ NoLegend()
  270. umap
  271. ggsave(umap, filename = paste0(resultsdirectory,"/plots/revisions_120D/UMAPs_120D/20250717_ControlMPSIIIA_CellTypes_E-ISplit_percent.mito_BIAD-id.pdf", sep = ""), height = 6, width = 10, dpi = 300, useDingbats = FALSE, limitsize = F)
  273. ## Cell types separated by each donor - Supp Fig 3A ##
  274. Sanfilippo_120D$BIAD_id <- factor(Sanfilippo_120D$BIAD_id, levels = c( "B04-cNP11","B13-cNP20","B12-cNP24","B06-cNP08","B05-cNP21","B11-cNP30", "B03-cNP13", "B02-cNP16", "B01-cNP06", "B10-cNP25"))
  276. biad_ids <- unique(Sanfilippo_120D$BIAD_id)[1:10]
  278. umap_plots <- list()
  280. for (id in biad_ids) {
  281. umap <- DimPlot(Sanfilippo_120D, reduction = "umap", group.by = "final_cell_type",
  282. cells = which(Sanfilippo_120D$BIAD_id == id), # Subset cells by BIAD_id
  283. cols = c("firebrick4", "darkorchid2", "orange", "forestgreen", "mediumblue","deeppink2"),
  284. pt.size = 0.5, alpha = 1.0, order = TRUE, label = FALSE, repel = FALSE) +
  285. ggtitle(id) +
  286. theme(legend.position = "bottom",
  287. legend.text = element_text(size = 10),
  288. legend.justification = "center",
  289. axis.text.x = element_text(size = 22),
  290. axis.text.y = element_text(size = 22),
  291. axis.title.x = element_text(size = 22),
  292. axis.title.y = element_text(size = 22),
  293. panel.border = element_rect(fill = NA, colour = "black", size = 0.5),
  294. plot.title = element_text(size = 22, hjust = 0.5)) +
  295. coord_cartesian(ylim = c(-10, 10)) +
  296. NoAxes() #+ NoLegend() # Uncomment NoLegend if desired
  298. umap_plots[[id]] <- umap
  299. }
  301. combined_plot <- wrap_plots(umap_plots, nrow = 2, ncol = 5)
  302. combined_plot
  304. ggsave(combined_plot, filename = paste0(resultsdirectory,"/plots/revisions_120D/UMAPs_120D/20250718_ControlMPSIIIA_BIAD-id_neurons-split_noaxes_percent.mito_BIAD-id.pdf", sep = ""), height = 10.5, width = 20, dpi = 300, useDingbats = FALSE, limitsize = F)
  306. tbl <- [email hidden] %>%
  307. group_by(orig.ident, BIAD_id, final_cell_type) %>%
  308. summarise(n = n()) %>%
  309. mutate(freq = (n / sum(n))*100)
  310. tbl
  312. write.csv(tbl, file = paste0(resultsdirectory, '/sheets/revisions_120D/20250604_120D_BIAD-id_nodonordownsample_cell_types.csv'))
  314. ## Cell cycle analysis - Supp Fig 3F ##
  315. umap_density <- [email hidden][,c("orig.ident","cell_cycle")]
  316. umap_density <- cbind(umap_density, Sanfilippo_120D@reductions$[email hidden])
  318. scatterPlot <- ggplot(umap_density[na.omit(umap_density$cell_cycle == "NonCycling"),], aes(x = umap_1, y = umap_2)) +
  319. geom_point(data = umap_density[na.omit(umap_density$cell_cycle != "NonCycling"),], color = "grey90", size = 0.5) +
  320. geom_pointdensity(adjust = 0.2, size = 0.5) +
  321. scale_colour_gradient(low = "ghostwhite", high = "firebrick1", breaks = c(0,45,90), limits=c(0,90), na.value = "black") + #
  322. facet_wrap(~orig.ident, labeller = labeller(orig.ident = c("Control_120D" = "Healthy\n n= 6464 (98%)", "Sanfilippo_120D" = "MPSIIIA\n n= 6348 (96%)"))) +
  323. theme_bw() +
  324. theme(plot.title = element_text(size = 20, hjust = 0.5), plot.tag.position = "top",plot.subtitle = element_text(size = 20, hjust = 0.5, face = "plain"), plot.title.position = "plot", legend.position = "bottom", legend.text = element_text(size = 16, vjust = 0.5), legend.title = element_blank(), axis.text.x = element_text(size = 20), axis.text.y = element_text(size = 20), axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20), panel.border = element_rect(fill=NA, colour = "black", size=0.5), panel.grid = element_blank(), strip.text.x = element_text(size=20),strip.background = element_blank()) + NoAxes()
  325. scatterPlot
  327. scatterPlot <- ggplot(umap_density[na.omit(umap_density$cell_cycle == "G1S"),], aes(x = umap_1, y = umap_2)) +
  328. geom_point(data = umap_density[na.omit(umap_density$cell_cycle != "G1S"),], color = "grey90", size = 0.5) +
  329. geom_pointdensity(adjust = 0.2, size = 0.5) +
  330. scale_colour_gradient(low = "ghostwhite", high = "green1", breaks = c(0,25), limits=c(0,25), na.value = "black") + #
  331. facet_wrap(~orig.ident, labeller = labeller(orig.ident = c("Control_120D" = "Healthy\n n= 91 (1%)", "Sanfilippo_120D" = "MPSIIIA\n n= 153 (2%)"))) +
  332. theme_bw() +
  333. theme(plot.title = element_text(size = 20, hjust = 0.5), plot.tag.position = "top",plot.subtitle = element_text(size = 20, hjust = 0.5, face = "plain"), plot.title.position = "plot", legend.position = "bottom", legend.text = element_text(size = 16, vjust = 0.5), legend.title = element_blank(), axis.text.x = element_text(size = 20), axis.text.y = element_text(size = 20), axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20), panel.border = element_rect(fill=NA, colour = "black", size=0.5), panel.grid = element_blank(), strip.text.x = element_text(size=20),strip.background = element_blank()) + NoAxes()
  334. scatterPlot
  336. scatterPlot <- ggplot(umap_density[na.omit(umap_density$cell_cycle == "G2M"),], aes(x = umap_1, y = umap_2)) +
  337. geom_point(data = umap_density[na.omit(umap_density$cell_cycle != "G2M"),], color = "grey90", size = 0.5) +
  338. geom_pointdensity(adjust = 0.2, size = 0.5) +
  339. scale_colour_gradient(low = "ghostwhite", high = "gold", breaks = c(0,16,32), limits=c(0,32), na.value = "black") + #
  340. facet_wrap(~orig.ident, labeller = labeller(orig.ident = c("Control_120D" = "Healthy\n n= 64 (1%)", "Sanfilippo_120D" = "MPSIIIA\n n= 118 (2%)"))) +
  341. theme_bw() +
  342. theme(plot.title = element_text(size = 20, hjust = 0.5), plot.tag.position = "top",plot.subtitle = element_text(size = 20, hjust = 0.5, face = "plain"), plot.title.position = "plot", legend.position = "bottom", legend.text = element_text(size = 16, vjust = 0.5), legend.title = element_blank(), axis.text.x = element_text(size = 20), axis.text.y = element_text(size = 20), axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20), panel.border = element_rect(fill=NA, colour = "black", size=0.5), panel.grid = element_blank(), strip.text.x = element_text(size=20),strip.background = element_blank()) + NoAxes()
  343. scatterPlot
  345. ggsave(scatterPlot, filename = paste0(resultsdirectory,"/plots/revisions_120D/UMAPs_120D/20250523_Control_MPSIIIA_G2M_Density_final.pdf", sep = ""), height = 6.5, width = 10, dpi = 300, useDingbats = FALSE, limitsize = F)
  347. tbl <- [email hidden] %>%
  348. group_by(final_cell_type, new_cell_types) %>%
  349. summarise(n = n()) %>%
  350. mutate(freq = (n / sum(n))*100)
  351. tbl
  353. ```
  355. 4. ***Now we isolate just the neurons and separate them into excitatory and inhibitory neurons***
  356. ```{Identify excitatory and inhibitory neurons}
  357. ## Isolate the 'Neurons' from broad_cell_types to further separate them into excitatory and inhibitory neurons
  358. table([email hidden][['final_cell_type']])
  359. Idents(Sanfilippo_120D) <- 'final_cell_type'
  360. Sanfilippo_120D_Neurons <- subset(Sanfilippo_120D, idents = c('Excitatory', 'Inhibitory'))
  361. table(Sanfilippo_120D_Neurons$orig.ident)
  362. Idents(Sanfilippo_120D_Neurons) <- "orig.ident"
  363. Sanfilippo_120D_Neurons <- subset(Sanfilippo_120D_Neurons, downsample = 1015)
  364. [email hidden] %>% group_by(orig.ident) %>% summarise(n = n())
  366. Sanfilippo_120D_Neurons$neuron_types_ident <- paste(Sanfilippo_120D_Neurons$orig.ident, Sanfilippo_120D_Neurons$final_cell_type, sep = "_")
  367. Sanfilippo_120D_Neurons$biad_id_type <- paste(Sanfilippo_120D_Neurons$BIAD_id, Sanfilippo_120D_Neurons$final_cell_type, sep = "_")
  369. ## Save out proportion of excitatory and inhibitory neurons
  370. tbl <- [email hidden] %>%
  371. group_by(orig.ident, final_cell_type) %>%
  372. summarise(n = n()) %>%
  373. mutate(freq = (n / sum(n))*100)
  374. tbl
  375. write.csv(tbl, file = paste0(resultsdirectory, '/sheets/revisions_120D/20250604_120D_BIAD-id_nodonordownsample_Split_Neuron-type-proportions.csv'))
  377. # proportions used in Supp Fig 20B #
  378. tbl <- [email hidden] %>%
  379. group_by(orig.ident, BIAD_id, final_cell_type) %>%
  380. summarise(n = n()) %>%
  381. mutate(freq = (n / sum(n))*100)
  382. tbl
  384. write.csv(tbl, file = paste0(resultsdirectory, '/sheets/[Sanfilippo] DE_120D/20250604_MPSIIIA-Healthy_NeuronType_BIAD-id_proportions.csv'))
  386. ## Prepare object for UMAP reduction
  387. Sanfilippo_120D_Neurons <- SCTransform(Sanfilippo_120D_Neurons, vars.to.regress = c("percent.mt", "BIAD_id"))
  389. ## Save out seurat object
  390. save(Sanfilippo_120D_Neurons_New, file = paste0(datadirectory, "/RData/120D/revisions/Combined_Control_MPSIIIA_120D_NOambientRNA_NeuronsOnly.RData")) # As above, save after running SCTransform to save time
  391. load(file = paste0(datadirectory, "/RData/120D/revisions/Combined_Control_MPSIIIA_120D_NOambientRNA_NeuronsOnly.RData"))
  393. ## Re-run PCA and clustering for visualisation of neuronal subtypes
  394. var_features <- Sanfilippo_120D_Neurons_New@assays[["SCT"]]@var.features
  395. var_features <- var_features[grep("^LINC|^LOC|^AC[0-9]*\\.[0-9]|^AP[0-9]*\\.[0-9]|^AL[0-9]*\\.[0-9]",var_features, invert = T)]
  397. Sanfilippo_120D_Neurons_New <- RunPCA(Sanfilippo_120D_Neurons_New, npcs = 100, assay = "SCT", features = var_features)
  398. Sanfilippo_120D_Neurons_New <- NormalizeData(Sanfilippo_120D_Neurons_New, assay = "RNA", normalization.method = "LogNormalize")
  400. min.pc <- significant_pcs(Sanfilippo_120D_Neurons_New)
  401. Sanfilippo_120D_Neurons_New <- FindNeighbors(Sanfilippo_120D_Neurons_New, dims = 1:min.pc, reduction = "pca")
  402. Sanfilippo_120D_Neurons_New <- FindClusters(Sanfilippo_120D_Neurons_New, resolution = 0.1)
  403. Sanfilippo_120D_Neurons_New <- RunUMAP(Sanfilippo_120D_Neurons_New, reduction = "pca", dims= 1:min.pc, n.components = 2L)
  405. ## Visualise spread of excitatory/inhibitory neurons - Fig 7A ##
  406. umap <- DimPlot(Sanfilippo_120D_Neurons_New, reduction = "umap", group.by = "nowakowski_celltypes", split.by = "orig.ident", pt.size = 1.25, alpha = 1.0, order = T, label = FALSE, repel = FALSE) + ##you can replace the group.by variable to any column on the metadata
  407. theme(legend.position = "bottom", legend.text = element_text(size = 10), legend.justification = "center",
  408. axis.text.x = element_text(size = 22), axis.text.y = element_text(size = 22),
  409. axis.title.x = element_text(size = 22), axis.title.y = element_text(size = 22),
  410. panel.border = element_rect(fill=NA, colour = "black", size=0.5), plot.title= element_blank())# +
  411. #coord_cartesian(ylim = c(-10, 10))
  412. umap <- umap + NoAxes() #+ NoLegend()
  413. umap
  415. ggsave(umap, filename = paste0(resultsdirectory,"/plots/revisions_120D/UMAPs_120D/20250529_Excit-Inhib_percent.mt_BIAD-id_nodonordownsample.pdf", sep = ""), height = 6, width = 10, dpi = 300, useDingbats = FALSE, limitsize = F)
  417. ## Confirm neuronal populations express correct excitatory or inhibitory markers - Fig 7A ##
  418. features <- c("NXPH1","GAD1","GAD2","SOX6","SLC6A1","SHANK2","GRIA2","GRM5","SOX5","CAMK2A") # excitatory and inhibitory markers
  419. Idents(Sanfilippo_120D_Neurons_New) <- 'neuron_types_ident'
  420. dotplot <- DotPlot(Sanfilippo_120D_Neurons_New, features = features, cols = c("blue", "firebrick2"), idents = c("Sanfilippo_120D_Excitatory", "Sanfilippo_120D_Inhibitory")) + theme_bw() +
  421. #labs(title = "Excitatory Neurons", subtitle = "[n = 26 cells]") +
  422. theme(plot.title = element_text(size = 20, hjust = 0.5), plot.tag.position = "top",plot.subtitle = element_text(size = 20, hjust = 0.5, face = "plain"), plot.title.position = "plot", legend.position = "right", legend.text = element_text(size = 16, vjust = 0.5), axis.text.x = element_blank(), axis.text.y = element_text(size = 20), axis.title.x = element_blank(), axis.title.y = element_text(size = 20), panel.border = element_rect(fill=NA, colour = "black", size=0.5), panel.grid = element_blank(), strip.text.x = element_text(size=20),strip.background = element_blank()) + coord_flip() #+ scale_y_reordered()
  424. dotplot$data$id <- factor(dotplot$data$id,
  425. levels = c("Sanfilippo_120D_Excitatory", "Sanfilippo_120D_Inhibitory"))
  427. dotplot
  428. ggsave(dotplot, filename = paste0(resultsdirectory,"/plots/revisions_120D/20250529_MPSIIIA_ExcVInhibNeurons_DotPlot_nodonordownsamp.pdf"), height = 5, width = 4.5, dpi = 300)
  430. ## Density plots separating excitatory and inhibitory neurons - Supp Fig 20C ##
  431. ## Excitatory neurons
  432. umap_density <- [email hidden][,c("orig.ident","final_cell_type")]
  433. umap_density <- cbind(umap_density, Sanfilippo_120D_Neurons_New@reductions$[email hidden])
  435. scatterPlot <- ggplot(umap_density[na.omit(umap_density$final_cell_type == "Excitatory"),], aes(x = umap_1, y = umap_2)) +
  436. geom_point(data = umap_density[na.omit(umap_density$final_cell_type != "Excitatory"),], color = "grey90", size = 1.25) +
  437. geom_pointdensity(adjust = 0.2, size = 1.25) +
  438. scale_colour_gradient(low = "ghostwhite", high = "deeppink2", breaks = c(0,24), limits=c(0,24), na.value = "black") + #
  439. facet_wrap(~orig.ident, labeller = labeller(orig.ident = c("Control_120D" = "Healthy\n n= 320 (31.5%)", "Sanfilippo_120D" = "MPSIIIA\n n= 505 (49.8%)"))) +
  440. theme_bw() +
  441. theme(plot.title = element_text(size = 20, hjust = 0.5), plot.tag.position = "top",plot.subtitle = element_text(size = 20, hjust = 0.5, face = "plain"), plot.title.position = "plot", legend.position = "bottom", legend.text = element_text(size = 16, vjust = 0.5), legend.title = element_blank(), axis.text.x = element_text(size = 20), axis.text.y = element_text(size = 20), axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20), panel.border = element_rect(fill=NA, colour = "black", size=0.5), panel.grid = element_blank(), strip.text.x = element_text(size=20),strip.background = element_blank()) + NoAxes()
  442. scatterPlot
  444. ggsave(scatterPlot, filename = paste0(resultsdirectory,"/plots/revisions_120D/UMAPs_120D/20250529_Control_MPSIIIA_ExcitatoryNeurons_Density_nodonordownsamp.pdf", sep = ""), height = 6, width = 10, dpi = 300, useDingbats = FALSE, limitsize = F)
  447. ## Inhibitory neurons
  448. scatterPlot <- ggplot(umap_density[na.omit(umap_density$final_cell_type == "Inhibitory"),], aes(x = umap_1, y = umap_2)) +
  449. geom_point(data = umap_density[na.omit(umap_density$final_cell_type != "Inhibitory"),], color = "grey90", size = 1.25) +
  450. geom_pointdensity(adjust = 0.2, size = 1.25) +
  451. scale_colour_gradient(low = "ghostwhite", high = "mediumblue", breaks = c(0,19), limits=c(0,19), na.value = "black") + #
  452. facet_wrap(~orig.ident, labeller = labeller(orig.ident = c("Control_120D" = "Healthy\n n= 695 (68.5%)", "Sanfilippo_120D" = "MPSIIIA\n n= 510 (50.2%)"))) +
  453. theme_bw() +
  454. theme(plot.title = element_text(size = 20, hjust = 0.5), plot.tag.position = "top",plot.subtitle = element_text(size = 20, hjust = 0.5, face = "plain"), plot.title.position = "plot", legend.position = "bottom", legend.text = element_text(size = 16, vjust = 0.5), legend.title = element_blank(), axis.text.x = element_text(size = 20), axis.text.y = element_text(size = 20), axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20), panel.border = element_rect(fill=NA, colour = "black", size=0.5), panel.grid = element_blank(), strip.text.x = element_text(size=20),strip.background = element_blank()) + NoAxes()
  455. scatterPlot
  457. ggsave(scatterPlot, filename = paste0(resultsdirectory,"/plots/revisions_120D/UMAPs_120D/20250529_Control_MPSIIIA_Inhib_Density_redone_nodonordownsamp.pdf", sep = ""), height = 6, width = 10, dpi = 300, useDingbats = FALSE, limitsize = F)
  459. ## Density plot looking at split of excitatory and inhibitory neurons by donor - Supp Fig 20E ##
  460. Sanfilippo_120D_Neurons_New$BIAD_id <- factor(Sanfilippo_120D_Neurons_New$BIAD_id, levels = c( "B04-cNP11","B13-cNP20","B12-cNP24","B06-cNP08","B05-cNP21","B11-cNP30", "B03-cNP13", "B02-cNP16", "B01-cNP06", "B10-cNP25"))
  461. umap_density <- [email hidden][,c("BIAD_id","final_cell_type")]
  462. umap_density <- cbind(umap_density, Sanfilippo_120D_Neurons_New@reductions$[email hidden])
  464. library(ggnewscale)
  465. scatterPlot <- ggplot(umap_density, aes(x = umap_1, y = umap_2)) +
  466. geom_point(data = umap_density[!(umap_density$final_cell_type %in% c("Excitatory", "Inhibitory")),], color = "grey80", size = 1.25) +
  467. geom_pointdensity(data = umap_density[umap_density$final_cell_type == "Excitatory",], adjust = 0.2, size = 1.25) +
  468. scale_colour_gradient(name = "Excitatory", low = "ghostwhite", high = 'deeppink2', breaks = c(0,24), limits=c(0,24), na.value = 'black') +
  469. new_scale_colour() +
  470. geom_pointdensity(data = umap_density[umap_density$final_cell_type == "Inhibitory",], adjust = 0.2, size = 1.25) +
  471. scale_colour_gradient(name = "Inhibitory", low = "ghostwhite", high = "mediumblue", breaks = c(0,19), limits=c(0,19), na.value = "black") +#breaks = c(0,24), limits = c(0,24),
  472. facet_wrap(~BIAD_id, labeller = labeller(BIAD_id = c(
  473. "B12-cNP24" = 'BIAD12 (2 clones)',
  474. "B06-cNP08" = 'BIAD06 (2 clones)',
  475. "B13-cNP20" = 'BIAD13 (2 clones)',
  476. "B04-cNP11" = 'BIAD04 (2 clones)',
  477. "B05-cNP21" = 'BIAD05 (1 clone)',
  478. "B11-cNP30" = 'BIAD11 (2 clones)',
  479. 'B10-cNP25' = 'BIAD10 (2 clones)',
  480. "B01-cNP06" = 'BIAD01 (2 clones)',
  481. "B03-cNP13" = 'BIAD03 (2 clones)',
  482. "B02-cNP16" = 'BIAD02 (2 clones)')), ncol = 5) +
  483. theme_bw() +
  484. theme(
  485. plot.title = element_text(size = 20, hjust = 0.5),
  486. plot.subtitle = element_text(size = 16, hjust = 0.5, face = "italic"),
  487. legend.position = "bottom",
  488. legend.text = element_text(size = 16, vjust = 0.5),
  489. axis.text.x = element_text(size = 20),
  490. axis.text.y = element_text(size = 20),
  491. axis.title.x = element_text(size = 20),
  492. axis.title.y = element_text(size = 20),
  493. panel.border = element_rect(fill = NA, colour = "black", size = 0.25),
  494. panel.grid = element_blank(),
  495. strip.text.x = element_text(size = 20),
  496. strip.background = element_blank()
  497. ) +
  498. coord_equal() + NoAxes()
  500. # Display the scatter plot
  501. scatterPlot
  502. ggsave(scatterPlot, filename = paste0(resultsdirectory,"/plots/revisions_120D/UMAPs_120D/20250529_ControlMPSIIIA_UMAP_DonorID_ExcitatoryInhibitory_Density_nodonordownsamp.pdf", sep = ""), height = 12, width = 14, dpi = 300, useDingbats = FALSE, limitsize = F)
  504. ```
  506. 5. We perform differential expression analysis on different neuronal subtypes
  507. ``` {DE Analysis and visulation for excitatory and inhibitory neurons}
  508. ## We conduct a differential expression to identify genes that are significantly expressed in Sanfilippo compared to Control.
  509. ## Excitatory neurons
  510. Idents(Sanfilippo_120D_Neurons_New) <- 'neuron_types_ident'
  511. DE_Control_MPSIIIA_EN <- FindMarkers(Sanfilippo_120D_Neurons_New, assay = 'RNA', slot = 'data', ident.1 = "Sanfilippo_120D_Excitatory", ident.2 = "Control_120D_Excitatory", logfc.threshold = 0, min.pct = 0.1, min.diff.pct = 0.1, test.use = "wilcox")
  512. DE_Control_MPSIIIA_EN$hgnc_symbol <- rownames(DE_Control_MPSIIIA_EN)
  513. DE_Control_MPSIIIA_EN <- significantly_expressed_genes(DE_Control_MPSIIIA_EN)
  515. ## Save out DE results
  516. write.csv(DE_Control_MPSIIIA_EN, file = paste0(resultsdirectory, "/sheets/revisions_120D/DE_120D/20250718_DE_Control_MPSIIIA_ExcitatoryNeurons_wilcox_nodonordownsample.csv"))
  518. # #Inhibitory neurons
  519. DE_Control_MPSIIIA_IN <- FindMarkers(Sanfilippo_120D_Neurons_New, assay = 'RNA', slot = 'data', ident.1 = "Sanfilippo_120D_Inhibitory", ident.2 = "Control_120D_Inhibitory", logfc.threshold = 0, min.pct = 0.1, min.diff.pct = 0.1, test.use = "wilcox")
  520. DE_Control_MPSIIIA_IN$hgnc_symbol <- rownames(DE_Control_MPSIIIA_IN)
  521. DE_Control_MPSIIIA_IN <- significantly_expressed_genes(DE_Control_MPSIIIA_IN)
  523. ## Save out DE results
  524. write.csv(DE_Control_MPSIIIA_IN, file = paste0(resultsdirectory, "/sheets/revisions_120D/DE_120D/20250718_DE_Control_MPSIIIA_InhibitoryNeurons_wilcox_nodonordownsample.csv"))
  526. ## Dot plots to visualise key genes up or down regulated in excitatory and inhibitory neurons (disease v healthy) - Figure 5B
  527. excitatory_genes <- c("PEG3", "AC069277.1", "XIST", "AC006115.2", "GRM7-AS3", "IQCJ-SCHIP1", "NELL2", "NEUROD6", "SATB2", "TAFA1") # top 5 up/down DEG
  528. inhibitory_genes <- c("PEG3", "AC069277.1", "MTRNR2L1", "XIST", "GRM7-AS3", "CRLF1", "PTGDS", "CHL1", "IFITM3", "NEUROD6") # top 5 up/down DEG
  530. excit_v_inhib <- DotPlot(Sanfilippo_120D_Neurons_New, features = inhibitory_genes, cols = c("blue", "firebrick2"), idents = c("Sanfilippo_120D_Inhibitory", "Control_120D_Inhibitory")) + theme_bw() +
  531. #labs(title = "Excitatory Neurons", subtitle = "[n = 26 cells]") +
  532. theme(plot.title = element_text(size = 20, hjust = 0.5), plot.tag.position = "top",plot.subtitle = element_text(size = 20, hjust = 0.5, face = "plain"), plot.title.position = "plot", legend.position = "right", legend.text = element_text(size = 16, vjust = 0.5), axis.text.x = element_blank(), axis.text.y = element_text(size = 20), axis.title.x = element_blank(), axis.title.y = element_text(size = 20), panel.border = element_rect(fill=NA, colour = "black", size=0.5), panel.grid = element_blank(), strip.text.x = element_text(size=20),strip.background = element_blank()) + coord_flip() #+ scale_y_reordered()
  534. excit_v_inhib$data$id <- factor(excit_v_inhib$data$id,
  535. levels = c("Sanfilippo_120D_Inhibitory", "Control_120D_Inhibitory"))
  537. excit_v_inhib
  538. ggsave(excit_v_inhib, filename = paste0(resultsdirectory,"/plots/revisions_120D/DE_120D/20250529_SanfilvControl_InhibitoryNeurons_DotPlot_nodonordownsample.pdf"), height = 5, width = 4.75, dpi = 300)
  540. ## Volcano Plot - Fig 7B (swap the data between excitatory or inhibitory DE results) ##
  541. VP <- ggplot(DE_Control_MPSIIIA_EN, aes(avg_log2FC, log10p_adj, label=delabel, color=diffexpressed)) + geom_point(alpha = 1, size = 3, fill = "grey90") + ## transform to change the order of the plots in facet_wrap()
  542. theme_light() +
  543. scale_fill_gradientn(colours = c(rainbow(100)), ## colors in the scale
  544. breaks = seq(0, 1, 0.2),
  545. limits = c(0, 1)) +
  546. scale_color_manual(values = c("grey60","blue","red")) +
  547. geom_text_repel(data=DE_Control_MPSIIIA_EN, size = 8, max.overlaps = 75, direction = "both") +
  548. theme(legend.position = "right", axis.text.x = element_text(size = 28), axis.text.y = element_text(size = 28), axis.title.x = element_text(size = 30), axis.title.y = element_text(size = 30), panel.border = element_rect(fill=NA, colour = "black", size=0.5)) +
  549. labs(y=expression('-Log'[10]*' P'[adj]), x=expression(atop('Log'[2]*' fold change', paste("[MPSIIIA/Control]")))) +
  550. scale_x_continuous(breaks = seq(-7, 7, by = 2)) +
  551. scale_y_continuous(breaks = seq(0, 90, by = 20)) +
  552. coord_cartesian(xlim = c(-7, 7), ylim = c(0, 80)) +
  553. geom_hline(yintercept = -log10(0.05), linetype = "dashed", col = "darkgrey", size = 1) + NoLegend()
  554. VP
  555. ggsave(VP, filename = paste0(resultsdirectory,"/plots/revisions_120D/DE_120D/20250529_VP_MPSIIIA_Control_DE_ExcitatoryNeurons_nodonordownsample.pdf"), height = 10, width = 10, dpi = 300)
  557. # Looking at expression of E/I genes by patient - Fig 7F #
  558. Idents(Sanfilippo_120D_Neurons_New) <- "BIAD_id"
  559. average_expression <- AverageExpression(Sanfilippo_120D_Neurons_New, assay = "RNA", slot = "counts", features = c("NLGN1", "CDKL5", "DISC1", "SHANK2", "DLG4", "GPHN", "NLGN2", "GABRA1", "GLRB", "GRIA1", "GRIA2", "GRIA3", "GRIA4", "HOMER1", "CAMK2A", "ARHGEF9", "KCC2", "GRIN1", "GRIN2A", "GRIN2B", "GABARAP", "GABRA2"), group.by = c("BIAD_id", "final_cell_type"))
  560. average_expression <- as.data.frame(average_expression)
  561. average_expression <- t(average_expression)
  562. write.csv(average_expression, file = paste0(resultsdirectory, "/sheets/revisions_120D/DE_120D/20250717_synapsegenes_e-i-neurons_biad-id_rawcounts_120d.csv"))
  563. ```
  565. 6.***We perform GSEA for excitatory and inhibitory neurons***
  566. ```{r GSEA for differentially expressed genes in excitatory and inhibitory neurons}
  567. ### Use this if you want to run fgsea analysis
  568. library(fgsea)
  569. library(data.table)
  571. ## Isolate avglog2FC values from previously ranked DE list for excitatory and inhibitory neurons
  572. upreg_EN <- DE_Control_MPSIIIA_EN[which(DE_Control_MPSIIIA_EN$avg_log2FC > 0 & DE_Control_MPSIIIA_EN$p_val_adj < 0.05),]
  573. downreg_EN <- DE_Control_MPSIIIA_EN[which(DE_Control_MPSIIIA_EN$avg_log2FC < 0 & DE_Control_MPSIIIA_EN$p_val_adj < 0.05),]
  574. dysreg_EN <- rbind(upreg_EN, downreg_EN)
  576. ranked_genes_EN <- dysreg_EN$avg_log2FC
  577. names(ranked_genes_EN) <- rownames(dysreg_EN)
  579. upreg_IN <- DE_Control_MPSIIIA_IN[which(DE_Control_MPSIIIA_IN$avg_log2FC > 0 & DE_Control_MPSIIIA_IN$p_val_adj < 0.05),]
  580. downreg_IN <- DE_Control_MPSIIIA_IN[which(DE_Control_MPSIIIA_IN$avg_log2FC < 0 & DE_Control_MPSIIIA_IN$p_val_adj < 0.05),]
  581. dysreg_IN <- rbind(upreg_IN, downreg_IN)
  583. ranked_genes_IN <- dysreg_IN$avg_log2FC
  584. names(ranked_genes_IN) <- rownames(dysreg_IN)
  586. Idents(Sanfilippo_120D_Neurons_New) <- 'neuron_types_ident'
  588. # Create a dataframe with avg log2FC values for a list of synaptic genes
  589. av_exp_df <- as.data.frame(AverageExpression(Sanfilippo_120D_Neurons_New, assays = "RNA", slot = "data", features = unique(unlist(synaptic_genes$hgnc_symbol))))
  590. av_exp_df$Logadj_Control_EN <- log2((av_exp_df$RNA.Control.120D.Excitatory))
  591. av_exp_df$Logadj_MPSIIIA_EN <- log2((av_exp_df$RNA.Sanfilippo.120D.Excitatory))
  592. av_exp_df$Logadj_Control_IN <- log2((av_exp_df$RNA.Control.120D.Inhibitory))
  593. av_exp_df$Logadj_MPSIIIA_IN <- log2((av_exp_df$RNA.Sanfilippo.120D.Inhibitory))
  595. av_exp_df$avg_log2FC_ENs <- av_exp_df$Logadj_MPSIIIA_EN - av_exp_df$Logadj_Control_EN
  596. av_exp_df$avg_log2FC_ENs[is.na(av_exp_df$avg_log2FC_ENs)] <- 0
  597. av_exp_df$avg_log2FC_ENs[is.infinite(av_exp_df$avg_log2FC_ENs)] <- 0
  599. av_exp_df$avg_log2FC_INs <- av_exp_df$Logadj_MPSIIIA_IN - av_exp_df$Logadj_Control_IN
  600. av_exp_df$avg_log2FC_INs[is.na(av_exp_df$avg_log2FC_INs)] <- 0
  601. av_exp_df$avg_log2FC_INs[is.infinite(av_exp_df$avg_log2FC_INs)] <- 0
  603. ranked_genes <- as.vector(av_exp_df[,10])
  604. names(ranked_genes) <- rownames(av_exp_df)
  605. ranked_genes <- sort(ranked_genes, decreasing = T)
  607. ## Load list of synaptic genes
  608. pathways <- list("DEGs" = dysreg_EN$hgnc_symbol)
  609. pathways_IN <- list("DEGs" = dysreg_IN$hgnc_symbol)
  611. gse_all_EN <- fgsea(pathways,
  612. stats = ranked_genes,
  613. minSize = 1,
  614. maxSize = 500,
  615. nPermSimple = 10000,
  616. gseaParam = 1,
  617. eps = 0)
  619. gse_all_IN <- fgsea(pathways_IN,
  620. stats = ranked_genes,
  621. minSize = 1,
  622. maxSize = 500,
  623. nPermSimple = 10000,
  624. gseaParam = 1,
  625. eps = 0)
  627. ## Visualise enrichment for excitatory and inhibitory neurons - Fig 7E ##
  628. pdf(paste0(resultsdirectory,"/plots/revisions_120D/DE_120D/20250623_SynapticGene_Enrichment_Plot_InhibitoryNeurons_nodonordownsample.pdf"), width=7.5, height=5)
  629. EP <- EnrichmentPlot(pathways$DEGs,
  630. ranked_genes, gseaParam = 2, ticksSize = 0.2)
  632. EP
  633. dev.off()
  635. test_list <- list(gse_all$leadingEdge)
  637. gse_all_EN <- as.data.frame(gse_all_EN)
  638. gse_all_EN <- gse_all_EN[1:7]
  639. gse_all_IN <- as.data.frame(gse_all_IN)
  640. gse_all_IN <- gse_all_IN[1:7]
  642. write.csv(gse_all_EN, file = paste0(resultsdirectory, '/sheets/revisions_120D/20250623_SynapticEnrich_ExcitatoryNeurons_nodonordownsamp.csv'))
  643. write.csv(gse_all_IN, file = paste0(resultsdirectory, '/sheets/revisions_120D/20250623_SynapticEnrich_InhibitoryNeurons_nodonordownsamp.csv'))
  645. # Perform GO Pathway analysis on significantly dysregulated genes
  646. library(msigdbr)
  647. library(org.Hs.eg.db)
  648. library(enrichplot)
  649. library(clusterProfiler)
  650. library(DOSE)
  651. library(simplifyEnrichment)
  652. library(GOSemSim)
  654. ## Create a vector of all gene names in count matrix to use as background set of genes
  655. all_genes <- as.matrix(Sanfilippo_120D_Neurons_New@assays[["RNA"]]@layers[["counts"]])
  656. rownames(all_genes) <- rownames(Sanfilippo_120D_Neurons_New@assays[['RNA']])
  657. all_genes <- row.names(all_genes)
  659. # Create vector of gene names that are significantly dysregulated for GO enrichment
  660. geneList <- dysreg_EN$hgnc_symbol
  661. geneList_IN <- dysreg_IN$hgnc_symbol
  663. test_enrichment <- enrichGO(gene = geneList, # make sure it's running for either the excitatory list or inhibitory
  664. OrgDb = 'org.Hs.eg.db',
  665. keyType = 'SYMBOL',
  666. ont = 'BP',
  667. universe = all_genes,
  668. pvalueCutoff = 0.05,
  669. pAdjustMethod = 'BH',
  670. qvalueCutoff = 0.1,
  671. minGSSize = 10,
  672. maxGSSize = 500,
  673. readable = TRUE)
  675. job::job({inhib_enrichment <- enrichGO(gene = geneList_IN, # to be able to run two enrichment analyses at once
  676. OrgDb = 'org.Hs.eg.db',
  677. keyType = 'SYMBOL',
  678. ont = 'BP',
  679. universe = all_genes,
  680. pvalueCutoff = 0.05,
  681. pAdjustMethod = 'BH',
  682. qvalueCutoff = 0.1,
  683. minGSSize = 10,
  684. maxGSSize = 500,
  685. readable = TRUE)})
  687. ## Isolate GO ID's for similarity clustering
  688. go_id_EN <- test_enrichment@result$ID
  689. job::job({go_mat <- GO_similarity(go_id_EN, ont = "BP")})
  691. go_id <- inhib_enrichment@result$ID
  692. job::job({go_mat_inhib <- GO_similarity(go_id, ont = "BP")})
  694. ## generate clusters using louvain clustering.
  695. go_cluster_louvain <- simplifyGO(go_mat_inhib, method = "louvain") ## NEED TO RUN THIS FOR INHIBITORY NEURONS AND FINISH CHORD PLOT TO UPDATE FIGURE
  697. test_enrichment@result$Cluster_louvain <- go_cluster_louvain$cluster[match(test_enrichment@result$ID, go_cluster_louvain$id)]
  698. inhib_enrichment@result$Cluster_louvain <- go_cluster_louvain$cluster[match(inhib_enrichment@result$ID, go_cluster_louvain$id)]
  700. enrichment_result <- as.data.frame(test_enrichment@result)
  701. enrichment_result <- as.data.frame(inhib_enrichment@result)
  703. most_common_word=function(x){
  705. #Split every word into single words for counting
  706. splitTest=unlist(strsplit(x," "))
  708. #Counting words
  709. count=table(splitTest)
  711. #Sorting to select only the highest value, which is the first one
  712. count=count[order(count, decreasing=TRUE)][1:50]
  714. #Return the desired character.
  715. #By changing this you can choose whether it show the number of times a word repeats
  716. return(names(count))
  717. }
  719. cluster_words <- vector("list", 12) # List to store results for each cluster
  721. # Loop through clusters 1 to 18
  722. for (i in 1:12) {
  723. # Subset descriptions for the current cluster
  724. cluster_descriptions <- test_enrichment@result[test_enrichment@result$Cluster_louvain == i,]
  726. # Apply most_common_word to the descriptions
  727. cluster_words[[i]] <- most_common_word(cluster_descriptions)
  728. }
  730. # Optionally, name the list elements for clarity
  731. names(cluster_words) <- paste0("Cluster_", 1:18)
  733. ## Add description based on pathways in the clusters and most common word function above
  734. # Excitatory
  735. test_enrichment@result$new_description <- ifelse(test_enrichment@result$Cluster_louvain == "1", "cellular differentiation and development", ifelse(test_enrichment@result$Cluster_louvain == "2", "ions and transport", ifelse(test_enrichment@result$Cluster_louvain == "3", "signaling pathways", ifelse(test_enrichment@result$Cluster_louvain == "4", "regulation and organisation", ifelse(test_enrichment@result$Cluster_louvain == "5", "regulation and organisation", ifelse(test_enrichment@result$Cluster_louvain == "6", "cellular metabolism", "other"
  736. ))))))
  738. # Inhibitory
  739. inhib_enrichment@result$new_description <- ifelse(inhib_enrichment@result$Cluster_louvain == "1", "cellular differentiation and development", ifelse(inhib_enrichment@result$Cluster_louvain == "2", "ions and transport", ifelse(inhib_enrichment@result$Cluster_louvain == "3", "signaling pathways", ifelse(inhib_enrichment@result$Cluster_louvain == "4", "regulation and organisation", ifelse(inhib_enrichment@result$Cluster_louvain == "5", "regulation and organisation", ifelse(inhib_enrichment@result$Cluster_louvain == "6", "cellular metabolism", "other"
  740. ))))))
  742. # Reorder by P.adjust value and save out to create table of top 5 GO IDs - Supp Fig 20E
  743. test_enrichment <- test_enrichment@result[order(test_enrichment@result$p.adjust),]
  744. inhib_enrichment <- inhib_enrichment@result[order(inhib_enrichment@result$p.adjust),]
  745. write.csv(test_enrichment, file = paste0(resultsdirectory, '/sheets/revisions_120D/20250604_Excitatory_neuron_DEG_PathwayAnalysis_nodonordownsample.csv'))
  747. enrichGO_list <- list()
  748. enrichGO_list$results <- data.frame("Category" = "BP","ID" = test_enrichment$Cluster_louvain, "Term" = test_enrichment$new_description, "Genes" = test_enrichment$geneID, "adj_pval" = test_enrichment$p.adjust) # swap enrichment result object between excitatory or inhibitory for cord plot visualisation
  749. enrichGO_list$results$Genes <- gsub("/",",",enrichGO_list$results$Genes)
  751. ## Make list of dysregulated genes in excitatory and inhibitory with their avglog2FC from DE analysis
  752. library(GOplot)
  753. enrichGO_list$geneList <- data.frame("ID" = rownames(dysreg_IN), "logFC" = dysreg_IN$avg_log2FC)
  754. circ <- circle_dat(enrichGO_list$results, enrichGO_list$geneList)
  755. enrichGO_list$process <- unique(circ$term)
  756. chord <- chord_dat(circ, enrichGO_list$geneList, enrichGO_list$process)
  757. chord <- chord[c(1:25, 751:755),] # Filter for the top 15 upregulated and top 15 downregulated genes
  758. chord <- chord[order(chord[,6], decreasing = T), ] # Order by decreasing logFC
  760. ## Order by different cluster descriptions
  761. chord <- chord[, c("cellular metabolism", "signaling pathways", "regulation and organisation", "cellular differentiation and development", "ions and transport", "logFC")] #
  763. ## Visualise chord plot - Fig 7D
  764. pdf(file = paste0(resultsdirectory,"/plots/revisions_120D/DE_120D/20250604_GOChord_top30DEgenes_upanddown_louvainclustering_inhibitoryneurons_nodonordownsample.pdf"), width = 14, height =14)
  765. cplot <- GOChord(chord,
  766. nlfc = 1,
  767. limit = c(0,0),
  768. space = 0.02,
  769. gene.space = 0.25,
  770. gene.order = "logFC",
  771. border.size = 0,
  772. ribbon.col = c("firebrick2", "orange","green3", "dodgerblue2", "deeppink"), #
  773. gene.size = 5,
  774. lfc.min = -7,
  775. lfc.max = 7,
  776. lfc.col = c("red", "grey90", "blue"))
  777. cplot
  778. dev.off()
  780. ```
  782. 7. ddressing reviewer concerns regarding BIAD03's high SOX2 iPSC expression and effect that has on phenotypes identified
  783. ```{Addressing reviewer concerns regarding BIAD03}
  784. ## Isolating excitatory neurons then performing FindMarkers on BIAD03 compared to other excitatory neurons
  785. Idents(Sanfilippo_120D_Neurons_New) <- 'neuron_types_ident'
  786. excitatory_only <- subset(Sanfilippo_120D_Neurons_New, idents = 'Sanfilippo_120D_Excitatory')
  787. excitatory_only <- NormalizeData(excitatory_only, assay = "RNA", normalization.method = "LogNormalize")
  789. # Perform DE on BIAD03 vs all MPSIIIA patients to understand unique population #
  790. Idents(excitatory_only) <- 'biad_id_type'
  791. BIAD03_unique <- FindMarkers(excitatory_only, assay = 'RNA', slot = 'data', ident.1 = "B03-cNP13_Excitatory", ident.2 = NULL, logfc.threshold = 0, min.pct = 0.1, min.diff.pct = 0.1, test.use = "wilcox")
  792. BIAD03_unique$hgnc_symbol <- rownames(BIAD03_unique)
  793. BIAD03_unique <- significantly_expressed_genes(BIAD03_unique)
  794. BIAD03_unique_upreg <- BIAD03_unique[which(BIAD03_unique$avg_log2FC > 0 &BIAD03_unique$p_val_adj < 0.05),] # find uniquely upregulated for GO analysis
  795. BIAD03_unique_downreg <- BIAD03_unique[which(BIAD03_unique$avg_log2FC < 0 &BIAD03_unique$p_val_adj < 0.05),] # find uniquely downregulated for GO analysis
  797. write.csv(BIAD03_unique, file = paste0(resultsdirectory, "/sheets/revisions_120D/DE_120D/20251203_DE_BIAD03-ENs_AllDonors_wilcox_nodonordownsample.csv"))
  799. # Volcano plot of uniquely dysregulated genes
  800. VP <- ggplot(BIAD03_unique, aes(avg_log2FC, log10p_adj, label=delabel, color=diffexpressed)) + geom_point(alpha = 1, size = 3, fill = "grey90") + ## transform to change the order of the plots in facet_wrap()
  801. theme_light() +
  802. scale_fill_gradientn(colours = c(rainbow(100)), ## colors in the scale
  803. breaks = seq(0, 1, 0.2),
  804. limits = c(0, 1)) +
  805. scale_color_manual(values = c("grey60","blue","red")) +
  806. geom_text_repel(data=BIAD03_unique, size = 8, max.overlaps = 75, direction = "both") +
  807. theme(legend.position = "right", axis.text.x = element_text(size = 28), axis.text.y = element_text(size = 28), axis.title.x = element_text(size = 30), axis.title.y = element_text(size = 30), panel.border = element_rect(fill=NA, colour = "black", size=0.5)) +
  808. labs(y=expression('-Log'[10]*' P'[adj]), x=expression(atop('Log'[2]*' fold change', paste("[MPSIIIA/Control]")))) +
  809. scale_x_continuous(breaks = seq(-7, 7, by = 2)) +
  810. scale_y_continuous(breaks = seq(0, 90, by = 20)) +
  811. coord_cartesian(xlim = c(-7, 7), ylim = c(0, 80)) +
  812. geom_hline(yintercept = -log10(0.05), linetype = "dashed", col = "darkgrey", size = 1) + NoLegend()
  813. VP
  814. ggsave(VP, filename = paste0(resultsdirectory,"/plots/revisions_120D/DE_120D/20251203_VP_BIAD03_ENs_GeneExpression_nodonordownsample.pdf"), height = 10, width = 10, dpi = 300)
  816. # pathways of DEGs for BIAD03 ENs
  817. library(msigdbr)
  818. library(org.Hs.eg.db)
  819. library(enrichplot)
  820. library(clusterProfiler)
  821. library(DOSE)
  822. library(simplifyEnrichment)
  823. library(GOSemSim)
  825. ## Create a vector of all gene names in count matrix to use as background set of genes
  826. all_genes <- as.matrix(excitatory_only@assays[["RNA"]]@layers[["counts"]])
  827. rownames(all_genes) <- rownames(excitatory_only@assays[['RNA']])
  828. all_genes <- row.names(all_genes)
  830. # Create vector of gene names that are significantly dysregulated for GO enrichment
  831. geneList <- BIAD03_unique_upreg$hgnc_symbol # swap this between upreg and downreg genes
  833. test_enrichment <- enrichGO(gene = geneList, # make sure it's running for either the excitatory list or inhibitory
  834. OrgDb = 'org.Hs.eg.db',
  835. keyType = 'SYMBOL',
  836. ont = 'BP',
  837. universe = all_genes,
  838. pvalueCutoff = 0.05,
  839. pAdjustMethod = 'BH',
  840. qvalueCutoff = 0.1,
  841. minGSSize = 10,
  842. maxGSSize = 500,
  843. readable = TRUE)
  845. go_id <- test_enrichment@result$ID
  847. GO_pathways <- dotplot(test_enrichment, showCategory = 20)
  848. GO_pathways
  850. ggsave(GO_pathways, filename = paste0(resultsdirectory,"/plots/revisions_120D/DE_120D/20251207_BIAD03_Pathways-upregDEGs.pdf"), height = 10, width = 7, dpi = 300)
  852. ## SOX2 expression per donor
  853. cortical_mark <- c("SOX2")
  855. control_matrix <- GetAssayData(object = Sanfilippo_120D,
  856. assay = "RNA",
  857. layer = "counts")
  858. control_matrix <- as.matrix(control_matrix)
  859. control_matrix <- control_matrix[cortical_mark,, drop = F]
  860. control_matrix <- as.data.frame(t(control_matrix))
  862. control_matrix$BIAD_id[match(rownames(control_matrix), rownames([email hidden]))] <- Sanfilippo_120D$BIAD_id[match(rownames(control_matrix),rownames([email hidden]))]
  864. control_matrix$condition[match(rownames(control_matrix), rownames([email hidden]))] <- Sanfilippo_120D$orig.ident[match(rownames(control_matrix),rownames([email hidden]))]
  866. MPSIIIA_matrix <- control_matrix[control_matrix$condition == "Sanfilippo_120D",]
  867. control_matrix <- control_matrix[control_matrix$condition == "Control_120D",]
  868. control_matrix_noexp <- control_matrix[rowSums(control_matrix[1]) == 0,]
  869. MPSIIIA_matrix_noexp <- MPSIIIA_matrix[rowSums(MPSIIIA_matrix[1]) == 0,]
  871. SOX2_gene_prop <- data.frame("Control_Number" = colSums(control_matrix > 0), "Control_Percent" = (colSums(control_matrix > 0)/nrow(control_matrix))*100,"MPSIIIA_Number" = colSums(MPSIIIA_matrix > 0), "MPSIIIA_Percent" = (colSums(MPSIIIA_matrix > 0)/nrow(MPSIIIA_matrix))*100)
  873. control_gene_prop_per_donor <- control_matrix |>
  874. group_by(BIAD_id) |>
  875. summarize(n = n(),
  876. across(-n, ~sum(. > 0)))
  878. control_gene_prop_per_donor$pct_expressing <- (control_gene_prop_per_donor$SOX2/control_gene_prop_per_donor$n)*100
  879. control_gene_prop_per_donor$condition <- "Healthy"
  880. control_gene_prop_per_donor$DIV <- "120D"
  882. MPSIIIA_gene_prop_per_donor <- MPSIIIA_matrix |>
  883. group_by(BIAD_id) |>
  884. summarize(n = n(),
  885. across(-n, ~sum(. > 0)))
  886. MPSIIIA_gene_prop_per_donor$pct_expressing <- (MPSIIIA_gene_prop_per_donor$SOX2/MPSIIIA_gene_prop_per_donor$n)*100
  887. MPSIIIA_gene_prop_per_donor$condition <- "MPSIIIA"
  888. MPSIIIA_gene_prop_per_donor$DIV <- "120D"
  890. av_exp_120D <- AggregateExpression(Sanfilippo_120D, features = "SOX2", group.by = c("BIAD_id"), assays = "RNA", slot = "counts")
  891. av_exp_120D <- as.data.frame(t(av_exp_120D$RNA))
  892. colnames(av_exp_120D) <- "SOX2_sum_exp"
  893. av_exp_120D$BIAD_id <- rownames(av_exp_120D)
  895. cortical_gene_prop_per_donor <- rbind(control_gene_prop_per_donor, MPSIIIA_gene_prop_per_donor)
  896. rownames(cortical_gene_prop_per_donor) <- cortical_gene_prop_per_donor$BIAD_id
  898. cortical_gene_prop_per_donor <- merge(cortical_gene_prop_per_donor, av_exp_120D, by = "BIAD_id")
  900. write.csv(cortical_gene_prop_per_donor, file = paste0(resultsdirectory, "/sheets/revisions_120D/Control_MPSIIIA_neurons_cortical_SOX2gene_prop_per_donor_nodownsample.csv"))
  901. ```
  903. 6. ***Look at proportion of neurons at 30DIV to address reviewer concerns about uneven proportions at 120D***
  904. ```{Assessing expression of key synaptic genes in 30DIV data}
  905. ### LOAD 30 DAY DATA ###
  906. ## Load required datasets
  907. load(file = paste0(datadirectory, "/RData/30D/Combined_Control_MPSIIIA_30D.RData"))
  909. tbl <- [email hidden] %>% group_by(orig.ident) %>% summarise(n = n())
  910. tbl
  911. Idents(Combined_Healthy_MPSIIIA) <- 'orig.ident'
  913. Combined_Healthy_MPSIIIA <- subset(Combined_Healthy_MPSIIIA, downsample = 4444)
  914. [email hidden] %>% group_by(orig.ident) %>% summarise(n = n())
  916. # Sort into excitatory and inhibitory neurons
  917. Combined_30D_neurons <- subset(Combined_Healthy_MPSIIIA, idents = c('Control_30D', 'Sanfilippo_30D'))
  918. table(Combined_30D_neurons$orig.ident)
  919. Idents(Combined_30D_neurons) <- "final_cell_type"
  920. Combined_30D_neurons <- subset(Combined_30D_neurons, idents = c('Excitatory', 'Inhibitory'))
  921. Idents(Combined_30D_neurons) <- 'orig.ident'
  922. Combined_30D_neurons <- subset(Combined_30D_neurons, downsample = 1952)
  924. # Donor proportions for excitatory and inhibitory neurons - Supp Fig 20A #
  925. tbl <- [email hidden] %>% group_by(orig.ident, BIAD_id, final_cell_type) %>% summarise(n = n())
  926. tbl
  928. write.csv(tbl, file = paste0(resultsdirectory, "/sheets/revisions_120D/DE_120D/20250612_30DIV_neuron-types_biad-id_nodwonsample.csv"))
  929. ```

Differential_Expression_Analysis_120D.Rmd at commit 04a5121, no license · at the source

Overview

Authors: Paris Mazzachi1,2, Ella McDonald1,2, Zarina Greenberg1,2, Alejandra Noreña Puerta1, Jenne Tran1, Manam Inushi De Silva1,2, Cade Christensen1, Robert Adams1,2, Sebastian Loskarn1, Helen Beard2, Michael Zabolocki1,2, Meera Elmasri1,2, Megan Maack3,4, Kristina L. Elvidge3,4, Mark R. Hutchinson5,6, Cara O’Neill7, Kim M. Hemsley2, Lisa Melton3, Nicholas Smith8,9, Cedric Bardy1,2,10
ORCID iDs: Paris Mazzachi, Ella McDonald, Zarina Greenberg, Alejandra Noreña Puerta, Manam Inushi De Silva, Cade Christensen, Robert Adams, Sebastian Loskarn, Michael Zabolocki, Megan Maack, Kristina L. Elvidge, Mark R. Hutchinson, Cara O’Neill, Lisa Melton, Nicholas Smith, Cedric Bardy
  1. Laboratory for Human Neurophysiology and Genetics, South Australian Health and Medical Research Institute (SAHMRI),Adelaide, SA Australia
  2. Flinders Health and Medical Research Institute, College of Medicine and Public Health, Flinders University,Adelaide, SA Australia
  3. Sanfilippo Children’s Foundation, Sydney, NSW Australia
  4. Childhood Dementia Initiative, Sydney, NSW Australia
  5. School of Biomedicine, Adelaide University,Adelaide, SA Australia
  6. Institute for Photonic and Advanced Sensing, Adelaide University,Adelaide, SA Australia
  7. Cure Sanfilippo Foundation, Columbia, SC USA
  8. Department of Neurology and Clinical Neurophysiology, Women’s and Children’s Health Network,Adelaide, SA Australia
  9. Paediatric Neurodegenerative Disease Research Group, Discipline of Paediatrics, College of Health, Adelaide University,Adelaide, SA Australia
  10. Brain Organoid Therapeutics, Adelaide, SA Australia
Institutions: South Australian Health and Medical Research Institute (Australia); Flinders University (Australia); Sanfilippo Children's Foundation (Australia); The University of Adelaide (Australia); Women's and Children's Health Network (Australia)
Journal: Nature communications, volume 17, issue 1, article 3161
Dates: received 12 September 2024; accepted 13 March 2026; published online 7 April 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41467-026-71112-9 · PMID 41946741 · PMCID PMC13057175 · OpenAlex W7151452559
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: intracellular / patch clamp (modality), human (organism), Alzheimer's / dementia (population), cellular / molecular (subfield)
Methods: Spectral & time-frequency, Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, Evoked potentials, Single-unit activity, calcium imaging
Keywords: Cellular neuroscience, Phenotypic screening, Patch clamp, Paediatric neurological disorders, Induced pluripotent stem cells
MeSH: Cerebral Cortex*, Dementia*, Induced Pluripotent Stem Cells*, Mucopolysaccharidosis III*, Synapses*, Action Potentials, Child, Humans, Nerve Net, Neurons (* major topic)
Topic: Lysosomal Storage Disorders Research (Physiology, Medicine), according to OpenAlex
Funding: National Health and Medical Research Council (MRFF MRF2024419); Australian Research Council (FT180100565)
Citations: cited by 1 paper (Europe PMC); 128 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

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

bardylab/MPSIIIA_snRNAseq_Synaptic_Paper_2026

License: none: the authors keep all their rights
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 04a512141b0b5beeb0fe7b925e8b3186311d53a9, 12 February 2026
Languages: R (3)
Size: 12 files, 3 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, 3 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: ggplot2 (3 files), Seurat (3 files), tidyverse (3 files), clusterProfiler (1 file), data.table (1 file), Harmony (1 file), patchwork (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
4 files
At the source: github.com/bardylab/MPSIIIA_snRNAseq_Synaptic_Paper_2026

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s41467-026-71112-9.

Tracing map

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

What the map holds:

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

No dataset and no data link were found in the paper.

Data availability statement

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

  • it says that the data are available on request

Read it in the paper: doi.org/10.1038/s41467-026-71112-9.

Versions

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

Version 1, 29 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 20 authors, 5 keywords, 10 MeSH terms, 2 funders, 128 references.

Cite

This paper

Mazzachi, P., McDonald, E., Greenberg, Z., Puerta, A. N., Tran, J., De Silva, M. I., Christensen, C., Adams, R., Loskarn, S., Beard, H., Zabolocki, M., Elmasri, M., Maack, M., Elvidge, K. L., Hutchinson, M. R., O’Neill, C., Hemsley, K. M., Melton, L., Smith, N., & Bardy, C. (2026). Modelling synaptic dysfunction in childhood dementia using human iPSC-derived cortical networks. Nature communications, 17(1), 3161. https://doi.org/10.1038/s41467-026-71112-9

BibTeX

@article{mazzachi2026modelling,
author = {Mazzachi, Paris and McDonald, Ella and Greenberg, Zarina and Puerta, Alejandra Noreña and Tran, Jenne and De Silva, Manam Inushi and Christensen, Cade and Adams, Robert and Loskarn, Sebastian and Beard, Helen and Zabolocki, Michael and Elmasri, Meera and Maack, Megan and Elvidge, Kristina L. and Hutchinson, Mark R. and O’Neill, Cara and Hemsley, Kim M. and Melton, Lisa and Smith, Nicholas and Bardy, Cedric},
title = {{Modelling synaptic dysfunction in childhood dementia using human iPSC-derived cortical networks}},
journal = {Nature communications},
year = {2026},
month = apr,
volume = {17},
number = {1},
pages = {3161},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-71112-9},
url = {https://doi.org/10.1038/s41467-026-71112-9},
pmid = {41946741},
pmcid = {PMC13057175}
}

RIS

TY - JOUR
AU - Mazzachi, Paris
AU - McDonald, Ella
AU - Greenberg, Zarina
AU - Puerta, Alejandra Noreña
AU - Tran, Jenne
AU - De Silva, Manam Inushi
AU - Christensen, Cade
AU - Adams, Robert
AU - Loskarn, Sebastian
AU - Beard, Helen
AU - Zabolocki, Michael
AU - Elmasri, Meera
AU - Maack, Megan
AU - Elvidge, Kristina L.
AU - Hutchinson, Mark R.
AU - O’Neill, Cara
AU - Hemsley, Kim M.
AU - Melton, Lisa
AU - Smith, Nicholas
AU - Bardy, Cedric
TI - Modelling synaptic dysfunction in childhood dementia using human iPSC-derived cortical networks
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/04/07
VL - 17
IS - 1
SP - 3161
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-71112-9
UR - https://doi.org/10.1038/s41467-026-71112-9
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-71112-9",
"type": "article-journal",
"title": "Modelling synaptic dysfunction in childhood dementia using human iPSC-derived cortical networks",
"container-title": "Nature communications",
"author": [
{
"family": "Mazzachi",
"given": "Paris"
},
{
"family": "McDonald",
"given": "Ella"
},
{
"family": "Greenberg",
"given": "Zarina"
},
{
"family": "Puerta",
"given": "Alejandra Noreña"
},
{
"family": "Tran",
"given": "Jenne"
},
{
"family": "De Silva",
"given": "Manam Inushi"
},
{
"family": "Christensen",
"given": "Cade"
},
{
"family": "Adams",
"given": "Robert"
},
{
"family": "Loskarn",
"given": "Sebastian"
},
{
"family": "Beard",
"given": "Helen"
},
{
"family": "Zabolocki",
"given": "Michael"
},
{
"family": "Elmasri",
"given": "Meera"
},
{
"family": "Maack",
"given": "Megan"
},
{
"family": "Elvidge",
"given": "Kristina L."
},
{
"family": "Hutchinson",
"given": "Mark R."
},
{
"family": "O’Neill",
"given": "Cara"
},
{
"family": "Hemsley",
"given": "Kim M."
},
{
"family": "Melton",
"given": "Lisa"
},
{
"family": "Smith",
"given": "Nicholas"
},
{
"family": "Bardy",
"given": "Cedric"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "3161",
"DOI": "10.1038/s41467-026-71112-9",
"PMID": "41946741",
"PMCID": "PMC13057175",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-71112-9",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
7
]
]
}
}

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

Similar papers

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

[1] doi:10.1038/s41467-026-76837-1 [code]
Drug screen and machine learning predict neuroprotective agents in a preclinical human model of childhood dementia.
Journal: Nature communications
In common: Harmony, clusterProfiler, Seurat, 3 other tools, Alzheimer's / dementia, 34 references, 13 authors
[2] doi:10.1093/bib/bbag175 [code]
Uncovering causal relationships in single-cell omic studies with causarray.
Journal: Briefings in bioinformatics
In common: Harmony, clusterProfiler, Seurat, 4 other tools, Alzheimer's / dementia, 2 references
[3] doi:10.1038/s41467-026-71595-6 [code]
A single-cell and spatial atlas of early human olfactory development.
Journal: Nature communications
In common: Harmony, Seurat, patchwork, 2 other tools, 4 references
[4] doi:10.1038/s41593-026-02367-0 [code]
A reproducible three-dimensional model of human brain tissue to investigate physiological and disease-associated microglia phenotypes.
Journal: Nature neuroscience
In common: Harmony, clusterProfiler, Seurat, 4 other tools, Alzheimer's / dementia, cellular / molecular, 2 references
[5] doi:10.1038/s41586-026-10877-x [code]
Human brain organoids record the passage of time over multiple years.
Journal: Nature
In common: Harmony, Seurat, data.table, 2 other tools, 4 references
[6] doi:10.1002/alz.71823 [code]
Cellular transcriptomic signatures underpinning the heterogeneity of depression in Alzheimer's disease.
Journal: Alzheimer's & dementia : the journal of the Alzheimer's Association
In common: Harmony, Seurat, data.table, 3 other tools, Alzheimer's / dementia, cellular / molecular, 2 references
[7] doi:10.1038/s41586-026-10214-2 [code]
Multidimensional profiling of heterogeneity in supratentorial ependymomas.
Journal: Nature
In common: Harmony, clusterProfiler, Seurat, 4 other tools, 2 references
[8] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: Harmony, clusterProfiler, Seurat, 4 other tools, cellular / molecular, 1 reference
[9] doi:10.1038/s41467-026-69944-6 [code]
Multi-modal dissection of cell-type specific TDP-43 pathology in the motor cortex.
Journal: Nature communications
In common: clusterProfiler, Seurat, data.table, 3 other tools, 2 references
[10] doi:10.1016/j.cell.2026.05.026 [code]
The critical role of the endogenous immune compartment after CAR T cell therapy in recurrent GBM.
Journal: Cell
In common: Harmony, clusterProfiler, Seurat, 4 other tools, 1 reference

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.