OSCR

Brain Transcriptomic Changes Associated With the Socially Regulated Behavioural Specialisation of Ant Queens.

Code ↔ Paper

9 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 9 matches
  1. [1] § Results › Effect of the Presence of Larvae on the Expression of the Genes of Interest (Experiment C) ↔ R script_Mol_Ecol_2026.R, lines 577 to 650 · score 1.00 · Pheromone binding protein, Lnig_g15862, gLnigGR81, Juvenile hormone esterase, Lnig_g01559, Lnig_g04508 · Report a mismatch
  2. [2] § Results › Effect of Separated Workers on the Expression of the Genes of Interest (Experiment B) ↔ R script_Mol_Ecol_2026.R, lines 577 to 650 · score 0.98 · Juvenile hormone esterase, Lnig_g01559, Lnig_g04508, Lnig_g11740, Lnig_g11948, Lnig_g13064 · Report a mismatch
  3. [3] § Methods › Gene Expression Analyses › Analysis of the Effect of Worker Presence on Queen Brain Gene Expression ↔ R script_Mol_Ecol_2026.R, lines 378 to 416 · score 0.83 · topGO, GO annotation, biological process, GO terms, Fisher, Ontology · Report a mismatch
  4. [4] § Methods › Gene Expression Analyses › Analysis of the Effect of Worker Presence on Queen Brain Gene Expression ↔ R script_Mol_Ecol_2026.R, lines 337 to 375 · score 0.81 · explanatory variables, full model, reduced model, DESeq2, LRT, batch · Report a mismatch
  5. [5] § Results › Effect of Worker Removal on the Expression of the Genes of Interest (Experiment A) ↔ R script_Mol_Ecol_2026.R, lines 786 to 827 · score 0.75 · Juvenile hormone esterase, Lnig_g17128, Lnig_g01448, Lnig_g08195, Kr h1, Insulin · Report a mismatch
  6. [6] § Methods › Sample Collection › Effect of the Experimental Addition of Workers and Separated Workers ↔ R script_Mol_Ecol_2026.R, lines 1407 to 1477 · score 0.68 · wire mesh, antennal contacts, separated workers, queens passed, touching · Report a mismatch
  7. [7] § Results › Effect of the Presence of Workers on Brain Gene Expression in Queens ↔ R script_Mol_Ecol_2026.R, lines 418 to 456 · score 0.55 · biological processes, molecular functions, GO term, genes · Report a mismatch
  8. [8] § Results › Effect of Separated Workers on the Expression of the Genes of Interest (Experiment B) ↔ R script_Mol_Ecol_2026.R, lines 1407 to 1477 · score 0.54 · wire mesh, antennal contacts, separated workers, behaviour, queen · Report a mismatch
  9. [9] § Results › Effect of the Presence of Workers on Brain Gene Expression in Queens ↔ R script_Mol_Ecol_2026.R, lines 418 to 456 · score 0.52 · biological processes, molecular functions, GO term, DEG, gene · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 1,888 lines · 59 KB · no license · 9 matches

  1. #R code for the statistical analyses of :
  3. #Transcriptomic changes associated with the socially regulated behavioral specialization of ant queens
  5. #2026
  7. #Hugo Le Lay, Maximilian F. Bolder, Thomas J. Colgan, Romain Libbrecht
  9. rm(list=ls())
  11. #Load packages-----
  13. library(viridis)
  14. library(readxl) # Import Excel files (.xls, .xlsx)
  15. library(car) # Tools for regression analysis and ANOVA
  16. library(lme4) # Linear and generalized linear mixed-effects models (LMM, GLMM)
  17. library(emmeans) # Estimated marginal means and post-hoc comparisons
  18. library(DHARMa) # Model diagnostics using simulated residuals
  19. library(tximport) # Import and summarize transcript-level RNA-seq quantification data
  20. library(DESeq2) # Differential gene expression analysis for RNA-seq data
  21. library(topGO) # Gene Ontology (GO) enrichment analysis
  22. library(tidyverse) # Collection of packages for data manipulation and visualization
  23. library(ggrepel) # Add non-overlapping text labels to ggplot2 plots
  24. library(ggpubr) # Easy-to-use functions for publication-ready plots
  25. library(smplot2) # Themes and functions for scientific data visualization
  26. library(MASS) # Advanced statistical methods and distributions
  27. library(dplyr) # Data manipulation (filter, mutate, summarise, etc.)
  28. library(stats) # Base R statistical functions (tests, linear models, etc.)
  29. library(ggplot2) # Data visualization using the grammar of graphics
  30. library(grid) # Low-level graphical system for fine plot control
  31. library(GO.db) # Gene Ontology annotation database
  32. library(dplyr) # Data manipulation (loaded again here)
  33. library(stringr) # String manipulation functions
  34. library(cowplot) # Plot arrangement and annotation for ggplot2
  35. library(svglite) # Export plots to SVG format
  36. library(multcomp) # Multiple comparisons and post-hoc tests
  37. library(patchwork) # Combine multiple ggplot2 plots into one layout
  38. library(openxlsx) # Export table in excel files
  40. #Build functions to produce plots------
  42. # Create a gene expression plot showing normalized counts
  43. # across worker conditions, with jittered points, mean ± SE,
  44. # formatted adjusted p-value in the subtitle, and layout control for multi-panel grids.
  45. make_gene_plot <- function(row, i, ncol_grid = 3) {
  46. counts <- plotCounts(
  47. dds50_Q1,
  48. gene = row$gene_id,
  49. intgroup = c("batch","workers"),
  50. returnData = TRUE,
  51. normalized = TRUE
  52. )
  54. adj_pval_txt <- fmt_adj_pval(padj_vals[row$padj_key])
  55. if (nzchar(row$extra_desc)) {
  56. subtitle_txt <- sprintf("(%s, %s, %s)",
  57. row$gene_id,
  58. row$extra_desc,
  59. adj_pval_txt)
  60. } else {
  61. subtitle_txt <- sprintf("(%s, %s)",
  62. row$gene_id,
  63. adj_pval_txt)
  64. }
  65. p <- ggplot(counts,
  66. aes(x = workers, y = count)) +
  67. geom_jitter(
  68. colour = jitter_pt_colour,
  69. fill = jitter_pt_fill,
  70. shape = jitter_pt_shape,
  71. size = jitter_pt_size,
  72. width = 0.1,
  73. height = 0
  74. ) +
  75. stat_summary(
  76. aes(colour = workers),
  77. fun.data = mean_se,
  78. geom = "pointrange",
  79. size = mean_se_pt_size,
  80. show.legend = FALSE
  81. ) +
  82. labs(
  83. title = row$display_name,
  84. subtitle = subtitle_txt,
  85. x = NULL,
  86. y = "Normalized read counts"
  87. ) +
  88. scale_color_manual(values = worker_cols) +
  89. scale_x_discrete(labels = c("Without \nworkers", "With \nworkers")) +
  90. theme_classic(base_size = 18) +
  91. theme(
  92. plot.title = element_text(size = title_txt_size, face = "bold.italic"),
  93. plot.subtitle = element_text(size = subtitle_txt_size),
  94. axis.text.x = element_text(size = x_axis_txt_size),
  95. axis.text.y = element_text(size = y_axis_txt_size),
  96. axis.title.y = element_text(size = axis_title_y_sz),
  97. panel.border = element_blank(),
  98. panel.grid.major = element_blank(),
  99. panel.grid.minor = element_blank(),
  100. axis.line = element_line(colour = "black")
  101. )
  103. col_index <- ((i - 1) %% ncol_grid) + 1
  104. if (col_index != 1) {
  105. p <- p + theme(axis.title.y = element_blank())
  106. }
  107. p
  108. }
  111. # Create a gene expression plot using dataset dds50_3.1,
  112. # displaying normalized counts across three worker conditions,
  113. # adding mean ± SE, dynamic y-axis scaling, optional compact letter display (CLD)
  114. # for statistical group comparisons, and grid-aware y-axis formatting.
  115. make_q3_panel <- function(row, i, ncol_grid = 3, cld_df = NULL) {
  117. counts <- plotCounts(
  118. dds50_3.1,
  119. gene = row$gene_id,
  120. intgroup = c("batch","workers"),
  121. returnData = TRUE,
  122. normalized = TRUE
  123. )
  126. padj_val <- padj_vec_3[row$padj_key]
  127. subtitle_txt <- sprintf("(%s%s, %s)",
  128. row$gene_id,
  129. ifelse(nzchar(row$extra_desc),
  130. paste0(", ", row$extra_desc), ""),
  131. fmt_adj_pval(padj_val))
  134. p <- ggplot(counts, aes(x = workers, y = count)) +
  135. geom_jitter(colour = "grey", fill = "grey",
  136. size = 3, width = 0.1, height = 0) +
  137. labs(title = row$display_name,
  138. subtitle = subtitle_txt,
  139. x = NULL,
  140. y = "Normalized read counts") +
  141. scale_color_manual(values = cb_palette) +
  142. scale_x_discrete(labels = c("Without \nworkers",
  143. "With \nworkers",
  144. "Workers \nremoved")) +
  145. theme_classic(base_size = 18) +
  146. theme(
  147. plot.title = element_text(size = 20, face = "bold.italic",
  148. colour = title_col[i]),
  149. plot.subtitle = element_text(size = 18),
  150. axis.text.x = element_text(size = 18),
  151. axis.text.y = element_text(size = 20),
  152. axis.title.y = element_text(size = 20),
  153. panel.border = element_blank(),
  154. panel.grid = element_blank(),
  155. axis.line = element_line(colour = "black")
  156. ) +
  157. stat_summary(aes(colour = workers),
  158. fun.data = mean_se,
  159. geom = "pointrange",
  160. size = 2,
  161. show.legend = FALSE)
  164. max_y <- max(counts$count, na.rm = TRUE)
  165. p <- p + scale_y_continuous(limits = c(0, max_y * 1.2))
  168. if (!is.null(cld_df) && nrow(cld_df) > 0) {
  169. ypos <- max_y * 1.2
  170. cld_df$ypos <- ypos
  171. p <- p + geom_text(data = cld_df,
  172. aes(x = workers, y = ypos, label = letter),
  173. size = 6, colour = "black")
  174. }
  177. if ((i - 1) %% ncol_grid != 0) {
  178. p <- p + theme(axis.title.y = element_blank())
  179. }
  180. p
  181. }
  184. # Create a gene expression plot using dataset dds50_3.2,
  185. # displaying normalized counts across three worker conditions,
  186. # adding mean ± SE, dynamic y-axis scaling, optional compact letter display (CLD)
  187. # for statistical group comparisons, and grid-aware y-axis formatting.
  188. make_q4_panel <- function(row, i, ncol_grid = 3, cld_df = NULL) {
  189. counts <- plotCounts(
  190. dds50_3.2,
  191. gene = row$gene_id,
  192. intgroup = c("batch","workers"),
  193. returnData = TRUE,
  194. normalized = TRUE
  195. )
  198. padj_val <- padj_vec_4[row$padj_key]
  199. subtitle_txt <- sprintf("(%s%s, %s)",
  200. row$gene_id,
  201. ifelse(nzchar(row$extra_desc),
  202. paste0(", ", row$extra_desc), ""),
  203. fmt_adj_pval(padj_val))
  206. p <- ggplot(counts, aes(x = workers, y = count)) +
  207. geom_jitter(colour = "grey", fill = "grey",
  208. size = 3, width = 0.1, height = 0) +
  209. labs(title = row$display_name,
  210. subtitle = subtitle_txt,
  211. x = NULL,
  212. y = "Normalized read counts") +
  213. scale_color_manual(values = cb_palette_4) +
  214. scale_x_discrete(labels = c("Without\nworkers",
  215. "With\nworkers",
  216. "Workers\nseparated")) +
  217. theme_classic(base_size = 18) +
  218. theme(
  219. plot.title = element_text(size = 20, face = "bold.italic",
  220. colour = title_col[i]),
  221. plot.subtitle = element_text(size = 18),
  222. axis.text.x = element_text(size = 18),
  223. axis.text.y = element_text(size = 20),
  224. axis.title.y = element_text(size = 20),
  225. panel.border = element_blank(),
  226. panel.grid = element_blank(),
  227. axis.line = element_line(colour = "black")
  228. ) +
  229. stat_summary(aes(colour = workers),
  230. fun.data = mean_se,
  231. geom = "pointrange",
  232. size = 2,
  233. show.legend = FALSE)
  236. max_y <- max(counts$count, na.rm = TRUE)
  237. p <- p + scale_y_continuous(limits = c(0, max_y * 1.2))
  240. if (!is.null(cld_df) && nrow(cld_df) > 0 && any(cld_df$letter != "")) {
  241. ypos <- max_y * 1.2
  242. cld_df$ypos <- ypos
  243. p <- p + geom_text(data = cld_df,
  244. aes(x = workers, y = ypos, label = letter),
  245. size = 6, colour = "black")
  246. }
  249. if ((i - 1) %% ncol_grid != 0) {
  250. p <- p + theme(axis.title.y = element_blank())
  251. }
  252. p
  253. }
  255. # Create a gene expression plot using dataset dds50_3.3,
  256. # comparing queens with and without larvae,
  257. # showing jittered normalized counts, mean ± SE,
  258. # formatted adjusted p-values, dynamic y-axis limits,
  259. # and grid-aware y-axis formatting.
  260. make_q5_panel <- function(row, i, ncol_grid = 3) {
  262. counts <- plotCounts(
  263. dds50_3.3,
  264. gene = row$gene_id,
  265. intgroup = c("batch","workers"),
  266. returnData = TRUE,
  267. normalized = TRUE
  268. )
  270. counts$workers <- factor(counts$workers,
  271. levels = c("with_larvae","without_larvae"))
  274. padj_val <- padj_vec_5[row$padj_key]
  275. subtitle_txt <- sprintf("(%s%s, %s)",
  276. row$gene_id,
  277. ifelse(nzchar(row$extra_desc),
  278. paste0(", ", row$extra_desc), ""),
  279. fmt_adj_pval(padj_val))
  282. p <- ggplot(counts, aes(x = workers, y = count)) +
  283. geom_jitter(colour = "grey", fill = "grey",
  284. size = 3, width = 0.1, height = 0) +
  285. stat_summary(aes(colour = workers),
  286. fun.data = mean_se,
  287. geom = "pointrange",
  288. size = 2,
  289. show.legend = FALSE) +
  290. labs(title = row$display_name,
  291. subtitle = subtitle_txt,
  292. x = NULL,
  293. y = "Normalized read counts") +
  294. scale_color_manual(values = cb_palette_5) +
  295. scale_x_discrete(labels = c(
  296. without_larvae = "Without\nlarvae",
  297. with_larvae = "With\nlarvae"
  298. )) +
  299. theme_classic(base_size = 18) +
  300. theme(
  301. plot.title = element_text(size = 20, face = "bold.italic",
  302. colour = title_col[i]),
  303. plot.subtitle = element_text(size = 18),
  304. axis.text.x = element_text( size = 18),
  305. axis.text.y = element_text(size = 20),
  306. axis.title.y = element_text(size = 20),
  307. panel.border = element_blank(),
  308. panel.grid = element_blank(),
  309. axis.line = element_line(colour = "black")
  310. )
  312. max_y <- max(counts$count, na.rm = TRUE)
  313. p <- p + scale_y_continuous(limits = c(0, max_y * 1.2))
  316. if ((i - 1) %% ncol_grid != 0) {
  317. p <- p + theme(axis.title.y = element_blank())
  318. }
  319. p
  320. }
  323. #P-value adjustment function to homogenize naming and threshold
  325. fmt_adj_pval <- function(p) {
  326. if (is.na(p)) {
  327. return("adjusted-pvalue = NA")
  328. }
  329. if (p < 0.0001) {
  330. return("adjusted-pvalue < 0.0001")
  331. }
  333. return(sprintf("adjusted-pvalue = %.2g", p))
  334. }
  337. #Set the working directory-----
  338. #Path to the RDS files
  339. setwd(choose.dir())
  341. #Part 1: Comparison of queens with and without workers-----
  343. ##1.1. Load count table and metadata information------
  345. txi_with_without_workers <- readRDS("txi_with_without_workers.rds")
  347. ##1.2. Differential expression analysis-----
  349. #Create dds with batch (i.e., experiment) and presence of workers in the full model
  350. dds_Q1=DESeqDataSetFromMatrix(countData = txi_with_without_workers$counts,
  351. colData= txi_with_without_workers$meta,
  352. design=~batch+workers)
  354. #To only keep genes expressed in more than half of the samples,
  355. #we only keep rows with at least 25 samples with a non-zero value
  356. dds50_Q1=dds_Q1[rowSums(counts(dds_Q1) >0)>= 25]
  358. #Analysis of the impact of the presence of workers, while controlling for batch,
  359. #by conducting a LRT comparison of the full model with a reduced model only containing batch as explanatory variable
  360. dds_test_Q1=DESeq(dds50_Q1, test="LRT", reduced=~batch)
  362. #Produce the result table
  363. res_Q1=results(dds_test_Q1, alpha=0.05)
  364. summary(res_Q1)
  366. #Extract the genes with adjusted-p<0.05
  367. DEG_Q1=subset(res_Q1, padj<0.05)
  368. rownames(DEG_Q1)
  369. nrow(DEG_Q1)
  370. #Note: positive log-fold = With workers > Without workers
  372. #Produce an Excel table with the output of the analysis for the DEG
  373. #DEG_Q1_Supplementary=DEG_Q1
  374. #DEG_Q1_Supplementary$Gene_names=rownames(DEG_Q1_Supplementary)
  375. #write.xlsx(DEG_Q1_Supplementary, file = "DEG_Q1.xlsx")
  378. ##1.3. Run the GO term analysis-----
  380. GO=read.table("GO_terms.txt")
  381. colnames(GO) = c("Transcripts", "GOterms")
  382. all(GO$Transcripts %in% txi_with_without_workers$tx2gene$transcript)
  384. # Merge GO annotations with transcript-to-gene mapping,
  385. # keeping all rows from GO (left join style)
  386. GO_merged = merge(GO, txi_with_without_workers$tx2gene, by.x = "Transcripts", by.y = "transcript", all.x = TRUE)
  387. colnames(GO_merged) = c("Transcripts", "GOterms", "Genes")
  389. # Create a list mapping each gene to its associated GO terms
  390. gene2GO = split(GO_merged$GOterms, GO_merged$Genes)
  391. length(gene2GO)
  392. head(gene2GO)
  394. #Create a factor where:
  395. #1 means the gene is in the DEG list
  396. #0 means the gene is in the background (all genes with GO annotation)
  398. all_genes = unique(txi_with_without_workers$tx2gene$gene)
  399. all_genes = all_genes[all_genes %in% rownames(res_Q1)]
  400. gene_list = factor(as.integer(all_genes %in% rownames(DEG_Q1)))
  401. names(gene_list) = all_genes
  403. #Create a topGO Object for Biological process
  404. GO_data_Q1_BP = new("topGOdata",
  405. ontology = "BP",
  406. allGenes = gene_list,
  407. annot = annFUN.gene2GO,
  408. gene2GO = gene2GO)
  410. #Run Fisher's Exact Test
  411. resultFisher_GO_Q1_BP = runTest(GO_data_Q1_BP, algorithm = "weight01", statistic = "fisher")
  412. print(resultFisher_GO_Q1_BP)
  414. #Retrieve Significant GO Terms represented by at least three genes in the list of DEG
  415. GO_results_Q1_BP = GenTable(GO_data_Q1_BP, resultFisher_GO_Q1_BP, topNodes = 28)
  416. print(subset(GO_results_Q1_BP, Significant >= 3))
  418. #Create a topGO Object for Molecular function
  419. GO_data_Q1_MF = new("topGOdata",
  420. ontology = "MF",
  421. allGenes = gene_list,
  422. annot = annFUN.gene2GO,
  423. gene2GO = gene2GO)
  425. #Run Fisher's Exact Test
  426. resultFisher_GO_Q1_MF = runTest(GO_data_Q1_MF, algorithm = "weight01", statistic = "fisher")
  427. print(resultFisher_GO_Q1_MF)
  429. #Retrieve Significant GO Terms represented by at least three genes in the list of DEG
  430. GO_results_Q1_MF = GenTable(GO_data_Q1_MF, resultFisher_GO_Q1_MF, topNodes = 9)
  431. print(GO_results_Q1_MF)
  432. print(subset(GO_results_Q1_MF, Significant >= 3))
  434. #Extraction of DEGs with significant GOs
  436. #Biological Process
  437. GO2genes_list_BP = genesInTerm(GO_data_Q1_BP)
  439. GO_signif_Q1_BP = subset(GO_results_Q1_BP, Significant >= 3)
  440. GO_signif_Q1_BP$Genes = sapply(GO_signif_Q1_BP$GO.ID, function(go_id) {
  441. genes = GO2genes_list_BP[[go_id]]
  442. paste(genes, collapse = ", ")
  443. genes = intersect(genes, rownames(DEG_Q1))
  444. })
  445. GO_signif_Q1_BP$Genes
  448. #Molecular function
  449. GO2genes_list_MF = genesInTerm(GO_data_Q1_MF)
  450. GO_signif_Q1_MF = subset(GO_results_Q1_MF, Significant >= 3)
  451. GO_signif_Q1_MF$Genes = sapply(GO_signif_Q1_MF$GO.ID, function(go_id) {
  452. genes = GO2genes_list_MF[[go_id]]
  453. paste(genes, collapse = ", ")
  454. genes = intersect(genes, rownames(DEG_Q1))
  455. })
  456. GO_signif_Q1_MF$Genes
  458. ##1.4. Produce Figure 1-----
  460. #Data preparation
  461. res_Q1_df <- as.data.frame(res_Q1)
  462. res_Q1_df <- res_Q1_df[!is.na(res_Q1_df$padj),]
  463. res_Q1_df$diffexpressed <- "No"
  464. res_Q1_df$diffexpressed[res_Q1_df$log2FoldChange > 0 & res_Q1_df$padj < 0.05] <- "Overexpressed with workers"
  465. res_Q1_df$diffexpressed[res_Q1_df$log2FoldChange < 0 & res_Q1_df$padj < 0.05] <- "Underexpressed with workers"
  466. res_Q1_df$Comparison <- ifelse(res_Q1_df$diffexpressed == "Overexpressed with workers", "Workers superior",
  467. ifelse(res_Q1_df$diffexpressed == "Underexpressed with workers", "Workers inferior", "No"))
  469. #Arrow Preparation
  470. y_arrow <- max(-log10(res_Q1_df$padj), na.rm = TRUE) * 1.05
  472. # Volcano plot
  473. Volcano_plot_final <- ggplot(res_Q1_df, aes(x = log2FoldChange, y = -log10(padj), col = Comparison)) +
  474. geom_point(size = 2, alpha = 0.8) +
  475. geom_hline(yintercept = -log10(0.05), color = "red", linetype = "dashed", size = 0.8) +
  476. geom_vline(xintercept = 0, color = "black", linetype = "dashed", size = 0.8) +
  477. annotate("segment", x = 0, xend = 3.5, y = y_arrow, yend = y_arrow,
  478. color = "#0072B2", size = 1) +
  479. annotate("text", x = 1.75, y = y_arrow + 0.4,
  480. label = "Overexpressed in the brain of queens \nprovided with workers", color = "#0072B2", size = 4) +
  481. annotate("segment", x = 0, xend = -3.5, y = y_arrow, yend = y_arrow,
  482. color = "darkgreen", size = 1) +
  483. annotate("text", x = -1.75, y = y_arrow + 0.4,
  484. label = "Overexpressed in the brain of queens \nkept without workers", color = "darkgreen", size = 4) +
  487. scale_color_manual(values = c("Workers superior" = "#0072B2",
  488. "Workers inferior" = "darkgreen",
  489. "No" = "grey70")) +
  492. scale_x_continuous(limits = c(-3.5, 3.5), breaks = seq(-3, 3, 1)) +
  493. labs(x = expression(log[2]~"fold change"), y = expression(-log[10]~"(adjusted p-value)")) + ##C1:Change p Value name as you wish
  495. theme_minimal(base_size = 12) +
  496. theme(
  497. legend.position = "none",
  498. axis.title = element_text(size = 14),
  499. axis.text = element_text(size = 10)
  500. )
  502. Volcano_plot_final
  504. #Barplot
  506. #Recombining the datasets
  507. GO_signif_Q1_BP$Category <- "Biological processes"
  508. GO_signif_Q1_MF$Category <- "Molecular functions"
  509. GO_signif_Q1_all <- rbind(GO_signif_Q1_BP, GO_signif_Q1_MF)
  510. GO_signif_Q1_all <- GO_signif_Q1_all %>%
  511. group_by(Category) %>%
  512. arrange(desc(Significant)) %>%
  513. ungroup() %>%
  514. mutate(Term = factor(Term, levels = rev(Term)))
  516. n_BP <- nrow(GO_signif_Q1_BP)
  517. n_all <- nrow(GO_signif_Q1_all)
  519. GO_signif_Q1_all <- GO_signif_Q1_all %>%
  520. rowwise() %>%
  521. mutate(
  522. FullTerm = Term(GOTERM[[GO.ID]]) %>% as.character()
  523. ) %>%
  524. ungroup()
  526. #Wrap for long names. Very useful
  527. GO_signif_Q1_all$FullTerm <- str_wrap(GO_signif_Q1_all$FullTerm, width = 30)
  529. GO_signif_Q1_all <- GO_signif_Q1_all %>%
  530. group_by(Category) %>%
  531. arrange(desc(Significant)) %>%
  532. ungroup() %>%
  533. mutate(FullTerm = factor(FullTerm, levels = rev(FullTerm)))
  535. GO_barplot <- ggplot(GO_signif_Q1_all, aes(x = Significant, y = FullTerm, fill = as.numeric(result1))) +
  536. geom_col() +
  537. scale_fill_gradient(low = "grey20", high = "grey90", name = "p-value") + #C1 Check p-value
  538. labs(x = "Number of genes in the list of DEG", y = NULL) +
  539. theme_classic() +
  540. theme(
  541. axis.text.y = element_text(size = 12),
  542. axis.text.x = element_text(size = 10),
  543. axis.title.x = element_text(size = 12)
  544. ) +
  545. geom_hline(yintercept = nrow(GO_signif_Q1_all %>% filter(Category == "Biological processes")) + 0.5, linewidth = 0.8) +
  546. annotate(
  547. "text",
  548. x = 9,
  549. y = nrow(GO_signif_Q1_all %>% filter(Category == "Biological processes")) - 0.1,
  550. label = "Biological processes",
  551. fontface = "bold",
  552. hjust = 0,
  553. size = 3
  554. ) +
  555. annotate(
  556. "text",
  557. x = 9,
  558. y = nrow(GO_signif_Q1_all %>% filter(Category == "Biological processes")) +1.1,
  559. label = "Molecular functions",
  560. fontface = "bold",
  561. hjust = 0,
  562. size = 3
  563. ) +
  564. scale_x_continuous(limits = c(0, max(GO_signif_Q1_all$Significant, 21)),breaks = seq(0,21,by = 3))
  567. GO_barplot
  569. Combined_plot <- ggarrange(
  570. Volcano_plot_final,
  571. GO_barplot,
  572. ncol = 2,
  573. widths = c(1.3, 1),
  574. labels = c("A", "B")
  575. )
  577. #Figure 1 Result
  578. Combined_plot
  580. #Export
  581. #ggsave(
  582. # filename = "Figure_1.svg",
  583. # plot = Combined_plot,
  584. # width = 14,
  585. # height = 7,
  586. # dpi = 600,
  587. # bg = "white",
  588. # device = "svg")
  591. #Part 2: Identification of genes of interests for the social regulation of queen brood care behavior-----
  592. ##2.1. List of genes of interest-----
  594. DEG_Q1_Genes_of_Interest=DEG_Q1
  595. DEG_Q1_Genes_of_Interest$Gene_name=rownames(DEG_Q1_Genes_of_Interest)
  597. Krh1=subset(DEG_Q1_Genes_of_Interest, Gene_name=="Lnig_g01448")
  598. dPDK1=subset(DEG_Q1_Genes_of_Interest, Gene_name=="Lnig_g11740")
  599. Insulin=subset(DEG_Q1_Genes_of_Interest, Gene_name=="Lnig_g17128")
  600. Juv_esterase=subset(DEG_Q1_Genes_of_Interest, Gene_name=="Lnig_g08195")
  601. E4=subset(DEG_Q1_Genes_of_Interest, Gene_name=="Lnig_g15862")
  602. Dsx=subset(DEG_Q1_Genes_of_Interest, Gene_name=="Lnig_g04508")
  603. Henna=subset(DEG_Q1_Genes_of_Interest, Gene_name=="Lnig_g13064")
  604. Prohormone=subset(DEG_Q1_Genes_of_Interest, Gene_name=="Lnig_g06899")
  605. GP9=subset(DEG_Q1_Genes_of_Interest, Gene_name=="gLnigOBP1")
  606. GR81= subset(DEG_Q1_Genes_of_Interest, Gene_name=="gLnigGR81")
  607. Tret1a=subset(DEG_Q1_Genes_of_Interest, Gene_name=="Lnig_g01559")
  608. Tret1b=subset(DEG_Q1_Genes_of_Interest, Gene_name=="Lnig_g11948")
  609. Tret1_2=subset(DEG_Q1_Genes_of_Interest, Gene_name=="Lnig_g10144")
  611. ##2.2. Produce Supplementary Figure 1 -----
  613. plot_info <- data.frame(
  614. gene_id = c("Lnig_g01448","Lnig_g11740","Lnig_g17128","Lnig_g08195",
  615. "Lnig_g15862","Lnig_g04508","Lnig_g13064",
  616. "Lnig_g06899","gLnigOBP1","gLnigGR81",
  617. "Lnig_g01559","Lnig_g11948","Lnig_g10144"),
  618. padj_key = c("Krh1","dPDK1","Insulin","JuvEsterase",
  619. "E4","Dsx","Henna",
  620. "Prohormone","GP9","GR81",
  621. "Tret1a","Tret1b","Tret1_2"),
  622. display_name = c("Kr-h1","PDK1","Insulin","Juvenile hormone esterase",
  623. "Esterase E4","Protein doublesex",
  624. "Protein henna","Prohormone-2",
  625. "Pheromone binding protein GP9",
  626. "GR81","Tret1a","Tret1b","Tret1-2 homolog"
  627. ),
  628. extra_desc = c("", "", "", "", "", "", "", "", "", "", "", "", ""),
  629. stringsAsFactors = FALSE
  630. )
  633. Genes_of_interest_Overexpressed_in_worker_presence= c("Lnig_g01448", "Lnig_g17128", "Lnig_g08195", "Lnig_g11948", "Lnig_g11740", "Lnig_g04508")
  634. Genes_of_interest_Overexpressed_in_worker_absence= c( "Lnig_g15862", "Lnig_g13064", "gLnigOBP1", "gLnigGR81", "Lnig_g01559", "Lnig_g10144", "Lnig_g06899")
  636. res_Q1_plot=as.data.frame(res_Q1)
  637. res_Q1_plot$gene=rownames(res_Q1_plot)
  638. res_Q1_plot$Interest=NA
  639. for (i in seq_len(nrow(res_Q1_plot))) {
  640. gene_name <- rownames(res_Q1_plot)[i]
  641. if(gene_name %in% Genes_of_interest_Overexpressed_in_worker_presence) {
  642. res_Q1_plot$Interest[i]="Genes of interest Overexpressed in worker presence"
  643. } else if (gene_name %in% Genes_of_interest_Overexpressed_in_worker_absence) {
  644. res_Q1_plot$Interest[i] = "Genes of interest Overexpressed in worker absence"
  645. }
  646. }
  648. res_Q1_plot <- res_Q1_plot %>%
  649. left_join(plot_info[, c("gene_id", "display_name")],
  650. by = c("gene" = "gene_id"))
  652. #Volcano plot with genes of interests
  654. Volcano_plot_Supplementary <- ggplot(data=res_Q1_plot, aes(x=log2FoldChange, y=-log10(padj), col=Interest)) +
  655. geom_point(size = 2, alpha = 0.8) +
  656. geom_hline(yintercept = -log10(0.05), color = "red", linetype = "dashed", size = 0.8) +
  657. geom_vline(xintercept = 0, color = "black", linetype = "dashed", size = 0.8) +
  658. geom_text_repel(data=res_Q1_plot[res_Q1_plot$padj < 0.05 & res_Q1_plot$Interest!="NA", ], aes(label=display_name),
  659. size=5, max.overlaps=15) +
  660. scale_color_manual(values=c("darkgreen","#0072B2"),na.value="lightgrey") +
  661. annotate("segment", x = 0, xend = 3.5, y = y_arrow, yend = y_arrow,
  662. color = "#0072B2", size = 1) +
  663. annotate("text", x = 1.75, y = y_arrow + 0.4,
  664. label = "Overexpressed in the brain of queens \nprovided with workers", color = "#0072B2", size = 4) +
  665. annotate("segment", x = 0, xend = -3.5, y = y_arrow, yend = y_arrow,
  666. color = "darkgreen", size = 1) +
  667. annotate("text", x = -1.75, y = y_arrow + 0.4,
  668. label = "Overexpressed in the brain of queens \nkept without workers", color = "darkgreen", size = 4) +
  671. scale_x_continuous(limits = c(-3.5, 3.5), breaks = seq(-3, 3, 1)) +
  672. labs(x = expression(log[2]~"fold change"), y = expression(-log[10]~"(adjusted p-value)")) +
  674. theme_minimal(base_size = 12) +
  675. theme(
  676. legend.position = "none",
  677. axis.title = element_text(size = 14),
  678. axis.text = element_text(size = 10)
  679. )
  681. Volcano_plot_Supplementary
  683. #ggsave(
  684. # filename = "Lasius_Supplementary_Figure_1.svg",
  685. # plot = Volcano_plot_Supplementary,
  686. # width = 12,
  687. # height = 7,
  688. # dpi = 600,
  689. # bg = "white",
  690. # device = "svg")
  692. ##2.3. Produce Figure 2-----
  694. padj_vals <- c(
  695. Krh1 = Krh1$padj,
  696. dPDK1 = dPDK1$padj,
  697. Insulin = Insulin$padj,
  698. JuvEsterase= Juv_esterase$padj,
  699. E4 = E4$padj,
  700. Dsx = Dsx$padj,
  701. GP9 = GP9$padj,
  702. GR81 = GR81$padj,
  703. Tret1a = Tret1a$padj,
  704. Tret1b = Tret1b$padj,
  705. Tret1_2 = Tret1_2$padj,
  706. Henna = Henna$padj,
  707. Prohormone = Prohormone$padj
  708. )
  711. worker_cols <- c("darkgreen","#0072B2")
  712. jitter_pt_size <- 3
  713. jitter_pt_colour <- "grey"
  714. jitter_pt_fill <- "grey"
  715. jitter_pt_shape <- 21
  716. mean_se_pt_size <- 2
  717. x_axis_txt_size <- 20
  718. y_axis_txt_size <- 20
  719. title_txt_size <- 20
  720. subtitle_txt_size<- 18
  721. axis_title_y_sz <- 20
  722. ncol_grid <- 3
  724. all_plots <- lapply(seq_len(nrow(plot_info)), function(i) {
  725. make_gene_plot(plot_info[i, ], i, ncol_grid = ncol_grid)
  726. })
  729. final_figure <- plot_grid(
  730. plotlist = all_plots,
  731. nrow = ceiling(length(all_plots) / ncol_grid),
  732. ncol = ncol_grid,
  733. labels = LETTERS[seq_along(all_plots)],
  734. label_size = 20,
  735. label_fontface = "bold",
  736. align = "hv"
  737. )
  739. final_figure
  741. #Export
  742. #ggsave(
  743. # filename = "Figure_2.svg",
  744. # plot = final_figure,
  745. # width = 50,
  746. # height = 75,
  747. # units = "cm",
  748. # device = "svg")
  750. #Part 3: Comparison of the genes of interest among queens with workers, without workers and with workers removed-----
  753. ##3.1. Load count table and metadata information------
  754. #Only the brain samples for experiment E11 are used in this analysis
  756. txi_with_without_removed_workers <- readRDS("txi_with_without_removed_workers.rds")
  758. ##3.2. Analysis of expression differences-----
  760. #Create dds
  761. dds_3.1=DESeqDataSetFromMatrix(countData = txi_with_without_removed_workers$counts,
  762. colData = txi_with_without_removed_workers$meta,
  763. design=~workers)
  765. #To only keep genes expressed in at least half of the samples,
  766. #we only keep rows at least 19 samples with a non-zero value
  767. dds50_3.1=dds_3.1[rowSums(counts(dds_3.1) >0)>= 19]
  769. #Normalized counts extraction
  770. dds_test_3.1=DESeq(dds50_3.1, test="LRT", reduced=~1)
  771. normalized_counts.3.1 <- counts(dds_test_3.1, normalized = TRUE)
  773. #Pivot longer the data frame
  774. normalized_counts_long.3.1 <- normalized_counts.3.1 %>%
  775. as.data.frame() %>%
  776. rownames_to_column(var = "gene") %>%
  777. pivot_longer(
  778. cols = -gene,
  779. names_to = "sample",
  780. values_to = "normalized_count"
  781. )
  783. normalized_counts_long.3.1 <- normalized_counts_long.3.1 %>%
  784. left_join(as.data.frame(txi_with_without_removed_workers$meta) %>% dplyr::select(sample, workers), by = "sample")
  786. #Comparison of the expression levels for each of the genes of interest
  788. ###Kr-h1
  789. #Subset Kr-h1 counts
  790. Krup.3.1=subset(normalized_counts_long.3.1, gene=="Lnig_g01448")
  792. #Model
  793. modK.3.1=lm(log(normalized_count+1)~workers, data=Krup.3.1)
  794. simulationOutput = simulateResiduals(fittedModel= modK.3.1, plot = T)
  796. #Results
  797. car::Anova(modK.3.1)
  798. emmeans(modK.3.1,pairwise~workers)
  801. ###PDK1
  802. #Subset the PDK1 counts
  803. dPDK1.3.1 =subset(normalized_counts_long.3.1, gene=="Lnig_g11740")
  805. #Model
  806. moddPDK1.3.1=lm(normalized_count~workers, data=dPDK1.3.1)
  807. simulationOutput = simulateResiduals(fittedModel= moddPDK1.3.1, plot = T)
  809. #Results
  810. car::Anova(moddPDK1.3.1)
  811. emmeans(moddPDK1.3.1,pairwise~workers)
  813. ###Insulin
  814. #Subset the Insulin counts
  815. Ins.3.1 =subset(normalized_counts_long.3.1, gene=="Lnig_g17128")
  817. #Model
  818. modIns.3.1=lm(normalized_count~workers, data=Ins.3.1)
  819. simulationOutput = simulateResiduals(fittedModel= modIns.3.1, plot = T)
  821. ##Results
  822. car::Anova(modIns.3.1)
  823. emmeans(modIns.3.1,pairwise~workers)
  825. ###Juvenile hormone esterase
  826. #Subset the Juvenile hormone esterase counts
  827. Carb.3.1 =subset(normalized_counts_long.3.1, gene=="Lnig_g08195")
  829. #Model
  830. modCarb.3.1=lm(normalized_count~workers, data=Carb.3.1)
  831. simulationOutput = simulateResiduals(fittedModel= modCarb.3.1, plot = T)
  833. #Results
  834. car::Anova(modCarb.3.1)
  835. emmeans(modCarb.3.1,pairwise~workers)
  838. ###Esterase E4
  839. #Subset the Esterase E4 counts
  840. E4.3.1 =subset(normalized_counts_long.3.1, gene=="Lnig_g15862")
  842. #Model
  843. modE4.3.1=lm(normalized_count~workers, data=E4.3.1)
  844. simulationOutput = simulateResiduals(fittedModel= modE4.3.1, plot = T)
  846. #Results
  847. car::Anova(modE4.3.1)
  848. emmeans(modE4.3.1,pairwise~workers)
  851. ###Protein doublesex
  852. #Subset the Doublesex counts
  853. Dbs.3.1 =subset(normalized_counts_long.3.1, gene=="Lnig_g04508")
  855. #Model
  856. modDbs.3.1=lm(normalized_count~workers, data=Dbs.3.1)
  857. simulationOutput = simulateResiduals(fittedModel= modDbs.3.1, plot = T)
  859. #Results
  860. car::Anova(modDbs.3.1)
  863. ###Protein henna
  864. #Subset the Tret1-2 homolog counts
  865. Henna.3.1=subset(normalized_counts_long.3.1, gene=="Lnig_g13064")
  867. #Model
  868. modHenna.3.1=lm(normalized_count~workers, data=Henna.3.1)
  869. simulationOutput = simulateResiduals(fittedModel= modHenna.3.1, plot = T)
  871. #Results
  872. car::Anova(modHenna.3.1)
  873. emmeans(modHenna.3.1,pairwise~workers)
  876. ###Prohormone-2
  877. #Subset the Prohormone-2 counts
  878. Pro.3.1 =subset(normalized_counts_long.3.1, gene=="Lnig_g06899")
  880. #Model
  881. modPro.3.1=lm(normalized_count~workers, data=Pro.3.1)
  882. simulationOutput = simulateResiduals(fittedModel= modPro.3.1, plot = T)
  884. #Results
  885. car::Anova(modPro.3.1)
  886. emmeans(modPro.3.1,pairwise~workers)
  889. ###Pheromone-binding protein Gp-9
  890. #Subset the GP9 counts
  891. GP9.3.1=subset(normalized_counts_long.3.1, gene=="gLnigOBP1")
  893. #Model
  894. modGP9.3.1=lm(normalized_count~workers, data=GP9.3.1)
  895. simulationOutput = simulateResiduals(fittedModel= modGP9.3.1, plot = T)
  897. #Results
  898. car::Anova(modGP9.3.1)
  899. emmeans(modGP9.3.1,pairwise~workers)
  901. ###GR81
  902. #Subset the GR81 counts
  903. GR81.3.1 =subset(normalized_counts_long.3.1, gene=="gLnigGR81")
  905. #Model
  906. modGR81.3.1=lm(normalized_count~workers, data=GR81.3.1)
  907. simulationOutput = simulateResiduals(fittedModel= modGR81.3.1, plot = T)
  909. #Results
  910. car::Anova(modGR81.3.1)
  911. emmeans(modGR81.3.1,pairwise~workers)
  914. ###Tret1a
  915. #Subset the Tret1a counts
  916. Tret1a.3.1 =subset(normalized_counts_long.3.1, gene=="Lnig_g01559")
  918. #Model
  919. modTret1a.3.1=lm(normalized_count~workers, data=Tret1a.3.1)
  920. simulationOutput = simulateResiduals(fittedModel= modTret1a.3.1, plot = T)
  922. #Results
  923. car::Anova(modTret1a.3.1)
  926. ###Tret1b
  927. #Subset the Tret1b counts
  928. Tret1b.3.1 =subset(normalized_counts_long.3.1, gene=="Lnig_g11948")
  930. #Model
  931. modTret1b.3.1=lm(normalized_count~workers, data=Tret1b.3.1)
  932. simulationOutput = simulateResiduals(fittedModel= modTret1b.3.1, plot = T)
  934. #Results
  935. car::Anova(modTret1b.3.1)
  936. emmeans(modTret1b.3.1,pairwise~workers)
  939. ###Tret1-2 homolog
  940. #Subset the Tret1-2 homolog counts
  941. Tret1_2.3.1=subset(normalized_counts_long.3.1, gene=="Lnig_g10144")
  943. #Model
  944. modTret1_2.3.1=lm(normalized_count~workers, data=Tret1_2.3.1)
  945. simulationOutput = simulateResiduals(fittedModel= modTret1_2.3.1, plot = T)
  947. #Results
  948. car::Anova(modTret1_2.3.1)
  949. emmeans(modTret1_2.3.1,pairwise~workers)
  951. ##3.3. Produce Supplementary Figure 2----------
  953. #data
  954. models <- list(
  955. Krh1 = modK.3.1,
  956. dPDK1 = moddPDK1.3.1,
  957. Insulin = modIns.3.1,
  958. JuvEsterase= modCarb.3.1,
  959. E4 = modE4.3.1,
  960. Dsx = modDbs.3.1,
  961. GP9 = modGP9.3.1,
  962. GR81 = modGR81.3.1,
  963. Tret1a = modTret1a.3.1,
  964. Tret1b = modTret1b.3.1,
  965. Tret1_2 = modTret1_2.3.1,
  966. Henna = modHenna.3.1,
  967. Prohormone = modPro.3.1
  968. )
  970. #Color of titles and means. I did it manually here. Recheck important
  971. title_pat <- gsub(" ", "", "brb brr rbb rrr b")
  972. title_vec <- strsplit(title_pat, "")[[1]]
  973. title_col <- ifelse(title_vec == "r", "darkgrey", "black")
  974. cb_palette <- c("darkgreen", "#0072B2", "#E69F00")
  978. padj_3.1 = list(
  979. Krh1 = car::Anova(modK.3.1)[["Pr(>F)"]][1],
  980. dPDK1 = car::Anova(moddPDK1.3.1)[["Pr(>F)"]][1],
  981. Insulin = car::Anova(modIns.3.1)[["Pr(>F)"]][1],
  982. JuvEsterase = car::Anova(modCarb.3.1)[["Pr(>F)"]][1],
  983. E4 = car::Anova(modE4.3.1)[["Pr(>F)"]][1],
  984. Dsx = car::Anova(modDbs.3.1)[["Pr(>F)"]][1],
  985. GP9 = car::Anova(modGP9.3.1)[["Pr(>F)"]][1],
  986. GR81 = car::Anova(modGR81.3.1)[["Pr(>F)"]][1],
  987. Tret1a = car::Anova(modTret1a.3.1)[["Pr(>F)"]][1],
  988. Tret1b = car::Anova(modTret1b.3.1)[["Pr(>F)"]][1],
  989. Tret1_2 = car::Anova(modTret1_2.3.1)[["Pr(>F)"]][1],
  990. Henna = car::Anova(modHenna.3.1)[["Pr(>F)"]][1],
  991. Prohormone = car::Anova(modPro.3.1)[["Pr(>F)"]][1])
  993. padj_vec_3 <- setNames(padj_3.1, names(padj_3.1))
  994. names(padj_vec_3)[names(padj_vec_3) == "Doublesex"] <- "Dsx"
  995. padj_vec_3 <- unlist(padj_vec_3, use.names = TRUE)
  998. #Letter build
  999. raw_cld_list <- vector("list", nrow(plot_info))
  1001. for (i in seq_len(nrow(plot_info))) {
  1002. padj_val <- padj_vec_3[plot_info$padj_key[i]]
  1003. if (!is.na(padj_val) && padj_val < 0.05) {
  1005. tryCatch({
  1006. emm <- emmeans(models[[plot_info$padj_key[i]]], pairwise ~ workers)
  1007. cld <- cld(emm)
  1010. cld$.group <- toupper(trimws(cld$.group))
  1013. num2letter <- function(s) {
  1014. s <- gsub(" ", "", s)
  1015. if (s == "") return("")
  1016. digits <- strsplit(s, "")[[1]]
  1017. letters <- sapply(digits,
  1018. function(d) {
  1019. if (d %in% as.character(1:26)) LETTERS[as.integer(d)] else d
  1020. })
  1021. paste0(letters, collapse = "")
  1022. }
  1023. cld$.group <- sapply(cld$.group, num2letter)
  1025. raw_cld_list[[i]] <- data.frame(
  1026. gene_id = plot_info$gene_id[i],
  1027. workers = cld$workers,
  1028. letter = cld$.group,
  1029. stringsAsFactors = FALSE
  1030. )
  1031. }, error = function(e) {
  1032. raw_cld_list[[i]] <<- NULL
  1033. message("CLD failed for ", plot_info$gene_id[i], ": ", e$message)
  1034. })
  1035. } else {
  1036. raw_cld_list[[i]] <- NULL
  1037. }
  1038. }
  1040. cld_list <- raw_cld_list
  1042. #Manual adjustments to make without workers group assignment always A
  1043. cld_list[[1]]$letter <- c("C", "A", "B")
  1044. cld_list[[5]]$letter <- c("B", "B", "A")
  1045. cld_list[[7]]$letter <- c("B", "B", "A")
  1046. cld_list[[8]]$letter <- c("B", "A", "A")
  1047. cld_list[[10]]$letter <- c("B", "AB", "A")
  1048. cld_list[[13]]$letter <- c("B", "A", "A")
  1051. if (length(title_col) < nrow(plot_info))
  1052. title_col <- rep(title_col, length.out = nrow(plot_info))
  1055. all_plots <- lapply(seq_len(nrow(plot_info)), function(i) {
  1056. make_q3_panel(plot_info[i, ], i,
  1057. ncol_grid = 3,
  1058. cld_df = cld_list[[i]])
  1059. })
  1062. final_figure <- plot_grid(
  1063. plotlist = all_plots,
  1064. nrow = 5,
  1065. ncol = 3,
  1066. labels = LETTERS,
  1067. label_size = 20,
  1068. label_fontface = "bold",
  1069. align = "hv"
  1070. )
  1072. print(final_figure)
  1074. #ggsave(
  1075. # filename = "Supplementary_Figure_2.svg",
  1076. # plot = final_figure,
  1077. # width = 50,
  1078. # height = 75,
  1079. # units = "cm",
  1080. # device = "svg"
  1081. #)
  1084. #Part 4: Comparison of the genes of interest among queens with workers, without workers and with workers separated-----
  1086. ##4.1. Load count table and metadata information------
  1087. #Only the brain samples for experiment E4 are used in this analysis
  1089. txi_with_without_separated_workers <- readRDS("txi_with_without_separated_workers.rds")
  1091. ##4.2. Analysis of expression differences-----
  1093. #Create dds
  1094. dds_3.2=DESeqDataSetFromMatrix(countData = txi_with_without_separated_workers$counts,
  1095. colData=txi_with_without_separated_workers$meta,
  1096. design=~workers)
  1099. #To only keep genes expressed in at least half of the samples, we only keep rows at least 19 samples with a non-zero value
  1100. dds50_3.2=dds_3.2[rowSums(counts(dds_3.2) >0)>= 19]
  1102. #Normalized counts extraction
  1103. dds_test_3.2=DESeq(dds50_3.2, test="LRT", reduced=~1)
  1104. normalized_counts.3.2 = counts(dds_test_3.2, normalized = TRUE)
  1107. #Pivot longer the dataframe
  1108. normalized_counts_long.3.2 = normalized_counts.3.2 %>%
  1109. as.data.frame() %>%
  1110. rownames_to_column(var = "gene") %>%
  1111. pivot_longer(
  1112. cols = -gene,
  1113. names_to = "sample",
  1114. values_to = "normalized_count"
  1115. )
  1117. normalized_counts_long.3.2 = normalized_counts_long.3.2 %>%
  1118. left_join(as.data.frame(txi_with_without_separated_workers$meta) %>% dplyr::select(sample, workers), by = "sample")
  1121. ###Kr-h1
  1122. #Subset kr-h1 counts
  1123. Krup.3.2=subset(normalized_counts_long.3.2, gene=="Lnig_g01448")
  1125. #Model
  1126. modK.3.2=lm(normalized_count~workers, data=Krup.3.2)
  1127. simulationOutput = simulateResiduals(fittedModel= modK.3.2, plot = T)
  1129. #Results
  1130. car::Anova(modK.3.2)
  1131. emmeans(modK.3.2,pairwise~workers)
  1133. ###PDK1
  1134. #Subset the PDK1 counts
  1135. dPDK1.3.2 =subset(normalized_counts_long.3.2, gene=="Lnig_g11740")
  1137. #Model
  1138. moddPDK1.3.2=lm(normalized_count~workers, data=dPDK1.3.2)
  1139. simulationOutput = simulateResiduals(fittedModel= moddPDK1.3.2, plot = T)
  1141. #Results
  1142. car::Anova(moddPDK1.3.2)
  1143. emmeans(moddPDK1.3.2,pairwise~workers)
  1145. ###Insulin
  1146. #Subset the Insulin counts
  1147. Ins.3.2 =subset(normalized_counts_long.3.2, gene=="Lnig_g17128")
  1149. #Model
  1150. modIns.3.2=lm(normalized_count~workers, data=Ins.3.2)
  1151. simulationOutput = simulateResiduals(fittedModel= modIns.3.2, plot = T)
  1153. ##Results
  1154. car::Anova(modIns.3.2)
  1155. emmeans(modIns.3.2,pairwise~workers)
  1158. ###Juvenile hormone esterase
  1159. #Subset the Juvenile hormone esterase counts
  1160. Carb.3.2 =subset(normalized_counts_long.3.2, gene=="Lnig_g08195")
  1162. #Model
  1163. modCarb.3.2=lm(normalized_count~workers, data=Carb.3.2)
  1164. simulationOutput = simulateResiduals(fittedModel= modCarb.3.2, plot = T)
  1166. #Results
  1167. car::Anova(modCarb.3.2)
  1168. emmeans(modCarb.3.2,pairwise~workers)
  1171. ###Esterase E4
  1172. #Subset the Esterase E4 counts
  1173. E4.3.2 =subset(normalized_counts_long.3.2, gene=="Lnig_g15862")
  1175. #Model
  1176. modE4.3.2=lm(normalized_count~workers, data=E4.3.2)
  1177. simulationOutput = simulateResiduals(fittedModel= modE4.3.2, plot = T)
  1179. #Results
  1180. car::Anova(modE4.3.2)
  1183. ###Protein doublesex
  1184. #Subset the Doublesex counts
  1185. Dbs.3.2 =subset(normalized_counts_long.3.2, gene=="Lnig_g04508")
  1187. #Model
  1188. modDbs.3.2=lm(normalized_count~workers, data=Dbs.3.2)
  1189. simulationOutput = simulateResiduals(fittedModel= modDbs.3.2, plot = T)
  1191. #Results
  1192. car::Anova(modDbs.3.2)
  1193. emmeans(modDbs.3.2,pairwise~workers)
  1196. ###Protein henna
  1197. #Subset the Protein henna counts
  1198. Henna.3.2 =subset(normalized_counts_long.3.2, gene=="Lnig_g13064")
  1200. #Model
  1201. modHenna.3.2=lm(log(normalized_count+1)~workers, data=Henna.3.2)
  1202. simulationOutput = simulateResiduals(fittedModel= modHenna.3.2, plot = T)
  1204. #Results
  1205. car::Anova(modHenna.3.2)
  1206. emmeans(modHenna.3.2,pairwise~workers)
  1209. ###Prohormone-2
  1210. #Subset the Prohormone-2 counts
  1211. Pro.3.2 =subset(normalized_counts_long.3.2, gene=="Lnig_g06899")
  1213. #Model
  1214. modPro.3.2=lm(normalized_count~workers, data=Pro.3.2)
  1215. simulationOutput = simulateResiduals(fittedModel= modPro.3.2, plot = T)
  1217. #Results
  1218. car::Anova(modPro.3.2)
  1219. emmeans(modPro.3.2,pairwise~workers)
  1222. ###Pheromone binding protein GP-9
  1223. #Subset the GP9 counts
  1224. GP9.3.2=subset(normalized_counts_long.3.2, gene=="gLnigOBP1")
  1226. #Model
  1227. modGP9.3.2=lm(normalized_count~workers, data=GP9.3.2)
  1228. simulationOutput = simulateResiduals(fittedModel= modGP9.3.2, plot = T)
  1230. #Results
  1231. car::Anova(modGP9.3.2)
  1234. ###GR81
  1235. #Subset the GR81 counts
  1236. GR81.3.2=subset(normalized_counts_long.3.2, gene=="gLnigGR81")
  1238. #Model
  1239. modGR81.3.2=lm(normalized_count~workers, data=GR81.3.2)
  1240. simulationOutput = simulateResiduals(fittedModel= modGR81.3.2, plot = T)
  1242. #Results
  1243. car::Anova(modGR81.3.2)
  1244. emmeans(modGR81.3.2,pairwise~workers)
  1247. ###Tret1a
  1248. #Subset the Tret1a counts
  1249. Tret1a.3.2 =subset(normalized_counts_long.3.2, gene=="Lnig_g01559")
  1251. #Model
  1252. modTret1a.3.2=lm(normalized_count~workers, data=Tret1a.3.2)
  1253. simulationOutput = simulateResiduals(fittedModel= modTret1a.3.2, plot = T)
  1255. #Results
  1256. car::Anova(modTret1a.3.2)
  1257. emmeans(modTret1a.3.2,pairwise~workers)
  1260. ###Tret1b
  1261. #Subset the Tret1b counts
  1262. Tret1b.3.2 =subset(normalized_counts_long.3.2, gene=="Lnig_g11948")
  1264. #Model
  1265. modTret1b.3.2=lm(log(normalized_count+1)~workers, data=Tret1b.3.2)
  1266. simulationOutput = simulateResiduals(fittedModel= modTret1b.3.2, plot = T)
  1268. #Results
  1269. car::Anova(modTret1b.3.2)
  1270. emmeans(modTret1b.3.2,pairwise~workers)
  1273. ###Tret1-2 homolog
  1274. #Subset the Tret1-2 homolog counts
  1275. Tret1_2.3.2=subset(normalized_counts_long.3.2, gene=="Lnig_g10144")
  1277. #Model
  1278. modTret1_2.3.2=lm(normalized_count~workers, data=Tret1_2.3.2)
  1279. simulationOutput = simulateResiduals(fittedModel= modTret1_2.3.2, plot = T)
  1281. #Results
  1282. car::Anova(modTret1_2.3.2)
  1283. emmeans(modTret1_2.3.2, pairwise~workers)
  1286. ##4.3. Produce Supplementary Figure 4-----
  1288. #Again manually. Important to check
  1289. title_pat <- gsub(" ", "", "bbb brb brr rbb b")
  1290. title_vec <- strsplit(title_pat, "")[[1]]
  1291. title_col <- ifelse(title_vec == "r", "darkgrey", "black")
  1292. title_col <- rep(title_col, length.out = 13)
  1293. cb_palette_4 <- c("darkgreen", "#0072B2", "purple")
  1296. models_4 <- list(
  1297. Krh1 = modK.3.2,
  1298. dPDK1 = moddPDK1.3.2,
  1299. Insulin = modIns.3.2,
  1300. JuvEsterase = modCarb.3.2,
  1301. E4 = modE4.3.2,
  1302. Dsx = modDbs.3.2,
  1303. Henna = modHenna.3.2,
  1304. Prohormone = modPro.3.2,
  1305. GP9 = modGP9.3.2,
  1306. GR81 = modGR81.3.2,
  1307. Tret1a = modTret1a.3.2,
  1308. Tret1b = modTret1b.3.2,
  1309. Tret1_2 = modTret1_2.3.2
  1310. )
  1312. padj_3.2 = list(
  1313. Krh1 = car::Anova(modK.3.2)[["Pr(>F)"]][1],
  1314. dPDK1 = car::Anova(moddPDK1.3.2)[["Pr(>F)"]][1],
  1315. Insulin = car::Anova(modIns.3.2)[["Pr(>F)"]][1],
  1316. JuvEsterase = car::Anova(modCarb.3.2)[["Pr(>F)"]][1],
  1317. E4 = car::Anova(modE4.3.2)[["Pr(>F)"]][1],
  1318. Dsx = car::Anova(modDbs.3.2)[["Pr(>F)"]][1],
  1319. Henna = car::Anova(modHenna.3.2)[["Pr(>F)"]][1],
  1320. Prohormone = car::Anova(modPro.3.2)[["Pr(>F)"]][1],
  1321. GP9 = car::Anova(modGP9.3.2)[["Pr(>F)"]][1],
  1322. GR81 = car::Anova(modGR81.3.2)[["Pr(>F)"]][1],
  1323. Tret1a = car::Anova(modTret1a.3.2)[["Pr(>F)"]][1],
  1324. Tret1b = car::Anova(modTret1b.3.2)[["Pr(>F)"]][1],
  1325. Tret1_2 = car::Anova(modTret1_2.3.2)[["Pr(>F)"]][1])
  1327. padj_vec_4 <- setNames(padj_3.2, names(padj_3.2))
  1328. padj_vec_4 <- unlist(padj_vec_4, use.names = TRUE)
  1331. #Letter build
  1332. raw_cld_list_4 <- vector("list", nrow(plot_info))
  1334. for (i in seq_len(nrow(plot_info))) {
  1335. padj_val <- padj_vec_4[plot_info$padj_key[i] ]
  1337. if (!is.na(padj_val) && padj_val < 0.05) {
  1338. cld_tmp <- tryCatch({
  1339. emm <- emmeans(models_4[[ plot_info$padj_key[i] ]], pairwise ~ workers)
  1340. cld <- cld(emm)
  1341. cld$.group <- toupper(trimws(cld$.group))
  1344. num2letter <- function(s) {
  1345. s <- gsub(" ", "", s)
  1346. if (s == "") return("")
  1347. digits <- strsplit(s, "")[[1]]
  1348. letters <- sapply(digits,
  1349. function(d) {
  1350. if (d %in% as.character(1:26))
  1351. LETTERS[as.integer(d)] else d
  1352. })
  1353. paste0(letters, collapse = "")
  1354. }
  1355. cld$.group <- sapply(cld$.group, num2letter)
  1357. data.frame(
  1358. gene_id = plot_info$gene_id[i],
  1359. workers = cld$workers,
  1360. letter = cld$.group,
  1361. stringsAsFactors = FALSE
  1362. )
  1363. }, error = function(e) {
  1364. message("CLD failed for ", plot_info$gene_id[i], ": ", e$message)
  1365. NULL
  1366. })
  1369. raw_cld_list_4[[i]] <- if (!is.null(cld_tmp)) cld_tmp else NULL
  1370. } else {
  1372. raw_cld_list_4[[i]] <- NULL
  1373. }
  1374. }
  1377. cld_list_4 <- raw_cld_list_4
  1378. #Manual adjustments to make without workers group assigment always A
  1379. cld_list_4[[7]]$letter <- c("B", "A", "A")
  1380. cld_list_4[[8]]$letter <- c("B", "AB", "A")
  1381. cld_list_4[[10]]$letter <- c("B", "AB", "A")
  1382. cld_list_4[[11]]$letter <- c("A", "B", "B")
  1383. cld_list_4[[13]]$letter <- c("B", "A", "A")
  1385. if (length(title_col) < nrow(plot_info))
  1386. title_col <- rep(title_col, length.out = nrow(plot_info))
  1389. all_plots_4 <- lapply(seq_len(nrow(plot_info)), function(i) {
  1390. make_q4_panel(plot_info[i, ], i,
  1391. ncol_grid = 3,
  1392. cld_df = cld_list_4[[i]])
  1393. })
  1395. final_figure_4 <- plot_grid(
  1396. plotlist = all_plots_4,
  1397. nrow = 5,
  1398. ncol = 3,
  1399. labels = LETTERS[1:13],
  1400. label_size = 20,
  1401. label_fontface = "bold",
  1402. align = "hv"
  1403. )
  1405. print(final_figure_4)
  1407. #ggsave(
  1408. # filename = "Supplementary_Figure_4.svg",
  1409. # plot = final_figure_4,
  1410. # width = 50,
  1411. # height = 75,
  1412. # units = "cm",
  1413. # device = "svg"
  1414. #)
  1416. ##4.4. Analysis of queen behavior and antennation through the wire mesh -----
  1418. #Import data
  1419. Antennation=read.table("data_queen_antennation_wire_mesh.txt", h=T, stringsAsFactors = T)
  1420. Antennation_subject=subset(as.data.frame(txi_with_without_separated_workers$meta), workers=="separated_workers")
  1421. Antennation_subject=sub("B$","",Antennation_subject$sample)
  1422. Antennation=subset(Antennation, Subject %in% Antennation_subject)
  1424. #Number of occurrence per behavior
  1425. Antennation_counts <- Antennation %>%
  1426. group_by(Subject , Behavior) %>%
  1427. summarise(nb_occurrences = n(), .groups = "drop")
  1430. All_queens <- data.frame(
  1431. Subject = unique(Antennation_subject)
  1432. )
  1434. #Number of each queen behavior
  1435. Antennation_all <- All_queens %>%
  1436. left_join(Antennation_counts, by = "Subject") %>%
  1437. mutate(nb_occurrences = ifelse(is.na(nb_occurrences), 0, nb_occurrences))
  1440. #Plot Number of antennation through the wire mesh for each queen
  1441. Antennation_Occurence <- ggplot(
  1442. Antennation_all,
  1443. aes(
  1444. x = Subject,
  1445. y = nb_occurrences,
  1446. fill = Behavior
  1447. )
  1448. ) +
  1449. geom_col() +
  1450. labs(
  1451. x = "Queen identity",
  1452. y = "Number of occurrences"
  1453. ) +
  1454. theme_classic(base_size = 18) +
  1455. scale_fill_viridis_d(
  1456. option = "E",
  1457. na.translate = FALSE,
  1458. name=NULL,
  1459. labels = c(
  1460. "Number of times that queens made antennal contact with workers" = "Number of times that queens made antennal\ncontact with workers",
  1461. "Number of times that queens passed their antennae through or above the wire mesh without touching a worker" = "Number of times that queens passed their antennae\nthrough or above the wire mesh without touching a worker"
  1462. )
  1463. ) +
  1464. theme(
  1465. plot.title = element_text(size = 20, face = "bold", hjust = 0.5),
  1466. strip.text = element_text(size = 16, face = "bold"),
  1467. panel.border = element_blank(),
  1468. panel.grid.major = element_blank(),
  1469. panel.grid.minor = element_blank(),
  1470. axis.line = element_line(colour = "black"),
  1472. legend.position = "top",
  1473. legend.justification = "center",
  1474. legend.direction = "horizontal"
  1475. )
  1477. print(Antennation_Occurence)
  1480. ###Duration of antennal contacts with workers
  1481. Duration=Antennation
  1482. Duration$Antennation=NA
  1484. #Creation of a new column with the duration of both antennation behaviors with the workers,
  1485. #and 0 for all other behaviors
  1486. Duration$Antennation <- ifelse(Duration$Behavior =="Number of times that queens made antennal contact with workers",
  1487. as.numeric(Duration$Duration),0)
  1489. #Addition of the E4_53 subject that did not perform any behavior
  1490. Duration <- Duration %>%
  1491. dplyr::bind_rows(data.frame(Subject = "E4_53",Antennation = 0)
  1492. )
  1494. Duration$Duration=as.numeric(Duration$Duration)
  1496. #Creation of a dataframe with the total duration of antennation for each queen
  1497. Duration_Both_Antennation <- All_queens %>%
  1498. left_join(
  1499. Duration %>%
  1500. group_by(Subject)%>%
  1501. summarise(
  1502. Total_duration = sum(Duration, na.rm = TRUE),
  1503. .groups = "drop"
  1504. ),
  1505. by = "Subject"
  1506. )
  1509. #Calculation of the mean and standard deviation of the d'antennation
  1510. Duration_mean=subset(Duration_Both_Antennation, Total_duration >0)
  1511. mean(Duration_mean$Total_duration)
  1512. sd(Duration_mean$Total_duration)
  1515. #Antennation duration plot
  1516. Duration_plot <- ggplot(
  1517. Duration_Both_Antennation,
  1518. aes(x = Subject, y = Total_duration, fill = Subject)
  1519. ) +
  1520. geom_col(width = 0.7, fill="#00204DFF") +
  1521. labs(
  1522. x = "Queen identity",
  1523. y = "Total duration of antennal contacts with workers (s)"
  1524. ) +
  1525. theme_classic(base_size = 18) +
  1526. theme(
  1527. axis.text.x = element_text(size = 14),
  1528. plot.title = element_text(color = "#00204DFF", size = 14, face = "bold"),
  1529. panel.border = element_blank(),
  1530. panel.grid.major = element_blank(),
  1531. panel.grid.minor = element_blank(),
  1532. axis.line = element_line(colour = "black"),
  1533. legend.position = "none"
  1534. )
  1536. print(Duration_plot)
  1538. #Creation of a panel with both plots
  1539. Antennation_2plots=ggarrange(Antennation_Occurence,Duration_plot, labels="AUTO")
  1540. Antennation_2plots
  1542. #Export
  1543. #ggsave(
  1544. # filename = "Supplementary_Figure_3.svg",
  1545. # plot = Antennation_2plots,
  1546. # width = 22,
  1547. # height = 10,
  1548. # dpi = 600,
  1549. # bg = "white",
  1550. # device = "svg")
  1552. #Part 5: Comparison of the genes of interest between queens with and without larvae-----
  1554. ##5.1. Load count table and metadata information------
  1556. txi_with_without_larvae <- readRDS("txi_with_without_larvae.rds")
  1558. ##5.2. Analysis of expression differences-----
  1560. #Create dds
  1561. dds_3.3=DESeqDataSetFromMatrix(countData = txi_with_without_larvae$counts,
  1562. colData=txi_with_without_larvae$meta,
  1563. design=~workers)
  1566. #To only keep genes expressed in at least half of the samples, we only keep rows with at least 13 samples with a non-zero value
  1567. dds50_3.3=dds_3.3[rowSums(counts(dds_3.3) >0)>= 13]
  1569. #Normalized counts extraction
  1570. dds_test_3.3=DESeq(dds50_3.3, test="LRT", reduced=~1)
  1571. normalized_counts.3.3 = counts(dds_test_3.3, normalized = TRUE)
  1573. #Pivot longer the dataframe
  1574. normalized_counts_long.3.3 = normalized_counts.3.3 %>%
  1575. as.data.frame() %>%
  1576. rownames_to_column(var = "gene") %>%
  1577. pivot_longer(
  1578. cols = -gene,
  1579. names_to = "sample",
  1580. values_to = "normalized_count"
  1581. )
  1583. normalized_counts_long.3.3 = normalized_counts_long.3.3 %>%
  1584. left_join(as.data.frame(txi_with_without_larvae$meta) %>% dplyr::select(sample, workers), by = "sample")
  1586. ###Kr-h1
  1587. #Subset Kr-h1 counts
  1588. Krup.3.3=subset(normalized_counts_long.3.3, gene=="Lnig_g01448")
  1590. #Model
  1591. modK.3.3=lm(normalized_count~workers, data=Krup.3.3)
  1592. simulationOutput = simulateResiduals(fittedModel= modK.3.3, plot = T)
  1594. #Results
  1595. car::Anova(modK.3.3)
  1596. summary(modK.3.3)
  1599. ###PDK1
  1600. #Subset the dPDK1 counts
  1601. dPDK1.3.3 =subset(normalized_counts_long.3.3, gene=="Lnig_g11740")
  1603. #Model
  1604. moddPDK1.3.3=lm(normalized_count~workers, data=dPDK1.3.3)
  1605. simulationOutput = simulateResiduals(fittedModel= moddPDK1.3.3, plot = T)
  1607. #Results
  1608. car::Anova(moddPDK1.3.3)
  1609. summary(moddPDK1.3.3)
  1611. ###Insulin
  1612. #Subset the Insulin counts
  1613. Ins.3.3 =subset(normalized_counts_long.3.3, gene=="Lnig_g17128")
  1615. #Model
  1616. modIns.3.3=lm(normalized_count~workers,data=Ins.3.3)
  1617. simulationOutput = simulateResiduals(fittedModel= modIns.3.3, plot = T)
  1619. ##Results
  1620. car::Anova(modIns.3.3)
  1621. summary(modIns.3.3)
  1624. ###Juvenile hormone esterase
  1625. #Subset the Juvenile hormone esterase counts
  1626. Carb.3.3 =subset(normalized_counts_long.3.3, gene=="Lnig_g08195")
  1628. #Model
  1629. modCarb.3.3=lm(normalized_count~workers, data=Carb.3.3)
  1630. simulationOutput = simulateResiduals(fittedModel= modCarb.3.3, plot = T)
  1632. #Results
  1633. car::Anova(modCarb.3.3)
  1634. summary(modCarb.3.3)
  1636. ###Esterase E4
  1637. #Subset the Esterase E4 counts
  1638. E4.3.3 =subset(normalized_counts_long.3.3, gene=="Lnig_g15862")
  1640. #Model
  1641. modE4.3.3=lm(normalized_count~workers, data=E4.3.3)
  1642. simulationOutput = simulateResiduals(fittedModel= modE4.3.3, plot = T)
  1644. #Results
  1645. car::Anova(modE4.3.3)
  1646. summary(modE4.3.3)
  1649. ###Doublesex
  1650. #Subset the Doublesex counts
  1651. Dbs.3.3 =subset(normalized_counts_long.3.3, gene=="Lnig_g04508")
  1653. #Model
  1654. modDbs.3.3=lm(normalized_count~workers, data=Dbs.3.3)
  1655. simulationOutput = simulateResiduals(fittedModel= modDbs.3.3, plot = T)
  1657. #Results
  1658. car::Anova(modDbs.3.3)
  1659. summary(modDbs.3.3)
  1662. ###Protein henna
  1663. #Subset the Protein henna counts
  1664. Henna.3.3=subset(normalized_counts_long.3.3, gene=="Lnig_g13064")
  1666. #Model
  1667. modHenna.3.3=lm(normalized_count~workers, data=Henna.3.3)
  1668. simulationOutput = simulateResiduals(fittedModel= modHenna.3.3, plot = T)
  1670. #Results
  1671. car::Anova(modHenna.3.3)
  1672. summary(modHenna.3.3)
  1675. ###Prohormone-2
  1676. #Subset the Prohormone-2 counts
  1677. Pro.3.3 =subset(normalized_counts_long.3.3, gene=="Lnig_g06899")
  1679. #Model
  1680. modPro.3.3=lm(normalized_count~workers, data=Pro.3.3)
  1681. simulationOutput = simulateResiduals(fittedModel= modPro.3.3, plot = T)
  1683. #Results
  1684. car::Anova(modPro.3.3)
  1685. summary(modPro.3.3)
  1687. ###Pheromone binding protein GP-9
  1688. #Subset the GP9 counts
  1689. GP9.3.3=subset(normalized_counts_long.3.3, gene=="gLnigOBP1")
  1691. #Model
  1692. modGP9.3.3=lm(normalized_count~workers, data=GP9.3.3)
  1693. simulationOutput = simulateResiduals(fittedModel= modGP9.3.3, plot = T)
  1695. #Results
  1696. car::Anova(modGP9.3.3)
  1697. summary(modGP9.3.3)
  1700. ### GR81
  1701. #Subset the GR81 counts
  1702. GR81.3.3=subset(normalized_counts_long.3.3, gene=="gLnigGR81")
  1704. #Model
  1705. modGR81.3.3=lm(normalized_count~workers, data=GR81.3.3)
  1706. simulationOutput = simulateResiduals(fittedModel= modGR81.3.3, plot = T)
  1708. #Results
  1709. car::Anova(modGR81.3.3)
  1710. summary(modGR81.3.3)
  1713. ###Tret1a
  1714. #Subset the Tret1a counts
  1715. Tret1a.3.3 =subset(normalized_counts_long.3.3, gene=="Lnig_g01559")
  1717. #Model
  1718. modTret1a.3.3=lm(normalized_count~workers,data=Tret1a.3.3)
  1719. simulationOutput = simulateResiduals(fittedModel= modTret1a.3.3, plot = T)
  1721. #Results
  1722. car::Anova(modTret1a.3.3)
  1723. summary(modTret1a.3.3)
  1726. ###Tret1b
  1727. #Subset the Tret1b counts
  1728. Tret1b.3.3 =subset(normalized_counts_long.3.3, gene=="Lnig_g11948")
  1730. #Model
  1731. modTret1b.3.3=lm(normalized_count~workers, data=Tret1b.3.3)
  1732. simulationOutput = simulateResiduals(fittedModel= modTret1b.3.3, plot = T)
  1734. #Results
  1735. car::Anova(modTret1b.3.3)
  1736. summary(modTret1b.3.3)
  1739. ###Tret1-2 homolog
  1740. #Subset the Tret1-2 counts
  1741. Tret1_2.3.3=subset(normalized_counts_long.3.3, gene=="Lnig_g10144")
  1743. #Model
  1744. modTret1_2.3.3=lm(normalized_count~workers, data=Tret1_2.3.3)
  1745. simulationOutput = simulateResiduals(fittedModel= modTret1_2.3.3, plot = T)
  1747. #Results
  1748. car::Anova(modTret1_2.3.3)
  1749. summary(modTret1_2.3.3)
  1751. ##5.3. Produce Supplementary Figure 5-----
  1753. #Manual title colors
  1754. title_pat <- gsub(" ", "", "bbb bbb bbb bbb b")
  1755. title_vec <- strsplit(title_pat, "")[[1]]
  1756. title_col <- ifelse(title_vec == "r", "darkgrey", "black")
  1757. title_col <- rep(title_col, length.out = 13)
  1758. cb_palette_5 <- c(without_larvae = "#46D3E9" ,with_larvae = "#D55E00")
  1760. models_5 <- list(
  1761. Krh1 = modK.3.3,
  1762. dPDK1 = moddPDK1.3.3,
  1763. Insulin = modIns.3.3,
  1764. JHE = modCarb.3.3,
  1765. E4 = modE4.3.3,
  1766. Dsx = modDbs.3.3,
  1767. GP9 = modGP9.3.3,
  1768. GR81 = modGR81.3.3,
  1769. Tret1a = modTret1a.3.3,
  1770. Tret1b = modTret1b.3.3,
  1771. Tret1_2 = modTret1_2.3.3,
  1772. Henna = modHenna.3.3,
  1773. Prohormone = modPro.3.3
  1774. )
  1776. padj_3.3 <- list(
  1777. Krh1 = car::Anova(modK.3.3) [["Pr(>F)"]][1],
  1778. dPDK1 = car::Anova(moddPDK1.3.3)[["Pr(>F)"]][1],
  1779. Insulin = car::Anova(modIns.3.3)[["Pr(>F)"]][1],
  1780. JuvEsterase = car::Anova(modCarb.3.3)[["Pr(>F)"]][1],
  1781. E4 = car::Anova(modE4.3.3)[["Pr(>F)"]][1],
  1782. Dsx = car::Anova(modDbs.3.3)[["Pr(>F)"]][1],
  1783. GP9 = car::Anova(modGP9.3.3)[["Pr(>F)"]][1],
  1784. GR81 = car::Anova(modGR81.3.3)[["Pr(>F)"]][1],
  1785. Tret1a = car::Anova(modTret1a.3.3)[["Pr(>F)"]][1],
  1786. Tret1b = car::Anova(modTret1b.3.3)[["Pr(>F)"]][1],
  1787. Tret1_2 = car::Anova(modTret1_2.3.3)[["Pr(>F)"]][1],
  1788. Henna = car::Anova(modHenna.3.3)[["Pr(>F)"]][1],
  1789. Prohormone = car::Anova(modPro.3.3)[["Pr(>F)"]][1]
  1790. )
  1792. padj_vec_5 <- setNames(padj_3.3, names(padj_3.3))
  1793. padj_vec_5 <- unlist(padj_vec_5, use.names = TRUE)
  1796. all_plots_5 <- lapply(seq_len(nrow(plot_info)), function(i) {
  1797. make_q5_panel(plot_info[i, ], i, ncol_grid = 3)
  1798. })
  1801. final_figure_5 <- plot_grid(
  1802. plotlist = all_plots_5,
  1803. nrow = 5,
  1804. ncol = 3,
  1805. labels = LETTERS[1:13],
  1806. label_size = 20,
  1807. label_fontface = "bold",
  1808. align = "hv"
  1809. )
  1811. print(final_figure_5)
  1813. #ggsave(
  1814. # filename = "Supplementary_Figure_5.svg",
  1815. # plot = final_figure_5,
  1816. # width = 50,
  1817. # height = 75,
  1818. # units = "cm",
  1819. # device = "svg"
  1820. #)
  1822. #Production of Figure 3-----
  1824. #Addition of y axis label to plot D
  1825. all_plots[[3]] <- all_plots[[3]] +
  1826. ggplot2::ylab("Normalized read counts") +
  1827. ggplot2::theme(axis.title.y = ggplot2::element_text(angle = 90,size = 20))
  1830. figure3_plots <- c(
  1831. all_plots[1], all_plots_4[1],
  1832. all_plots_5[1],all_plots[3], all_plots_4[3],
  1833. all_plots_5[3],all_plots[4], all_plots_4[4],
  1834. all_plots_5[4],all_plots[13], all_plots_4[13],
  1835. all_plots_5[13]
  1836. )
  1840. keep_y <- seq(1, length(figure3_plots), by = 3)
  1842. figure3_plots_clean <- lapply(seq_along(figure3_plots), function(i) {
  1843. if (!(i %in% keep_y)) {
  1844. figure3_plots[[i]] +
  1845. ggplot2::theme(axis.title.y = ggplot2::element_blank())
  1846. } else {
  1847. figure3_plots[[i]]
  1848. }
  1849. })
  1852. figure3_grid <- cowplot::plot_grid(
  1853. plotlist = figure3_plots_clean,
  1854. nrow = 5,
  1855. ncol = 3,
  1856. labels = LETTERS,
  1857. label_size = 20,
  1858. label_fontface = "bold",
  1859. align = "hv"
  1860. )
  1862. col_titles <- plot_grid(
  1863. ggdraw() + draw_label("Experiment A", fontface = "bold", size = 30, x = 0.55),
  1864. ggdraw() + draw_label("Experiment B", fontface = "bold", size = 30, x = 0.55),
  1865. ggdraw() + draw_label("Experiment C", fontface = "bold", size = 30, x = 0.55),
  1866. ncol = 3
  1867. )
  1870. figure3 <- plot_grid(
  1871. col_titles,
  1872. figure3_grid,
  1873. ncol = 1,
  1874. rel_heights = c(0.05, 0.95)
  1875. )
  1878. figure3
  1881. #ggsave(
  1882. # filename = "Figure_3.svg",
  1883. # plot = figure3,
  1884. # width = 50,
  1885. # height = 75,
  1886. #units = "cm",
  1887. # device = "svg"
  1888. #)

R script_Mol_Ecol_2026.R, no license · at the source

Overview

Authors: Hugo E Le Lay1, Maximilian F Bolder2, Thomas J Colgan2, Romain Libbrecht1,2
ORCID iDs: Hugo E Le Lay, Maximilian F Bolder, Thomas J Colgan, Romain Libbrecht
  1. Insect Biology Research Institute, UMR 7261, CNRS, University of Tours, Tours, France
  2. Institute of Organismic and Molecular Evolution, Johannes Gutenberg University, Mainz, Germany
Institutions: Université de Tours (France); Centre National de la Recherche Scientifique (France); Institut de Recherche sur la Biologie de l'Insecte UMR 7261 (France); Johannes Gutenberg University Mainz (Germany)
Journal: Molecular ecology, volume 35, issue 21, article e70614
Dates: received 7 May 2026; accepted 15 October 2026; published online 1 November 2026; in print November 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1111/mec.70614 · PMID 42968400 · PMCID PMC13724878 · OpenAlex W7225107947
Open access: hybrid, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: genetics / omics (modality), other (organism), cellular / molecular (subfield)
Methods: Statistics, Preprocessing
Keywords: colony foundation, division of labour, gene expression, social behaviour, social insects
MeSH: Ants*, Behavior, Animal*, Brain*, Social Behavior*, Transcriptome*, Animals, Female, Gene Expression Profiling, Juvenile Hormones (* major topic)
Topic: Insect and Arachnid Ecology and Behavior (Insect Science, Agricultural and Biological Sciences), according to OpenAlex
Funding: Deutsche Forschungsgemeinschaft (GRK2526/1-407023052); Agence Nationale de la Recherche (ANR-23-CE02-0004)
Citations: not cited yet (Europe PMC); 130 references in the paper

Abstract

Insect societies are characterised by a reproductive division of labour between queens that monopolise egg production and workers that perform all non‐reproductive tasks. The queen specialisation in egg laying is typically considered to be constitutive and independent of environmental conditions, but this common assumption was recently challenged by the finding that the presence of workers both induces and maintains the queen specialisation in several species of ants. However, the molecular changes that translate worker presence into queen behavioural specialisation remain unknown. To investigate the transcriptomic basis of the socially regulated specialisation of ant queens, we produced 93 individual brain transcriptomes from queens of the black garden ant Lasius niger that were experimentally maintained in different social environments. Our transcriptomic analyses revealed that 123 genes were differentially expressed between the brains of queens maintained with and without workers. These genes were enriched for processes related to brain energy metabolism and neuronal plasticity and signalling. By investigating the response of selected genes of interest to further manipulations of the social environment, we identified four candidate genes that fully mirrored the queen behavioural changes. Specifically, their expression changes in response to worker presence were reversible and they did not respond to mere cues of the worker presence. These candidate genes included the transcription factor Krüppel‐homologue 1 and other genes in the juvenile hormone and insulin signalling pathways, both known to regulate multiple aspects of division of labour in social insects. By producing a robust, well‐replicated and well‐controlled investigation of the transcriptomic basis of the behavioural specialisation of ant queens, this study provides a better understanding of the emergence and maintenance of reproductive division of labour in insect societies.

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

supp:PMC13724878/MEC-35-e70614-s003.zip

License: none: the authors keep all their rights
State: the link answers, verified on 3 November 2026
Evidence: files inventoried
Languages: R (1)
Size: 1 file, 1 script
Software Heritage: not checked
Found in: the supplementary material
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: car (1 file), cowplot (1 file), DESeq2 (1 file), emmeans (1 file), ggplot2 (1 file), ggpubr (1 file), lme4 (1 file), multcomp (1 file), patchwork (1 file), tidyverse (1 file)
Availability: 1 check, the latest on 3 November 2026: the link answers (HTTP 200)
  • 3 November 2026: the link answers (HTTP 200)
1 file
At the source: pmc-oa-opendata.s3.amazonaws.com/PMC13724878.1/MEC-35-e70614-s003.zip

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

Tracing map

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

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 1 script, each with its path and the digest of its content;
  • 9 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 code and input data for all analyses are provided in Data S1. The raw RNA‐seq sequences used in this study are publicly available from the NCBI BioProject PRJNA1548310.

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, 11 November 2026

  • Publisher: n/a → Wiley

Version 1, 3 November 2026: the first record

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

Cite

This paper

Le Lay, H. E., Bolder, M. F., Colgan, T. J., & Libbrecht, R. (2026). Brain Transcriptomic Changes Associated With the Socially Regulated Behavioural Specialisation of Ant Queens. Molecular ecology, 35(21), e70614. https://doi.org/10.1111/mec.70614

BibTeX

@article{lelay2026brain,
author = {Le Lay, Hugo E and Bolder, Maximilian F and Colgan, Thomas J and Libbrecht, Romain},
title = {{Brain Transcriptomic Changes Associated With the Socially Regulated Behavioural Specialisation of Ant Queens}},
journal = {Molecular ecology},
year = {2026},
month = nov,
volume = {35},
number = {21},
pages = {e70614},
publisher = {Wiley},
issn = {0962-1083},
doi = {10.1111/mec.70614},
url = {https://doi.org/10.1111/mec.70614},
pmid = {42968400},
pmcid = {PMC13724878}
}

RIS

TY - JOUR
AU - Le Lay, Hugo E
AU - Bolder, Maximilian F
AU - Colgan, Thomas J
AU - Libbrecht, Romain
TI - Brain Transcriptomic Changes Associated With the Socially Regulated Behavioural Specialisation of Ant Queens
T2 - Molecular ecology
J2 - Mol Ecol
PY - 2026
DA - 2026/11/01
VL - 35
IS - 21
SP - e70614
SN - 0962-1083
PB - Wiley
DO - 10.1111/mec.70614
UR - https://doi.org/10.1111/mec.70614
LA - en
ER -

CSL-JSON

{
"id": "10.1111/mec.70614",
"type": "article-journal",
"title": "Brain Transcriptomic Changes Associated With the Socially Regulated Behavioural Specialisation of Ant Queens",
"container-title": "Molecular ecology",
"author": [
{
"family": "Le Lay",
"given": "Hugo E"
},
{
"family": "Bolder",
"given": "Maximilian F"
},
{
"family": "Colgan",
"given": "Thomas J"
},
{
"family": "Libbrecht",
"given": "Romain"
}
],
"container-title-short": "Mol Ecol",
"volume": "35",
"issue": "21",
"page": "e70614",
"DOI": "10.1111/mec.70614",
"PMID": "42968400",
"PMCID": "PMC13724878",
"ISSN": "0962-1083",
"publisher": "Wiley",
"URL": "https://doi.org/10.1111/mec.70614",
"language": "en",
"issued": {
"date-parts": [
[
2026,
11,
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.1111/mec.70360 [code]
Unmated Queens Show Worker-Like Behaviour and Gene Expression in Polygynous Colonies of the Ant Stigmatomma pallipes.
Journal: Molecular ecology
In common: other, genetics / omics, cellular / molecular, 16 references, 2 authors
[2] doi:10.1038/s42003-023-05273-4 [code]
Task-specific odorant receptor expression in worker antennae indicates that sensory filters regulate division of labor in ants.
Journal: Communications biology
In common: other, genetics / omics, cellular / molecular, 16 references, author Romain Libbrecht
[3] doi:10.1371/journal.pbio.3002763 [code]
An evolutionarily conserved pathway mediated by neuroparsin-A regulates reproductive plasticity in ants.
Journal: PLoS biology
In common: other, cellular / molecular, 12 references
[4] doi:10.1242/jeb.246785
Effects of microRNA-305 knockdown on brain gene expression associated with division of labor in honey bee colonies (Apis mellifera).
Journal: The Journal of experimental biology
In common: other, genetics / omics, cellular / molecular, 9 references
[5] doi:10.1038/s41586-026-10747-6 [code]
Ancient feeding-related neuropeptides regulate alloparenting in ants.
Journal: Nature
In common: patchwork, ggpubr, ggplot2, 1 other tool, other, 7 references
[6] doi:10.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: multcomp, emmeans, car, 7 other tools, genetics / omics, 2 references
[7] doi:10.1038/s43587-025-00816-2 [code]
Interleukin-12 signaling drives Alzheimer's disease pathology through disrupting neuronal and oligodendrocyte homeostasis.
Journal: Nature aging
In common: multcomp, car, DESeq2, 6 other tools, cellular / molecular, 1 reference
[8] doi:10.1111/gbb.70007 [code]
The brain atlas of a subsocial bee reflects that of eusocial Hymenoptera.
Journal: Genes, brain, and behavior
In common: patchwork, cowplot, ggplot2, 1 other tool, other, genetics / omics, 4 references
[9] doi:10.1371/journal.pone.0353990 [code]
Positive mood enhances accessibility of unrelated concepts in the first language but not in the foreign language.
Journal: PloS one
In common: multcomp, emmeans, car, 5 other tools, 1 reference
[10] doi:10.1002/aur.70042 [code]
Multisensory Integration of Naturalistic Speech and Gestures in Autistic Adults.
Journal: Autism research : official journal of the International Society for Autism Research
In common: multcomp, emmeans, car, 6 other tools

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.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1111/mec.70614.svg)](https://openscicode.org/badge/doi/10.1111/mec.70614)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1111/mec.70614.svg
   :target: https://openscicode.org/badge/doi/10.1111/mec.70614
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1111/mec.70614"><img src="https://openscicode.org/badge/doi/10.1111/mec.70614.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge 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: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

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