OSCR

Prenatal acoustic communication triggers adaptive vascular programming in the developing avian brain.

Code ↔ Paper

14 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 14 matches · 1 of them tie a paragraph to a whole file, not to given lines: a weak match, whose lines are not tinted
  1. [1] § MATERIALS AND METHODS › Isoform switch analysis ↔ scripts/splicing_isoform_switch_pipeline.R, lines 131–173 · score 0.89 · analyzeAlternativeSplicing, single isoform genes, isoform switching, splice isoform, DEXSeq, FASTA
  2. [2] § MATERIALS AND METHODS › Weighted gene co-expression network analysis ↔ scripts/WGCNA_Analysis.R, lines 389–440 · score 0.82 · adjacency matrix, Kamada Kawai, kWithin, hub genes, green module, ggraph
  3. [3] § MATERIALS AND METHODS › Weighted gene co-expression network analysis ↔ scripts/WGCNA_Analysis.R, lines 85–124 · score 0.78 · soft thresholding, scale free topology, Networking construction, fit, connectivity, power
  4. [4] § MATERIALS AND METHODS › Weighted gene co-expression network analysis ↔ scripts/WGCNA_Analysis.R, lines 196–260 · score 0.75 · module trait relationships, module assignments, module eigengenes, Student, heatmap, enrichment
  5. [5] § MATERIALS AND METHODS › RNA extraction, sequencing and statistical analysis ↔ scripts/WGCNA_Analysis.R, lines 42–83 · score 0.73 · variance stabilization transformed, Quality control, DESeq2, VST, Sex, genes
  6. [6] § RESULTS › Heat call exposure alters isoform usage in contractile genes ↔ scripts/splicing_isoform_switch_pipeline.R, lines 131–173 · score 0.68 · single isoform gene, isoform switch, alternative splice, gene expression, filtering, Seq
  7. [7] § RESULTS › Vascular gene expression changes localize to cerebrovascular cell types ↔ scripts/WGCNA_Analysis.R, lines 442–528 · score 0.67 · HYPOMAP cell, strong enrichment, FDR corrected, module genes, hypergeometric, Heatmap
  8. [8] § MATERIALS AND METHODS › RNA extraction, sequencing and statistical analysis ↔ scripts/WGCNA_Analysis.R, lines 442–528 · score 0.66 · background gene, gene symbols, hypergeometric, overlapping, FDR, enrichment
  9. [9] § RESULTS › Co-expression network analysis reveals coordinated vascular gene regulation ↔ scripts/WGCNA_Analysis.R, lines 287–331 · score 0.64 · intramodular connectivity, kWithin, hub genes, module genes, adjacency, WGCNA
  10. [10] § RESULTS › Co-expression network analysis reveals coordinated vascular gene regulation ↔ scripts/WGCNA_Analysis.R, lines 287–331 · score 0.62 · Intramodular connectivity, kWithin, hub genes, GO, adjacency, modules
  11. [11] § MATERIALS AND METHODS › RNA extraction, sequencing and statistical analysis ↔ scripts/global_rnaseq_deseq2_pipeline.R, lines 305–380 · score 0.61 · Principal component, RNA Seq, PCA, females, VST, variance
  12. [12] § MATERIALS AND METHODS › Deconvolution analysis ↔ scripts/hypomap_run_cibersortx_hires.sh, the whole file · a weak match · score 0.59 · custom signature matrix, CIBERSORTx, mixture, HYPOMAP
  13. [13] § MATERIALS AND METHODS › RNA extraction, sequencing and statistical analysis ↔ scripts/global_rnaseq_deseq2_pipeline.R, lines 159–223 · score 0.58 · likelihood ratio, DESeq2, shrinkage, apeglm, model, sex
  14. [14] § RESULTS › Co-expression network analysis reveals coordinated vascular gene regulation ↔ scripts/WGCNA_Analysis.R, lines 196–260 · score 0.54 · Module trait relationships, module eigengenes, Student, Heatmap, WGCNA, network

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 · 528 lines · 17 KB · MIT · 9 matches

  1. # Setup and Data Preparation
  2. # Set working directory (change as appropriate)
  3. setwd("your/working/directory")
  5. # Load required packages
  6. library(WGCNA)
  7. library(rtracklayer)
  8. library(dplyr)
  9. library(ggplot2)
  10. library(clusterProfiler)
  11. library(pheatmap)
  12. library(DESeq2) # For VST normalization
  14. # Avoid factor conversion issues
  15. options(stringsAsFactors = FALSE)
  17. # Increase memory limit for large datasets
  18. options(future.globals.maxSize = 250000 * 1024^2)
  20. # Step 1: Prepare gene length data from GFF annotation
  21. my_tags <- c("Name", "Dbxref", "gene")
  22. my_columns <- c("seqid", "start", "end", "strand", "type")
  23. my_filter <- list(type = "gene")
  25. dat <- readGFF("GCF_003957565.2_bTaeGut1.4.pri_genomic.gff",
  26. tags = my_tags,
  27. columns = my_columns,
  28. filter = my_filter)
  30. dat$gene_name <- dat$gene
  31. dat$interval_start <- dat$start
  32. dat$interval_stop <- dat$end
  34. dat_1 <- as.data.frame(dat)
  35. all_data_from_dat <- dplyr::select(dat_1, gene_name, start, end)
  37. # Calculate gene length as end - start - 1
  38. all_data_from_dat$geneLength <- (all_data_from_dat$end - all_data_from_dat$start) - 1
  40. gene_name_length <- dplyr::select(all_data_from_dat, gene_name, geneLength)
  42. # Step 2: Read and process raw gene counts data
  43. data0_csv <- read.csv("WGCNA_gene_count_7_June_2023.csv", header = TRUE)
  45. # Rename first column for merging
  46. colnames(data0_csv)[1] <- "gene_name"
  48. # Merge gene length data with counts (geneLength retained for reference)
  49. data0 <- merge(data0_csv, gene_name_length, by = "gene_name", all.x = TRUE)
  51. # Set rownames as gene names and remove gene_name column
  52. colnames(data0)[1] <- ""
  53. rownames(data0) <- data0[, 1]
  54. data0 <- data0[, -1]
  56. # Convert numeric columns to integer counts (raw counts)
  57. data0[] <- lapply(data0, function(x) if(is.numeric(x)) as.integer(x) else x)
  59. # Step 3: Load sample metadata
  60. sample_metadata <- read.csv(file = "WGCNA_Heat_Call_Experimental_Design.csv")
  61. colnames(sample_metadata) <- c('sample_ID', 'Sex', 'Condition')
  62. rownames(sample_metadata) <- sample_metadata$sample_ID
  64. # Normalization and Quality Control
  65. # Step 4: Normalize raw counts using DESeq2 VST (Variance Stabilizing Transformation)
  67. # Remove geneLength column before normalization
  68. counts_for_norm <- data0[, !colnames(data0) %in% "geneLength"]
  70. # Ensure sample names match exactly between counts and metadata
  71. stopifnot(setequal(colnames(counts_for_norm), rownames(sample_metadata)))
  73. # Create DESeq2 dataset (design ~1 for normalization only)
  74. dds <- DESeqDataSetFromMatrix(countData = counts_for_norm,
  75. colData = sample_metadata,
  76. design = ~ 1)
  78. # Apply variance stabilizing transformation (blind = TRUE)
  79. vsd <- vst(dds, blind = TRUE)
  81. # Extract normalized expression matrix and transpose for WGCNA
  82. datExpr <- assay(vsd)
  83. datExpr <- t(datExpr)
  85. # Step 5: Sample Quality Control with Outlier Detection
  86. sampleTree <- hclust(dist(datExpr), method = "average")
  88. pdf("sample_clustering_outliers.pdf", width = 12, height = 6)
  89. plot(sampleTree, main = "Sample Clustering for Outlier Detection")
  90. dev.off()
  92. # Optional: Remove outliers identified visually
  93. # outlier_samples <- c("Sample1", "Sample2")
  94. # datExpr <- datExpr[!rownames(datExpr) %in% outlier_samples, ]
  96. # Network Construction and Module Detection
  97. # Step 6: Soft Threshold Diagnostics for Network Construction
  98. powers <- c(1:20, seq(22, 30, 2))
  100. # Disable parallel processing to avoid forking issues
  101. disableWGCNAThreads()
  103. sft <- pickSoftThreshold(
  104. datExpr,
  105. powerVector = powers,
  106. verbose = 5,
  107. networkType = "signed",
  108. blockSize = ncol(datExpr)
  109. )
  111. pdf("soft_threshold_diagnostics.pdf", width = 10, height = 5)
  112. par(mfrow = c(1, 2))
  113. plot(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2],
  114. xlab = "Soft Threshold (power)", ylab = "Scale Free Topology Fit (R²)",
  115. main = "Scale Independence", type = "n")
  116. text(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2],
  117. labels = powers, col = "red")
  118. abline(h = 0.85, col = "red", lty = 2)
  120. plot(sft$fitIndices[, 1], sft$fitIndices[, 5],
  121. xlab = "Soft Threshold (power)", ylab = "Mean Connectivity",
  122. main = "Mean Connectivity", type = "n")
  123. text(sft$fitIndices[, 1], sft$fitIndices[, 5], labels = powers, col = "red")
  124. dev.off()
  126. # Select optimal power
  127. if (is.na(sft$powerEstimate)) {
  128. optimal_power <- sft$fitIndices$Power[which.max(sft$fitIndices$SFT.R.sq)]
  129. message("No power reached R² > 0.85. Using max R² at power: ", optimal_power)
  130. } else {
  131. optimal_power <- sft$powerEstimate
  132. }
  134. save.image("06_July_2025_Pre_Parameter.Rdata")
  136. # Step 7: Parameter Sensitivity Analysis for Module Detection
  137. parameter_grid <- expand.grid(
  138. minModuleSize = c(20, 30, 40),
  139. mergeCutHeight = c(0.15, 0.25, 0.35)
  140. )
  142. module_results <- list()
  144. for (i in 1:nrow(parameter_grid)) {
  145. params <- parameter_grid[i, ]
  147. assign("cor", WGCNA::cor, envir = .GlobalEnv)
  149. net <- blockwiseModules(
  150. datExpr,
  151. power = optimal_power,
  152. TOMType = "signed",
  153. minModuleSize = params$minModuleSize,
  154. mergeCutHeight = params$mergeCutHeight,
  155. numericLabels = TRUE,
  156. saveTOMs = FALSE,
  157. verbose = 3
  158. )
  160. assign("cor", stats::cor, envir = .GlobalEnv)
  162. module_results[[i]] <- list(
  163. params = params,
  164. n_modules = length(unique(net$colors)),
  165. module_sizes = table(net$colors)
  166. )
  168. mergedColors <- labels2colors(net$colors)
  169. pdf(paste0("module_dendrogram_minSize", params$minModuleSize,
  170. "_mergeCut", params$mergeCutHeight, ".pdf"),
  171. width = 10, height = 6)
  172. plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
  173. "Module colors", dendroLabels = FALSE)
  174. dev.off()
  175. }
  177. parameter_summary <- do.call(rbind, lapply(module_results, function(x) {
  178. data.frame(
  179. minModuleSize = x$params$minModuleSize,
  180. mergeCutHeight = x$params$mergeCutHeight,
  181. n_modules = x$n_modules,
  182. min_module_size = min(x$module_sizes),
  183. max_module_size = max(x$module_sizes)
  184. )
  185. }))
  187. print(parameter_summary)
  189. final_params <- list(
  190. minModuleSize = 30,
  191. mergeCutHeight = 0.25
  192. )
  194. save.image("06_July_2025_pre_network_WGCNA.RData")
  196. # Step 8: Final Network Construction
  197. enableWGCNAThreads(nThreads = 20)
  199. assign("cor", WGCNA::cor, envir = .GlobalEnv)
  201. net <- blockwiseModules(
  202. datExpr,
  203. power = optimal_power,
  204. TOMType = "signed",
  205. minModuleSize = final_params$minModuleSize,
  206. mergeCutHeight = final_params$mergeCutHeight,
  207. numericLabels = TRUE,
  208. saveTOMs = TRUE,
  209. saveTOMFileBase = "network_TOM",
  210. verbose = 3,
  211. maxBlockSize = 10000
  212. )
  214. assign("cor", stats::cor, envir = .GlobalEnv)
  216. moduleColors <- labels2colors(net$colors)
  217. write.csv(data.frame(Gene = colnames(datExpr), Module = moduleColors),
  218. "final_module_assignments.csv", row.names = FALSE)
  220. pdf("final_module_dendrogram.pdf", width = 15, height = 10)
  221. plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]],
  222. "Module colors", dendroLabels = FALSE)
  223. dev.off()
  225. save.image("06_July_2025_post_network_WGCNA.RData")
  227. #Module-Trait Relationships
  228. # Step 9: Module-Trait Relationships
  229. traits <- sample_metadata[rownames(datExpr), c("Sex", "Condition")]
  230. traits$Sex <- as.numeric(factor(traits$Sex))
  231. traits$Condition <- as.numeric(factor(traits$Condition))
  233. MEs <- moduleEigengenes(datExpr, moduleColors)$eigengenes
  235. moduleTraitCor <- cor(MEs, traits, use = "p")
  236. moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nrow(datExpr))
  238. pdf("module_trait_relationships.pdf", width = 10, height = 8)
  239. textMatrix <- paste(signif(moduleTraitCor, 2), "\n(",
  240. signif(moduleTraitPvalue, 1), ")", sep = "")
  241. dim(textMatrix) <- dim(moduleTraitCor)
  243. par(mar = c(6, 8.5, 3, 3))
  244. labeledHeatmap(Matrix = moduleTraitCor,
  245. xLabels = names(traits),
  246. yLabels = names(MEs),
  247. ySymbols = names(MEs),
  248. colorLabels = FALSE,
  249. colors = blueWhiteRed(50),
  250. textMatrix = textMatrix,
  251. setStdMargins = FALSE,
  252. cex.text = 0.5,
  253. zlim = c(-1, 1),
  254. main = "Module-Trait Relationships")
  255. dev.off()
  257. # Gene Ontology Enrichment
  258. library(org.Gg.eg.db, lib.loc = "/project/cugbf/software/R/4.4.1/library")
  260. module_genes <- split(colnames(datExpr), moduleColors)
  262. perform_GO_analysis <- function(module_name, genes) {
  263. if (length(genes) < 10) {
  264. message("Skipping GO for ", module_name, " (only ", length(genes), " genes)")
  265. return(NULL)
  266. }
  268. genes_entrez <- bitr(genes, fromType = "SYMBOL", toType = "ENTREZID",
  269. OrgDb = org.Gg.eg.db, drop = TRUE)
  271. if (nrow(genes_entrez) < 5) {
  272. message("Skipping GO for ", module_name, " (only ", nrow(genes_entrez), " ENTREZ IDs)")
  273. return(NULL)
  274. }
  276. enrichGO(
  277. gene = genes_entrez$ENTREZID,
  278. OrgDb = org.Gg.eg.db,
  279. keyType = "ENTREZID",
  280. ont = "BP",
  281. pAdjustMethod = "BH",
  282. pvalueCutoff = 0.05,
  283. qvalueCutoff = 0.2
  284. )
  285. }
  287. go_results <- lapply(names(module_genes), function(module) {
  288. result <- tryCatch({
  289. perform_GO_analysis(module, module_genes[[module]])
  290. }, error = function(e) {
  291. message("GO failed for ", module, ": ", conditionMessage(e))
  292. return(NULL)
  293. })
  295. if (!is.null(result) && nrow(result) > 0) {
  296. pdf(paste0("GO_enrichment_", module, ".pdf"), width = 10, height = 8)
  297. print(barplot(result, showCategory = 20, title = paste("GO Enrichment -", module)))
  298. dev.off()
  299. }
  301. return(result)
  302. })
  303. names(go_results) <- names(module_genes)
  305. saveRDS(go_results, "06_July_2025_go_enrichment_results.rds")
  307. # Hub Gene Identification and Validation
  308. ADJ <- adjacency(datExpr, power = optimal_power, type = "signed")
  309. kIN <- intramodularConnectivity(ADJ, moduleColors, scaleByMax = TRUE)
  311. hub_genes <- lapply(unique(moduleColors), function(mod) {
  312. modGenes <- (moduleColors == mod)
  313. kIN_mod <- kIN[modGenes, ]
  314. head(kIN_mod[order(kIN_mod$kWithin, decreasing = TRUE), ], 10)
  315. })
  316. names(hub_genes) <- unique(moduleColors)
  318. saveRDS(hub_genes, "hub_genes_per_module.rds")
  320. pdf("hub_gene_validation.pdf", width = 12, height = 8)
  321. for (mod in names(hub_genes)) {
  322. top_gene <- rownames(hub_genes[[mod]])[1]
  323. plot(datExpr[, top_gene],
  324. MEs[, paste0("ME", mod)],
  325. main = paste("Hub Gene Validation -", mod),
  326. xlab = top_gene,
  327. ylab = "Module Eigengene")
  328. }
  329. dev.off()
  331. save.image("07_july_2025_WGCNA_analysis_final.RData")
  333. # Network Visualization
  334. library(igraph)
  335. library(ggraph)
  336. library(tidygraph)
  337. library(dplyr)
  338. library(ggplot2)
  340. module_of_interest <- "green"
  342. module_genes <- colnames(datExpr)[moduleColors == module_of_interest]
  344. ADJ <- adjacency(datExpr, power = optimal_power, type = "signed")
  346. module_adj <- ADJ[module_genes, module_genes]
  348. edge_threshold <- 0.1
  349. module_adj[module_adj < edge_threshold] <- 0
  351. graph <- graph_from_adjacency_matrix(module_adj, mode = "undirected", weighted = TRUE, diag = FALSE)
  353. kIN_module <- kIN[moduleColors == module_of_interest, ]
  355. top_hubs <- rownames(kIN_module)[order(kIN_module$kWithin, decreasing = TRUE)[1:10]]
  357. V(graph)$hub <- V(graph)$name %in% top_hubs
  359. V(graph)$color <- ifelse(V(graph)$hub, "red", "green")
  361. V(graph)$size <- ifelse(V(graph)$hub, 8, 4)
  363. E(graph)$width <- E(graph)$weight * 5
  365. p <- ggraph(graph, layout = "fr") +
  366. geom_edge_link(aes(width = weight), alpha = 0.3, color = "grey50") +
  367. geom_node_point(aes(color = color, size = size)) +
  368. geom_node_text(aes(label = ifelse(hub, name, "")), repel = TRUE, size = 4, color = "black", fontface = "bold") +
  369. scale_color_identity() +
  370. scale_size_identity() +
  371. theme_void() +
  372. ggtitle("WGCNA Green Module Network with Hub Genes Highlighted") +
  373. theme(plot.title = element_text(hjust = 0.5, size = 18, face = "bold"))
  375. print(p)
  377. library(WGCNA)
  378. exportNetworkToCytoscape(
  379. adjMat = module_adj,
  380. edgeFile = "CytoscapeInput-edges-green.txt",
  381. nodeFile = "CytoscapeInput-nodes-green.txt",
  382. weighted = TRUE,
  383. threshold = 0.2,
  384. nodeNames = module_genes,
  385. nodeAttr = moduleColors[moduleColors == module_of_interest],
  386. includeColNames = TRUE
  387. )
  389. # Less busy plot with top 5 hubs and their neighbors
  390. top_hubs <- rownames(kIN_module)[order(kIN_module$kWithin, decreasing = TRUE)[1:5]]
  392. edge_threshold <- 0.4
  394. neighbors <- unique(unlist(lapply(top_hubs, function(hub) {
  395. hub_edges <- module_adj[hub, ]
  396. top_neighbors <- names(sort(hub_edges, decreasing = TRUE)[1:5])
  397. top_neighbors[hub_edges[top_neighbors] > edge_threshold]
  398. })))
  399. genes_to_plot <- unique(c(top_hubs, neighbors))
  401. sub_adj <- module_adj[genes_to_plot, genes_to_plot]
  402. sub_graph <- graph_from_adjacency_matrix(sub_adj, mode = "undirected", weighted = TRUE, diag = FALSE)
  404. V(sub_graph)$hub <- V(sub_graph)$name %in% top_hubs
  405. V(sub_graph)$color <- ifelse(V(sub_graph)$hub, "red", "green")
  406. V(sub_graph)$size <- ifelse(V(sub_graph)$hub, 8, 4)
  408. p_sub <- ggraph(sub_graph, layout = "fr", niter = 4000, area = 40000, repulserad = 20000) +
  409. geom_edge_link(aes(width = weight), alpha = 0.3, color = "grey50") +
  410. geom_node_point(aes(color = color, size = size)) +
  411. geom_node_text(aes(label = ifelse(hub, name, "")), repel = TRUE, size = 4, color = "black", fontface = "bold", max.overlaps = Inf) +
  412. scale_color_identity() +
  413. scale_size_identity() +
  414. theme_void() +
  415. ggtitle("Green Module: Hub Genes and Their Top Neighbors") +
  416. theme(plot.title = element_text(hjust = 0.5, size = 18, face = "bold"))
  418. print(p_sub)
  420. hub_color <- "#0072B2"
  421. neighbor_color <- "#D55E00"
  423. V(sub_graph)$color <- ifelse(V(sub_graph)$hub, hub_color, neighbor_color)
  424. V(sub_graph)$size <- ifelse(V(sub_graph)$hub, 10, 5)
  425. V(sub_graph)$frame.color <- "black"
  427. p_kk <- ggraph(sub_graph, layout = "kk") +
  428. geom_edge_link(aes(width = weight), alpha = 0.25, color = "grey70") +
  429. geom_node_point(aes(color = color, size = size), shape = 21, stroke = 1.2) +
  430. geom_node_text(aes(label = ifelse(hub, name, "")), repel = TRUE, size = 4, color = "black", fontface = "bold", bg.color = "white", bg.r = 0.15, max.overlaps = Inf) +
  431. scale_color_identity() +
  432. scale_size_identity() +
  433. theme_void() +
  434. ggtitle("Green Module: Hub Genes and Their Top Neighbors (Kamada-Kawai)") +
  435. theme(plot.title = element_text(hjust = 0.5, size = 18, face = "bold"), legend.position = "none")
  437. ggsave("green_module_hub_subnetwork_nature.pdf", p_kk, width = 12, height = 12, dpi = 600)
  438. ggsave("green_module_hub_subnetwork_nature.png", p_kk, width = 12, height = 12, dpi = 600)
  440. print(p_kk)
  442. # Cell-Type Enrichment Analysis
  443. library(dplyr)
  444. library(clusterProfiler)
  445. library(org.Hs.eg.db)
  446. library(pheatmap)
  447. library(RColorBrewer)
  449. # Load cell-type marker genes (converted to gene symbols)
  450. top_markers_converted <- read.csv("hypomap_cell_markers_gene_id.csv")
  452. # Create a named list of marker genes per cell type
  453. cell_markers_list <- split(top_markers_converted$gene, top_markers_converted$cluster)
  455. # Load module assignments
  456. module_gene_df <- read.csv("final_module_assignments.csv") # columns: Gene, Module
  458. # Create a named list of genes per module
  459. module_genes <- split(module_gene_df$Gene, module_gene_df$Module)
  461. # Define the background gene universe (all genes assigned to modules)
  462. all_genes <- unique(module_gene_df$Gene)
  464. # Perform hypergeometric enrichment test with FDR correction
  465. enrichment_pvals <- matrix(NA, nrow = length(module_genes), ncol = length(cell_markers_list))
  466. rownames(enrichment_pvals) <- names(module_genes)
  467. colnames(enrichment_pvals) <- names(cell_markers_list)
  469. for (mod in names(module_genes)) {
  470. for (ct in names(cell_markers_list)) {
  471. overlap <- intersect(module_genes[[mod]], cell_markers_list[[ct]])
  472. enrichment_pvals[mod, ct] <- phyper(
  473. length(overlap) - 1,
  474. length(cell_markers_list[[ct]]),
  475. length(all_genes) - length(cell_markers_list[[ct]]),
  476. length(module_genes[[mod]]),
  477. lower.tail = FALSE
  478. )
  479. }
  480. }
  482. # Multiple testing correction (FDR)
  483. enrichment_fdr <- matrix(
  484. p.adjust(as.vector(enrichment_pvals), method = "BH"),
  485. nrow = nrow(enrichment_pvals),
  486. dimnames = dimnames(enrichment_pvals)
  487. )
  489. # Annotate FDR values with asterisks for significance
  490. annotations <- matrix("", nrow = nrow(enrichment_fdr), ncol = ncol(enrichment_fdr))
  491. annotations[enrichment_fdr < 0.001] <- "***"
  492. annotations[enrichment_fdr < 0.01 & enrichment_fdr >= 0.001] <- "**"
  493. annotations[enrichment_fdr < 0.05 & enrichment_fdr >= 0.01] <- "*"
  495. display_numbers <- matrix("", nrow = nrow(enrichment_fdr), ncol = ncol(enrichment_fdr))
  496. for (i in seq_len(nrow(enrichment_fdr))) {
  497. for (j in seq_len(ncol(enrichment_fdr))) {
  498. if (enrichment_fdr[i, j] < 0.001) {
  499. display_numbers[i, j] <- paste0("<0.001", annotations[i, j])
  500. } else {
  501. display_numbers[i, j] <- paste0(sprintf("%.3f", enrichment_fdr[i, j]), annotations[i, j])
  502. }
  503. }
  504. }
  506. # Define color palette: low FDR (strong enrichment) = dark red, high FDR = white
  507. color_palette <- colorRampPalette(c("darkred", "white"))(100)
  509. # Plot publication-quality heatmap
  510. pheatmap(
  511. enrichment_fdr,
  512. color = color_palette,
  513. cluster_rows = TRUE,
  514. cluster_cols = TRUE,
  515. fontsize_row = 10,
  516. fontsize_col = 10,
  517. angle_col = 45,
  518. treeheight_row = 0,
  519. treeheight_col = 0,
  520. legend = TRUE,
  521. display_numbers = display_numbers,
  522. number_color = "black",
  523. main = "Module Enrichment in Cell Types (FDR)",
  524. legend_title = "FDR values",
  525. width = 12,
  526. height = 10,
  527. filename = "module_cell_type_enrichment_FDR_reversed.pdf"
  528. )

WGCNA_Analysis.R at commit e7f779e, under MIT · at the source

Overview

Authors: Prakrit Subba1, Mylene M. Mariette2,3, Katerina A. Palios4, Michael G. Emmerson4, Elisabetta Versace4, Katherine L. Buchanan3, David F. Clayton5, Julia M. George1
ORCID iDs: Prakrit Subba, Mylene M. Mariette, Katerina A. Palios, Elisabetta Versace, Katherine L. Buchanan, David F. Clayton, Julia M. George
  1. Department of Biological Sciences, Clemson University, Clemson, SC 29634, USA
  2. Doñana Biological Station, EBD-CSIC, Seville 41092, Spain
  3. School of Life and Environmental Sciences, Deakin University, Waurn Ponds, VIC 3288, Australia
  4. School of Biological and Behavioral Sciences, Queen Mary University of London, London E1 4NS, UK
  5. Department of Genetics and Biochemistry, Clemson University, Clemson, SC 29634, USA
Institutions: Clemson University (United States); Deakin University (Australia); Estación Biológica de Doñana (Spain); Queen Mary University of London (United Kingdom)
Journal: The Journal of experimental biology, volume 229, issue 11, article jeb252287
Dates: received 1 February 2026; accepted 24 April 2026; published online 15 June 2026; in print June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1242/jeb.252287 · PMID 42276015 · PMCID PMC13327539 · OpenAlex W7164324788
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: other (organism), developmental (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Connectivity, Spectral & time-frequency
Keywords: Thermal physiology, Acoustic communication, Developmental plasticity, Developmental programming, Blood–brain barrier, Zebra finch
MeSH: Brain*, Finches*, Hypothalamus*, Vocalization, Animal*, Animals, Embryo, Nonmammalian, Female, Gene Expression Regulation, Developmental (* major topic)
Topic: Animal Vocal Communication and Behavior (Developmental Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: Australian Research Council (DP180101207, DE170100824, FT140100131, DP210101238); RCUK | Biotechnology and Biological Sciences Research Council (BBSRC) (BB/S003223/1); Ministerio de Ciencia e Innovación (RYC2019-028066-I, PID2021-128494NA-I00); Clemson University
Citations: not cited yet (Europe PMC); 107 references in the paper

Abstract

Developmental plasticity allows organisms to adjust their phenotypes to match environmental conditions, but how sensory cues program specific physiological systems remains poorly understood. In Australian zebra finches, incubating parents emit heat calls during extreme temperatures, and embryos exposed to these acoustic signals develop enhanced thermal tolerance and altered growth trajectories as adults, a striking example of anticipatory programming. We hypothesized that heat call exposure alters embryonic hypothalamic gene expression, given this brain region's central role in integrating environmental signals and regulating metabolism, thermoregulation and growth. We exposed zebra finch embryos to playback of parental heat calls or control calls during late incubation and used RNA-sequencing of hypothalamic tissue to identify transcriptional responses. Contrary to predictions of widespread neuroendocrine reprogramming, heat call exposure produced targeted changes: robust downregulation of genes regulating vascular smooth muscle contraction and cytoskeletal dynamics, with coordinated isoform switching. Cell-type analyses revealed these molecular changes localized to vascular endothelial cells, smooth muscle cells and ependymal cells, the cellular components that control cerebral blood flow and regulate the brain's vascular barrier. Gene expression patterns suggest increased vascular plasticity that may protect against heat-induced cellular damage. Remarkably, these adaptive modifications occurred in response to an acoustic signal alone, without thermal exposure. Our results provide transcriptional evidence that prenatal acoustic cues may program cerebrovascular function through cell type-specific gene regulation, providing a novel mechanism for sensory-mediated developmental plasticity. This targeted vascular programming may represent a conserved strategy for anticipatory adaptation to predictable thermal challenges across endothermic vertebrates.

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

praxsubba/Heat_Call_Reprogramming_Bulk_RNA_Seq

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: e7f779eab452a321a8016f665965747842a64171, 18 February 2026
Languages: R (4), Shell (1)
Size: 16 files, 5 scripts
Software Heritage: not archived
Found in: the end of the paper
Holds: README, license file, documentation
Not found: CITATION.cff, environment file, tests, continuous integration
Tools: DESeq2 (3 files), ggplot2 (3 files), tidyverse (3 files), clusterProfiler (1 file), igraph (1 file), pheatmap (1 file), WGCNA (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
7 files
At the source: github.com/praxsubba/Heat_Call_Reprogramming_Bulk_RNA_Seq

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;
  • 5 scripts, each with its path and the digest of its content;
  • 14 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

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

  • Publisher: n/a → The Company of Biologists

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 8 authors, 6 keywords, 8 MeSH terms, 4 funders, 106 references.

Cite

This paper

Subba, P., Mariette, M. M., Palios, K. A., Emmerson, M. G., Versace, E., Buchanan, K. L., Clayton, D. F., & George, J. M. (2026). Prenatal acoustic communication triggers adaptive vascular programming in the developing avian brain. The Journal of experimental biology, 229(11), jeb252287. https://doi.org/10.1242/jeb.252287

BibTeX

@article{subba2026prenatal,
author = {Subba, Prakrit and Mariette, Mylene M. and Palios, Katerina A. and Emmerson, Michael G. and Versace, Elisabetta and Buchanan, Katherine L. and Clayton, David F. and George, Julia M.},
title = {{Prenatal acoustic communication triggers adaptive vascular programming in the developing avian brain}},
journal = {The Journal of experimental biology},
year = {2026},
month = jun,
volume = {229},
number = {11},
pages = {jeb252287},
publisher = {The Company of Biologists},
issn = {0022-0949},
doi = {10.1242/jeb.252287},
url = {https://doi.org/10.1242/jeb.252287},
pmid = {42276015},
pmcid = {PMC13327539}
}

RIS

TY - JOUR
AU - Subba, Prakrit
AU - Mariette, Mylene M.
AU - Palios, Katerina A.
AU - Emmerson, Michael G.
AU - Versace, Elisabetta
AU - Buchanan, Katherine L.
AU - Clayton, David F.
AU - George, Julia M.
TI - Prenatal acoustic communication triggers adaptive vascular programming in the developing avian brain
T2 - The Journal of experimental biology
J2 - J Exp Biol
PY - 2026
DA - 2026/06/11
VL - 229
IS - 11
SP - jeb252287
SN - 0022-0949
PB - The Company of Biologists
DO - 10.1242/jeb.252287
UR - https://doi.org/10.1242/jeb.252287
LA - en
ER -

CSL-JSON

{
"id": "10.1242/jeb.252287",
"type": "article-journal",
"title": "Prenatal acoustic communication triggers adaptive vascular programming in the developing avian brain",
"container-title": "The Journal of experimental biology",
"author": [
{
"family": "Subba",
"given": "Prakrit"
},
{
"family": "Mariette",
"given": "Mylene M."
},
{
"family": "Palios",
"given": "Katerina A."
},
{
"family": "Emmerson",
"given": "Michael G."
},
{
"family": "Versace",
"given": "Elisabetta"
},
{
"family": "Buchanan",
"given": "Katherine L."
},
{
"family": "Clayton",
"given": "David F."
},
{
"family": "George",
"given": "Julia M."
}
],
"container-title-short": "J Exp Biol",
"volume": "229",
"issue": "11",
"page": "jeb252287",
"DOI": "10.1242/jeb.252287",
"PMID": "42276015",
"PMCID": "PMC13327539",
"ISSN": "0022-0949",
"publisher": "The Company of Biologists",
"URL": "https://doi.org/10.1242/jeb.252287",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
11
]
]
}
}

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

Similar papers

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

[1] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: WGCNA, igraph, DESeq2, 4 other tools, 3 references
[2] doi:10.1038/s41380-026-03629-w [code]
Maternal fasting during early gestation induces epigenetic alterations and schizophrenia-related phenotypes.
Journal: Molecular psychiatry
In common: WGCNA, igraph, DESeq2, 4 other tools, 3 references
[3] doi:10.1038/s41467-026-72598-z [code]
Functional impact of genetic background on variable expressivity in neurodevelopmental disorders.
Journal: Nature communications
In common: WGCNA, DESeq2, clusterProfiler, 3 other tools, 4 references
[4] doi:10.1038/s44318-026-00806-z [code]
Interspecific diversity in the neuronal composition of the mammalian cortex arises from heterochrony in neurogenesis.
Journal: The EMBO journal
In common: WGCNA, igraph, DESeq2, 4 other tools, developmental, 1 reference
[5] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: WGCNA, igraph, DESeq2, 4 other tools, 1 reference
[6] doi:10.1038/s41467-026-76675-1 [code]
Long-read proteogenomic atlas of human neuronal differentiation reveals isoform diversity informing neurodevelopmental risk mechanisms.
Journal: Nature communications
In common: WGCNA, igraph, DESeq2, 3 other tools, developmental, 2 references
[7] doi:10.1038/s41593-026-02367-0 [code]
A reproducible three-dimensional model of human brain tissue to investigate physiological and disease-associated microglia phenotypes.
Journal: Nature neuroscience
In common: WGCNA, igraph, DESeq2, 4 other tools, 1 reference
[8] 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: WGCNA, igraph, DESeq2, 4 other tools, 1 reference
[9] doi:10.1016/j.isci.2026.115196 [code]
Transcriptional and cellular maturation of the chick spinal cord in the context of distinct neuromuscular circuits.
Journal: iScience
In common: WGCNA, igraph, DESeq2, 3 other tools, other, 2 references
[10] doi:10.1038/s41586-026-10629-x [code]
Whole-genome duplication shaped cell-type evolution in the vertebrate brain.
Journal: Nature
In common: igraph, DESeq2, clusterProfiler, 3 other tools, other, 2 references

Contribute

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

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

Request its removal

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

Discussion, reproductions, activity

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

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

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