OSCR

Transcriptional Response to Chronic Long-Access Fentanyl Self-Administration in Rat Habenula and Amygdala.

Code ↔ Paper

12 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 12 matches
  1. [1] § Materials and Methods Summary › Experiment Summary and RNA‐Seq Data Generation ↔ code/09_SRA/03_metadata.R, lines 1–55 · score 0.98 · inside sound attenuating, operant chambers housed, bilateral tissue punches, Ribosomal RNA depletion, Cayman Chemical, coronal slabs
  2. [2] § Results › Shared and Unique Gene Expression Changes in Hb and Amyg Following Chronic Fentanyl Intake ↔ code/06_GO_KEGG/01_GO_KEGG_Analyses.R, lines 548–592 · score 0.97 · gastric acid secretion, ECM receptor interaction, Parkinson disease, actin cytoskeleton, oxidative phosphorylation, glutamatergic synapse
  3. [3] § Results › Shared and Unique Gene Expression Changes in Hb and Amyg Following Chronic Fentanyl Intake ↔ code/05_DEA/01_Modeling.R, lines 477–532 · score 0.96 · Atp5mc2, Atp6v0e2, Atp6v1e1, Slc39a10, Col4a3, KEGG term
  4. [4] § Results › Shared and Unique Gene Expression Changes in Hb and Amyg Following Chronic Fentanyl Intake ↔ code/05_DEA/01_Modeling.R, lines 477–532 · score 0.91 · Cox6b1, Scn1a, Col9a3, Gria4, Grik1, Kcnj10
  5. [5] § Results › Shared and Unique Gene Expression Changes in Hb and Amyg Following Chronic Fentanyl Intake ↔ code/06_GO_KEGG/01_GO_KEGG_Analyses.R, lines 548–592 · score 0.87 · postsynaptic membrane potential, Col9a3, voltage gated, cation channels, Epha4, Gsn
  6. [6] § Materials and Methods Summary › Experiment Summary and RNA‐Seq Data Generation ↔ code/09_SRA/03_metadata.R, lines 1–55 · score 0.76 · bulk RNA sequencing, jugular catheters, house light, illumination, Drug, daily
  7. [7] § Materials and Methods Summary › Differential Gene Expression (DGE) ↔ code/06_GO_KEGG/01_GO_KEGG_Analyses.R, lines 73–188 · score 0.75 · enriched GO term, cellular components, Biological processes, molecular functions, enrichment, Amyg
  8. [8] § Results › Enrichment of Genes Associated With Fentanyl Intake in Human and Rodent Hb and Amyg Cell Types ↔ code/08_GSEA/01_enrich_DEGs_vs_cell_type_markers.R, lines 896–934 · score 0.65 · fine cell, DRD1, LAMP5, SATB2, nostrin, S14
  9. [9] § Results › Enrichment of Genes Associated With Fentanyl Intake in Human and Rodent Hb and Amyg Cell Types ↔ code/08_GSEA/01_enrich_DEGs_vs_cell_type_markers.R, lines 1295–1361 · score 0.63 · LHb.2, LHb.6, LHb.7, astrocytes, enriched, orthologous
  10. [10] § Materials and Methods Summary › Experiment Summary and RNA‐Seq Data Generation ↔ code/03_Data_preparation/01_build_objects.R, lines 47–85 · score 0.62 · tissue punches, library preparation, sequenced, RNA, brains, Amyg
  11. [11] § Materials and Methods Summary › RNA‐Seq Data Processing and Quality Control ↔ code/04_EDA/01_QCA.R, lines 2–87 · score 0.59 · quality control, expressed genes, QC metrics, RNA extraction, brain regions, filtered
  12. [12] § Materials and Methods Summary › Experiment Summary and RNA‐Seq Data Generation ↔ code/04_EDA/03_Explore_gene_level_effects.R, lines 2–55 · score 0.54 · linear regression, infusion slopes, mg, hour, RNA, intake

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,005 lines · 57 KB · MIT · 3 matches

  2. library(dplyr)
  3. library(tidyr)
  4. library(tibble)
  5. library(purrr)
  6. library(here)
  7. library(SummarizedExperiment)
  8. library(clusterProfiler)
  9. library(org.Rn.eg.db)
  10. library(cowplot)
  11. library(ggplot2)
  12. library(biomaRt)
  13. library(sessioninfo)
  16. ####################### Functional Enrichment Analysis #######################
  18. load(here('processed-data/05_DEA/de_genes_Substance_habenula.Rdata'), verbose = TRUE)
  19. load(here('processed-data/05_DEA/de_genes_Substance_amygdala.Rdata'), verbose = TRUE)
  20. load(here('processed-data/05_DEA/results_Substance_uncorr_vars_habenula.Rdata'), verbose = TRUE)
  21. load(here('processed-data/05_DEA/results_Substance_uncorr_vars_amygdala.Rdata'), verbose = TRUE)
  24. ## Groups of DEGs
  26. ########################
  27. ## Habenula DEGs
  28. ########################
  29. ## Up and down habenula DEGs from model with uncorrelated sample variables
  30. up_hab <- de_genes_habenula[which(de_genes_habenula$logFC>0),]
  31. down_hab <- de_genes_habenula[which(de_genes_habenula$logFC<0),]
  33. ########################
  34. ## Amygdala DEGs
  35. ########################
  36. ## Up and down amygdala DEGs from model with uncorrelated sample variables
  37. up_amy <- de_genes_amygdala[which(de_genes_amygdala$logFC>0),]
  38. down_amy <- de_genes_amygdala[which(de_genes_amygdala$logFC<0),]
  40. #########################################
  41. ## Up/Down unique/shared in Hb/Amyg
  42. #########################################
  43. only_up_hab <- de_genes_habenula[which(!de_genes_habenula$ensemblID %in% de_genes_amygdala$ensemblID & de_genes_habenula$logFC>0),]
  44. only_down_hab <- de_genes_habenula[which(!de_genes_habenula$ensemblID %in% de_genes_amygdala$ensemblID & de_genes_habenula$logFC<0),]
  46. only_up_amy <- de_genes_amygdala[which(!de_genes_amygdala$ensemblID %in% de_genes_habenula$ensemblID & de_genes_amygdala$logFC>0),]
  47. only_down_amy <- de_genes_amygdala[which(!de_genes_amygdala$ensemblID %in% de_genes_habenula$ensemblID & de_genes_amygdala$logFC<0),]
  49. shared_hab_amy <- inner_join(de_genes_habenula, de_genes_amygdala, by = colnames(de_genes_amygdala)[1:9], suffix = c(".hb", ".amyg"))
  51. shared_up_hab_up_amy <- shared_hab_amy[shared_hab_amy$logFC.hb>0 & shared_hab_amy$logFC.amyg>0, ]
  52. shared_up_hab_down_amy <- shared_hab_amy[shared_hab_amy$logFC.hb>0 & shared_hab_amy$logFC.amyg<0, ]
  53. shared_down_hab_up_amy <- shared_hab_amy[shared_hab_amy$logFC.hb<0 & shared_hab_amy$logFC.amyg>0, ]
  54. shared_down_hab_down_amy <- shared_hab_amy[shared_hab_amy$logFC.hb<0 & shared_hab_amy$logFC.amyg<0, ]
  57. ## Retrieve valid Entrez IDs
  58. only_up_hab_genes <- only_up_hab %>% dplyr::filter(!is.na(EntrezID) & !is.null(EntrezID) & EntrezID != "NULL" & EntrezID != "") %>% pull(EntrezID) %>% unique()
  59. only_down_hab_genes <- only_down_hab %>% dplyr::filter(!is.na(EntrezID) & !is.null(EntrezID) & EntrezID != "NULL" & EntrezID != "") %>% pull(EntrezID) %>% unique()
  60. only_up_amy_genes <- only_up_amy %>% dplyr::filter(!is.na(EntrezID) & !is.null(EntrezID) & EntrezID != "NULL" & EntrezID != "") %>% pull(EntrezID) %>% unique()
  61. only_down_amy_genes <- only_down_amy %>% dplyr::filter(!is.na(EntrezID) & !is.null(EntrezID) & EntrezID != "NULL" & EntrezID != "") %>% pull(EntrezID) %>% unique()
  62. shared_up_hab_up_amy_genes <- shared_up_hab_up_amy %>% dplyr::filter(!is.na(EntrezID) & !is.null(EntrezID) & EntrezID != "NULL" & EntrezID != "") %>% pull(EntrezID) %>% unique()
  63. shared_up_hab_down_amy_genes <- shared_up_hab_down_amy %>% dplyr::filter(!is.na(EntrezID) & !is.null(EntrezID) & EntrezID != "NULL" & EntrezID != "") %>% pull(EntrezID) %>% unique()
  64. shared_down_hab_up_amy_genes <- shared_down_hab_up_amy %>% dplyr::filter(!is.na(EntrezID) & !is.null(EntrezID) & EntrezID != "NULL" & EntrezID != "") %>% pull(EntrezID) %>% unique()
  65. shared_down_hab_down_amy_genes <- shared_down_hab_down_amy %>% dplyr::filter(!is.na(EntrezID) & !is.null(EntrezID) & EntrezID != "NULL" & EntrezID != "") %>% pull(EntrezID) %>% unique()
  68. ## Background genes (all genes assessed for DGE) -- same genes in Hb and Amyg DGE
  69. geneUniverse <- results_Substance_uncorr_vars_amygdala[[1]] %>%
  70. dplyr::filter(!is.na(EntrezID) & !is.null(EntrezID) & EntrezID != "NULL" & EntrezID != "") %>% pull(EntrezID) %>% unique()
  73. ## Function to find enriched GO and KEGG terms
  74. GO_KEGG<- function(sigGeneList, geneUniverse, name){
  76. ## GO terms
  77. ## Obtain biological processes
  78. goBP_Adj <- compareCluster(
  79. sigGeneList,
  80. fun = "enrichGO",
  81. universe = geneUniverse,
  82. OrgDb = org.Rn.eg.db,
  83. ont = "BP",
  84. pAdjustMethod = "BH",
  85. qvalueCutoff = 0.05,
  86. readable = TRUE
  87. )
  89. ## Save
  90. if (!is.null(goBP_Adj)){
  91. p1 <- dotplot(goBP_Adj, title="GO Enrichment Analysis: Biological processes")
  92. goBP_Adj <- as.data.frame(goBP_Adj)
  93. goBP_Adj$geneID <- sapply(goBP_Adj$geneID, function(row){gsub("/", ", ", row)})
  94. }
  97. ## Obtain molecular functions
  98. goMF_Adj <- compareCluster(
  99. sigGeneList,
  100. fun = "enrichGO",
  101. universe = geneUniverse,
  102. OrgDb = org.Rn.eg.db,
  103. ont = "MF",
  104. pAdjustMethod = "BH",
  105. qvalueCutoff = 0.05,
  106. readable = TRUE
  107. )
  109. if (!is.null(goMF_Adj)){
  110. p2 <- dotplot(goMF_Adj, title="GO Enrichment Analysis: Molecular function")
  111. goMF_Adj <- as.data.frame(goMF_Adj)
  112. goMF_Adj$geneID <- sapply(goMF_Adj$geneID, function(row){gsub("/", ", ", row)})
  113. }
  116. ## Obtain cellular components
  117. goCC_Adj <- compareCluster(
  118. sigGeneList,
  119. fun = "enrichGO",
  120. universe = geneUniverse,
  121. OrgDb = org.Rn.eg.db,
  122. ont = "CC",
  123. pAdjustMethod = "BH",
  124. qvalueCutoff = 0.05,
  125. readable = TRUE
  126. )
  128. if (!is.null(goCC_Adj)){
  129. p3 <- dotplot(goCC_Adj, title="GO Enrichment Analysis: Cellular components")
  130. goCC_Adj <- as.data.frame(goCC_Adj)
  131. goCC_Adj$geneID <- sapply(goCC_Adj$geneID, function(row){gsub("/", ", ", row)})
  132. }
  135. ## KEGG terms
  136. kegg_Adj <- compareCluster(
  137. sigGeneList,
  138. fun = "enrichKEGG",
  139. organism = 'rat',
  140. universe = geneUniverse,
  141. pAdjustMethod = "BH",
  142. qvalueCutoff = 0.05
  143. )
  145. if (!is.null(kegg_Adj)){
  146. p4 <- dotplot(kegg_Adj, title="KEGG Enrichment Analysis")
  147. ## Add symbols
  148. kegg_Adj <- as.data.frame(kegg_Adj)
  149. genes <-sapply(kegg_Adj$geneID, function(term_genes){unlist(strsplit(term_genes, "/"))})
  150. names(genes) <- kegg_Adj$ID
  151. mart = useMart("ensembl", dataset = "rnorvegicus_gene_ensembl")
  152. term_symbols <- lapply(genes, function(term_genes){
  154. symbols <- getBM(attributes = c("entrezgene_id", "external_gene_name", "ensembl_gene_id"),
  155. filters = "entrezgene_id",
  156. values = term_genes,
  157. mart = mart)
  158. symbols <- apply(symbols, 1, function(gene){if(!is.na(gene["external_gene_name"])){gene["external_gene_name"]}
  159. else if(!is.na(gene["ensembl_gene_id"])){gene["ensembl_gene_id"]}
  160. else{gene["entrezgene_id"]}})
  162. paste(symbols, collapse = ", ")
  164. })
  165. kegg_Adj$geneID <- do.call(rbind, term_symbols)
  166. }
  168. ## Plots
  169. if(name != "Hb_and_Amyg_Up_and_Down_unique_and_shared_DEGs"){
  170. h = 10
  171. w = 14
  172. } else{
  173. h = 35
  174. w = 49
  175. }
  176. plot_grid(p1, p2, p3, p4, ncol=2, align = 'vh')
  177. ggsave(paste("plots/06_GO_KEGG/GO_KEGG_", name, ".pdf", sep=""), height = h, width = w)
  179. ## Save results
  180. goList <- list(
  181. BP = goBP_Adj,
  182. MF = goMF_Adj,
  183. CC = goCC_Adj,
  184. KEGG = kegg_Adj
  185. )
  187. return(goList)
  188. }
  190. #-------------------------------------------------------------------------------
  191. ## 1. Analysis for all DEGs from each brain region
  193. ######################
  194. # Habenula
  195. ######################
  196. sigGeneList <- list("All"= unique(de_genes_habenula[which(!is.na(de_genes_habenula$EntrezID) & !de_genes_habenula$EntrezID=='NULL' & !de_genes_habenula$EntrezID==''), 'EntrezID']))
  197. goList_habenula_all_DEGs <- GO_KEGG(sigGeneList, geneUniverse, 'habenula_all_DEGs')
  198. save(goList_habenula_all_DEGs, file="processed-data/06_GO_KEGG/goList_habenula_all_DEGs.Rdata")
  200. ######################
  201. # Amygdala
  202. ######################
  203. sigGeneList <- list("All"= unique(de_genes_amygdala[which(!is.na(de_genes_amygdala$EntrezID) & !de_genes_amygdala$EntrezID=='NULL' & !de_genes_amygdala$EntrezID==''), 'EntrezID']))
  204. goList_amygdala_all_DEGs<-GO_KEGG(sigGeneList, geneUniverse, 'amygdala_all_DEGs')
  205. save(goList_amygdala_all_DEGs, file="processed-data/06_GO_KEGG/goList_amygdala_all_DEGs.Rdata")
  207. #-------------------------------------------------------------------------------
  208. ## 2. Analysis for up- and down-regulated DEGs from each brain region
  210. ######################
  211. # Habenula
  212. ######################
  214. ## List of DEG sets
  215. sigGeneList <- list("Up"=up_hab[which(!is.na(up_hab$EntrezID) & !up_hab$EntrezID=='NULL' & !up_hab$EntrezID==''), 'EntrezID'],
  216. "Down"=down_hab[which(!is.na(down_hab$EntrezID) & !down_hab$EntrezID=='NULL' & !down_hab$EntrezID==''), 'EntrezID'])
  218. goList_habenula_up_down_DEGs<-GO_KEGG(sigGeneList, geneUniverse, 'habenula_up_down_DEGs')
  219. save(goList_habenula_up_down_DEGs, file="processed-data/06_GO_KEGG/goList_habenula_up_down_DEGs.Rdata")
  221. ## Merge
  222. go_kegg_results_hab <- rbind(cbind(goList_habenula_up_down_DEGs$BP, Ontology = "BP"),
  223. cbind(goList_habenula_up_down_DEGs$MF, Ontology = "MF"),
  224. cbind(goList_habenula_up_down_DEGs$CC, Ontology = "CC"),
  225. cbind(goList_habenula_up_down_DEGs$KEGG[colnames(goList_habenula_up_down_DEGs$BP)], Ontology = "KEGG"))
  226. go_kegg_results_hab$DEGs_set <- go_kegg_results_hab$Cluster
  227. go_kegg_results_hab$Cluster <- NULL
  229. go_kegg_results_hab <- go_kegg_results_hab[, c("Ontology", "DEGs_set", "ID", "Description", "Count", "GeneRatio",
  230. "BgRatio", "FoldEnrichment", "pvalue", "p.adjust", "geneID")]
  231. go_kegg_results_hab <- go_kegg_results_hab[order(go_kegg_results_hab$Ontology, go_kegg_results_hab$DEGs_set, go_kegg_results_hab$p.adjust), ]
  234. write.table(go_kegg_results_hab, "processed-data/Supplementary_Tables/TableS8_GO_KEGG_results_hab.tsv", row.names = FALSE, col.names = TRUE, sep = '\t')
  236. ######################
  237. # Amygdala
  238. ######################
  240. sigGeneList <- list("Up"=up_amy[which(!is.na(up_amy$EntrezID) & !up_amy$EntrezID=='NULL' & !up_amy$EntrezID==''), 'EntrezID'],
  241. "Down"=down_amy[which(!is.na(down_amy$EntrezID) & !down_amy$EntrezID=='NULL' & !down_amy$EntrezID==''), 'EntrezID'])
  243. goList_amygdala_up_down_DEGs<-GO_KEGG(sigGeneList, geneUniverse, 'amygdala_up_down_DEGs')
  244. save(goList_amygdala_up_down_DEGs, file="processed-data/06_GO_KEGG/goList_amygdala_up_down_DEGs.Rdata")
  246. go_kegg_results_amy <- rbind(cbind(goList_amygdala_up_down_DEGs$BP, Ontology = "BP"),
  247. cbind(goList_amygdala_up_down_DEGs$MF, Ontology = "MF"),
  248. cbind(goList_amygdala_up_down_DEGs$CC, Ontology = "CC"),
  249. cbind(goList_amygdala_up_down_DEGs$KEGG[colnames(goList_amygdala_up_down_DEGs$BP)], Ontology = "KEGG"))
  250. go_kegg_results_amy$DEGs_set <- go_kegg_results_amy$Cluster
  251. go_kegg_results_amy$Cluster <- NULL
  253. go_kegg_results_amy <- go_kegg_results_amy[, c("Ontology", "DEGs_set", "ID", "Description", "Count", "GeneRatio",
  254. "BgRatio", "FoldEnrichment", "pvalue", "p.adjust", "geneID")]
  255. go_kegg_results_amy <- go_kegg_results_amy[order(go_kegg_results_amy$Ontology, go_kegg_results_amy$DEGs_set, go_kegg_results_amy$p.adjust), ]
  256. write.table(go_kegg_results_amy, "processed-data/Supplementary_Tables/TableS9_GO_KEGG_results_amy.tsv", row.names = FALSE, col.names = TRUE, sep = '\t')
  258. #-------------------------------------------------------------------------------
  259. ## 3. Analysis for up/down DEGs unique/shared in Hb and Amyg
  261. sigGeneList <- list("Unique in Hb - Up" = only_up_hab_genes,
  262. "Unique in Hb - Down" = only_down_hab_genes,
  263. "Unique in Amyg - Up" = only_up_amy_genes,
  264. "Unique in Amyg - Down" = only_down_amy_genes,
  265. "Shared: Up in Hb, Up in Amyg" = shared_up_hab_up_amy_genes,
  266. "Shared: Up in Hb, Down in Amyg" = shared_up_hab_down_amy_genes,
  267. "Shared: Down in Hb, Up in Amyg" = shared_down_hab_up_amy_genes,
  268. "Shared: Down in Hb, Down in Amyg" = shared_down_hab_down_amy_genes)
  270. goList_hb_and_amyg_DEGs <- GO_KEGG(sigGeneList, geneUniverse,
  271. 'Hb_and_Amyg_Up_and_Down_unique_and_shared_DEGs')
  272. save(goList_hb_and_amyg_DEGs, file="processed-data/06_GO_KEGG/goList_Hb_and_Amyg_Up_and_Down_unique_and_shared_DEGs.Rdata")
  274. ## Merge
  275. go_kegg_results <- rbind(cbind(goList_hb_and_amyg_DEGs$BP, Ontology = "BP"),
  276. cbind(goList_hb_and_amyg_DEGs$MF, Ontology = "MF"),
  277. cbind(goList_hb_and_amyg_DEGs$CC, Ontology = "CC"),
  278. cbind(goList_hb_and_amyg_DEGs$KEGG[colnames(goList_hb_and_amyg_DEGs$BP)], Ontology = "KEGG"))
  280. go_kegg_results$DEGs_set <- go_kegg_results$Cluster
  281. go_kegg_results$Cluster <- NULL
  283. go_kegg_results <- go_kegg_results[, c("Ontology", "DEGs_set", "ID", "Description", "Count", "GeneRatio",
  284. "BgRatio", "FoldEnrichment", "pvalue", "p.adjust", "geneID")]
  285. go_kegg_results <- go_kegg_results[order(go_kegg_results$Ontology, go_kegg_results$DEGs_set, go_kegg_results$p.adjust), ]
  287. write.table(go_kegg_results, "processed-data/Supplementary_Tables/TableS10_GO_KEGG_results_hab_vs_amyg.tsv", row.names = FALSE, col.names = TRUE, sep = '\t')
  289. ## ------
  290. ## Heatmap with GO & KEGG enrichment results for Hb vs Amyg
  292. ## Find enriched terms of interest (list provided by Kristen and Robin)
  294. ## For Hb:
  295. go_kegg_results %>% dplyr::filter(Ontology == "BP", Description == "regulation of synapse structure or activity") %>% .[,1:4]
  296. # Ontology DEGs_set ID Description
  297. # 1 BP Unique in Hb - Up GO:0050803 regulation of synapse structure or activity
  298. # 2 BP Unique in Amyg - Up GO:0050803 regulation of synapse structure or activity
  300. go_kegg_results %>% dplyr::filter(Ontology == "BP", Description == "regulation of presynaptic membrane potential") %>% .[,1:4]
  301. # Ontology DEGs_set ID Description
  302. # 1 BP Unique in Hb - Up GO:0099505 regulation of presynaptic membrane potential
  303. # 2 BP Shared: Up in Hb, Up in Amyg GO:0099505 regulation of presynaptic membrane potential
  305. go_kegg_results %>% dplyr::filter(Ontology == "BP", Description == "regulation of postsynaptic membrane potential") %>% .[,1:4]
  306. # Ontology DEGs_set ID Description
  307. # 1 BP Unique in Hb - Up GO:0060078 regulation of postsynaptic membrane potential
  309. go_kegg_results %>% dplyr::filter(Ontology == "BP", Description == "extracellular matrix organization") %>% .[,1:4]
  310. # Ontology DEGs_set ID Description
  311. # 1 BP Unique in Hb - Down GO:0030198 extracellular matrix organization
  312. # 2 BP Unique in Amyg - Down GO:0030198 extracellular matrix organization
  313. # 3 BP Shared: Down in Hb, Down in Amyg GO:0030198 extracellular matrix organization
  315. go_kegg_results %>% dplyr::filter(Ontology == "BP", Description == "oligodendrocyte differentiation") %>% .[,1:4]
  316. # Ontology DEGs_set ID Description
  317. # 1 BP Unique in Amyg - Down GO:0048709 oligodendrocyte differentiation
  318. # 2 BP Shared: Down in Hb, Down in Amyg GO:0048709 oligodendrocyte differentiation
  321. go_kegg_results %>% dplyr::filter(Ontology == "CC", Description == "presynaptic membrane") %>% .[,1:4]
  322. # Ontology DEGs_set ID Description
  323. # 1 CC Unique in Hb - Up GO:0042734 presynaptic membrane
  324. # 2 CC Shared: Up in Hb, Up in Amyg GO:0042734 presynaptic membrane
  326. go_kegg_results %>% dplyr::filter(Ontology == "CC", Description == "postsynaptic membrane") %>% .[,1:4]
  327. # Ontology DEGs_set ID Description
  328. # 1 CC Unique in Hb - Up GO:0045211 postsynaptic membrane
  329. # 2 CC Shared: Up in Hb, Up in Amyg GO:0045211 postsynaptic membrane
  331. go_kegg_results %>% dplyr::filter(Ontology == "CC", Description == "collagen-containing extracellular matrix") %>% .[,1:4]
  332. # Ontology DEGs_set ID Description
  333. # 1 CC Unique in Hb - Down GO:0062023 collagen-containing extracellular matrix
  334. # 2 CC Unique in Amyg - Down GO:0062023 collagen-containing extracellular matrix
  335. # 3 CC Shared: Down in Hb, Down in Amyg GO:0062023 collagen-containing extracellular matrix
  337. go_kegg_results %>% dplyr::filter(Ontology == "CC", Description == "myelin sheath") %>% .[,1:4]
  338. # Ontology DEGs_set ID Description
  339. # 1 CC Unique in Amyg - Down GO:0043209 myelin sheath
  340. # 2 CC Shared: Down in Hb, Down in Amyg GO:0043209 myelin sheath
  341. #
  343. go_kegg_results %>% dplyr::filter(Ontology == "MF", Description == "voltage-gated monoatomic cation channel activity") %>% .[,1:4]
  344. # Ontology DEGs_set ID Description
  345. # 1 MF Unique in Hb - Up GO:0022843 voltage-gated monoatomic cation channel activity
  346. # 2 MF Shared: Up in Hb, Up in Amyg GO:0022843 voltage-gated monoatomic cation channel activity
  348. go_kegg_results %>% dplyr::filter(Ontology == "MF", Description == "metal ion transmembrane transporter activity") %>% .[,1:4]
  349. # Ontology DEGs_set ID Description
  350. # 1 MF Unique in Hb - Up GO:0046873 metal ion transmembrane transporter activity
  353. go_kegg_results %>% dplyr::filter(Ontology == "KEGG", Description == "Morphine addiction") %>% .[,1:4]
  354. # Ontology DEGs_set ID Description
  355. # 1 KEGG Unique in Hb - Up rno05032 Morphine addiction
  357. go_kegg_results %>% dplyr::filter(Ontology == "KEGG", Description == "Glutamatergic synapse") %>% .[,1:4]
  358. # Ontology DEGs_set ID Description
  359. # 1 KEGG Unique in Hb - Up rno04724 Glutamatergic synapse
  360. # 2 KEGG Unique in Amyg - Up rno04724 Glutamatergic synapse
  361. # 3 KEGG Shared: Up in Hb, Down in Amyg rno04724 Glutamatergic synapse
  363. go_kegg_results %>% dplyr::filter(Ontology == "KEGG", Description == "Gastric acid secretion") %>% .[,1:4]
  364. # Ontology DEGs_set ID Description
  365. # 1 KEGG Unique in Hb - Up rno04971 Gastric acid secretion
  366. # 2 KEGG Unique in Hb - Down rno04971 Gastric acid secretion
  367. # 3 KEGG Shared: Up in Hb, Down in Amyg rno04971 Gastric acid secretion
  369. go_kegg_results %>% dplyr::filter(Ontology == "KEGG", Description == "ECM-receptor interaction") %>% .[,1:4]
  370. # Ontology DEGs_set ID Description
  371. # 1 KEGG Unique in Hb - Down rno04512 ECM-receptor interaction
  372. # 2 KEGG Unique in Amyg - Down rno04512 ECM-receptor interaction
  373. # 3 KEGG Shared: Up in Hb, Down in Amyg rno04512 ECM-receptor interaction
  376. ## For Amyg:
  377. go_kegg_results %>% dplyr::filter(Ontology == "BP", Description == "aerobic respiration") %>% .[,1:4]
  378. # Ontology DEGs_set ID Description
  379. # 1 BP Unique in Amyg - Up GO:0009060 aerobic respiration
  381. go_kegg_results %>% dplyr::filter(Ontology == "BP", Description == "vesicle-mediated transport in synapse") %>% .[,1:4]
  382. # Ontology DEGs_set ID Description
  383. # 1 BP Unique in Amyg - Up GO:0099003 vesicle-mediated transport in synapse
  385. go_kegg_results %>% dplyr::filter(Ontology == "BP", Description == "glial cell differentiation") %>% .[,1:4]
  386. # Ontology DEGs_set ID Description
  387. # 1 BP Unique in Amyg - Down GO:0010001 glial cell differentiation
  388. # 2 BP Shared: Down in Hb, Down in Amyg GO:0010001 glial cell differentiation
  390. go_kegg_results %>% dplyr::filter(Ontology == "BP", Description == "extracellular matrix organization") %>% .[,1:4]
  391. # Ontology DEGs_set ID Description
  392. # 1 BP Unique in Hb - Down GO:0030198 extracellular matrix organization
  393. # 2 BP Unique in Amyg - Down GO:0030198 extracellular matrix organization
  394. # 3 BP Shared: Down in Hb, Down in Amyg GO:0030198 extracellular matrix organization
  397. go_kegg_results %>% dplyr::filter(Ontology == "CC", Description == "mitochondrial inner membrane") %>% .[,1:4]
  398. # Ontology DEGs_set ID Description
  399. # 1 CC Unique in Amyg - Up GO:0005743 mitochondrial inner membrane
  401. go_kegg_results %>% dplyr::filter(Ontology == "CC", Description == "exocytic vesicle") %>% .[,1:4]
  402. # Ontology DEGs_set ID Description
  403. # 1 CC Unique in Amyg - Up GO:0070382 exocytic vesicle
  405. go_kegg_results %>% dplyr::filter(Ontology == "CC", Description == "extracellular matrix") %>% .[,1:4]
  406. # Ontology DEGs_set ID Description
  407. # 1 CC Unique in Hb - Down GO:0031012 extracellular matrix
  408. # 2 CC Unique in Amyg - Down GO:0031012 extracellular matrix
  410. go_kegg_results %>% dplyr::filter(Ontology == "CC", Description == "myelin sheath") %>% .[,1:4]
  411. # Ontology DEGs_set ID Description
  412. # 1 CC Unique in Amyg - Down GO:0043209 myelin sheath
  413. # 2 CC Shared: Down in Hb, Down in Amyg GO:0043209 myelin sheath
  416. go_kegg_results %>% dplyr::filter(Ontology == "MF", Description == "structural constituent of ribosome") %>% .[,1:4]
  417. # Ontology DEGs_set ID Description
  418. # 1 MF Unique in Amyg - Up GO:0003735 structural constituent of ribosome
  420. go_kegg_results %>% dplyr::filter(Ontology == "MF", Description == "proton-transporting ATP synthase activity, rotational mechanism") %>% .[,1:4]
  421. # Ontology DEGs_set ID Description
  422. # 1 MF Unique in Amyg - Up GO:0046933 proton-transporting ATP synthase activity, rotational mechanism
  424. go_kegg_results %>% dplyr::filter(Ontology == "MF", Description == "cytoskeletal motor activity") %>% .[,1:4]
  425. # Ontology DEGs_set ID Description
  426. # 1 MF Unique in Amyg - Down GO:0003774 cytoskeletal motor activity
  428. go_kegg_results %>% dplyr::filter(Ontology == "MF", Description == "extracellular matrix structural constituent") %>% .[,1:4]
  429. # Ontology DEGs_set ID Description
  430. # 1 MF Unique in Amyg - Down GO:0005201 extracellular matrix structural constituent
  431. # 2 MF Shared: Up in Hb, Down in Amyg GO:0005201 extracellular matrix structural constituent
  434. go_kegg_results %>% dplyr::filter(Ontology == "KEGG", Description == "Oxidative phosphorylation") %>% .[,1:4]
  435. # Ontology DEGs_set ID Description
  436. # 1 KEGG Unique in Amyg - Up rno00190 Oxidative phosphorylation
  438. go_kegg_results %>% dplyr::filter(Ontology == "KEGG", Description == "Parkinson disease") %>% .[,1:4]
  439. # Ontology DEGs_set ID Description
  440. # 1 KEGG Unique in Amyg - Up rno05012 Parkinson disease
  442. go_kegg_results %>% dplyr::filter(Ontology == "KEGG", Description == "ECM-receptor interaction") %>% .[,1:4]
  443. # Ontology DEGs_set ID Description
  444. # 1 KEGG Unique in Hb - Down rno04512 ECM-receptor interaction
  445. # 2 KEGG Unique in Amyg - Down rno04512 ECM-receptor interaction
  446. # 3 KEGG Shared: Up in Hb, Down in Amyg rno04512 ECM-receptor interaction
  448. go_kegg_results %>% dplyr::filter(Ontology == "KEGG", Description == "Regulation of actin cytoskeleton") %>% .[,1:4]
  449. # Ontology DEGs_set ID Description
  450. # 1 KEGG Unique in Amyg - Up rno04810 Regulation of actin cytoskeleton
  451. # 2 KEGG Unique in Amyg - Down rno04810 Regulation of actin cytoskeleton
  454. enriched_terms <- list("only_up_hab" = c("regulation of synapse structure or activity", "regulation of presynaptic membrane potential",
  455. "regulation of postsynaptic membrane potential",
  456. "presynaptic membrane", "postsynaptic membrane", "voltage-gated monoatomic cation channel activity",
  457. "metal ion transmembrane transporter activity", "Morphine addiction", "Glutamatergic synapse",
  458. "Gastric acid secretion"),
  459. "only_down_hab" = c("extracellular matrix organization", "collagen-containing extracellular matrix", "Gastric acid secretion",
  460. "ECM-receptor interaction", "extracellular matrix"),
  461. "only_up_amy" = c("regulation of synapse structure or activity", "Glutamatergic synapse", "aerobic respiration",
  462. "vesicle-mediated transport in synapse", "mitochondrial inner membrane", "exocytic vesicle",
  463. "structural constituent of ribosome", "proton-transporting ATP synthase activity, rotational mechanism",
  464. "Oxidative phosphorylation", "Parkinson disease", "Regulation of actin cytoskeleton"),
  465. "only_down_amy" = c("extracellular matrix organization", "oligodendrocyte differentiation",
  466. "collagen-containing extracellular matrix", "myelin sheath", "ECM-receptor interaction",
  467. "glial cell differentiation", "extracellular matrix", "cytoskeletal motor activity",
  468. "extracellular matrix structural constituent", "Regulation of actin cytoskeleton"),
  469. "shared_up_hab_up_amy" = c("regulation of presynaptic membrane potential", "presynaptic membrane",
  470. "postsynaptic membrane", "voltage-gated monoatomic cation channel activity"),
  471. "shared_up_hab_down_amy" = c("Glutamatergic synapse", "Gastric acid secretion", "ECM-receptor interaction",
  472. "extracellular matrix structural constituent"),
  473. "shared_down_hab_up_amy" = c(),
  474. "shared_down_hab_down_amy" = c("extracellular matrix organization", "oligodendrocyte differentiation",
  475. "collagen-containing extracellular matrix", "myelin sheath", "glial cell differentiation")
  476. )
  478. ## Annotate DEGs of interest in each term of interest (Kristen and Robin list) in specific group(s) (up/down unique/shared in Hb/Amyg)
  480. only_up_hab_genes_2_show = list()
  481. only_down_hab_genes_2_show = list()
  482. only_up_amy_genes_2_show = list()
  483. only_down_amy_genes_2_show = list()
  485. shared_up_hab_up_amy_genes_2_show = list("regulation of presynaptic membrane potential" = c("Kcnj3", "Kcnc2", "Kcnj9", "Kctd16"),
  486. "presynaptic membrane" = c("Kcnj3", "Kcnc2", "Kcnj9", "Kctd16"),
  487. "postsynaptic membrane" = c("Kcnc2", "LRRTM1", "Epha4", "Lrrtm2"),
  488. "voltage-gated monoatomic cation channel activity" = c("Kcnj3", "Kcnc2", "Kcnj9"))
  490. shared_up_hab_down_amy_genes_2_show = list("Glutamatergic synapse",
  491. "Gastric acid secretion",
  492. "extracellular matrix structural constituent" = c("Col4a3"),
  493. "ECM-receptor interaction" = c("Col4a3"))
  495. shared_down_hab_down_amy_genes_2_show = list("extracellular matrix organization" = c("Tgfbi", "Antxr1", "Col9a3", "Loxl4", "Sox9"),
  496. "oligodendrocyte differentiation" = c("Cnp", "Sox8", "Gsn", "Sox9", "Opalin"),
  497. "collagen-containing extracellular matrix" = c("Tgfbi", "Col9a3", "Loxl4", "Fgfr2"),
  498. "myelin sheath" = c("Cnp", "Tubb4a", "Gsn", "Tspan2"),
  499. "glial cell differentiation" = c("Cnp", "Gsn", "Sox9", "Opalin", "Tspan2"))
  503. ## Show additional top 5 most signif DEGs in each group
  504. only_up_hab_genes_additional_top = list("regulation of synapse structure or activity" = c(),
  505. "regulation of presynaptic membrane potential" = c(),
  506. "regulation of postsynaptic membrane potential" = c(),
  507. "presynaptic membrane" = c(),
  508. "postsynaptic membrane" = c(),
  509. "voltage-gated monoatomic cation channel activity" = c(),
  510. "metal ion transmembrane transporter activity" = c(),
  511. "Morphine addiction" = c(),
  512. "Glutamatergic synapse" = c(),
  513. "Gastric acid secretion" = c())
  516. only_down_hab_genes_additional_top = list("extracellular matrix organization" = c(),
  517. "collagen-containing extracellular matrix" = c(),
  518. "Gastric acid secretion" = c(),
  519. "ECM-receptor interaction" = c(),
  520. "extracellular matrix" = c(),
  521. "ECM-receptor interaction" = c())
  523. only_up_amy_genes_additional_top = list("regulation of synapse structure or activity" = c(),
  524. "Glutamatergic synapse" = c(),
  525. "aerobic respiration" = c(),
  526. "vesicle-mediated transport in synapse" = c(),
  527. "mitochondrial inner membrane" = c(),
  528. "exocytic vesicle" = c(),
  529. "structural constituent of ribosome" = c(),
  530. "proton-transporting ATP synthase activity, rotational mechanism" = c(),
  531. "Oxidative phosphorylation" = c(),
  532. "Parkinson disease" = c(),
  533. "Regulation of actin cytoskeleton" = c())
  535. only_down_amy_genes_2_show = list("extracellular matrix organization",
  536. "oligodendrocyte differentiation",
  537. "collagen-containing extracellular matrix",
  538. "myelin sheath",
  539. "ECM-receptor interaction",
  540. "glial cell differentiation",
  541. "extracellular matrix",
  542. "cytoskeletal motor activity",
  543. "extracellular matrix structural constituent",
  544. "ECM-receptor interaction",
  545. "Regulation of actin cytoskeleton")
  548. ## Extract DEGs of interest terms across DEG groups (include interest DEGs from Kristen and Robin list)
  549. genes_2_show_x_term_x_group <- list(list(), list(), list(), list(), list(), list(), list(), list())
  550. names(genes_2_show_x_term_x_group) <- unique(go_kegg_results$DEGs_set)
  552. ## Specific DEGs to highlight
  553. interest_genes <- unique(c("Ptpn3", "Kcnc2", "Kcnj9", "Kcnj3", "Kctd16",
  554. "LRRTM1", "Lrrtm1", "Epha4", "Lrrn3", "Fam107a", "LRRTM2", "Lrrtm2",
  555. "Loxl4", "Antxr1", "Tgfbi", "Col9a3", "Sox9",
  556. "Cnp", "Gsn", "Sox8", "Sox9", "Opalin",
  557. "Kcnc2", "Npy1r", "Kcnj9", "Epha4", "Kcnj3",
  558. "Kcnc2", "LRRTM1", "Epha4", "Kctd16", "Cdh10",
  559. "Loxl4", "Col9a3", "Tgfbi", "Fgfr2",
  560. "Cnp", "Gsn", "Tspan2", "Tubb4a",
  561. "Kcnc2", "Kcnj9", "Kcnj3",
  562. "Slc13a5", "Kcnc2", "Kcnj9", "Kcnj3",
  563. "Col9a3",
  564. "Ap3s1",
  565. "Col4a3"))
  567. ## Terms of interest
  568. terms <- c("regulation of presynaptic membrane potential", "regulation of postsynaptic membrane potential",
  569. "regulation of synapse structure or activity", "extracellular matrix organization", "oligodendrocyte differentiation",
  570. "presynaptic membrane", "postsynaptic membrane", "collagen-containing extracellular matrix", "myelin sheath",
  571. "voltage-gated monoatomic cation channel activity", "metal ion transmembrane transporter activity",
  572. "Morphine addiction", "Glutamatergic synapse", "Gastric acid secretion", "ECM-receptor interaction",
  573. "aerobic respiration", "vesicle-mediated transport in synapse", "glial cell differentiation", "extracellular matrix organization",
  574. "mitochondrial inner membrane", "exocytic vesicle", "extracellular matrix", "structural constituent of ribosome",
  575. "proton-transporting ATP synthase activity, rotational mechanism", "cytoskeletal motor activity",
  576. "extracellular matrix structural constituent", "Oxidative phosphorylation", "Parkinson disease",
  577. "Regulation of actin cytoskeleton")
  579. BP_terms <- c("regulation of presynaptic membrane potential", "regulation of postsynaptic membrane potential",
  580. "regulation of synapse structure or activity", "extracellular matrix organization",
  581. "oligodendrocyte differentiation", "aerobic respiration", "vesicle-mediated transport in synapse",
  582. "glial cell differentiation", "extracellular matrix organization")
  583. CC_terms <- c("presynaptic membrane", "postsynaptic membrane", "collagen-containing extracellular matrix", "myelin sheath", "mitochondrial inner membrane", "exocytic vesicle", "extracellular matrix")
  584. MF_terms <- c("voltage-gated monoatomic cation channel activity", "metal ion transmembrane transporter activity",
  585. "structural constituent of ribosome",
  586. "proton-transporting ATP synthase activity, rotational mechanism",
  587. "cytoskeletal motor activity", "extracellular matrix structural constituent")
  588. KEGG_terms <- c("Morphine addiction", "Glutamatergic synapse", "Gastric acid secretion", "ECM-receptor interaction",
  589. "Oxidative phosphorylation", "Parkinson disease",
  590. "Regulation of actin cytoskeleton")
  592. for(term in terms){
  594. ## Extract groups of DEGs where term is enriched
  595. DEG_groups_with_term <- go_kegg_results %>% dplyr::filter(Description == term) %>% pull(DEGs_set)
  597. for(group in DEG_groups_with_term){
  599. ## Genes in term and group
  600. intersection_genes <- strsplit(go_kegg_results %>% dplyr::filter(Description == term & DEGs_set == group) %>% pull(geneID), ", ") %>% unlist
  601. ## Subset to the specific genes of interest
  602. genes_of_interest_2_show <- intersect(intersection_genes, interest_genes)
  604. if(length(genes_of_interest_2_show) >= 3){
  605. genes_2_show_x_term_x_group[[group]][[term]] = genes_of_interest_2_show
  606. }
  608. else{
  609. ## Add top n most signif genes in group and term
  610. n = 3 - length(genes_of_interest_2_show)
  612. if(group == "Unique in Hb - Up"){
  613. top3_intersection_genes <- only_up_hab %>% dplyr::filter(Symbol %in% intersection_genes, !Symbol %in% interest_genes) %>% arrange(adj.P.Val) %>% slice(1:n) %>% pull(Symbol)
  614. }
  615. else if(group == "Unique in Hb - Down"){
  616. top3_intersection_genes <- only_down_hab %>% dplyr::filter(Symbol %in% intersection_genes, !Symbol %in% interest_genes) %>% arrange(adj.P.Val) %>% slice(1:n) %>% pull(Symbol)
  617. }
  618. else if(group == "Unique in Amyg - Up"){
  619. top3_intersection_genes <- only_up_amy %>% dplyr::filter(Symbol %in% intersection_genes, !Symbol %in% interest_genes) %>% arrange(adj.P.Val) %>% slice(1:n) %>% pull(Symbol)
  620. }
  621. else if(group == "Unique in Amyg - Down"){
  622. top3_intersection_genes <- only_down_amy %>% dplyr::filter(Symbol %in% intersection_genes, !Symbol %in% interest_genes) %>% arrange(adj.P.Val) %>% slice(1:n) %>% pull(Symbol)
  623. }
  624. else if(group == "Shared: Up in Hb, Up in Amyg"){
  625. top3_intersection_genes <- shared_up_hab_up_amy %>% dplyr::filter(Symbol %in% intersection_genes, !Symbol %in% interest_genes) %>%
  626. mutate(min_p = pmin(adj.P.Val.hb, adj.P.Val.amyg)) %>% arrange(min_p) %>% slice(1:n) %>% pull(Symbol)
  627. }
  628. else if(group == "Shared: Up in Hb, Down in Amyg"){
  629. top3_intersection_genes <- shared_up_hab_down_amy %>% dplyr::filter(Symbol %in% intersection_genes, !Symbol %in% interest_genes) %>%
  630. mutate(min_p = pmin(adj.P.Val.hb, adj.P.Val.amyg)) %>% arrange(min_p) %>% slice(1:n) %>% pull(Symbol)
  631. }
  632. else if(group == "Shared: Down in Hb, Up in Amyg"){
  633. top3_intersection_genes <- shared_down_hab_up_amy %>% dplyr::filter(Symbol %in% intersection_genes, !Symbol %in% interest_genes) %>%
  634. mutate(min_p = pmin(adj.P.Val.hb, adj.P.Val.amyg)) %>% arrange(min_p) %>% slice(1:n) %>% pull(Symbol)
  635. }
  636. else if(group == "Shared: Down in Hb, Down in Amyg"){
  637. top3_intersection_genes <- shared_down_hab_down_amy %>% dplyr::filter(Symbol %in% intersection_genes, !Symbol %in% interest_genes) %>%
  638. mutate(min_p = pmin(adj.P.Val.hb, adj.P.Val.amyg)) %>% arrange(min_p) %>% slice(1:n) %>% pull(Symbol)
  639. }
  641. genes_2_show_x_term_x_group[[group]][[term]] = c(genes_of_interest_2_show, top3_intersection_genes)
  642. }
  644. }
  645. }
  647. l <- sapply(names(genes_2_show_x_term_x_group), function(set){
  648. imap_dfr(genes_2_show_x_term_x_group[[set]], ~ enframe(.x, value = "Symbol") %>% mutate(term = .y, DEGs_set = set))
  649. })
  651. l <- do.call(rbind, l)
  654. ## Add logFC of selected genes in Hb and Amyg
  655. df <- l %>% left_join(results_Substance_uncorr_vars_habenula[[1]][, c("Symbol", "logFC", "t", "P.Value", "adj.P.Val")], by = "Symbol", multiple = "any") %>% left_join(results_Substance_uncorr_vars_amygdala[[1]][, c("Symbol", "logFC", "t", "P.Value", "adj.P.Val")], by = "Symbol", multiple = "any", suffix = c(".Hb", ".Amy"))
  657. df$DEGs_set <- factor(df$DEGs_set, levels = unique(df$DEGs_set))
  659. ## Order terms by alp order x group
  660. df <- df %>% arrange(DEGs_set, term) %>% as.data.frame()
  662. num_genes_x_term_x_group <- df %>% group_by(DEGs_set, term) %>% summarise(count = length(unique(Symbol)))
  663. term_indices <- map2(.x = c(1, head(cumsum(num_genes_x_term_x_group$count), -1)+1),
  664. .y = cumsum(num_genes_x_term_x_group$count),
  665. ~(1:dim(df)[1])[.x:.y])
  666. names(term_indices) <- num_genes_x_term_x_group$term
  668. la = rowAnnotation(Group = df$DEGs_set)
  669. ra = rowAnnotation(term = anno_block(align_to = term_indices,
  670. panel_fun = function(index, nm){
  671. grid.text(nm, rot = 0, just = "left", name = "term",
  672. gp = gpar(fontsize = 7), x = 3.2)}))
  674. h <- Heatmap(as.matrix(df[, c("logFC.Hb", "logFC.Amy")]),
  675. name = "logFC",
  676. border = T,
  677. row_labels = df$Symbol,
  678. row_names_side = "right",
  679. column_labels = c("logFC in Hb", "logFC in Amyg"),
  680. column_names_rot = 90,
  681. column_names_centered = F,
  682. column_names_gp = gpar(fontsize = 7),
  683. row_names_gp = gpar(fontsize = 5, fontface = "italic"),
  684. row_split = df[, c("DEGs_set", "term")],
  685. gap = unit(0.75, "mm"),
  686. left_annotation = la,
  687. right_annotation = ra,
  688. row_title_gp = gpar(fontsize = 0),
  689. cluster_rows = FALSE,
  690. cluster_columns = FALSE,
  691. height = unit(35, "cm"),
  692. heatmap_width = unit(5, "cm"))
  695. pdf(file = paste0("plots/06_GO_KEGG/GO_KEGG_heatmap_Hb_vs_Amyg.pdf"), height = 15, width = 10)
  696. draw(h, heatmap_legend_side = "left")
  697. dev.off()
  700. ## Tile plot
  701. df_wide <- l %>% dplyr::select(Symbol, DEGs_set, term) %>%
  702. mutate(present = 1) %>%
  703. pivot_wider(
  704. names_from = term,
  705. values_from = present,
  706. values_fill = 0
  707. )
  709. df_longer <- df_wide %>% pivot_longer(cols = setdiff(colnames(df_wide), c("Symbol", "DEGs_set")))
  711. df_longer$DEGs_set <- factor(df_longer$DEGs_set, levels = unique(l$DEGs_set))
  712. df_longer <- df_longer %>% arrange(DEGs_set, name)
  713. df_longer$Symbol <- factor(df_longer$Symbol, levels = unique(df_longer$Symbol))
  714. df_longer$name <- factor(df_longer$name, levels = rev(unique(df_longer$name)))
  716. df_longer <- df_longer %>%
  717. mutate("Ontology" = case_when(name %in% BP_terms ~ "BP",
  718. name %in% CC_terms ~ "CC",
  719. name %in% MF_terms ~ "MF",
  720. name %in% KEGG_terms ~ "KEGG"))
  722. df_longer$Ontology <- factor(df_longer$Ontology, levels = c("BP", "CC", "MF", "KEGG"))
  724. ## Add logFC in Hb and Amyg
  725. df_logFCs <- results_Substance_uncorr_vars_habenula[[1]] %>%
  726. subset(Symbol %in% df_longer$Symbol) %>%
  727. dplyr::select(Symbol, logFC) %>%
  728. left_join(subset(results_Substance_uncorr_vars_amygdala[[1]],
  729. results_Substance_uncorr_vars_amygdala[[1]]$Symbol %in% df_longer$Symbol)[, c("Symbol", "logFC")],
  730. by = "Symbol", multiple = "any", suffix = c(".Hb", ".Amyg"))
  732. df_logFCs$logFC.Hb <- df_logFCs$logFC.Hb - mean(df_logFCs$logFC.Hb)
  733. df_logFCs$logFC.Amyg <- df_logFCs$logFC.Amyg - mean(df_logFCs$logFC.Amyg)
  734. df_logFCs <- df_logFCs[-which(duplicated(df_logFCs$Symbol)),]
  736. df_logFCs <- pivot_longer(df_logFCs, cols = c("logFC.Hb", "logFC.Amyg")) %>%
  737. left_join(unique(df_longer[, c("Symbol", "DEGs_set")]), by = "Symbol") %>%
  738. mutate(Ontology = case_when(name == "logFC.Hb" ~ "LogFC in Hb",
  739. name == "logFC.Amyg" ~ "LogFC in Amyg"))
  742. df_longer <- rbind(df_longer, df_logFCs[, colnames(df_longer)])
  744. noKEGGgenes <- df_longer %>%
  745. dplyr::filter(Ontology %in% c("BP", "CC" ,"MF")) %>%
  746. dplyr::group_by(Symbol) %>%
  747. summarise(sum(value)) %>%
  748. .[.[,2]> 0, ] %>%
  749. dplyr::select(Symbol) %>% unlist()
  751. df_longer_noKEGG <- df_longer %>% dplyr::filter(Symbol %in% noKEGGgenes)
  752. toplot <- df_longer_noKEGG %>% dplyr::filter(Ontology %in% c("BP", "CC" ,"MF"))
  753. bottomplot <- df_longer_noKEGG %>% dplyr::filter(Ontology %in% c("LogFC in Hb", "LogFC in Amyg"))
  755. p1 <- ggplot(toplot, aes(x = Symbol, y = name, fill = value)) +
  756. geom_tile(color = "gray80", linewidth = 0.005) +
  757. facet_grid(rows = vars(Ontology), cols = vars(DEGs_set),
  758. scales = "free", space = "free") +
  759. scale_fill_gradient(low = "white", high = "gray40") +
  760. guides(fill = "none") +
  761. theme_bw() +
  762. labs(y = "Enriched term") +
  763. theme(axis.text.x = element_text(size = 6, angle = 90, hjust = 1, face = 3),
  764. axis.text.y = element_text(size = 7),
  765. strip.text.x = element_blank(),
  766. strip.background = element_rect(fill="white", color = "white"),
  767. panel.spacing = unit(0.1, "lines"))
  770. p2 <- ggplot(bottomplot, aes(x = Symbol, y = name, fill = value)) +
  771. geom_tile(color = "gray80", linewidth = 0.00) +
  772. facet_grid(cols = vars(DEGs_set),
  773. scale = "free", space = "free") +
  774. scale_x_discrete(expand=c(0,0)) +
  775. scale_y_discrete(expand=c(0,0)) +
  776. scale_fill_gradient2(low = "blue3",
  777. mid = "gray90",
  778. high = "red3") +
  779. theme(axis.title.x = element_blank(),
  780. axis.text.x = element_blank(),
  781. axis.ticks.x = element_blank(),
  782. axis.text.y = element_text(size = 7),
  783. strip.text.x = element_text(size = 8, angle = 90, hjust = 0, face = 3),
  784. strip.background = element_rect(fill="white", color = "white"),
  785. panel.spacing = unit(0.1, "lines"))
  787. plot_grid(p2, p1, nrow = 2, align = "v", axis = "tblr", rel_heights = c(0.7, 1))
  789. ggsave("plots/06_GO_KEGG/GO_tile_Hb_vs_Amyg.pdf", height = 6, width = 10)
  792. ## Supp table for KEGG
  793. KEGGgenes <- df_longer %>%
  794. dplyr::filter(Ontology %in% c("KEGG")) %>%
  795. dplyr::group_by(Symbol) %>%
  796. summarise(sum(value)) %>%
  797. .[.[,2]> 0, ] %>%
  798. dplyr::select(Symbol) %>% unlist()
  800. df_longer_KEGGgenes <- df_longer %>% dplyr::filter(Symbol %in% KEGGgenes)
  801. toplot <- df_longer_KEGGgenes %>% dplyr::filter(Ontology %in% c("KEGG"))
  802. bottomplot <- df_longer_KEGGgenes %>% dplyr::filter(Ontology %in% c("LogFC in Hb", "LogFC in Amyg"))
  804. p1 <- ggplot(toplot, aes(x = Symbol, y = name, fill = value)) +
  805. geom_tile(color = "gray80", linewidth = 0.005) +
  806. facet_grid(rows = vars(Ontology), cols = vars(DEGs_set),
  807. scales = "free", space = "free") +
  808. scale_fill_gradient(low = "white", high = "gray40") +
  809. guides(fill = "none") +
  810. theme_bw() +
  811. labs(y = "Enriched term") +
  812. theme(axis.text.x = element_text(size = 6, angle = 90, hjust = 1, face = 3),
  813. axis.text.y = element_text(size = 7),
  814. strip.text.x = element_blank(),
  815. strip.background = element_rect(fill="white", color = "white"),
  816. panel.spacing = unit(0.1, "lines"))
  818. p2 <- ggplot(bottomplot, aes(x = Symbol, y = name, fill = value)) +
  819. geom_tile(color = "gray80", linewidth = 0.00) +
  820. facet_grid(cols = vars(DEGs_set),
  821. scale = "free", space = "free") +
  822. scale_x_discrete(expand=c(0,0)) +
  823. scale_y_discrete(expand=c(0,0)) +
  824. scale_fill_gradient2(low = "blue3",
  825. mid = "gray90",
  826. high = "red3") +
  827. theme(axis.title.x = element_blank(),
  828. axis.text.x = element_blank(),
  829. axis.ticks.x = element_blank(),
  830. axis.text.y = element_text(size = 7),
  831. strip.text.x = element_text(size = 8, angle = 90, hjust = 0, face = 3),
  832. strip.background = element_rect(fill="white", color = "white"),
  833. panel.spacing = unit(0.1, "lines"))
  835. plot_grid(p2, p1, nrow = 2, align = "v", axis = "tblr", rel_heights = c(0.9, 1))
  837. ggsave("plots/06_GO_KEGG/KEGG_tile_Hb_vs_Amyg.pdf", height = 5, width = 6)
  846. ## Reproducibility information
  848. options(width = 120)
  849. session_info()
  851. # ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
  852. # setting value
  853. # version R version 4.3.2 (2023-10-31)
  854. # os macOS Monterey 12.5.1
  855. # system aarch64, darwin20
  856. # ui RStudio
  857. # language (EN)
  858. # collate en_US.UTF-8
  859. # ctype en_US.UTF-8
  860. # tz America/Mexico_City
  861. # date 2024-04-23
  862. # rstudio 2023.12.1+402 Ocean Storm (desktop)
  863. # pandoc 3.1.1 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)
  864. #
  865. # ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
  866. # package * version date (UTC) lib source
  867. # abind 1.4-5 2016-07-21 [1] CRAN (R 4.3.0)
  868. # AnnotationDbi * 1.64.1 2023-11-02 [1] Bioconductor
  869. # AnnotationHub 3.10.0 2023-10-26 [1] Bioconductor
  870. # ape 5.7-1 2023-03-13 [1] CRAN (R 4.3.0)
  871. # aplot 0.2.2 2023-10-06 [1] CRAN (R 4.3.1)
  872. # Biobase * 2.62.0 2023-10-26 [1] Bioconductor
  873. # BiocFileCache 2.10.1 2023-10-26 [1] Bioconductor
  874. # BiocGenerics * 0.48.1 2023-11-02 [1] Bioconductor
  875. # BiocManager 1.30.22 2023-08-08 [1] CRAN (R 4.3.0)
  876. # BiocParallel 1.36.0 2023-10-26 [1] Bioconductor
  877. # BiocVersion 3.18.1 2023-11-18 [1] Bioconductor 3.18 (R 4.3.2)
  878. # Biostrings 2.70.2 2024-01-30 [1] Bioconductor 3.18 (R 4.3.2)
  879. # bit 4.0.5 2022-11-15 [1] CRAN (R 4.3.0)
  880. # bit64 4.0.5 2020-08-30 [1] CRAN (R 4.3.0)
  881. # bitops 1.0-7 2021-04-24 [1] CRAN (R 4.3.0)
  882. # blob 1.2.4 2023-03-17 [1] CRAN (R 4.3.0)
  883. # cachem 1.0.8 2023-05-01 [1] CRAN (R 4.3.0)
  884. # cli 3.6.2 2023-12-11 [1] CRAN (R 4.3.1)
  885. # clusterProfiler * 4.10.0 2023-11-06 [1] Bioconductor
  886. # codetools 0.2-19 2023-02-01 [1] CRAN (R 4.3.2)
  887. # colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.3.0)
  888. # cowplot 1.1.3 2024-01-22 [1] CRAN (R 4.3.1)
  889. # crayon 1.5.2 2022-09-29 [1] CRAN (R 4.3.0)
  890. # curl 5.2.1 2024-03-01 [1] CRAN (R 4.3.1)
  891. # data.table 1.15.2 2024-02-29 [1] CRAN (R 4.3.1)
  892. # DBI 1.2.2 2024-02-16 [1] CRAN (R 4.3.2)
  893. # dbplyr 2.4.0 2023-10-26 [1] CRAN (R 4.3.1)
  894. # DelayedArray 0.28.0 2023-11-06 [1] Bioconductor
  895. # digest 0.6.34 2024-01-11 [1] CRAN (R 4.3.1)
  896. # DOSE 3.28.2 2023-12-12 [1] Bioconductor 3.18 (R 4.3.2)
  897. # dplyr 1.1.4 2023-11-17 [1] CRAN (R 4.3.1)
  898. # ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.3.0)
  899. # enrichplot 1.22.0 2023-11-06 [1] Bioconductor
  900. # evaluate 0.23 2023-11-01 [1] CRAN (R 4.3.1)
  901. # fansi 1.0.6 2023-12-08 [1] CRAN (R 4.3.1)
  902. # farver 2.1.1 2022-07-06 [1] CRAN (R 4.3.0)
  903. # fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.3.0)
  904. # fastmatch 1.1-4 2023-08-18 [1] CRAN (R 4.3.0)
  905. # fgsea 1.28.0 2023-10-26 [1] Bioconductor
  906. # filelock 1.0.3 2023-12-11 [1] CRAN (R 4.3.1)
  907. # fs 1.6.3 2023-07-20 [1] CRAN (R 4.3.0)
  908. # generics 0.1.3 2022-07-05 [1] CRAN (R 4.3.0)
  909. # GenomeInfoDb * 1.38.6 2024-02-10 [1] Bioconductor 3.18 (R 4.3.2)
  910. # GenomeInfoDbData 1.2.11 2024-02-17 [1] Bioconductor
  911. # GenomicRanges * 1.54.1 2023-10-30 [1] Bioconductor
  912. # ggforce 0.4.2 2024-02-19 [1] CRAN (R 4.3.1)
  913. # ggfun 0.1.4 2024-01-19 [1] CRAN (R 4.3.1)
  914. # ggplot2 3.5.0 2024-02-23 [1] CRAN (R 4.3.1)
  915. # ggplotify 0.1.2 2023-08-09 [1] CRAN (R 4.3.0)
  916. # ggraph 2.2.0 2024-02-27 [1] CRAN (R 4.3.1)
  917. # ggrepel 0.9.5 2024-01-10 [1] CRAN (R 4.3.1)
  918. # ggtree 3.10.1 2024-02-27 [1] Bioconductor 3.18 (R 4.3.2)
  919. # glue 1.7.0 2024-01-09 [1] CRAN (R 4.3.1)
  920. # GO.db 3.18.0 2024-02-17 [1] Bioconductor
  921. # GOSemSim 2.28.1 2024-01-20 [1] Bioconductor 3.18 (R 4.3.2)
  922. # graphlayouts 1.1.0 2024-01-19 [1] CRAN (R 4.3.1)
  923. # gridExtra 2.3 2017-09-09 [1] CRAN (R 4.3.0)
  924. # gridGraphics 0.5-1 2020-12-13 [1] CRAN (R 4.3.0)
  925. # gson 0.1.0 2023-03-07 [1] CRAN (R 4.3.0)
  926. # gtable 0.3.4 2023-08-21 [1] CRAN (R 4.3.0)
  927. # HDO.db 0.99.1 2023-05-28 [1] Bioconductor
  928. # here * 1.0.1 2020-12-13 [1] CRAN (R 4.3.0)
  929. # htmltools 0.5.7 2023-11-03 [1] CRAN (R 4.3.1)
  930. # httpuv 1.6.14 2024-01-26 [1] CRAN (R 4.3.1)
  931. # httr 1.4.7 2023-08-15 [1] CRAN (R 4.3.0)
  932. # igraph 2.0.2 2024-02-17 [1] CRAN (R 4.3.1)
  933. # interactiveDisplayBase 1.40.0 2023-10-26 [1] Bioconductor
  934. # IRanges * 2.36.0 2023-10-26 [1] Bioconductor
  935. # jsonlite 1.8.8 2023-12-04 [1] CRAN (R 4.3.1)
  936. # KEGGREST 1.42.0 2023-10-26 [1] Bioconductor
  937. # knitr 1.45 2023-10-30 [1] CRAN (R 4.3.1)
  938. # labeling 0.4.3 2023-08-29 [1] CRAN (R 4.3.0)
  939. # later 1.3.2 2023-12-06 [1] CRAN (R 4.3.1)
  940. # lattice 0.22-5 2023-10-24 [1] CRAN (R 4.3.1)
  941. # lazyeval 0.2.2 2019-03-15 [1] CRAN (R 4.3.0)
  942. # lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.3.1)
  943. # magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.3.0)
  944. # MASS 7.3-60.0.1 2024-01-13 [1] CRAN (R 4.3.1)
  945. # Matrix 1.6-5 2024-01-11 [1] CRAN (R 4.3.1)
  946. # MatrixGenerics * 1.14.0 2023-10-26 [1] Bioconductor
  947. # matrixStats * 1.2.0 2023-12-11 [1] CRAN (R 4.3.1)
  948. # memoise 2.0.1 2021-11-26 [1] CRAN (R 4.3.0)
  949. # mime 0.12 2021-09-28 [1] CRAN (R 4.3.0)
  950. # munsell 0.5.0 2018-06-12 [1] CRAN (R 4.3.0)
  951. # nlme 3.1-164 2023-11-27 [1] CRAN (R 4.3.1)
  952. # org.Rn.eg.db * 3.18.0 2024-02-17 [1] Bioconductor
  953. # patchwork 1.2.0 2024-01-08 [1] CRAN (R 4.3.1)
  954. # pillar 1.9.0 2023-03-22 [1] CRAN (R 4.3.0)
  955. # pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.0)
  956. # plyr 1.8.9 2023-10-02 [1] CRAN (R 4.3.1)
  957. # png 0.1-8 2022-11-29 [1] CRAN (R 4.3.0)
  958. # polyclip 1.10-6 2023-09-27 [1] CRAN (R 4.3.1)
  959. # promises 1.2.1 2023-08-10 [1] CRAN (R 4.3.0)
  960. # purrr 1.0.2 2023-08-10 [1] CRAN (R 4.3.0)
  961. # qvalue 2.34.0 2023-10-26 [1] Bioconductor
  962. # R6 2.5.1 2021-08-19 [1] CRAN (R 4.3.0)
  963. # rappdirs 0.3.3 2021-01-31 [1] CRAN (R 4.3.0)
  964. # RColorBrewer 1.1-3 2022-04-03 [1] CRAN (R 4.3.0)
  965. # Rcpp 1.0.12 2024-01-09 [1] CRAN (R 4.3.1)
  966. # RCurl 1.98-1.14 2024-01-09 [1] CRAN (R 4.3.1)
  967. # reshape2 1.4.4 2020-04-09 [1] CRAN (R 4.3.0)
  968. # rlang 1.1.3 2024-01-10 [1] CRAN (R 4.3.1)
  969. # rmarkdown 2.26 2024-03-05 [1] CRAN (R 4.3.1)
  970. # rprojroot 2.0.4 2023-11-05 [1] CRAN (R 4.3.1)
  971. # RSQLite 2.3.5 2024-01-21 [1] CRAN (R 4.3.1)
  972. # rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.3.0)
  973. # S4Arrays 1.2.0 2023-10-26 [1] Bioconductor
  974. # S4Vectors * 0.40.2 2023-11-25 [1] Bioconductor 3.18 (R 4.3.2)
  975. # scales 1.3.0 2023-11-28 [1] CRAN (R 4.3.1)
  976. # scatterpie 0.2.1 2023-06-07 [1] CRAN (R 4.3.0)
  977. # sessioninfo * 1.2.2 2021-12-06 [1] CRAN (R 4.3.0)
  978. # shadowtext 0.1.3 2024-01-19 [1] CRAN (R 4.3.1)
  979. # shiny 1.8.0 2023-11-17 [1] CRAN (R 4.3.1)
  980. # SparseArray 1.2.4 2024-02-10 [1] Bioconductor 3.18 (R 4.3.2)
  981. # stringi 1.8.3 2023-12-11 [1] CRAN (R 4.3.1)
  982. # stringr 1.5.1 2023-11-14 [1] CRAN (R 4.3.1)
  983. # SummarizedExperiment * 1.32.0 2023-11-06 [1] Bioconductor
  984. # tibble 3.2.1 2023-03-20 [1] CRAN (R 4.3.0)
  985. # tidygraph 1.3.1 2024-01-30 [1] CRAN (R 4.3.1)
  986. # tidyr 1.3.1 2024-01-24 [1] CRAN (R 4.3.1)
  987. # tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.3.0)
  988. # tidytree 0.4.6 2023-12-12 [1] CRAN (R 4.3.1)
  989. # treeio 1.26.0 2023-11-06 [1] Bioconductor
  990. # tweenr 2.0.3 2024-02-26 [1] CRAN (R 4.3.1)
  991. # utf8 1.2.4 2023-10-22 [1] CRAN (R 4.3.1)
  992. # vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.3.1)
  993. # viridis 0.6.5 2024-01-29 [1] CRAN (R 4.3.1)
  994. # viridisLite 0.4.2 2023-05-02 [1] CRAN (R 4.3.0)
  995. # withr 3.0.0 2024-01-16 [1] CRAN (R 4.3.1)
  996. # xfun 0.42 2024-02-08 [1] CRAN (R 4.3.1)
  997. # xtable 1.8-4 2019-04-21 [1] CRAN (R 4.3.0)
  998. # XVector 0.42.0 2023-10-26 [1] Bioconductor
  999. # yaml 2.3.8 2023-12-11 [1] CRAN (R 4.3.1)
  1000. # yulab.utils 0.1.4 2024-01-28 [1] CRAN (R 4.3.1)
  1001. # zlibbioc 1.48.0 2023-10-26 [1] Bioconductor
  1002. #
  1003. # [1] /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library
  1004. #
  1005. # ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

01_GO_KEGG_Analyses.R at commit 4284eba, under MIT · at the source

Overview

Authors: Robin Magnard1, Daianna Gonzalez‐Padilla2, Ege A Yalcinbas2,3, Emma Chaloux‐Pinette4,5, Nicholas J Eagles2, Michael S Totty6, Patricia H Janak1,4,7, Leonardo Collado‐Torres2,6,8, Kristen R Maynard2,3,4
ORCID iDs: Robin Magnard, Daianna Gonzalez‐Padilla, Ege A Yalcinbas, Emma Chaloux‐Pinette, Nicholas J Eagles, Michael S Totty, Patricia H Janak, Leonardo Collado‐Torres, Kristen R Maynard
  1. Department of Psychological and Brain Sciences, Krieger School of Arts and Sciences, Johns Hopkins University, Baltimore, Maryland, USA
  2. Lieber Institute for Brain Development, Johns Hopkins Medical Campus, Baltimore, Maryland, USA
  3. Department of Psychiatry and Behavioral Sciences, Johns Hopkins School of Medicine, Baltimore, Maryland, USA
  4. The Solomon H. Snyder Department of Neuroscience, Johns Hopkins School of Medicine, Baltimore, Maryland, USA
  5. Medical Scientist Training Program, Johns Hopkins University School of Medicine, Baltimore, Maryland, USA
  6. Department of Biostatistics, Johns Hopkins Bloomberg School of Public Health, Baltimore, Maryland, USA
  7. Kavli Neuroscience Discovery Institute, Johns Hopkins University, Baltimore, Maryland, USA
  8. Center for Computational Biology, Johns Hopkins University, Baltimore, Maryland, USA
Institutions: Johns Hopkins University (United States); Lieber Institute for Brain Development (United States); Johns Hopkins Medicine (United States)
Journal: Addiction biology, volume 31, issue 7, article e70179
Dates: received 25 January 2026; accepted 24 June 2026; published online 14 July 2026; in print July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1111/adb.70179 · PMID 42444546 · PMCID PMC13366401 · OpenAlex W4417361577
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.11.25.690517
Status: code verified
Categories: genetics / omics (modality), human (organism), rat (organism), other condition (population), pain (population), cellular / molecular (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions
MeSH: Amygdala*, Analgesics, Opioid*, Fentanyl*, Habenula*, Transcriptome*, Animals, Gene Expression Profiling, Male, Rats, Rats, Sprague-Dawley, Self Administration (* major topic)
Topic: Neurogenesis and neuroplasticity mechanisms (Developmental Neuroscience, Neuroscience), according to OpenAlex
Funding: NIDA NIH HHS (R01 DA035943, R21 DA060407); NIH HHS (R01DA035943, T32MH015330, R21DA060407); National Institutes of Health (T32MH015330, R01DA035943, R21DA060407); NIMH NIH HHS (T32 MH015330); Lieber Institute for Brain Development
Citations: cited by 1 paper (Europe PMC); 104 references in the paper

Abstract

Fentanyl is a potent synthetic opioid associated with overdose. However, little is known about fentanyl‐induced molecular adaptations in the habenula and amygdala, two brain regions implicated in opioid use and withdrawal. We performed bulk RNA‐sequencing in the rat habenula and amygdala to identify transcriptomic changes associated with fentanyl intake. Male rats self‐administered intravenous saline or fentanyl over 22–24 days. Ninety minutes following the final session, brains were collected for transcriptomic profiling. In Hb, we identified 453 differentially expressed genes (DEGs) between saline and fentanyl rats, with upregulated genes associated with synaptic transmission and ionic conductance. In the amygdala, we identified 3041 fentanyl‐associated DEGs with upregulated genes implicated in metabolic and vesicular functions. Downregulated genes in both regions were enriched for extracellular matrix functions. Integration of DEGs with single‐cell RNA‐sequencing data from rodents and humans revealed that fentanyl DEGs were enriched in specific habenula and amygdala cell type markers. Furthermore, fentanyl downregulated DEGs in the amygdala were enriched in genes associated with the risk for substance use disorders. Together, we define how fentanyl intake alters transcriptional programs in the rat habenula and amygdala, and we link these changes to specific human cell types and risk genes for neuropsychiatric disorders and addiction.

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

Repositories

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

LieberInstitute/fentanyl_rat_hb_amy

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 4284ebac39e53fe23aea2a23812dc71441c14e10, 15 July 2026
Languages: R (17), Shell (9)
Size: 489 files, 26 scripts
Software Heritage: not archived
Found in: “Data Availability Statement”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: cowplot (9 files), ggplot2 (9 files), tidyverse (9 files), ComplexHeatmap (3 files), edgeR (2 files), pheatmap (2 files), circlize (1 file), clusterProfiler (1 file), data.table (1 file), limma (1 file), lme4 (1 file), Nextflow (1 file), reshape2 (1 file), rstatix (1 file), Seurat (1 file), SingleCellExperiment (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
27 files
At the source: github.com/LieberInstitute/fentanyl_rat_hb_amy

Zenodo 17573970

License: CC-BY-4.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
At the source:

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:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 26 scripts, each with its path and the digest of its content;
  • 12 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 source FASTQ files are publicly available from the NCBI Sequence Read Archive BioProject PRJNA1179901. All analysis code is available at https://github.com/LieberInstitute/fentanyl_rat_hb_amy [104].

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 2, 28 September 2026

  • Authors: added Ege A Yalcinbas (0000-0002-9480-7192); Emma Chaloux‐Pinette (0000-0002-6832-6862); Nicholas J Eagles (0000-0002-9808-5254); Michael S Totty (0000-0002-9292-8556); Patricia H Janak (0000-0002-3333-9049); Kristen R Maynard (0000-0003-0031-8468); removed Ege A Yalcinbas; Emma Chaloux‐Pinette; Nicholas J Eagles; Michael S Totty; Patricia H Janak; Kristen R Maynard

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 9 authors, 11 MeSH terms, 5 funders, 103 references.

Cite

This paper

Magnard, R., Gonzalez‐Padilla, D., Yalcinbas, E. A., Chaloux‐Pinette, E., Eagles, N. J., Totty, M. S., Janak, P. H., Collado‐Torres, L., & Maynard, K. R. (2026). Transcriptional Response to Chronic Long-Access Fentanyl Self-Administration in Rat Habenula and Amygdala. Addiction biology, 31(7), e70179. https://doi.org/10.1111/adb.70179

BibTeX

@article{magnard2026transcriptional,
author = {Magnard, Robin and Gonzalez‐Padilla, Daianna and Yalcinbas, Ege A and Chaloux‐Pinette, Emma and Eagles, Nicholas J and Totty, Michael S and Janak, Patricia H and Collado‐Torres, Leonardo and Maynard, Kristen R},
title = {{Transcriptional Response to Chronic Long-Access Fentanyl Self-Administration in Rat Habenula and Amygdala}},
journal = {Addiction biology},
year = {2026},
month = jul,
volume = {31},
number = {7},
pages = {e70179},
publisher = {Wiley},
issn = {1355-6215},
doi = {10.1111/adb.70179},
url = {https://doi.org/10.1111/adb.70179},
pmid = {42444546},
pmcid = {PMC13366401}
}

RIS

TY - JOUR
AU - Magnard, Robin
AU - Gonzalez‐Padilla, Daianna
AU - Yalcinbas, Ege A
AU - Chaloux‐Pinette, Emma
AU - Eagles, Nicholas J
AU - Totty, Michael S
AU - Janak, Patricia H
AU - Collado‐Torres, Leonardo
AU - Maynard, Kristen R
TI - Transcriptional Response to Chronic Long-Access Fentanyl Self-Administration in Rat Habenula and Amygdala
T2 - Addiction biology
J2 - Addict Biol
PY - 2026
DA - 2026/07/01
VL - 31
IS - 7
SP - e70179
SN - 1355-6215
PB - Wiley
DO - 10.1111/adb.70179
UR - https://doi.org/10.1111/adb.70179
LA - en
ER -

CSL-JSON

{
"id": "10.1111/adb.70179",
"type": "article-journal",
"title": "Transcriptional Response to Chronic Long-Access Fentanyl Self-Administration in Rat Habenula and Amygdala",
"container-title": "Addiction biology",
"author": [
{
"family": "Magnard",
"given": "Robin"
},
{
"family": "Gonzalez‐Padilla",
"given": "Daianna"
},
{
"family": "Yalcinbas",
"given": "Ege A"
},
{
"family": "Chaloux‐Pinette",
"given": "Emma"
},
{
"family": "Eagles",
"given": "Nicholas J"
},
{
"family": "Totty",
"given": "Michael S"
},
{
"family": "Janak",
"given": "Patricia H"
},
{
"family": "Collado‐Torres",
"given": "Leonardo"
},
{
"family": "Maynard",
"given": "Kristen R"
}
],
"container-title-short": "Addict Biol",
"volume": "31",
"issue": "7",
"page": "e70179",
"DOI": "10.1111/adb.70179",
"PMID": "42444546",
"PMCID": "PMC13366401",
"ISSN": "1355-6215",
"publisher": "Wiley",
"URL": "https://doi.org/10.1111/adb.70179",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
1
]
]
}
}

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

Similar papers

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

[1] doi:10.1016/j.xcrm.2026.102682 [code]
TET CpG sequence-context-specific DNA demethylation shapes progression of IDH-mutant gliomas.
Journal: Cell reports. Medicine
In common: Nextflow, edgeR, limma, 11 other tools, genetics / omics, other condition
[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: SingleCellExperiment, edgeR, limma, 11 other tools, genetics / omics, other condition, cellular / molecular
[3] doi:10.1002/imt2.70163 [code]
Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: SingleCellExperiment, edgeR, limma, 11 other tools, genetics / omics, cellular / molecular
[4] doi:10.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: edgeR, limma, rstatix, 10 other tools, genetics / omics, other condition, 1 reference
[5] doi:10.1038/s41467-026-71542-5 [code]
Astrocyte fatty acid metabolism as a driver of risk for major depressive disorder.
Journal: Nature communications
In common: SingleCellExperiment, edgeR, limma, 7 other tools, genetics / omics, cellular / molecular, 4 references
[6] 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: SingleCellExperiment, edgeR, limma, 10 other tools, cellular / molecular
[7] doi:10.1038/s41467-026-75723-0 [code]
Spatial transcriptomics reveals distinct cell type dynamics following opioid dependence in female mice with the common human μ-opioid receptor variant Oprm1 A118G.
Journal: Nature communications
In common: clusterProfiler, pheatmap, Seurat, 3 other tools, pain, genetics / omics, other condition, 1 other category, 6 references
[8] doi:10.1038/s41586-026-10214-2 [code]
Multidimensional profiling of heterogeneity in supratentorial ependymomas.
Journal: Nature
In common: SingleCellExperiment, edgeR, circlize, 9 other tools, genetics / omics, other condition
[9] doi:10.1038/s42003-026-10034-0 [code]
Region- and cell type-specific changes in gene expression in the cerebellum after classical fear conditioning.
Journal: Communications biology
In common: SingleCellExperiment, limma, rstatix, 9 other tools, genetics / omics, cellular / molecular
[10] doi:10.1038/s41467-026-73305-8 [code]
Comparative analysis of the cellular landscape in mammalian striatum.
Journal: Nature communications
In common: SingleCellExperiment, edgeR, clusterProfiler, 8 other tools, genetics / omics, cellular / molecular, 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.