OSCR

Ischemic stroke triggers brain-wide synaptic remodeling within four hours.

Code ↔ Paper

1 match 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 1 match
  1. [1] § Methods › Quantification and statistical analysis › Statistic analysis of proteomic and transcriptomic data. ↔ Dream_Synaptosome_DE_code.R, lines 43–91 · score 0.76 · linear mixed models, Rattus norvegicus, Salmon, Dream, variancePartition, gene

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 · 1,237 lines · 53 KB · CC-BY-4.0 · 1 match

  1. library(tximeta)
  2. library(DESeq2)
  3. library(tidyverse)
  4. library(readr)
  5. library(apeglm)
  6. library(ashr)
  7. library (EnhancedVolcano)
  8. library(gprofiler2)
  9. library(dplyr)
  10. library(IHW)
  11. library("variancePartition")
  12. library(edgeR)
  13. library(tximport)
  14. library(patchwork)
  15. library(cowplot)
  17. ################################################################################
  18. #################################################################### Synaptosome
  19. pathSyn <- "/BMK_DATA_20230515093157_1_synaptosome/transcripts_quant"
  21. directoryIDSyn <- c("BE096-02T0001_good_quant", "BE096-02T0002_good_quant", "BE096-02T0003_good_quant",
  22. "BE096-02T0004_good_quant", "BE096-02T0005_good_quant", "BE096-02T0006_good_quant",
  23. "BE096-02T0007_good_quant", "BE096-02T0008_good_quant", "BE096-02T0009_good_quant",
  24. "BE096-02T0010_good_quant", "BE096-02T0011_good_quant", "BE096-02T0012_good_quant")
  26. filesSyn <- file.path(pathSyn, directoryIDSyn, "quant.sf")
  28. files <- c(filesSyn)
  29. file.exists(files)
  30. rep <- c("CA1", "CA2", "CA3", "4z1", "4z2", "4z3", "4c1", "4c2", "4c3", "4p1", "4p2", "4p3")
  31. region <- c("CA", "CA", "CA", "4z", "4z", "4z", "4c", "4c", "4c", "4p", "4p", "4p")
  32. data <- c("synaptosome","synaptosome","synaptosome","synaptosome","synaptosome","synaptosome","synaptosome","synaptosome","synaptosome","synaptosome","synaptosome","synaptosome")
  33. rat <- c("1", "2", "3", "4", "5", "6", "4", "5", "6", "4", "5", "6")
  34. condition <- c("control", "control", "control", "stroke", "stroke","stroke","stroke", "stroke", "stroke","stroke","stroke", "stroke")
  36. coldata_syn <- data.frame(files, data=rep, names=rep, region=region,rat = rat, condition=condition, stringsAsFactors=FALSE)
  37. coldata_syn$region <- as.factor(coldata_syn$region)
  38. coldata_syn$data <- as.factor(coldata_syn$data)
  39. coldata_syn$names <- as.factor(coldata_syn$names)
  40. coldata_syn$rat <- as.factor(coldata_syn$rat)
  41. coldata_syn$condition <- as.factor(coldata_syn$condition)
  43. # Build taxome file
  44. indexDir <- file.path("/Resources/Salmon_indexes/Rattus_norvegicus/salmon_index")
  45. fastaFTP <- c("/Resources/Salmon_indexes/Rattus_norvegicus")
  46. gtfPath <- file.path("/Resources/Rattus_norvegicus.mRatBN7.2.111.gtf.gz")
  48. makeLinkedTxome(indexDir=indexDir,
  49. source="LocalEnsembl",
  50. organism="Rattus norvegicus",
  51. release="111",
  52. genome="mRatBN7.2",
  53. fasta=fastaFTP,
  54. gtf=gtfPath,
  55. write=FALSE)
  56. # Syn
  57. se_syn <- tximeta(coldata_syn)
  58. edb <- retrieveDb(se_syn)
  59. k <- keys(edb, keytype = "TXNAME")
  60. tx2gene <- AnnotationDbi::select(edb, k, "GENEID", "TXNAME")
  62. txi <- tximport(files, type = "salmon", tx2gene = tx2gene, countsFromAbundance = "lengthScaledTPM", ignoreTxVersion = T)
  65. y <- DGEList(txi$counts)
  66. design <- model.matrix(~region, data = coldata_syn)
  68. isexpr <- filterByExpr(y, design)
  69. y <- y[isexpr, ]
  70. y <- calcNormFactors(y)
  72. param <- SnowParam(4, "SOCK", progressbar = TRUE)
  73. form <- ~ region
  75. # estimate weights using linear mixed model of dream
  76. vobjDream <- voomWithDreamWeights(y, form, coldata_syn, BPPARAM = param)
  77. write_tsv(as.data.frame(vobjDream$E) %>%
  78. mutate(gene = rownames(vobjDream$E)) %>%
  79. rename( CA1 = Sample1, CA2 = Sample2, CA3 = Sample3, "4z1" = Sample4,
  80. "4z2" = Sample5, "4z3" = Sample6, "4c1" = Sample7, "4c2" = Sample8,
  81. "4c3" = Sample9, "4p1" = Sample10, "4p2" = Sample11, "4p3" = Sample12),
  82. file = "q/VariancePartitioning/Syn_voomWithDreamWeights-normalised.txt")
  84. fitmm <- dream(vobjDream, form, coldata_syn, ddf = "Kenward-Roger")
  85. fitmm <- eBayes(fitmm)
  87. L <- makeContrastsDream(form, coldata_syn,
  88. contrasts = c(
  89. compare4c_4z = "region4c - region4z",
  90. compare4c_4p = "region4c - region4p",
  91. compare4z_4p = "region4z - region4p"))
  93. # fit dream model with contrasts
  94. fit <- dream(vobjDream, form, coldata_syn, L)
  95. fit <- eBayes(fit)
  97. compare4c_4z = topTable(fit, coef = "compare4c_4z", number = 47000)
  98. compare4c_4z$input = rownames(compare4c_4z)
  99. compare4c_4zGenes = gprofiler2::gconvert(rownames(compare4c_4z), organism = "rnorvegicus")
  100. compare4c_4z = left_join(compare4c_4z, compare4c_4zGenes, by = "input")
  102. compare4c_4p = topTable(fit, coef = "compare4c_4p", number = 47000)
  103. compare4c_4p$input = rownames(compare4c_4p)
  104. compare4c_4pGenes = gprofiler2::gconvert(rownames(compare4c_4p), organism = "rnorvegicus")
  105. compare4c_4p = left_join(compare4c_4p, compare4c_4pGenes, by = "input")
  107. compare4z_4p = topTable(fit, coef = "compare4z_4p", number = 47000)
  108. compare4z_4p$input = rownames(compare4z_4p)
  109. compare4z_4pGenes = gprofiler2::gconvert(rownames(compare4z_4p), organism = "rnorvegicus")
  110. compare4z_4p = left_join(compare4z_4p, compare4z_4pGenes, by = "input")
  112. volcano_res1 = EnhancedVolcano(compare4z_4p,
  113. lab = compare4z_4p$name,
  114. x = 'logFC',
  115. y = 'adj.P.Val',
  116. pCutoff = 0.05,
  117. FCcutoff = 1,
  118. xlim = c(-11, 11),
  119. ylim = c(0, 6),
  120. pointSize = 1.5,
  121. labSize = 3,
  122. title = '4z vs 4p in Syn',
  123. # subtitle = 'Differential expression',
  124. # caption = 'FC cutoff, 1; p-value cutoff, 10e-4',
  125. legendPosition = "right",
  126. legendLabSize = 14,
  127. col = c('grey30', '#E9E9E9', "#0072B2", '#D55E00'),
  128. colAlpha = 0.5,
  129. #drawConnectors = TRUE,
  130. hline = c(10e-8),
  131. widthConnectors = 0.5)
  132. volcano_res2 = EnhancedVolcano(compare4c_4z,
  133. lab = compare4c_4z$name,
  134. x = 'logFC',
  135. y = 'adj.P.Val',
  136. pCutoff = 0.05,
  137. FCcutoff = 1,
  138. xlim = c(-11, 11),
  139. ylim = c(0, 6),
  140. pointSize = 1.5,
  141. labSize = 3,
  142. title = '4c vs 4z in Syn',
  143. # subtitle = 'Differential expression',
  144. # caption = 'FC cutoff, 1; p-value cutoff, 10e-4',
  145. legendPosition = "right",
  146. legendLabSize = 14,
  147. col = c('grey30', '#E9E9E9', "#0072B2", '#D55E00'),
  148. colAlpha = 0.5,
  149. #drawConnectors = TRUE,
  150. hline = c(10e-8),
  151. widthConnectors = 0.5)
  152. volcano_res3 = EnhancedVolcano(compare4c_4p,
  153. lab = compare4c_4p$name,
  154. x = 'logFC',
  155. y = 'adj.P.Val',
  156. pCutoff = 0.05,
  157. FCcutoff = 1,
  158. xlim = c(-11, 11),
  159. ylim = c(0, 6),
  160. pointSize = 1.5,
  161. labSize = 3,
  162. title = '4c vs 4p in Syn',
  163. # subtitle = 'Differential expression',
  164. # caption = 'FC cutoff, 1; p-value cutoff, 10e-4',
  165. legendPosition = "right",
  166. legendLabSize = 14,
  167. col = c('grey30', '#E9E9E9', "#0072B2", '#D55E00'),
  168. colAlpha = 0.5,
  169. #drawConnectors = TRUE,
  170. hline = c(10e-8),
  171. widthConnectors = 0.5)
  173. down_genes_compare4z_4p <- compare4z_4p %>% filter(adj.P.Val <= 0.05 & (logFC < (-1)))
  174. up_genes_compare4z_4p <- compare4z_4p %>% filter(adj.P.Val <= 0.05 & (logFC > 1))
  176. gost_down_genes_compare4z_4p = gprofiler2::gost(down_genes_compare4z_4p$name, organism = "rnorvegicus")
  177. gost_up_genes_compare4z_4p = gprofiler2::gost(up_genes_compare4z_4p$name, organism = "rnorvegicus")
  179. down_genes_compare4c_4z <- compare4c_4z %>% filter(adj.P.Val <= 0.05 & (logFC < (-1)))
  180. up_genes_compare4c_4z <- compare4c_4z %>% filter(adj.P.Val <= 0.05 & (logFC > 1))
  182. gost_down_genes_compare4c_4z = gprofiler2::gost(down_genes_compare4c_4z$name, organism = "rnorvegicus")
  183. gost_up_genes_compare4c_4z = gprofiler2::gost(up_genes_compare4c_4z$name, organism = "rnorvegicus")
  185. down_genes_compare4c_4p <- compare4c_4p %>% filter(adj.P.Val <= 0.05 & (logFC < (-1)))
  186. up_genes_compare4c_4p <- compare4c_4p %>% filter(adj.P.Val <= 0.05 & (logFC > 1))
  188. gost_down_genes_compare4c_4p = gprofiler2::gost(down_genes_compare4c_4p$name, organism = "rnorvegicus")
  189. gost_up_genes_compare4c_4p = gprofiler2::gost(up_genes_compare4c_4p$name, organism = "rnorvegicus")
  192. ############ Each versus CA
  193. compare4c_CA = topTable(fit, coef = "region4c", number = 47000)
  194. compare4c_CA$input = rownames(compare4c_CA)
  195. compare4c_CAGenes = gprofiler2::gconvert(rownames(compare4c_CA), organism = "rnorvegicus")
  196. compare4c_CA = left_join(compare4c_CA, compare4c_CAGenes, by = "input")
  198. compare4p_CA = topTable(fit, coef = "region4p", number = 47000)
  199. compare4p_CA$input = rownames(compare4p_CA)
  200. compare4p_CAGenes <- gprofiler2::gconvert(rownames(compare4p_CA), organism = "rnorvegicus")
  201. compare4p_CA = left_join(compare4p_CA, compare4p_CAGenes, by = "input")
  203. compare4z_CA = topTable(fit, coef = "region4z", number = 47000)
  204. compare4z_CA$input = rownames(compare4z_CA)
  205. compare4z_CAGenes = gprofiler2::gconvert(rownames(compare4z_CA), organism = "rnorvegicus")
  206. compare4z_CA = left_join(compare4z_CA, compare4z_CAGenes, by = "input")
  208. volcano_res4 = EnhancedVolcano(compare4z_CA,
  209. lab = compare4z_CA$name,
  210. x = 'logFC',
  211. y = 'adj.P.Val',
  212. pCutoff = 0.05,
  213. FCcutoff = 1,
  214. xlim = c(-11, 11),
  215. ylim = c(0, 6),
  216. pointSize = 1.5,
  217. labSize = 3,
  218. title = '4z vs CA in Synaptosome',
  219. # subtitle = 'Differential expression',
  220. # caption = 'FC cutoff, 1; p-value cutoff, 10e-4',
  221. legendPosition = "right",
  222. legendLabSize = 14,
  223. col = c('grey30', '#E9E9E9', "#0072B2", '#D55E00'),
  224. colAlpha = 0.5,
  225. #drawConnectors = TRUE,
  226. hline = c(10e-8),
  227. widthConnectors = 0.5)
  228. volcano_res5 = EnhancedVolcano(compare4c_CA,
  229. lab = compare4c_CA$name,
  230. x = 'logFC',
  231. y = 'adj.P.Val',
  232. pCutoff = 0.05,
  233. FCcutoff = 1,
  234. xlim = c(-11, 11),
  235. ylim = c(0, 6),
  236. pointSize = 1.5,
  237. labSize = 3,
  238. title = '4c vs CA in Synaptosome',
  239. # subtitle = 'Differential expression',
  240. # caption = 'FC cutoff, 1; p-value cutoff, 10e-4',
  241. legendPosition = "right",
  242. legendLabSize = 14,
  243. col = c('grey30', '#E9E9E9', "#0072B2", '#D55E00'),
  244. colAlpha = 0.5,
  245. #drawConnectors = TRUE,
  246. hline = c(10e-8),
  247. widthConnectors = 0.5)
  248. volcano_res6 = EnhancedVolcano(compare4p_CA,
  249. lab = compare4p_CA$name,
  250. x = 'logFC',
  251. y = 'adj.P.Val',
  252. pCutoff = 0.05,
  253. FCcutoff = 1,
  254. xlim = c(-11, 11),
  255. ylim = c(0, 6),
  256. pointSize = 1.5,
  257. labSize = 3,
  258. title = '4p vs CA in Synaptosome',
  259. # subtitle = 'Differential expression',
  260. # caption = 'FC cutoff, 1; p-value cutoff, 10e-4',
  261. legendPosition = "right",
  262. legendLabSize = 14,
  263. col = c('grey30', '#E9E9E9', "#0072B2", '#D55E00'),
  264. colAlpha = 0.5,
  265. #drawConnectors = TRUE,
  266. hline = c(10e-8),
  267. widthConnectors = 0.5)
  269. # Filter significant genes in 4z vs CA
  270. up_genes_compare4z_CA <- compare4z_CA %>% filter(adj.P.Val <= 0.05 & (logFC > 1))
  271. down_genes_compare4z_CA <- compare4z_CA %>% filter(adj.P.Val <= 0.05 & (logFC < (-1)))
  273. gost_up_genes_compare4z_CA = gprofiler2::gost(up_genes_compare4z_CA$name, organism = "rnorvegicus", evcodes = T)
  274. gost_down_genes_compare4z_CA = gprofiler2::gost(down_genes_compare4z_CA$name, organism = "rnorvegicus", evcodes = T)
  275. gost_up_genes_compare4z_CA = gprofiler2::gost(up_genes_compare4z_CA$name, organism = "rnorvegicus")
  277. # Filter significant genes in 4c vs CA
  278. up_genes_compare4c_CA <- compare4c_CA %>% filter(adj.P.Val <= 0.05 & (logFC > 1))
  279. down_genes_compare4c_CA <- compare4c_CA %>% filter(adj.P.Val <= 0.05 & (logFC < (-1)))
  282. gost_up_genes_compare4c_CA = gprofiler2::gost(up_genes_compare4c_CA$name, organism = "rnorvegicus")
  283. gost_down_genes_compare4c_CA = gprofiler2::gost(down_genes_compare4c_CA$name, organism = "rnorvegicus")
  285. x = as.data.frame(up_genes_compare4c_CA$name)
  287. compare4c_CA_gost = gostplot(gost_down_genes_compare4c_CA, interactive = F)
  288. compare4c_CA_publ = publish_gostplot(compare4c_CA_gost)
  289. compare4c_CA_table = publish_gosttable(gost_down_genes_compare4c_CA, highlight_terms = gost_down_genes_compare4c_CA$result[c(1:20),],
  290. use_colors = TRUE,
  291. show_columns = c("source", "term_name", "term_size", "intersection_size"),
  292. filename = NULL)
  294. gost_up_genes_compare4c_CA = gprofiler2::gost(up_genes_compare4c_CA$name, organism = "rnorvegicus")
  296. compare4c_CA_up_gost = gostplot(gost_up_genes_compare4c_CA, interactive = F)
  297. compare4c_CA_up_publ = publish_gostplot(compare4c_CA_up_gost)
  298. compare4c_CA_up_table = publish_gosttable(gost_up_genes_compare4c_CA, highlight_terms = gost_up_genes_compare4c_CA$result[c(1:20),],
  299. use_colors = TRUE,
  300. show_columns = c("source", "term_name", "term_size", "intersection_size"),
  301. filename = NULL) +
  302. labs(title = "Upregulated genes 4c versus CA")
  305. # Filter significant genes in 4p vs CA
  306. up_genes_compare4p_CA <- compare4p_CA %>% filter(adj.P.Val <= 0.05 & (logFC > 1))
  307. down_genes_compare4p_CA <- compare4p_CA %>% filter(adj.P.Val <= 0.05 & (logFC < (-1)))
  310. gost_up_genes_compare4p_CA = gprofiler2::gost(up_genes_compare4p_CA$name, organism = "rnorvegicus", evcodes = T)
  311. gost_down_genes_compare4p_CA = gprofiler2::gost(down_genes_compare4p_CA$name, organism = "rnorvegicus", evcodes = T)
  313. x = as.data.frame(up_genes_compare4p_CA$name)
  315. compare4p_CA_gost = gostplot(gost_down_genes_compare4p_CA, interactive = F)
  316. compare4p_CA_publ = publish_gostplot(compare4p_CA_gost)
  317. compare4p_CA_table = publish_gosttable(gost_down_genes_compare4p_CA, highlight_terms = gost_down_genes_compare4p_CA$result[c(1:20),],
  318. use_colors = TRUE,
  319. show_columns = c("source", "term_name", "term_size", "intersection_size"),
  320. filename = NULL)
  322. gost_up_genes_compare4p_CA = gprofiler2::gost(up_genes_compare4p_CA$name, organism = "rnorvegicus")
  324. compare4p_CA_up_gost = gostplot(gost_up_genes_compare4p_CA, interactive = F)
  325. compare4p_CA_up_publ = publish_gostplot(compare4p_CA_up_gost)
  326. compare4p_CA_up_table = publish_gosttable(gost_up_genes_compare4p_CA, highlight_terms = gost_up_genes_compare4p_CA$result[c(1:20),],
  327. use_colors = TRUE,
  328. show_columns = c("source", "term_name", "term_size", "intersection_size"),
  329. filename = NULL) +
  330. labs(title = "Upregulated genes 4p versus CA")
  332. ################################################################################
  333. ########## Plot to file
  335. #file_path <- "/VariancePartitioning/VariancePartition_Syn.pdf"
  337. pdf(file_path, width = 14, height = 12)
  339. plot(pca$x, cex = 4.0, col = c("#2f4f4f",
  340. "#228b22",
  341. "#7f0000",
  342. "#00008b"), pch = 18)
  343. legend(
  344. "topright",
  345. bty = "n",
  346. c("CA", "4z", "4c","4p"),
  347. fill = c("#2f4f4f",
  348. "#228b22",
  349. "#7f0000",
  350. "#00008b"),
  351. cex = 3.0)
  352. plot(tree)
  353. plotVarPart(vp)
  355. plotContrasts(L)
  357. plotPercentBars(vp[1:100, ])
  358. heatmap(syn_matrix)
  362. volcano_res1
  363. compare4z_4p_publ
  364. cat("Downreg in 4z versus 4p")
  365. compare4z_4p_table
  367. volcano_res2
  368. compare4c_4z_publ
  369. cat("Downreg in 4c versus 4z")
  370. compare4c_4z_table
  373. volcano_res3
  374. compare4c_4p_publ
  375. cat("Downreg in 4c versus 4p")
  376. compare4c_4p_table
  379. volcano_res4
  381. compare4z_CA_publ
  382. cat("Downreg in 4z versus CA")
  383. compare4z_CA_table
  384. cat("Upreg in 4z versus CA")
  385. compare4z_CA_up_publ
  386. compare4z_CA_up_table
  389. volcano_res5
  391. compare4c_CA_publ
  392. cat("Downreg in 4c versus CA")
  393. compare4c_CA_table
  394. cat("Upreg in 4c versus CA")
  395. compare4c_CA_up_publ
  396. compare4c_CA_up_table
  398. volcano_res6
  400. compare4p_CA_publ
  401. cat("Downreg in 4p versus CA")
  402. compare4p_CA_table
  403. cat("Upreg in 4p versus CA")
  404. compare4p_CA_up_publ
  405. compare4p_CA_up_table
  407. dev.off()
  409. ################################################################################
  410. ###################### Gene set enrichment analysis ############################
  411. library(GSVA)
  412. library(fgsea)
  413. library(GSEABase)
  414. library(readxl)
  415. library(dendextend)
  416. library(gplots)
  418. # Path to your Excel file with genes of interest
  419. xlsx_file <- "/Other_resources/mmc2.xlsx"
  420. # Add novel genes (Synapsed enriched novel)
  421. xlsx_file2 <- "/Other_resources/mmc2.xlsx"
  422. xlsx_file3 <- "/Other_resources/GeneSets_Synaptic_Presynaptic.xlsx"
  424. # Read the first/second sheet of the Excel file, and convert names
  425. # to match the expression matrix
  426. df = read_excel(xlsx_file, sheet = 1) %>% dplyr::select(GeneName, SynapseType)
  427. df2 = read_excel(xlsx_file, sheet = 2) %>% dplyr::select("hgnc symbol") %>%
  428. mutate(SynapseType = "Oostrum2023_novel") %>% rename(GeneName = "hgnc symbol")
  430. df3 = read_excel(xlsx_file3, sheet = 1) %>% pivot_longer(cols = c(Presynaptic, Postsynaptic),
  431. names_to = "SynapseType", values_to = "GeneName")
  433. df = rbind(df, df2, df3)
  435. df_ens = gprofiler2::gconvert(df$GeneName, organism = "rnorvegicus")
  436. df = left_join(df, df_ens, by = c("GeneName"="input")) %>% dplyr::select(target, SynapseType) %>% unique()
  438. # Define gene sets
  439. gene_set <- list()
  441. # Loop through unique entries in SynapseType
  442. unique_entries <- unique(df$SynapseType)
  443. for(entry in unique_entries) {
  444. # Extract values from GeneName corresponding to current SynapseType
  445. values <- df$target[df$SynapseType == entry]
  446. # Assign values to the list with name as current SynapseType
  447. gene_set[[entry]] <- values
  448. }
  449. # Convert to GeneSetCollection
  450. gene_set_list <- lapply(names(gene_set), function(x) GeneSet(setName = x, geneIds = gene_set[[x]]))
  452. ################################################################################
  453. gsc <- GeneSetCollection(gene_set_list)
  454. ################################################################################
  456. # Perform GSVA on voomWithDream normaalized counts (log2-counts per million (logCPM))
  457. colnames(vobjDream$E) = rep
  459. vobjDream_E = as.data.frame(vobjDream$E)
  460. vobjDream_E$genes = rownames(vobjDream_E)
  461. converted = gprofiler2::gconvert(rownames(vobjDream_E), organism = "rnorvegicus")
  462. vobjDream_E = left_join(vobjDream_E, converted, by = c("genes" = "input"))
  463. vobjDream_E = vobjDream_E[,c("name","CA1","CA2","CA3","4z1","4z2","4z3","4c1",
  464. "4c2","4c3","4p1","4p2","4p3")]
  465. vobjDream_E2 = as.matrix(vobjDream_E[,2:13])
  466. rownames(vobjDream_E2) = vobjDream_E$name
  467. vobjDream_E = as.data.frame(vobjDream$E)
  468. ################################################################################
  469. unique_genes_Oostrum23 <- unique(unlist(gene_set))
  470. unique_genes_Oostrum23_symbols = gprofiler2::gconvert(unique_genes_Oostrum23, organism = "rnorvegicus")
  471. Oostrum23_geneExpr = vobjDream_E2[rownames(vobjDream_E2) %in% unique_genes_Oostrum23_symbols$name,]
  472. ################################################################################
  473. gsva_res <- gsva(vobjDream$E, gsc, method = "gsva", kcdf = "Poisson", mx.diff = TRUE, verbose = FALSE)
  475. # Create design matrix
  476. design <- model.matrix(~ region, data = coldata_syn)
  478. # fit dream model with contrasts
  479. fit2 <- dream(gsva_res, form, coldata_syn, L)
  480. fit2 <- eBayes(fit2)
  482. # Get differential analysis results
  483. results4p_CA <- topTable(fit2, coef="region4p", adjust.method="BH", number=Inf)
  484. results4c_CA <- topTable(fit2, coef="region4c", adjust.method="BH", number=Inf)
  485. results4z_CA <- topTable(fit2, coef="region4z", adjust.method="BH", number=Inf)
  487. pca2 = prcomp(gsva_res)
  488. pca3 = prcomp(t(gsva_res))
  489. pca4 = prcomp(t(Oostrum23_geneExpr))
  491. ################################################################################
  492. ########## Correlation between terms
  494. # Filter gene sets to include only genes present in the gene expression matrix
  495. gene_sets <- lapply(gene_set, function(genes) {
  496. genes[genes %in% rownames(vobjDream_E)]
  497. })
  499. # Aggregate expression values for each gene set
  500. aggregate_expression <- sapply(gene_sets, function(genes) {
  501. if (length(genes) > 0) {
  502. colMeans(vobjDream_E[genes, , drop = FALSE])
  503. } else {
  504. rep(NA, ncol(vobjDream_E)) # Handle gene sets with no genes present in the data
  505. }
  506. })
  508. # Remove any columns with NA values (gene sets with no genes present in the data)
  509. aggregate_expression <- aggregate_expression[, colSums(is.na(aggregate_expression)) == 0]
  511. # Compute the correlation matrix for the aggregated expression values of the gene sets
  512. correlation_matrix <- cor(aggregate_expression)
  514. # Print the correlation matrix
  515. print(correlation_matrix)
  517. ################################################################################
  518. ########## Plot to file
  519. file_path <- "/GeneSet_Enrichments/Heaatmap_GeneSetOostrum2023_and_fromOleg_Syn.pdf"
  521. pdf(file_path, width = 18, height = 12)
  522. library(RColorBrewer)
  524. # Define color palette for heatmap
  525. colorLegend <- c("#0072B2", "#D55E00", "#CC79A7","#F0E442")
  527. #0072B2 (blue)
  528. #CC79A7 (pinky purple)
  529. #D55E00 (orange)
  530. #F0E442 (yellow)
  532. names(colorLegend) <- c("CA", "4z", "4c", "4p")
  533. sample.color.map <- colorLegend[coldata_syn$region]
  534. coldata_syn$rep <- rep
  535. names(sample.color.map) <- coldata_syn$rep
  536. sample.color.map["4z1"] = "#D55E00"
  537. sample.color.map["4z2"] = "#D55E00"
  538. sample.color.map["4z3"] = "#D55E00"
  540. sample.color.map["4p1"] = "#F0E442"
  541. sample.color.map["4p2"] = "#F0E442"
  542. sample.color.map["4p3"] = "#F0E442"
  544. sample.color.map["4c1"] = "#CC79A7"
  545. sample.color.map["4c2"] = "#CC79A7"
  546. sample.color.map["4c3"] = "#CC79A7"
  548. heatmap_colors <- colorRampPalette(c("lightblue", "white", "#FF457D"))(n = 1000)
  550. # Cluster the samples and gene terms
  551. geneSetClustering <- hclust(dist(gsva_res, method = "euclidean"), method = "complete")
  552. sampleClustering <- hclust(dist(t(gsva_res), method = "euclidean"), method = "complete")
  554. # Convert hclust objects to dendrograms
  555. geneSetDendrogram <- as.dendrogram(geneSetClustering)
  556. sampleDendrogram <- as.dendrogram(sampleClustering)
  558. sampleDendrogram <- rotate(sampleDendrogram, order = c("4p3", "4p1", "4p2",
  559. "4c2", "4c1", "4c3",
  560. "CA1", "CA2", "CA3",
  561. "4z2", "4z3", "4z1"))
  564. # Function to rotate a specific branch of a dendrogram given a path
  565. rotate_branch <- function(dend, path) {
  566. if (length(path) == 0) return(rev(dend))
  567. if (is.leaf(dend)) return(dend)
  568. branch_index <- path[1]
  569. dend[[branch_index]] <- rotate_branch(dend[[branch_index]], path[-1])
  570. dend
  571. }
  573. # Path to the [[2]][[2]][[1]] branch
  574. branch_path <- c(2, 2, 1)
  576. # Sample dendrogram (for demonstration purposes)
  577. # Create a sample dendrogram if needed
  578. # sampleDendrogram <- as.dendrogram(hclust(dist(USArrests), "ave"))
  580. # Rotate the specified branch
  581. sampleDendrogram <- rotate_branch(sampleDendrogram, branch_path)
  583. # Create heatmap with adjusted parameters
  584. heatmap(as.matrix(gsva_res),
  585. ColSideColors = sample.color.map,
  586. xlab = " ",
  587. ylab = " ",
  588. margins = c(5, 5),
  589. labRow = gsub("_", " ", rownames(gsva_res)), # Remove only "_"
  590. labCol = colnames(gsva_res),
  591. scale = "row",
  592. Colv = sampleDendrogram,
  593. Rowv = geneSetDendrogram,
  594. col = heatmap_colors,
  595. cexRow = 0.8, # 2 for gene set enrichment but 0.8 for cell type enrichment
  596. cexCol = 2,
  597. cex.lab = 2)
  599. # Add legend for heatmap colors
  600. min_value <- as.numeric(min(gsva_res))
  601. max_value <- as.numeric(max(gsva_res))
  602. mid_value <- (min_value + max_value) / 2
  603. min_value <- round(min_value, 3)
  604. max_value <- round(max_value, 3)
  605. mid_value <- round(mid_value, 3)
  607. # Create legend labels
  608. legend_labels <- c(min_value, mid_value, max_value)
  610. # Add legend for regions/reps
  611. legend("topleft", legend = names(colorLegend), fill = colorLegend, title = "Sections", cex = 1.5, bty = "n", inset = 0.05)
  613. # Legend plot for values
  614. legend("left", legend = legend_labels, fill = colorRampPalette(c("lightblue", "white", "#FF457D"))(n = 3), title = "Enrichment values", cex = 1.5, bty = "n", inset = 0.05)
  616. ################################################################################
  617. ########################## Samples PCA based on gene set/term enrichment scoring
  619. pca_data <- as.data.frame(pca3$x)
  620. pca_data$sample <- rownames(pca_data)
  623. # Assign colors to samples
  624. pca_data$color <- c("#0072B2","#0072B2","#0072B2", # blue
  625. "#D55E00","#D55E00","#D55E00", #orange
  626. "#CC79A7","#CC79A7","#CC79A7", #pinky purple
  627. "#F0E442","#F0E442","#F0E442") #yellow
  628. custom_colors <- c("#0072B2", "#CC79A7","#D55E00", "#F0E442")
  629. # Plot the PCA
  630. ggplot(pca_data, aes(x = PC1, y = PC2, label = sample, color = color)) +
  631. geom_point(size = 5) + # Adjust point size
  632. geom_text(vjust = 1.5, hjust = 1.5, size = 5) + # Adjust font size
  633. labs(title = "PCA of Gene Set Enrichment Scores (Samples)",
  634. x = "Principal Component 1",
  635. y = "Principal Component 2") +
  636. scale_color_manual(values = custom_colors) +
  637. theme_minimal() +
  638. theme(aspect.ratio = 1, # Make the plot square
  639. plot.title = element_text(size = 16), # Adjust title font size
  640. axis.title = element_text(size = 14), # Adjust axis titles font size
  641. axis.text = element_text(size = 12), # Adjust axis text font size
  642. legend.position = "none") # Remove legend
  644. ################################################################################
  645. ##################### Terms
  647. pca_data2 <- as.data.frame(pca2$x)
  648. pca_data2$sample <- rownames(pca_data2)
  650. base_colors <- brewer.pal(8, "Set1") # Use "Set1" for a vibrant color set
  651. custom_colors2 <- colorRampPalette(base_colors)(18)
  653. # Create a named vector for custom colors
  654. names(custom_colors2) <- rownames(pca_data2)
  657. # Plot the PCA
  658. ggplot(pca_data2, aes(x = PC1, y = PC2, label = sample, color = sample)) +
  659. geom_point(size = 5) + # Adjust point size
  660. geom_text_repel(aes(label = sample), size = 5, nudge_x = 0.05, direction = "y", hjust = 0) + # Adjust text label position
  661. labs(title = "PCA of Gene Set Enrichment Scores",
  662. x = "Principal Component 1",
  663. y = "Principal Component 2") +
  664. scale_color_manual(values = custom_colors2) +
  665. theme_minimal() +
  666. theme(aspect.ratio = 1, # Make the plot square
  667. plot.title = element_text(size = 16), # Adjust title font size
  668. axis.title = element_text(size = 14), # Adjust axis titles font size
  669. axis.text = element_text(size = 12), # Adjust axis text font size
  670. legend.position = "none") # Remove legend
  672. ################################################################################
  673. ############################ Gene expression fo teh genes in the terms/gene sets
  674. pca_data4 <- as.data.frame(pca4$x)
  675. pca_data4$sample <- rownames(pca_data4)
  677. # Assign colors to samples
  678. pca_data4$color <- c("#0072B2","#0072B2","#0072B2", # blue
  679. "#D55E00","#D55E00","#D55E00", #orange
  680. "#CC79A7","#CC79A7","#CC79A7", #pinky purple
  681. "#F0E442","#F0E442","#F0E442") #yellow
  683. custom_colors4 <- c("#0072B2","#0072B2","#0072B2", # blue
  684. "#D55E00","#D55E00","#D55E00", #orange
  685. "#CC79A7","#CC79A7","#CC79A7", #pinky purple
  686. "#F0E442","#F0E442","#F0E442") #yellow
  687. names(custom_colors4) <- pca_data4$sample
  689. # Plot the PCA
  690. ggplot(pca_data4, aes(x = PC1, y = PC2, label = sample, color = sample)) +
  691. geom_point(size = 5) + # Adjust point size
  692. geom_text_repel(aes(label = sample), size = 5, nudge_x = 0.05, direction = "y", hjust = 0) + # Adjust text label position
  693. labs(title = "PCA of Gene Set Expression",
  694. x = "Principal Component 1",
  695. y = "Principal Component 2") +
  696. scale_color_manual(values = custom_colors4) +
  697. theme_minimal() +
  698. theme(aspect.ratio = 1, # Make the plot square
  699. plot.title = element_text(size = 16), # Adjust title font size
  700. axis.title = element_text(size = 14), # Adjust axis titles font size
  701. axis.text = element_text(size = 12), # Adjust axis text font size
  702. legend.position = "none") # Remove legend
  705. heatmap.2(Oostrum23_geneExpr,
  706. trace = "none",
  707. col = heatmap_colors,
  708. margins = c(10, 10),
  709. key = TRUE,
  710. keysize = 1,
  711. density.info = "none",
  712. dendrogram = "none",
  713. main = "Heatmap of Gene Sets Gene Expression",
  714. Colv = FALSE,
  715. Rowv = FALSE)
  717. # Same thing as above, different format:
  718. heatmap(Oostrum23_geneExpr,
  719. col = heatmap_colors, Colv = NULL, Rowv = NULL)
  721. # Add legend for heatmap colors
  722. min_value <- as.numeric(min(Oostrum23_geneExpr))
  723. max_value <- as.numeric(max(Oostrum23_geneExpr))
  724. mid_value <- (min_value + max_value) / 2
  725. min_value <- round(min_value, 3)
  726. max_value <- round(max_value, 3)
  727. mid_value <- round(mid_value, 3)
  729. # Create legend labels
  730. legend_labels <- c(min_value, mid_value, max_value)
  732. # Legend plot for values
  733. legend("left", legend = legend_labels,
  734. fill = colorRampPalette(c("lightblue", "white", "#FF457D"))(n = 3),
  735. title = expression("Log"[2]*"CPM values"),
  736. cex = 1.5, bty = "n", inset = 0.05)
  738. ################## Correlation of terms expression
  739. # Define your color palette
  740. colorLegend <- c("#0072B2", "#D55E00", "#CC79A7", "#F0E442")
  741. names(colorLegend) <- c("CA", "4z", "4c", "4p")
  743. # Generate heatmap colors
  744. heatmap_colors <- colorRampPalette(c("lightblue", "white", "#FF457D"))(n = 1000)
  746. # Define the color breaks for the heatmap
  747. breaks <- seq(min(correlation_matrix, na.rm = TRUE), max(correlation_matrix, na.rm = TRUE), length.out = length(heatmap_colors) + 1)
  749. # Create the heatmap
  750. heatmap.2(
  751. correlation_matrix,
  752. col = heatmap_colors,
  753. breaks = breaks,
  754. trace = "none",
  755. margins = c(20, 20),
  756. main = "Correlation Matrix of Gene Sets (expression aggregated across samples)",
  757. key = TRUE,
  758. keysize = 1.5,
  759. dendrogram = "none",
  760. Colv = FALSE,
  761. Rowv = FALSE
  762. )
  764. ###################
  765. correlation_matrix2 = cor(t(gsva_res))
  767. # Plot the correlation heatmap
  768. heatmap.2(correlation_matrix2,
  769. trace = "none",
  770. col = bluered(256),
  771. margins = c(8, 8),
  772. main = "Correlation Matrix of Gene Sets (enrichment scoring)",
  773. density.info = "none")
  774. ###################
  775. # Ensure your gene sets contain only genes present in the gene expression matrix
  776. gene_sets <- lapply(gene_set, function(genes) {
  777. genes[genes %in% rownames(vobjDream_E)]
  778. })
  780. # Initialize a matrix to store average expressions per gene set per sample
  781. avg_expression_matrix <- matrix(NA, nrow = length(gene_sets), ncol = ncol(vobjDream_E))
  782. rownames(avg_expression_matrix) <- names(gene_sets)
  783. colnames(avg_expression_matrix) <- colnames(vobjDream_E)
  785. # Calculate the average expression for each gene set per sample
  786. for (i in seq_along(gene_sets)) {
  787. genes <- gene_sets[[i]]
  788. if (length(genes) > 0) {
  789. avg_expression_matrix[i, ] <- colMeans(vobjDream_E[genes, , drop = FALSE])
  790. }
  791. }
  793. # Compute the correlation matrix for the averaged expression values per sample
  794. correlation_matrix <- cor(t(avg_expression_matrix), use = "complete.obs")
  796. # Print the correlation matrix
  797. print(correlation_matrix)
  799. # Plot the correlation matrix as a heatmap
  800. heatmap.2(correlation_matrix,
  801. trace = "none",
  802. col = colorRampPalette(brewer.pal(9, "Blues"))(256),
  803. margins = c(8, 8),
  804. main = "Correlation Matrix of Gene Sets (gene expression average per sample)",
  805. density.info = "none",
  806. key.title = "Correlation",
  807. key.xlab = "Correlation",
  808. cexRow = 0.7,
  809. cexCol = 0.7)
  812. dev.off()
  815. ############## Plotting main figure 4p vs CA
  816. # Define colorblind-friendly colors manually
  817. PCA_colors <- c("#000000", # Black (for CA)
  818. "#000000", # Black (for CA)
  819. "#000000", # Black (for CA)
  820. "#D55E00", # Orange (for 4z)
  821. "#D55E00", # Orange (for 4z)
  822. "#D55E00", # Orange (for 4z)
  823. "#0072B2", # Sky Blue (for 4c)
  824. "#0072B2", # Sky Blue (for 4c)
  825. "#0072B2", # Sky Blue (for 4c)
  826. "#CC79A7", # Magenta (for 4p)
  827. "#CC79A7", # Magenta (for 4p)
  828. "#CC79A7") # Magenta (for 4p)
  830. pca_all = prcomp(t(syn_matrix))
  831. pca_data_all <- as.data.frame(pca_all$x)
  832. pca_data_all$sample <- rownames(pca_data_all)
  834. names(PCA_colors) <- pca_data_all$sample
  836. # Plot the PCA
  837. PCA_4p_CA_1 = ggplot(pca_data_all, aes(x = PC1, y = PC2, label = sample, color = sample)) +
  838. geom_point(size = 10) + # Adjust point size
  839. geom_text_repel(aes(label = sample), size = 15, nudge_x = 0.05, direction = "y", hjust = 0) + # Adjust text label position
  840. labs(x = "Principal Component 1",
  841. y = "Principal Component 2") +
  842. scale_color_manual(values = PCA_colors) +
  843. theme_minimal() +
  844. theme(
  845. aspect.ratio = 1, # Make the plot square
  846. plot.title = element_text(size = 35), # Adjust title font size
  847. axis.title = element_text(size = 40), # Adjust axis titles font size
  848. axis.text = element_text(size = 35), # Adjust axis text font size
  849. legend.position = "none", # Remove legend
  850. panel.grid.major = element_line(color = "grey90"), # Adjust major grid lines
  851. panel.grid.minor = element_line(color = "grey60"), # Adjust minor grid lines
  852. panel.border = element_rect(color = "black", fill = NA, size = 1) # Add a black border around the plot
  853. )
  855. pdf("/VariancePartitioning/4p_vs_CA_PCA_VariancePartition_Syn.pdf",
  856. width = 13.5, height = 12.5, bg = "transparent")
  857. print(PCA_4p_CA_1)
  858. dev.off() # Close the device
  861. volcano_syn_4p_CA = EnhancedVolcano(compare4p_CA,
  862. lab = compare4p_CA$name,
  863. subtitle = NULL,
  864. x = 'logFC',
  865. y = 'adj.P.Val',
  866. pCutoff = 0.05,
  867. FCcutoff = 1,
  868. xlim = c(-9, 11),
  869. ylim = c(0, 6),
  870. pointSize = 1.5,
  871. axisLabSize = 35,
  872. labSize = 10,
  873. title = NULL,
  874. legendPosition = "top",
  875. legendLabSize = 35,
  876. legendLabels = c("NS", expression(Log[2] ~ FC), "p-value", expression("p-value" ~ and
  877. ~ log[2] ~ FC)),
  878. col = c('grey30', '#E9E9E9', "#0072B2", '#D55E00'),
  879. colAlpha = 0.7,
  880. #drawConnectors = TRUE,
  881. widthConnectors = 0.5)
  884. pdf("/VariancePartitioning/4p_vs_CA_Volcano_VariancePartition_Syn.pdf",
  885. width = 17.5, height = 13.5, bg = "transparent")
  886. print(volcano_syn_4p_CA)
  887. dev.off() # Close the device
  890. #Prepare the data frame for up genes
  891. Syn_4p_CA_up_genes <- gost_up_genes_compare4p_CA$result %>%
  892. filter(source %in% c("GO:BP", "GO:CC", "GO:MF", "KEGG")) %>%
  893. mutate(Ratio = intersection_size / term_size) %>%
  894. group_by(source) %>%
  895. slice_min(p_value, n = 7) %>%
  896. ungroup() %>%
  897. mutate(term_name = factor(term_name, levels = unique(term_name))) # Maintain original order
  900. # Plot for down genes GO
  901. Bubble_4p_CA_1 <- ggplot(Syn_4p_CA_up_genes, aes(x = Ratio, y = term_name)) +
  902. geom_point(aes(size = intersection_size, color = -log10(p_value)), alpha = 0.7) +
  903. scale_color_gradient(
  904. low = "#0072B2",
  905. high = "#D55E00",
  906. breaks = c(min(-log10(Syn_4p_CA_up_genes$p_value)),
  907. max(-log10(Syn_4p_CA_up_genes$p_value))),
  908. labels = c(
  909. paste0(round(min(-log10(Syn_4p_CA_up_genes$p_value)), 2)),
  910. paste0(round(max(-log10(Syn_4p_CA_up_genes$p_value)), 2))
  911. )
  912. ) + # Color based on p-value with custom breaks
  913. scale_size(range = c(5, 10)) + # Adjust bubble sizes
  914. labs(size = "Intersection Size\n", color = "-log10(p-value)\n", y = NULL) + # Axis and legend labels
  915. theme_minimal() + # Clean theme
  916. theme(
  917. axis.text.y = element_text(size = 35), # Adjust y-axis text size
  918. axis.text.x = element_text(size = 30, hjust = 1), # Adjust x-axis text size and angle
  919. axis.title.x = element_text(size = 35), # Adjust x-axis title size
  920. legend.title = element_text(size = 30), # Font size of legend title
  921. legend.text = element_text(size = 25),
  922. legend.spacing.y = unit(2, "cm") # Increase spacing between legend title and text
  923. )
  925. pdf("/VariancePartitioning/4p_vs_CA_BubbleUp_VariancePartition_Syn.pdf",
  926. width = 22.5, height = 13.5, bg = "transparent")
  927. print(Bubble_4p_CA_1)
  928. dev.off() # Close the device
  931. # Apply the updated code for down genes without text wrapping < sorta better
  932. Syn_4p_CA_down_genes <- gost_down_genes_compare4p_CA$result %>%
  933. filter(source %in% c("GO:BP", "GO:CC", "GO:MF", "KEGG")) %>%
  934. mutate(Ratio = intersection_size / term_size) %>%
  935. group_by(source) %>%
  936. slice_min(p_value, n = 7) %>%
  937. ungroup() %>%
  938. mutate(term_name = factor(term_name, levels = unique(term_name))) # Maintain original order
  940. # Plot for down genes GO
  941. Bubble_4p_CA_2 <- ggplot(Syn_4p_CA_down_genes, aes(x = Ratio, y = term_name)) +
  942. geom_point(aes(size = intersection_size, color = -log10(p_value)), alpha = 0.7) +
  943. scale_color_gradient(
  944. low = "#0072B2",
  945. high = "#D55E00",
  946. breaks = c(min(-log10(Syn_4p_CA_down_genes$p_value)),
  947. max(-log10(Syn_4p_CA_down_genes$p_value))),
  948. labels = c(
  949. paste0(round(min(-log10(Syn_4p_CA_down_genes$p_value)), 2)),
  950. paste0(round(max(-log10(Syn_4p_CA_down_genes$p_value)), 2))
  951. )
  952. ) + # Color based on p-value with custom breaks
  953. scale_size(range = c(5, 10)) + # Adjust bubble sizes
  954. labs(size = "Intersection Size\n", color = "-log10(p-value)\n", y = NULL) + # Axis and legend labels
  955. theme_minimal() + # Clean theme
  956. theme(
  957. axis.text.y = element_text(size = 35), # Adjust y-axis text size
  958. axis.text.x = element_text(size = 30, hjust = 1), # Adjust x-axis text size and angle
  959. axis.title.x = element_text(size = 35), # Adjust x-axis title size
  960. legend.title = element_text(size = 30), # Font size of legend title
  961. legend.text = element_text(size = 25),
  962. legend.spacing.y = unit(2, "cm") # Increase spacing between legend title and text
  963. )
  965. pdf("/VariancePartitioning/4p_vs_CA_BubbleDown_VariancePartition_Syn.pdf",
  966. width = 22.5, height = 13.5, bg = "transparent")
  967. print(Bubble_4p_CA_2)
  968. dev.off() # Close the device
  970. # Create bubble plot of gene set enrichments
  971. gsva_res2 <- as.data.frame(gsva_res)
  973. # Convert row names to a column
  974. gsva_res2 <- gsva_res2 %>%
  975. tibble::rownames_to_column(var = "geneset")
  978. # Convert matrix to long format and add value column
  979. long_data <- gsva_res2 %>%
  980. pivot_longer(
  981. cols = -geneset, # Exclude the geneset column from pivoting
  982. names_to = "replicate",
  983. values_to = "enrichment_score"
  984. ) %>%
  985. filter(!is.na(enrichment_score)) %>% # Remove NA values if any
  986. mutate(value = enrichment_score, # Example: scaling enrichment score
  987. label = sprintf("%.2f", value)) # Format values to 2 decimal places
  989. # Create heatmap of gene set enrichments
  990. Genesets_4p_CA_heatmap = ggplot(long_data, aes(x = replicate, y = geneset, fill = enrichment_score)) +
  991. geom_tile(color = "white") + # Create heatmap tiles with white borders between cells
  992. scale_fill_gradient2(low = "#0072B2", mid = "white", high = "#D55E00", midpoint = 0, limits = c(min(long_data$enrichment_score), max(long_data$enrichment_score))) + # Color gradient with white midpoint
  993. geom_text(aes(label = label), color = "white", size = 10, vjust = 0.5, fontface = "bold") + # Add text labels inside tiles
  994. labs(
  995. x = "Replicates",
  996. y = "Gene Sets",
  997. fill = "GSVA\nEnrichment Score" # Use newline character (\n) to split the legend title
  998. ) +
  999. theme_minimal() + # Use minimal theme for a clean look
  1000. theme(
  1001. axis.text.x = element_text(size = 25), # Increase x-axis text size
  1002. axis.text.y = element_text(size = 35), # Increase y-axis text size
  1003. axis.title.x = element_text(size = 35), # Increase x-axis title size
  1004. axis.title.y = element_text(size = 35), # Increase y-axis title size
  1005. plot.title = element_text(size = 25, face = "bold"), # Increase plot title size and make it bold
  1006. legend.title = element_text(size = 30), # Font size of legend title
  1007. legend.text = element_text(size = 30),
  1008. legend.position = "right" # Position the legend on the right
  1009. )
  1011. pdf("/VariancePartitioning/4p_vs_CA_GeneSetsGSVAHeatmap_VariancePartition_Syn.pdf",
  1012. width = 24, height = 16, bg = "transparent")
  1013. #print(Genesets_4p_CA)
  1014. print(Genesets_4p_CA_heatmap)
  1015. dev.off() # Close the device
  1018. ################################################################################
  1019. ############## Plotting supplementary figure 4z vs CA
  1021. volcano_syn_4z_CA = EnhancedVolcano(compare4z_CA,
  1022. lab = compare4z_CA$name,
  1023. subtitle = NULL,
  1024. x = 'logFC',
  1025. y = 'adj.P.Val',
  1026. pCutoff = 0.05,
  1027. FCcutoff = 1,
  1028. xlim = c(-9, 11),
  1029. ylim = c(0, 6),
  1030. pointSize = 1.5,
  1031. axisLabSize = 35,
  1032. labSize = 10,
  1033. title = NULL,
  1034. legendPosition = "top",
  1035. legendLabSize = 35,
  1036. legendLabels = c("NS", expression(Log[2] ~ FC), "p-value", expression("p-value" ~ and
  1037. ~ log[2] ~ FC)),
  1038. col = c('grey30', '#E9E9E9', "#0072B2", '#D55E00'),
  1039. colAlpha = 0.7,
  1040. drawConnectors = TRUE,
  1041. arrowheads = FALSE,
  1042. max.overlaps = 15,
  1043. widthConnectors = 0.5)
  1045. pdf("/VariancePartitioning/4z_vs_CA_Volcano_VariancePartition_Syn.pdf",
  1046. width = 17.5, height = 13.5, bg = "transparent")
  1047. print(volcano_syn_4z_CA)
  1048. dev.off() # Close the device
  1050. # Apply the updated code for up genes with additional filtering and calculations
  1051. Syn_4z_CA_up_genes <- gost_up_genes_compare4z_CA$result %>%
  1052. filter(source %in% c("GO:BP", "GO:CC", "GO:MF", "KEGG")) %>%
  1053. mutate(Ratio = intersection_size / term_size) %>%
  1054. group_by(source) %>%
  1055. slice_min(p_value, n = 12) %>%
  1056. ungroup() %>%
  1057. mutate(term_name = factor(term_name, levels = unique(term_name))) # Maintain original order
  1059. # Plot for down genes GO
  1060. Bubble_4z_CA_1 <- ggplot(Syn_4z_CA_up_genes, aes(x = Ratio, y = term_name)) +
  1061. geom_point(aes(size = intersection_size, color = -log10(p_value)), alpha = 0.7) +
  1062. scale_color_gradient(
  1063. low = "#0072B2",
  1064. high = "#D55E00",
  1065. breaks = c(min(-log10(Syn_4z_CA_up_genes$p_value)),
  1066. max(-log10(Syn_4z_CA_up_genes$p_value))),
  1067. labels = c(
  1068. paste0(round(min(-log10(Syn_4z_CA_up_genes$p_value)), 2)),
  1069. paste0(round(max(-log10(Syn_4z_CA_up_genes$p_value)), 2))
  1070. )
  1071. ) + # Color based on p-value with custom breaks
  1072. scale_size(range = c(5, 10)) + # Adjust bubble sizes
  1073. labs(size = "Intersection Size\n", color = "-log10(p-value)\n", y = NULL) + # Axis and legend labels
  1074. theme_minimal() + # Clean theme
  1075. theme(
  1076. axis.text.y = element_text(size = 35), # Adjust y-axis text size
  1077. axis.text.x = element_text(size = 30, hjust = 1), # Adjust x-axis text size and angle
  1078. axis.title.x = element_text(size = 35), # Adjust x-axis title size
  1079. legend.title = element_text(size = 30), # Font size of legend title
  1080. legend.text = element_text(size = 25),
  1081. legend.spacing.y = unit(2, "cm") # Increase spacing between legend title and text
  1082. )
  1084. pdf("/VariancePartitioning/4z_vs_CA_BubbleUp_VariancePartition_Syn.pdf",
  1085. width = 22.5, height = 13.5, bg = "transparent")
  1086. print(Bubble_4z_CA_1)
  1087. dev.off() # Close the device
  1090. # Apply the updated code for down genes without text wrapping
  1091. Syn_4z_CA_down_genes <- gost_down_genes_compare4z_CA$result %>%
  1092. filter(source %in% c("GO:BP", "GO:CC", "GO:MF", "KEGG")) %>%
  1093. mutate(Ratio = intersection_size / term_size) %>%
  1094. group_by(source) %>%
  1095. slice_min(p_value, n = 7) %>%
  1096. ungroup() %>%
  1097. mutate(term_name = factor(term_name, levels = unique(term_name))) # Maintain original order
  1099. # Plot for down genes GO
  1100. Bubble_4z_CA_2 <- ggplot(Syn_4z_CA_down_genes, aes(x = Ratio, y = term_name)) +
  1101. geom_point(aes(size = intersection_size, color = -log10(p_value)), alpha = 0.7) +
  1102. scale_color_gradient(
  1103. low = "#0072B2",
  1104. high = "#D55E00",
  1105. breaks = c(min(-log10(Syn_4z_CA_down_genes$p_value)),
  1106. max(-log10(Syn_4z_CA_down_genes$p_value))),
  1107. labels = c(
  1108. paste0(round(min(-log10(Syn_4z_CA_down_genes$p_value)), 2)),
  1109. paste0(round(max(-log10(Syn_4z_CA_down_genes$p_value)), 2))
  1110. )
  1111. ) + # Color based on p-value with custom breaks
  1112. scale_size(range = c(5, 10)) + # Adjust bubble sizes
  1113. labs(size = "Intersection Size\n", color = "-log10(p-value)\n", y = NULL) + # Axis and legend labels
  1114. theme_minimal() + # Clean theme
  1115. theme(
  1116. axis.text.y = element_text(size = 35), # Adjust y-axis text size
  1117. axis.text.x = element_text(size = 30, hjust = 1), # Adjust x-axis text size and angle
  1118. axis.title.x = element_text(size = 35), # Adjust x-axis title size
  1119. legend.title = element_text(size = 30), # Font size of legend title
  1120. legend.text = element_text(size = 25),
  1121. legend.spacing.y = unit(2, "cm") # Increase spacing between legend title and text
  1122. )
  1124. pdf("/VariancePartitioning/4z_vs_CA_BubbleDown_VariancePartition_Syn.pdf",
  1125. width = 22.5, height = 10, bg = "transparent")
  1126. print(Bubble_4z_CA_2) # Print the plot to the PNG device
  1127. dev.off() # Close the device
  1129. ################################################################################
  1130. ############## Plotting supplementary figure 4c vs CA
  1131. volcano_syn_4c_CA = EnhancedVolcano(compare4c_CA,
  1132. lab = compare4c_CA$name,
  1133. subtitle = NULL,
  1134. x = 'logFC',
  1135. y = 'adj.P.Val',
  1136. pCutoff = 0.05,
  1137. FCcutoff = 1,
  1138. xlim = c(-9, 11),
  1139. ylim = c(0, 6),
  1140. pointSize = 1.5,
  1141. axisLabSize = 35,
  1142. labSize = 10,
  1143. title = NULL,
  1144. legendPosition = "top",
  1145. legendLabSize = 35,
  1146. legendLabels = c("NS", expression(Log[2] ~ FC), "p-value", expression("p-value" ~ and
  1147. ~ log[2] ~ FC)),
  1148. col = c('grey30', '#E9E9E9', "#0072B2", '#D55E00'),
  1149. colAlpha = 0.7,
  1150. #drawConnectors = TRUE,
  1151. widthConnectors = 0.5)
  1153. pdf("/VariancePartitioning/4c_vs_CA_Volcano_VariancePartition_Syn.pdf",
  1154. width = 17.5, height = 13.5, bg = "transparent")
  1155. print(volcano_syn_4c_CA) # Print the plot to the PNG device
  1156. dev.off() # Close the device
  1158. # Apply the updated code for up genes with additional filtering and calculations
  1159. Syn_4c_CA_up_genes <- gost_up_genes_compare4c_CA$result %>%
  1160. filter(source %in% c("GO:BP", "GO:CC", "GO:MF", "KEGG")) %>%
  1161. mutate(Ratio = intersection_size / term_size) %>%
  1162. group_by(source) %>%
  1163. slice_min(p_value, n = 7) %>%
  1164. ungroup() %>%
  1165. mutate(term_name = factor(term_name, levels = unique(term_name))) # Maintain original order
  1167. # Plot for down genes GO
  1168. Bubble_4c_CA_1 <- ggplot(Syn_4c_CA_up_genes, aes(x = Ratio, y = term_name)) +
  1169. geom_point(aes(size = intersection_size, color = -log10(p_value)), alpha = 0.7) +
  1170. scale_color_gradient(
  1171. low = "#0072B2",
  1172. high = "#D55E00",
  1173. breaks = c(min(-log10(Syn_4c_CA_up_genes$p_value)),
  1174. max(-log10(Syn_4c_CA_up_genes$p_value))),
  1175. labels = c(
  1176. paste0(round(min(-log10(Syn_4c_CA_up_genes$p_value)), 2)),
  1177. paste0(round(max(-log10(Syn_4c_CA_up_genes$p_value)), 2))
  1178. )
  1179. ) + # Color based on p-value with custom breaks
  1180. scale_size(range = c(5, 10)) + # Adjust bubble sizes
  1181. labs(size = "Intersection Size\n", color = "-log10(p-value)\n", y = NULL) + # Axis and legend labels
  1182. theme_minimal() + # Clean theme
  1183. theme(
  1184. axis.text.y = element_text(size = 35), # Adjust y-axis text size
  1185. axis.text.x = element_text(size = 30, hjust = 1), # Adjust x-axis text size and angle
  1186. axis.title.x = element_text(size = 35), # Adjust x-axis title size
  1187. legend.title = element_text(size = 30), # Font size of legend title
  1188. legend.text = element_text(size = 25),
  1189. legend.spacing.y = unit(2, "cm") # Increase spacing between legend title and text
  1190. )
  1192. pdf("/VariancePartitioning/4c_vs_CA_BubbleUp_VariancePartition_Syn.pdf",
  1193. width = 22.5, height = 13.5, bg = "transparent")
  1194. print(Bubble_4c_CA_1)
  1195. dev.off() # Close the device
  1198. # Apply the updated code for down genes without text wrapping
  1199. Syn_4c_CA_down_genes <- gost_down_genes_compare4c_CA$result %>%
  1200. filter(source %in% c("GO:BP", "GO:CC", "GO:MF", "KEGG")) %>%
  1201. mutate(Ratio = intersection_size / term_size) %>%
  1202. group_by(source) %>%
  1203. slice_min(p_value, n = 7) %>%
  1204. ungroup() %>%
  1205. mutate(term_name = factor(term_name, levels = unique(term_name))) # Maintain original order
  1207. # Plot for down genes GO
  1208. Bubble_4c_CA_2 <- ggplot(Syn_4c_CA_down_genes, aes(x = Ratio, y = term_name)) +
  1209. geom_point(aes(size = intersection_size, color = -log10(p_value)), alpha = 0.7) +
  1210. scale_color_gradient(
  1211. low = "#0072B2",
  1212. high = "#D55E00",
  1213. breaks = c(min(-log10(Syn_4c_CA_down_genes$p_value)),
  1214. max(-log10(Syn_4c_CA_down_genes$p_value))),
  1215. labels = c(
  1216. paste0(round(min(-log10(Syn_4c_CA_down_genes$p_value)), 2)),
  1217. paste0(round(max(-log10(Syn_4c_CA_down_genes$p_value)), 2))
  1218. )
  1219. ) + # Color based on p-value with custom breaks
  1220. scale_size(range = c(5, 10)) + # Adjust bubble sizes
  1221. labs(size = "Intersection Size\n", color = "-log10(p-value)\n", y = NULL) + # Axis and legend labels
  1222. theme_minimal() + # Clean theme
  1223. theme(
  1224. axis.text.y = element_text(size = 35), # Adjust y-axis text size
  1225. axis.text.x = element_text(size = 30, hjust = 1), # Adjust x-axis text size and angle
  1226. axis.title.x = element_text(size = 35), # Adjust x-axis title size
  1227. legend.title = element_text(size = 30), # Font size of legend title
  1228. legend.text = element_text(size = 25),
  1229. legend.spacing.y = unit(2, "cm") # Increase spacing between legend title and text
  1230. )
  1232. pdf("/VariancePartitioning/4c_vs_CA_BubbleDown_VariancePartition_Syn.pdf",
  1233. width = 22.5, height = 13.5, bg = "transparent")
  1234. print(Bubble_4c_CA_2)
  1235. dev.off() # Close the device

Dream_Synaptosome_DE_code.R, under CC-BY-4.0 · at the source

Overview

Authors: Huanhuan Chen1, Ye Wei1, Luminiţa Ruje2,3, Fusheng Du1, Zhendong Feng4, Qi Wan5, Mikhail Spivakov2,3, Oleg O. Glebov1,6
ORCID iDs: Zhendong Feng, Oleg O. Glebov
  1. Institute of Neuroregeneration and Neurorehabilitation, Qingdao University, Qingdao, Shandong, China
  2. MRC Laboratory of Medical Sciences, London, United Kingdom
  3. Institute of Clinical Sciences, Imperial College, London, United Kingdom
  4. Department of Neurobiology, School of Basic Medical Sciences, National Institute on Drug Dependence, Peking University, Beijing, China
  5. Faculty of Life and Health, Shenzhen University of Advanced Technology, Shenzhen, Guangdong, China
  6. Department of Psychological Medicine, Institute of Psychiatry, Psychology & Neuroscience, King’s College London, London, United Kingdom
Institutions: Qingdao University (China); MRC London Institute of Medical Sciences (United Kingdom); Imperial College London (United Kingdom); Peking University (China); Shenzhen University of Advanced Technology (China); King's College London (United Kingdom)
Journal: PLoS biology, volume 24, issue 3, article e3003608
Dates: received 26 August 2025; accepted 6 January 2026; published online 2 March 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1371/journal.pbio.3003608 · PMID 41770794 · PMCID PMC12981561 · OpenAlex W7133243098
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), histology / microscopy (modality), rat (organism), stroke (population), cellular / molecular (subfield)
Methods: Spectral & time-frequency, Statistics, Smoothing, state filtering, decompositions, Connectivity, Machine learning, fMRI & imaging
MeSH: Brain*, Brain Ischemia*, Ischemic Stroke*, Neuronal Plasticity*, Stroke*, Synapses*, Animals, Infarction, Middle Cerebral Artery, Male, Proteomics, Rats, Rats, Sprague-Dawley, Receptors, N-Methyl-D-Aspartate (* major topic)
Journal subjects: Medicine and Health Sciences, Medical Conditions, Cerebrovascular Diseases, Stroke, Ischemic Stroke, Neurology, Vascular Medicine, Biology and Life Sciences, Anatomy, Nervous System, Synapses, Synaptosomes, Physiology, Electrophysiology, Neurophysiology, Neuroscience, Genetics, Gene Expression, Research and Analysis Methods, Specimen Preparation and Treatment, Staining, Immunostaining, Cellular Neuroscience, Synaptic Plasticity, Developmental Neuroscience, Ischemia
Topic: Neuroscience and Neuropharmacology Research (Cellular and Molecular Neuroscience, Neuroscience), according to OpenAlex
Funding: Lewy Body Society (OOG2019/2020); National Natural Science Foundation of China (National Science Foundation of China) (32070772); UK Medical Research Council (MC-A652-5QA20)
Citations: cited by 2 papers (Europe PMC); 114 references in the paper

Abstract

Physiological mechanisms of the key hyperacute (0–24 hours) stage of stroke are poorly understood, hampering the development of new therapies. Synaptic plasticity has been strongly implicated in early stages of neurodegenerative and neurodevelopmental disorders, yet its relevance in early stroke remains unclear. Here, we describe the emergence of distinct region-specific forms of synaptic remodeling following middle cerebral artery occlusion in rats, arising within the critical 4-hour period. Synapses within the severely ischemic core region were rapidly lost, while those in the mildly ischemic penumbra, albeit largely structurally intact, were functionally diminished. In contrast, the contralateral cortex exhibited increased synaptic staining and synaptic vesicle cycling. Systemic pharmacological blockade of NMDA-type glutamate receptors abolished contralateral synaptic increase and exacerbated synaptic decline in the penumbra. Proteomic and transcriptomic analyses showed that cross-brain synaptic plasticity is independent of local gene expression and revealed metabolic rearrangement and synaptic downregulation in the penumbra. These findings identify brain-wide synaptic rebalancing as a potential mechanism for rapid functional compensation in hyperacute stroke, highlighting the extent of brain response to acute perturbation.

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

Repository

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

Zenodo 17987265

License: CC-BY-4.0
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Languages: R (1)
Size: 1 file, 1 script
Software Heritage: not checked
Found in: “Data Availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: cowplot (1 file), DESeq2 (1 file), edgeR (1 file), patchwork (1 file), tidyverse (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers (HTTP 200)
  • 30 September 2026: the link answers (HTTP 200)
1 file
At the source: zenodo.org/records/17987265

The paper's code and data availability statement is in the Data section.

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;
  • 1 script, each with its path and the digest of its content;
  • 1 match 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

Data links

Data Availability

All relevant data are contained within the paper and its Supporting information files, or freely available online. Proteomics data were analyzed by Oe Biotech using a proprietary pipeline. Mass spectrometry proteomics data is available at ProteomeXchange (https://proteomecentral.proteomexchange.org, dataset identifier PXD058834). Custom code for gene expression analysis has been uploaded to Zenodo (https://zenodo.org/records/17987265). RNAseq data is available at GEO (https://www.ncbi.nlm.nih.gov/geo/, accession number GSE283465). Numerical data is presented in S1 Table.

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

Versions

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

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 8 authors, 13 MeSH terms, 3 funders, 113 references.

Cite

This paper

Chen, H., Wei, Y., Ruje, L., Du, F., Feng, Z., Wan, Q., Spivakov, M., & Glebov, O. O. (2026). Ischemic stroke triggers brain-wide synaptic remodeling within four hours. PLoS biology, 24(3), e3003608. https://doi.org/10.1371/journal.pbio.3003608

BibTeX

@article{chen2026ischemic,
author = {Chen, Huanhuan and Wei, Ye and Ruje, Luminiţa and Du, Fusheng and Feng, Zhendong and Wan, Qi and Spivakov, Mikhail and Glebov, Oleg O.},
title = {{Ischemic stroke triggers brain-wide synaptic remodeling within four hours}},
journal = {PLoS biology},
year = {2026},
month = mar,
volume = {24},
number = {3},
pages = {e3003608},
publisher = {PLOS},
issn = {1544-9173},
doi = {10.1371/journal.pbio.3003608},
url = {https://doi.org/10.1371/journal.pbio.3003608},
pmid = {41770794},
pmcid = {PMC12981561}
}

RIS

TY - JOUR
AU - Chen, Huanhuan
AU - Wei, Ye
AU - Ruje, Luminiţa
AU - Du, Fusheng
AU - Feng, Zhendong
AU - Wan, Qi
AU - Spivakov, Mikhail
AU - Glebov, Oleg O.
TI - Ischemic stroke triggers brain-wide synaptic remodeling within four hours
T2 - PLoS biology
J2 - PLoS Biol
PY - 2026
DA - 2026/03/02
VL - 24
IS - 3
SP - e3003608
SN - 1544-9173
PB - PLOS
DO - 10.1371/journal.pbio.3003608
UR - https://doi.org/10.1371/journal.pbio.3003608
LA - en
ER -

CSL-JSON

{
"id": "10.1371/journal.pbio.3003608",
"type": "article-journal",
"title": "Ischemic stroke triggers brain-wide synaptic remodeling within four hours",
"container-title": "PLoS biology",
"author": [
{
"family": "Chen",
"given": "Huanhuan"
},
{
"family": "Wei",
"given": "Ye"
},
{
"family": "Ruje",
"given": "Luminiţa"
},
{
"family": "Du",
"given": "Fusheng"
},
{
"family": "Feng",
"given": "Zhendong"
},
{
"family": "Wan",
"given": "Qi"
},
{
"family": "Spivakov",
"given": "Mikhail"
},
{
"family": "Glebov",
"given": "Oleg O."
}
],
"container-title-short": "PLoS Biol",
"volume": "24",
"issue": "3",
"page": "e3003608",
"DOI": "10.1371/journal.pbio.3003608",
"PMID": "41770794",
"PMCID": "PMC12981561",
"ISSN": "1544-9173",
"publisher": "PLOS",
"URL": "https://doi.org/10.1371/journal.pbio.3003608",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
2
]
]
}
}

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.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: edgeR, cowplot, patchwork, 1 other tool, histology / microscopy, genetics / omics, cellular / molecular, 2 references
[2] 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: edgeR, DESeq2, cowplot, 2 other tools, genetics / omics, cellular / molecular, 1 reference
[3] doi:10.1002/advs.202521254 [code]
Persistently Increased Expression of PKMzeta and Unbiased Gene Expression Profiles Identify Hippocampal Molecular Traces of a Long-Term Active Place Avoidance Memory and "Shadow" Proteins.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: DESeq2, cowplot, tidyverse, rat, histology / microscopy, genetics / omics, 1 other category, 2 references
[4] doi:10.1186/s12974-026-03885-1 [code]
Shared transcriptomic signatures in perilesional and contralesional cortex after ischemic stroke.
Journal: Journal of neuroinflammation
In common: cowplot, tidyverse, stroke, genetics / omics, cellular / molecular, 3 references
[5] doi:10.1002/imt2.70163 [code]
Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: edgeR, DESeq2, cowplot, 2 other tools, genetics / omics, cellular / molecular
[6] doi:10.3389/fnmol.2026.1844705 [code]
Risperidone regulates the expression of schizophrenia-related genes in the forebrain of adult male mice.
Journal: Frontiers in molecular neuroscience
In common: edgeR, DESeq2, cowplot, 2 other tools, genetics / omics, cellular / molecular
[7] doi:10.1038/s41467-026-73305-8 [code]
Comparative analysis of the cellular landscape in mammalian striatum.
Journal: Nature communications
In common: edgeR, DESeq2, cowplot, 2 other tools, genetics / omics, cellular / molecular
[8] doi:10.1038/s41593-026-02300-5 [code]
Integrated single-cell and spatial transcriptomic profiling in ALS uncovers peripheral-to-central immune infiltration and reprogramming.
Journal: Nature neuroscience
In common: edgeR, DESeq2, cowplot, 2 other tools, genetics / omics, cellular / molecular
[9] doi:10.1016/j.celrep.2026.117073 [code]
Single-cell epigenomics uncovers heterochromatin instability and transcription factor dysfunction during mouse brain aging.
Journal: Cell reports
In common: edgeR, DESeq2, cowplot, 2 other tools, genetics / omics, cellular / molecular
[10] 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: edgeR, DESeq2, cowplot, 2 other tools, genetics / omics, cellular / molecular

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.