OSCR

Transcriptional predictors of rescue behaviour in ants.

Code ↔ Paper

10 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 10 matches
  1. [1] § RESULTS › Gene expression differences per tissue ↔ Script_Rescue.R, lines 568–614 · score 0.97 · gene_13149, gene_1607, gene_3426, gene_3661, gene_3662, gene_5298
  2. [2] § RESULTS › Gene expression differences across tissues ↔ Script_Rescue.R, lines 662–731 · score 0.90 · gene_13318, gene_13632, gene_14429, gene_1997, gene_7603, gene_14458
  3. [3] § MATERIALS AND METHODS › RNA extraction and sequencing ↔ Script_Rescue.R, lines 618–660 · score 0.84 · antennal lobes, optic lobes, mushroom bodies, brain tissue, RNA, dissect
  4. [4] § MATERIALS AND METHODS › Analysis of gene expression ↔ Script_Rescue.R, lines 662–731 · score 0.77 · variance stabilizing transformation, reduced model, Colony ID, DESeq2, variable, PCA
  5. [5] § MATERIALS AND METHODS › RNA extraction and sequencing ↔ Script_Rescue.R, lines 618–660 · score 0.74 · mRNA, DESeq2, brain tissues, sequencing, antennae, CC
  6. [6] § RESULTS › Non-annotated DEGs and GO-term enrichment ↔ Script_Rescue.R, lines 1074–1122 · score 0.72 · cellular components, GO enrichment, GO terms, molecular functions, biological, upregulated
  7. [7] § RESULTS › Non-annotated DEGs and GO-term enrichment ↔ Script_Rescue.R, lines 1074–1122 · score 0.72 · InterPro, GO term, gene_1243, gene_9900, Pfam, protein
  8. [8] § RESULTS › Gene expression differences per tissue ↔ Script_Rescue.R, lines 847–906 · score 0.65 · Venn diagram, variance stabilized transformed, log2 fold change, PCA, LRT, tissues
  9. [9] § MATERIALS AND METHODS › Behavioural experiment ↔ Script_Rescue.R, lines 500–566 · score 0.64 · travelled distance, simulateResiduals, outer circle, glmmTMB, latency, anova
  10. [10] § MATERIALS AND METHODS › Sensitivity analysis ↔ Script_Rescue.R, lines 1–72 · score 0.57 · individual ID, sensitivity, interaction, glmmTMB, predictors, rescue behaviour

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 1,169 lines · 54 KB · no license · 10 matches

  1. ### Transcriptional predictors of rescue behaviour in ants
  3. ### Jaimes-Nino, Luisa Maria1, Bar, Adi2, Scharf, Inon2, Foitzik, Susanne1
  5. ### 1 Institute of Organismic and Molecular Evolution, Johannes Gutenberg University Mainz, Germany
  6. ### 2 School of Zoology, George S Wise Faculty of Life Sciences, Tel Aviv University, Tel Aviv, Israel
  9. #### Brain size rescuers vs non-rescuers ####
  10. library(glmmTMB)
  11. library(emmeans)
  12. library(DHARMa)
  13. library(openxlsx)
  14. library(ggplot2)
  15. library(ggbeeswarm)
  16. library(patchwork)
  18. size_rescue<- read.xlsx("~/Documents/04_Catag_rescue/MS_rescue/New_subm/JExpBio/Revision_JExpB_rescue/Supl_material_File2.xlsx", sheet = 2)
  21. ggplot(data=size_rescue, aes(x=LengthB, y =WidthB, color = Rescue))+
  22. geom_point()+
  23. geom_smooth(method=lm)+
  24. ylab("Width (mm)")+
  25. xlab("Length (mm)")+
  26. scale_color_manual(values=c("blue3", "red3"))+
  27. theme_minimal()
  28. sizeL.model <- glmmTMB(data = size_rescue, LengthB ~ Rescue + (1|Individual_id))
  29. sizeL.model.wo <- glmmTMB(data = size_rescue, LengthB ~ 1 + (1|Individual_id))
  30. summary(sizeL.model)
  31. car::Anova(sizeL.model, type = "III")
  32. anova(sizeL.model, sizeL.model.wo )
  33. ##### Sensitivity analysis#####
  34. HR_track<- read.xlsx("~/Documents/04_Catag_rescue/MS_rescue/New_subm/JExpBio/Revision_JExpB_rescue/Supl_material_File2.xlsx", sheet = 6)
  36. HR_track$Average_Speed <- as.numeric(HR_track$Average_Speed)
  37. HR_track$Average_Speed_i<- as.numeric(HR_track$Average_Speed_i)
  38. HR_track$Average_Speed_o<- as.numeric(HR_track$Average_Speed_o)
  39. table(HR_track$Resolution, HR_track$Rescue)
  41. HR_track$ID <- paste(HR_track$Video, HR_track$Individual, sep = "_")
  43. speed_o_lm <- glmmTMB(data= HR_track, log(Average_Speed_o) ~ Resolution * Rescue + (1|ID))
  45. speed_o_lm_woint <- glmmTMB(data= HR_track, log(Average_Speed_o) ~ Resolution + Rescue + (1|ID))
  47. anova(speed_o_lm_woint, speed_o_lm) #Testing interaction term
  48. # Data: HR_track
  49. # Models:
  50. # speed_o_lm_woint: log(Average_Speed_o) ~ Resolution + Rescue + (1 | ID), zi=~0, disp=~1
  51. # speed_o_lm: log(Average_Speed_o) ~ Resolution * Rescue + (1 | ID), zi=~0, disp=~1
  52. # Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
  53. # speed_o_lm_woint 5 70.560 79.481 -30.280 60.560
  54. # speed_o_lm 6 70.987 81.692 -29.494 58.987 1.5733 1 0.2097
  55. speed_o_lm_woRes <- glmmTMB(data= HR_track, log(Average_Speed_o) ~ Resolution + (1|ID))
  56. anova(speed_o_lm_woint, speed_o_lm_woRes) # Testing Rescue
  57. # Data: HR_track
  58. # Models:
  59. # speed_o_lm_woRes: log(Average_Speed_o) ~ Resolution + (1 | ID), zi=~0, disp=~1
  60. # speed_o_lm_woint: log(Average_Speed_o) ~ Resolution + Rescue + (1 | ID), zi=~0, disp=~1
  61. # Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
  62. # speed_o_lm_woRes 4 79.903 87.040 -35.951 71.903
  63. # speed_o_lm_woint 5 70.560 79.481 -30.280 60.560 11.343 1 0.0007575 ***
  64. speed_o_lm_woReso <- glmmTMB(data= HR_track, log(Average_Speed_o) ~ Rescue + (1|ID))
  65. anova(speed_o_lm_woint, speed_o_lm_woReso) # Testing Resolution
  66. # Data: HR_track
  67. # Models:
  68. # speed_o_lm_woReso: log(Average_Speed_o) ~ Rescue + (1 | ID), zi=~0, disp=~1
  69. # speed_o_lm_woint: log(Average_Speed_o) ~ Resolution + Rescue + (1 | ID), zi=~0, disp=~1
  70. # Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
  71. # speed_o_lm_woReso 4 72.896 80.033 -32.448 64.896
  72. # speed_o_lm_woint 5 70.560 79.481 -30.280 60.560 4.336 1 0.03731 *
  74. # Estimated marginal means
  75. emm.resc <- emmeans(speed_o_lm_woint, ~ Resolution , type = "response")
  76. # Resolution response SE df asymp.LCL asymp.UCL
  77. # HR 1.163 0.1670 Inf 0.878 1.541
  78. # LR 0.816 0.0711 Inf 0.688 0.968
  80. # Resolution Rescue response SE df asymp.LCL asymp.UCL
  81. # HR No 1.534 0.2520 Inf 1.111 2.116
  82. # LR No 1.076 0.1370 Inf 0.839 1.380
  83. # HR Yes 0.882 0.1410 Inf 0.644 1.207
  84. # LR Yes 0.619 0.0634 Inf 0.506 0.756
  86. simulateResiduals(speed_o_lm_woint, plot = T)
  88. summary(speed_o_lm_woint)
  92. HR_track$xpos <- interaction(HR_track$Resolution, HR_track$Rescue,
  93. sep = "_")
  94. HR_track$xpos <- factor(HR_track$xpos , levels = c("HR_No", "HR_Yes", "LR_No", "LR_Yes"))
  96. Resol_outer <- ggplot(HR_track, aes(x = xpos, y = log(Average_Speed_o), colour = Rescue)) +
  97. geom_beeswarm() +
  98. theme_minimal()+
  99. stat_summary(fun.y= median, fun.ymin=median, fun.ymax=median, geom="crossbar", width=0.3, linewidth =0.3)+
  100. #annotate("text", color = "black",label = "n.s", x = 1.45, y = 1.65) +
  101. scale_color_manual(values=c("blue3", "red3"))+
  102. theme_classic()+
  103. ylab("log (Average speed outer circle)")+
  104. xlab("Video resolution and Rescue")+
  105. theme(legend.position="none",
  106. legend.text=element_text(size=12),
  107. axis.title = element_text(size = 12),
  108. axis.text = element_text(size = 12))
  111. Resol_in <- ggplot(HR_track, aes(x = xpos, y = log(Average_Speed_i), colour = Rescue)) +
  112. geom_beeswarm() +
  113. scale_color_manual(values = c("blue3", "red3")) +
  114. theme_minimal()+
  115. stat_summary(fun.y= median, fun.ymin=median, fun.ymax=median, geom="crossbar", width=0.3, linewidth =0.3)+
  116. #annotate("text", color = "black",label = "n.s", x = 1.45, y = 1.65) +
  117. scale_color_manual(values=c("blue3", "red3"))+
  118. theme_classic()+
  119. ylab("log (Average speed inner circle)")+
  120. xlab("Video resolution and Rescue")+
  121. theme(legend.position="none",
  122. legend.text=element_text(size=12),
  123. axis.title = element_text(size = 12),
  124. axis.text = element_text(size = 12))
  128. speed_i_lm <- glmmTMB(data= HR_track, log(Average_Speed_i) ~ Resolution * Rescue + (1|ID))
  129. speed_i_lm_wo <- glmmTMB(data= HR_track, log(Average_Speed_i) ~ Resolution + Rescue + (1|ID))
  130. anova(speed_i_lm_wo, speed_i_lm)
  131. # Data: HR_track
  132. # Models:
  133. # speed_i_lm_wo: log(Average_Speed_i) ~ Resolution + Rescue + (1 | ID), zi=~0, disp=~1
  134. # speed_i_lm: log(Average_Speed_i) ~ Resolution * Rescue + (1 | ID), zi=~0, disp=~1
  135. # Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
  136. # speed_i_lm_wo 5 76.712 84.767 -33.356 66.712
  137. # speed_i_lm 6 77.810 87.475 -32.905 65.810 0.9023 1 0.3422
  138. simulateResiduals(speed_i_lm_wo, plot = T)
  139. speed_i_lm_wo_Resc <- glmmTMB(data= HR_track, log(Average_Speed_i) ~ Resolution + (1|ID))
  140. anova(speed_i_lm_wo, speed_i_lm_wo_Resc)
  141. # Data: HR_track
  142. # Models:
  143. # speed_i_lm_wo_Resc: log(Average_Speed_i) ~ Resolution + (1 | ID), zi=~0, disp=~1
  144. # speed_i_lm_wo: log(Average_Speed_i) ~ Resolution + Rescue + (1 | ID), zi=~0, disp=~1
  145. # Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
  146. # speed_i_lm_wo_Resc 4 80.911 87.355 -36.455 72.911
  147. # speed_i_lm_wo 5 76.712 84.767 -33.356 66.712 6.1989 1 0.01278 *
  148. speed_i_lm_wo_Resol <- glmmTMB(data= HR_track, log(Average_Speed_i) ~ Rescue + (1|ID))
  149. anova(speed_i_lm_wo, speed_i_lm_wo_Resol)
  150. # Data: HR_track
  151. # Models:
  152. # speed_i_lm_wo_Resol: log(Average_Speed_i) ~ Rescue + (1 | ID), zi=~0, disp=~1
  153. # speed_i_lm_wo: log(Average_Speed_i) ~ Resolution + Rescue + (1 | ID), zi=~0, disp=~1
  154. # Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
  155. # speed_i_lm_wo_Resol 4 81.222 87.666 -36.611 73.222
  156. # speed_i_lm_wo 5 76.712 84.767 -33.356 66.712 6.5106 1 0.01072 *
  157. summary(speed_i_lm)
  159. # Estimated marginal means
  160. emm.resc <- emmeans(speed_i_lm, ~ Resolution * Rescue, type = "response")
  161. # Resolution Rescue response SE df asymp.LCL asymp.UCL
  162. # HR No 1.018 0.2840 Inf 0.589 1.759
  163. # LR No 0.735 0.1760 Inf 0.460 1.177
  164. # HR Yes 0.618 0.1270 Inf 0.413 0.924
  165. # LR Yes 0.334 0.0492 Inf 0.250 0.446
  167. # Average speed
  168. HR_track$Rescue <- factor(HR_track$Rescue, levels=c("Yes", "No"))
  169. speed_lm <- glmmTMB(data= HR_track, log(Average_Speed) ~ Resolution * Rescue + (1|ID))
  170. speed_lm_wo <- glmmTMB(data= HR_track, log(Average_Speed) ~ Resolution + Rescue + (1|ID))
  171. anova(speed_lm_wo,speed_lm)
  172. # Data: HR_track
  173. # Models:
  174. # speed_lm_wo: log(Average_Speed) ~ Resolution + Rescue + (1 | ID), zi=~0, disp=~1
  175. # speed_lm: log(Average_Speed) ~ Resolution * Rescue + (1 | ID), zi=~0, disp=~1
  176. # Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
  177. # speed_lm_wo 5 36.722 45.972 -13.361 26.722
  178. # speed_lm 6 37.221 48.322 -12.611 25.221 1.5005 1 0.2206
  179. speed_lm_wo_Resc <- glmmTMB(data= HR_track, log(Average_Speed) ~ Resolution + (1|ID))
  180. anova(speed_lm_wo, speed_lm_wo_Resc) # Testing Rescue
  181. # Data: HR_track
  182. # Models:
  183. # speed_lm_wo_Resc: log(Average_Speed) ~ Resolution + (1 | ID), zi=~0, disp=~1
  184. # speed_lm_wo: log(Average_Speed) ~ Resolution + Rescue + (1 | ID), zi=~0, disp=~1
  185. # Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
  186. # speed_lm_wo_Resc 4 34.750 42.150 -13.375 26.750
  187. # speed_lm_wo 5 36.722 45.972 -13.361 26.722 0.0279 1 0.8672
  188. speed_lm_wo_Resol <- glmmTMB(data= HR_track, log(Average_Speed) ~ Rescue + (1|ID))
  189. anova(speed_lm_wo, speed_lm_wo_Resol) # Testing Resolution
  190. # Data: HR_track
  191. # Models:
  192. # speed_lm_wo_Resol: log(Average_Speed) ~ Rescue + (1 | ID), zi=~0, disp=~1
  193. # speed_lm_wo: log(Average_Speed) ~ Resolution + Rescue + (1 | ID), zi=~0, disp=~1
  194. # Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
  195. # speed_lm_wo_Resol 4 49.707 57.108 -20.854 41.707
  196. # speed_lm_wo 5 36.722 45.972 -13.361 26.722 14.986 1 0.0001083 ***
  197. emm.resc <- emmeans(speed_lm_wo, ~ Resolution , type = "response")
  198. # Resolution response SE df asymp.LCL asymp.UCL
  199. # HR 1.20 0.0934 Inf 1.030 1.398
  200. # LR 0.88 0.0563 Inf 0.777 0.998
  201. Resol_general <- ggplot(HR_track, aes(x = xpos, y = log(Average_Speed), colour = Rescue)) +
  202. geom_beeswarm() +
  203. scale_color_manual(values = c("blue3", "red3")) +
  204. theme_minimal()+
  205. stat_summary(fun.y= median, fun.ymin=median, fun.ymax=median, geom="crossbar", width=0.3, linewidth =0.3)+
  206. #annotate("text", color = "black",label = "n.s", x = 1.45, y = 1.65) +
  207. scale_color_manual(values=c("blue3", "red3"))+
  208. theme_classic()+
  209. ylab("log (Average speed)")+
  210. xlab("Video resolution and Rescue")+
  211. theme(legend.position="none",
  212. legend.text=element_text(size=12),
  213. axis.title = element_text(size = 12),
  214. axis.text = element_text(size = 12))
  216. simulateResiduals(speed_lm_wo, plot = T)
  217. summary(speed_lm_wo)
  218. Resol_general + Resol_outer + Resol_in + plot_annotation(tag_levels = 'A')
  219. library(sjPlot)
  220. library(sjmisc)
  221. plot_model(speed_o_lm, type = "pred", terms = c("Resolution", "Rescue"))
  224. emm.resc <- emmeans(speed_lm_wo, ~ Resolution * Rescue, type = "response")
  225. # Resolution Rescue response SE df asymp.LCL asymp.UCL
  226. # HR No 1.266 0.1290 Inf 1.036 1.55
  227. # LR No 0.868 0.0835 Inf 0.719 1.05
  228. # HR Yes 1.297 0.1220 Inf 1.079 1.56
  229. # LR Yes 0.890 0.0764 Inf 0.752 1.05
  230. pairs(emm.resc, by = "Resolution")
  231. #### General tracking analysis ####
  232. library(patchwork)
  233. library(ggplot2)
  234. library(ggbeeswarm)
  235. library(dplyr)
  237. general_track<- read.xlsx("~/Documents/04_Catag_rescue/MS_rescue/New_subm/JExpBio/Revision_JExpB_rescue/Supl_material_File_2.xlsx", sheet = 4)
  239. gen_rescuers <- subset(general_track, Rescue == "Yes")
  240. gen_nonrescuers <- subset(general_track, Rescue == "No")
  241. gen_rescue_track <- rbind(gen_rescuers, gen_nonrescuers)
  242. gen_rescue_track$Average_Speed <- as.numeric(gen_rescue_track$Average_Speed)
  244. gen_plotA <- ggplot(data=gen_rescue_track, aes(x=Rescue, y=Average_Speed, color = Rescue))+
  245. geom_beeswarm(priority = "none")+
  246. stat_summary(fun.y= mean, fun.ymin=mean, fun.ymax=mean, geom="crossbar", width=0.3, linewidth =0.3)+
  247. annotate("text", color = "black",label = "n.s", x = 1.45, y = 1.65) +
  248. scale_color_manual(values=c("blue3", "red3"))+
  249. theme_classic()+
  250. ylab("Speed (cm/sec)")+
  251. xlab("Rescue")+
  252. theme(legend.position="none")
  254. speed_gen_lm <- glmmTMB(data= gen_rescue_track, Average_Speed ~ Rescue + (1|Video))
  255. simulateResiduals(speed_gen_lm, plot = T)
  256. summary(speed_gen_lm)
  257. # Family: gaussian ( identity )
  258. # Formula: Average_Speed ~ Rescue + (1 | Video)
  259. # Data: gen_rescue_track
  260. #
  261. # AIC BIC logLik -2*log(L) df.resid
  262. # 26.6 32.9 -9.3 18.6 31
  263. #
  264. # Random effects:
  265. #
  266. # Conditional model:
  267. # Groups Name Variance Std.Dev.
  268. # Video (Intercept) 0.03561 0.1887
  269. # Residual 0.07880 0.2807
  270. # Number of obs: 35, groups: Video, 7
  271. #
  272. # Dispersion estimate for gaussian family (sigma^2): 0.0788
  273. #
  274. # Conditional model:
  275. # Estimate Std. Error z value Pr(>|z|)
  276. # (Intercept) 0.93005 0.10361 8.977 <2e-16 ***
  277. # RescueYes 0.01454 0.09967 0.146 0.884
  278. speed_gen_lm_wo <- glmmTMB(data= gen_rescue_track, Average_Speed ~ 1 + (1|Video))
  279. anova(speed_gen_lm, speed_gen_lm_wo)
  280. # Data: gen_rescue_track
  281. # Models:
  282. # speed_gen_lm_wo: Average_Speed ~ 1 + (1 | Video), zi=~0, disp=~1
  283. # speed_gen_lm: Average_Speed ~ Rescue + (1 | Video), zi=~0, disp=~1
  284. # Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
  285. # speed_gen_lm_wo 3 24.652 29.318 -9.3260 18.652
  286. # speed_gen_lm 4 26.631 32.852 -9.3154 18.631 0.0213 1 0.884
  289. circle_track<- read.xlsx("~/Documents/04_Catag_rescue/MS_rescue/New_subm/JExpBio/Revision_JExpB_rescue/Supl_material_File2.xlsx", sheet = 5)
  290. colnames(circle_track)
  291. rescuers <- subset(circle_track, Rescue == "Yes")
  292. nonrescuers <- subset(circle_track, Rescue == "No")
  294. rescue_track <- rbind(rescuers, nonrescuers)
  295. rescue_track$Latency_i<- as.numeric(rescue_track$Latency_i)
  296. rescue_track$Latency_o<- as.numeric(rescue_track$Latency_o)
  297. rescue_track$Average_Speed_i<- as.numeric(rescue_track$Average_Speed_i)
  298. rescue_track$Average_Speed_o<- as.numeric(rescue_track$Average_Speed_o)
  300. table( rescue_track$Rescue)
  301. #rescue_track$Rescue <- factor(rescue_track$Rescue, levels = c("Yes", "No"))
  304. inner_plot <- ggplot(data=rescue_track, aes(x=Rescue, y=Average_Speed_i, colour = Rescue))+
  305. geom_beeswarm(priority = "none")+
  306. stat_summary(fun.y= mean, fun.ymin=mean, fun.ymax=mean, geom="crossbar", width=0.3, linewidth =0.3)+
  307. scale_color_manual(values=c("blue3", "red3"))+
  308. theme_classic()+
  309. annotate("text", color = "black",label = "**", x = 1.45, y = 3, size = 7) +
  310. ylab("Speed inner circle (cm/sec)")+
  311. xlab("Rescue")+
  312. theme(legend.position="none")
  314. speed_i_lm <- glmmTMB(data= rescue_track, log(Average_Speed_i) ~ Rescue + (1|Video))
  315. simulateResiduals(speed_i_lm, plot = T)
  316. summary(speed_i_lm)
  317. # Family: gaussian ( identity )
  318. # Formula: log(Average_Speed_i) ~ Rescue + (1 | Video)
  319. # Data: rescue_track
  320. #
  321. # AIC BIC logLik -2*log(L) df.resid
  322. # 48.0 53.2 -20.0 40.0 23
  324. # Random effects:
  325. #
  326. # Conditional model:
  327. # Groups Name Variance Std.Dev.
  328. # Video (Intercept) 0.3268 0.5717
  329. # Residual 0.1430 0.3782
  330. # Number of obs: 27, groups: Video, 7
  331. #
  332. # Dispersion estimate for gaussian family (sigma^2): 0.143
  333. #
  334. # Conditional model:
  335. # Estimate Std. Error z value Pr(>|z|)
  336. # (Intercept) -0.5094 0.2649 -1.923 0.05446 .
  337. # RescueYes -0.5776 0.1776 -3.252 0.00114 **
  338. speed_i_lm_wo <- glmmTMB(data= rescue_track, log(Average_Speed_i) ~ 1 + (1|Video))
  339. anova(speed_i_lm, speed_i_lm_wo)
  340. # Data: rescue_track
  341. # Models:
  342. # speed_i_lm_wo: log(Average_Speed_i) ~ 1 + (1 | Video), zi=~0, disp=~1
  343. # speed_i_lm: log(Average_Speed_i) ~ Rescue + (1 | Video), zi=~0, disp=~1
  344. # Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
  345. # speed_i_lm_wo 3 54.751 58.639 -24.376 48.751
  346. # speed_i_lm 4 47.970 53.153 -19.985 39.970 8.7817 1 0.003043 **
  348. outer_plot <- ggplot(data=rescue_track, aes(x=Rescue, y=Average_Speed_o, colour = Rescue))+
  349. geom_beeswarm(priority = "none")+
  350. stat_summary(fun.y= mean, fun.ymin=mean, fun.ymax=mean, geom="crossbar", width=0.3, linewidth =0.3)+
  351. annotate("text", color = "black",label = "***", x = 1.45, y = 3, size = 7) +
  352. scale_color_manual(values=c("blue3", "red3"))+
  353. theme_classic()+
  354. ylab("Speed outer circle (cm/sec)")+
  355. xlab("Rescue")+
  356. theme(legend.position="none")
  358. speed_o_lm <- glmmTMB(data= rescue_track, log(Average_Speed_o) ~ Rescue + (1|Video))
  359. simulateResiduals(speed_o_lm, plot = T)
  360. summary(speed_o_lm)
  361. # Family: gaussian ( identity )
  362. # Formula: log(Average_Speed_o) ~ Rescue + (1 | Video)
  363. # Data: rescue_track
  364. #
  365. # AIC BIC logLik -2*log(L) df.resid
  366. # 52.5 58.2 -22.3 44.5 27
  367. #
  368. # Random effects:
  369. #
  370. # Conditional model:
  371. # Groups Name Variance Std.Dev.
  372. # Video (Intercept) 0.06873 0.2622
  373. # Residual 0.19971 0.4469
  374. # Number of obs: 31, groups: Video, 7
  375. #
  376. # Dispersion estimate for gaussian family (sigma^2): 0.2
  377. #
  378. # Conditional model:
  379. # Estimate Std. Error z value Pr(>|z|)
  380. # (Intercept) 0.1223 0.1726 0.709 0.478
  381. # RescueYes -0.6902 0.1773 -3.893 9.88e-05 ***
  382. coefficients(speed_o_lm)
  383. exp(0.2329)
  384. exp(-0.6902)
  385. speed_o_lm_wo <- glmmTMB(data= rescue_track, log(Average_Speed_o) ~ 1 + (1|Video))
  386. anova(speed_o_lm, speed_o_lm_wo)
  387. # Data: rescue_track
  388. # Models:
  389. # speed_o_lm_wo: log(Average_Speed_o) ~ 1 + (1 | Video), zi=~0, disp=~1
  390. # speed_o_lm: log(Average_Speed_o) ~ Rescue + (1 | Video), zi=~0, disp=~1
  391. # Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
  392. # speed_o_lm_wo 3 63.059 67.361 -28.529 57.059
  393. # speed_o_lm 4 52.503 58.239 -22.251 44.503 12.556 1 0.0003949 ***
  394. speed_plotC <- gen_plotA+ outer_plot + inner_plot + plot_annotation(tag_levels = 'A')
  396. inner_plot <- ggplot(data=rescue_track, aes(x=Rescue, y=Exploration_relative_value_i, colour = Rescue))+
  397. geom_beeswarm(priority = "none")+
  398. stat_summary(fun.y= mean, fun.ymin=mean, fun.ymax=mean, geom="crossbar", width=0.3, linewidth =0.3)+
  399. scale_color_manual(values=c("blue3", "red3"))+
  400. annotate("text", color = "black",label = "***", x = 1.45, y = 1, size = 7) +
  401. theme_classic()+
  402. ylab("Relative exploration inner circle")+
  403. xlab("Rescue")+
  404. theme(legend.position="none")
  406. explor_i_lm <- glmmTMB(data= rescue_track, Exploration_relative_value_i ~ Rescue + (1|Video))
  407. simulateResiduals(explor_i_lm, plot = T)
  408. summary(explor_i_lm)
  409. # Conditional model:
  410. # Estimate Std. Error z value Pr(>|z|)
  411. # (Intercept) 0.30077 0.07447 4.039 5.38e-05 ***
  412. # RescueYes 0.48710 0.08848 5.505 3.68e-08 ***
  413. explor_i_lm_wo <- glmmTMB(data= rescue_track, Exploration_relative_value_i ~ 1 + (1|Video))
  414. anova(explor_i_lm, explor_i_lm_wo)
  415. # Data: rescue_track
  416. # Models:
  417. # explor_i_lm_wo: Exploration_relative_value_i ~ 1 + (1 | Video), zi=~0, disp=~1
  418. # explor_i_lm: Exploration_relative_value_i ~ Rescue + (1 | Video), zi=~0, disp=~1
  419. # Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
  420. # explor_i_lm_wo 3 34.014 38.680 -14.007 28.014
  421. # explor_i_lm 4 14.640 20.861 -3.320 6.640 21.374 1 3.778e-06 ***
  423. table(rescue_track$Rescue)
  424. outer_plot <- ggplot(data=rescue_track, aes(x=Rescue, y=Exploration_relative_value_o, colour = Rescue))+
  425. geom_beeswarm(priority = "none")+
  426. stat_summary(fun.y= mean, fun.ymin=mean, fun.ymax=mean, geom="crossbar", width=0.3, linewidth =0.3)+
  427. scale_color_manual(values=c("blue3", "red3"))+
  428. theme_classic()+
  429. annotate("text", color = "black",label = "**", x = 1.45, y = 1, size = 7) +
  430. ylab("Relative exploration outer circle")+
  431. xlab("Rescue")+
  432. theme(legend.position="none")
  433. explor_o_lm <- glmmTMB(data= rescue_track, Exploration_relative_value_o ~ Rescue + (1|Video))
  434. simulateResiduals(explor_o_lm, plot = T)
  435. summary(explor_o_lm)
  437. # Conditional model:
  438. # Estimate Std. Error z value Pr(>|z|)
  439. # (Intercept) 0.31087 0.08519 3.649 0.000263 ***
  440. # RescueYes 0.23932 0.08581 2.789 0.005285 **
  441. explor_o_lm_wo <- glmmTMB(data= rescue_track, Exploration_relative_value_o ~ 1 + (1|Video))
  442. anova(explor_o_lm, explor_o_lm_wo)
  443. # Data: rescue_track
  444. # Models:
  445. # explor_o_lm_wo: Exploration_relative_value_o ~ 1 + (1 | Video), zi=~0, disp=~1
  446. # explor_o_lm: Exploration_relative_value_o ~ Rescue + (1 | Video), zi=~0, disp=~1
  447. # Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
  448. # explor_o_lm_wo 3 20.172 24.838 -7.0861 14.172
  449. # explor_o_lm 4 15.302 21.523 -3.6510 7.302 6.8703 1 0.008764 **
  451. generalExpl_plot <- ggplot(data=gen_rescue_track, aes(x=Rescue, y=Exploration_relative_value, colour = Rescue))+
  452. geom_beeswarm(priority = "none")+
  453. stat_summary(fun.y= mean, fun.ymin=mean, fun.ymax=mean, geom="crossbar", width=0.3, linewidth =0.3)+
  454. scale_color_manual(values=c("blue3", "red3"))+
  455. theme_classic()+
  456. annotate("text", color = "black",label = "n.s", x = 1.45, y = 1) +
  457. ylab("Relative exploration arena")+
  458. xlab("Rescue")+
  459. theme(legend.position="none")
  462. explor_lm <- glmmTMB(data= gen_rescue_track, Exploration_relative_value ~ Rescue + (1|Video))
  463. simulateResiduals(explor_lm, plot = T)
  464. summary(explor_lm)
  465. # Conditional model:
  466. # Estimate Std. Error z value Pr(>|z|)
  467. # (Intercept) 0.41451 0.06695 6.191 5.97e-10 ***
  468. # RescueYes 0.03151 0.06690 0.471 0.638
  469. explor_lm_wo <- glmmTMB(data= gen_rescue_track, Exploration_relative_value ~ 1 + (1|Video))
  470. anova(explor_lm, explor_lm_wo)
  471. # Data: gen_rescue_track
  472. # Models:
  473. # explor_lm_wo: Exploration_relative_value ~ 1 + (1 | Video), zi=~0, disp=~1
  474. # explor_lm: Exploration_relative_value ~ Rescue + (1 | Video), zi=~0, disp=~1
  475. # Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
  476. # explor_lm_wo 3 -3.8369 0.8292 4.9184 -9.8369
  477. # explor_lm 4 -2.0586 4.1628 5.0293 -10.0586 0.2218 1 0.6377
  478. Exploration_plot_D <- generalExpl_plot + outer_plot + inner_plot + plot_annotation(tag_levels = 'A')
  482. inner_lat_plot <- ggplot(data=rescue_track, aes(x=Rescue, y=Latency_i, colour = Rescue))+
  483. geom_beeswarm(priority = "none")+
  484. stat_summary(fun.y= mean, fun.ymin=mean, fun.ymax=mean, geom="crossbar", width=0.3, linewidth =0.3)+
  485. scale_color_manual(values=c("blue3", "red3"))+
  486. annotate("text", color = "black",label = "**", x = 1.45, y = 850, size = 7) +
  487. theme_classic()+
  488. ylab("Latency to reach inner circle (s)")+
  489. xlab("Rescue")+
  490. theme(legend.position="none")
  492. latenc_i_lm <- glmmTMB(data= rescue_track, Latency_i ~ Rescue + (1|Video))
  493. simulateResiduals(latenc_i_lm, plot = T)
  494. summary(latenc_i_lm)
  495. # Conditional model:
  496. # Estimate Std. Error z value Pr(>|z|)
  497. # (Intercept) 594.66 89.15 6.671 2.55e-11 ***
  498. # RescueYes -254.03 90.10 -2.819 0.00481 **
  500. latenc_i_lm_wo <- glmmTMB(data= rescue_track, Latency_i ~ 1 + (1|Video))
  501. anova(latenc_i_lm, latenc_i_lm_wo)
  502. # Data: rescue_track
  503. # Models:
  504. # latenc_i_lm_wo: Latency_i ~ 1 + (1 | Video), zi=~0, disp=~1
  505. # latenc_i_lm: Latency_i ~ Rescue + (1 | Video), zi=~0, disp=~1
  506. # Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
  507. # latenc_i_lm_wo 3 379.91 383.80 -186.96 373.91
  508. # latenc_i_lm 4 375.33 380.52 -183.67 367.33 6.5789 1 0.01032 *
  510. outer_lat_plot <- ggplot(data=rescue_track, aes(x=Rescue, y=Latency_o, colour = Rescue))+
  511. geom_beeswarm(priority = "none")+
  512. stat_summary(fun.y= mean, fun.ymin=mean, fun.ymax=mean, geom="crossbar", width=0.3, linewidth =0.3)+
  513. scale_color_manual(values=c("blue3", "red3"))+
  514. annotate("text", color = "black",label = "**", x = 1.45, y = 850, size = 7) +
  515. theme_classic()+
  516. ylab("Latency to reach outer circle (s)")+
  517. xlab("Rescue")+
  518. theme(legend.position="none")
  520. latenc_o_lm <- glmmTMB(data= rescue_track, Latency_o ~ Rescue + (1|Video))
  521. simulateResiduals(latenc_o_lm, plot = T)
  522. summary(latenc_o_lm)
  523. # Conditional model:
  524. # Estimate Std. Error z value Pr(>|z|)
  525. # (Intercept) 484.45 75.43 6.422 1.34e-10 ***
  526. # RescueYes -189.42 72.78 -2.603 0.00925 **
  527. inner_lat_plot + outer_lat_plot
  528. latenc_o_lm_wo <- glmmTMB(data= rescue_track, Latency_o ~ 1 + (1|Video))
  529. anova(latenc_o_lm, latenc_o_lm_wo)
  530. # Data: rescue_track
  531. # Models:
  532. # latenc_o_lm_wo: Latency_o ~ 1 + (1 | Video), zi=~0, disp=~1
  533. # latenc_o_lm: Latency_o ~ Rescue + (1 | Video), zi=~0, disp=~1
  534. # Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
  535. # latenc_o_lm_wo 3 431.25 435.55 -212.62 425.25
  536. # latenc_o_lm 4 427.30 433.04 -209.65 419.30 5.9477 1 0.01474 *
  538. gen_rescue_track<- gen_rescue_track %>%
  539. mutate(across(c(3:14), as.numeric))
  540. travDist_plot <- ggplot(data=gen_rescue_track, aes(x=Rescue, y=Traveled_Dist, colour = Rescue))+
  541. geom_beeswarm(priority = "none")+
  542. stat_summary(fun.y= mean, fun.ymin=mean, fun.ymax=mean, geom="crossbar", width=0.3, linewidth =0.3)+
  543. scale_color_manual(values=c("blue3", "red3"))+
  544. annotate("text", color = "black",label = "n.s", x = 1.45, y = 1150) +
  545. theme_classic()+
  546. ylab("Travelled distance (cm)")+
  547. xlab("Rescue")+
  548. theme(legend.position="none")
  549. dist_lm <- glmmTMB(data= gen_rescue_track, log(Traveled_Dist) ~ Rescue + (1|Video))
  550. simulateResiduals(dist_lm, plot = T)
  551. summary(dist_lm)
  552. # Conditional model:
  553. # Intercept) 4.9692 0.2533 19.619 <2e-16 ***
  554. # RescueYes 0.1112 0.3078 0.361 0.718
  555. ((generalExpl_plot + outer_plot + inner_plot) /
  556. (travDist_plot + outer_lat_plot + inner_lat_plot ) )+ plot_annotation(tag_levels = 'A')
  557. dist_lm_wo <- glmmTMB(data= gen_rescue_track, log(Traveled_Dist) ~ 1 + (1|Video))
  558. simulateResiduals(dist_lm, plot = T)
  559. anova(dist_lm, dist_lm_wo)
  560. # Data: gen_rescue_track
  561. # Models:
  562. # dist_lm_wo: log(Traveled_Dist) ~ 1 + (1 | Video), zi=~0, disp=~1
  563. # dist_lm: log(Traveled_Dist) ~ Rescue + (1 | Video), zi=~0, disp=~1
  564. # Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
  565. # dist_lm_wo 3 99.291 103.96 -46.645 93.291
  566. # dist_lm 4 101.160 107.38 -46.580 93.160 0.1305 1 0.7179
  568. ######################## GOterms #######################
  569. library(openxlsx)
  570. library(ggplot2)
  571. library(Rgraphviz)
  572. library(topGO)
  573. library("RColorBrewer")
  574. library(tagcloud)
  577. #### MBs rescuers vs non-rescuers ####
  578. MB_up_genes <- as.data.frame(c("gene_3661", "gene_3662", "gene_3426", "gene_1607", "gene_2555", "gene_1243" , "gene_930" , "gene_5298" , "gene_6089", "gene_5245", "gene_13149", "gene_9895", "gene_9900", "gene_8499", "gene_7997"))
  580. colnames(MB_up_genes) <- c("genes")
  582. #### GOterms Upregulated when prone to rescue
  583. #we start with an enrichment analysis
  584. #load the universe reference file you created above using bash
  585. universe <- readMappings("~/Documents/Catagl_cognition/From_Eyal/Cnig_gn3.1.fasta.transcripts.fasta.trinotate.report.topGOinput.txt")
  586. # add newly found annotations:
  587. universe$gene_3426 <- c("GO:0008061", "GO:0005576")
  588. universe$gene_5298 <- c("GO:0005549")
  589. universe$gene_1607 <- c("GO:0035176","GO:0005549", "GO:0005615")
  590. universe$gene_9900 <- c("GO:0006955","GO:0005164", "GO:0005515", "GO:0016020", "GO:0048018", "GO:0005615")
  591. universe$gene_5301 <- c("GO:0035176","GO:0005549", "GO:0005615")
  593. # generate list of the protein IDs
  594. universe_genelist <- names(universe)
  596. # generate a reference list of the protein IDs
  597. #testset_genelist_idown <- readLines(LD_down_genes)
  598. #determine matches between both datasets
  599. gene_list_iup <- factor(as.integer(universe_genelist %in% MB_up_genes$genes))
  600. names(gene_list_iup) <- universe_genelist
  601. #set the ontology type you want to investigate, we first choose "Biological Processes"=BP, '93Cellular Components'94=CC, '93Molecular Function'94=MF
  602. ontology_type = "BP" # or CC or MF
  603. GO_data_iup <- new("topGOdata", description="GO_Enrichment", ontology=ontology_type, allGenes=gene_list_iup, annot = annFUN.gene2GO, gene2GO=universe)
  605. result_topGO_w01_iup <- runTest(GO_data_iup, algorithm = "weight01", statistic = "fisher")
  606. result_table_w01_iup_BP <- GenTable(GO_data_iup, Fisher = result_topGO_w01_iup, orderBy = "Fisher", ranksOf = "Fisher", topNodes = numSigGenes(GO_data_iup))
  607. result_table_w01_iup_BP
  608. #### Graphics ####
  609. result_table_w01_iup_BP_sig <- result_table_w01_iup_BP[which(result_table_w01_iup_BP$Fisher<0.05),]
  610. setwd("~/Documents/Catag_rescue/")
  611. #write.csv(result_table_w01_idown_BP_sig, file=paste0("topGO_result_weight01_significant_OL_Corrprop_BP_down_cooks.csv")
  612. colors <- colorRampPalette( brewer.pal( 12, "Paired" ) )( length(result_table_w01_iup_BP_sig$Term))
  613. #Create the tag cloud
  614. tagcloud(strmultline(result_table_w01_iup_BP_sig$Term), weights = - log(as.numeric(result_table_w01_iup_BP_sig$Fisher)),col=colors)
  618. ## DESeq2 - mRNA sequencing data from Cat niger. 7 samples of rescuers and non-rescuers. The brain tissue was sequenced by separating the optic lobes (OL), the mushroom bodies (MB), the rest of the brain tissue including antennal lobes and central complex (CC), and the antennae (AA).
  620. tissues <- read.xlsx("~/Documents/Catag_rescue/Sequencing_rescue/Gene_counts/Rescue_rawcounts_perGene.xlsx", colNames = T)
  621. tissue_counts <- tissues[-c(1:4),]
  622. colnames(tissue_counts)[1] <- c("Gene_ID")
  623. rownames(tissue_counts) <- tissue_counts$Gene_ID
  624. tissue_counts <- tissue_counts[,-c(1)]
  625. colnames(tissue_counts)
  626. # [1] "Ant_A_727" "Ant_B_724" "Ant_C_723" "Ant_D_727" "Ant_E_718" "Ant_F_718" "Ant_G_720" "Ant_H_720"
  627. # [9] "Ant_I_723" "Ant_J_724" "Ant_K_725" "Ant_L_725" "Ant_M_726" "Ant_N_726" "CC_A_727" "CC_B_724"
  628. # [17] "CC_C_723" "CC_D_727" "CC_F_718" "CC_G_720" "CC_H_720" "CC_I_723" "CC_J_724" "CC_K_725"
  629. # [25] "CC_L_725" "CC_M_726" "CC_N_726" "MB_A_727" "MB_C_723" "MB_D_727" "MB_E_718" "MB_F_718"
  630. # [33] "MB_G_720" "MB_H_720" "MB_I_723" "MB_J_724" "MB_K_725" "MB_L_725" "MB_M_726" "MB_N_726"
  631. # [41] "OL_A_727" "OL_B_724" "OL_C_723" "OL_E_718" "OL_F_718" "OL_G_720" "OL_H_720" "OL_I_723"
  632. # [49] "OL_J_724" "OL_K_725" "OL_L_725" "OL_M_726" "OL_N_726"
  634. # First check differences among tissues to highlight any outliers, or if the dissection failed.
  636. design_rescue<- read.xlsx("~/Documents/Catag_rescue/Samples_rescue.xlsx", sheet = 6)
  637. head(design_rescue)
  638. # Sample_ID Tissue Pool Colony_ID Rescue
  639. # 1 Ant_A_727 Ant A 727 No
  640. # 2 Ant_B_724 Ant B 724 Yes
  641. # 3 Ant_C_723 Ant C 723 Yes
  642. # 4 Ant_D_727 Ant D 727 Yes
  643. # 5 Ant_E_718 Ant E 718 No
  644. # 6 Ant_F_718 Ant F 718 Yes
  645. design_rescue[,1]
  646. # [1] "Ant_A_727" "Ant_B_724" "Ant_C_723" "Ant_D_727" "Ant_E_718" "Ant_F_718" "Ant_G_720" "Ant_H_720"
  647. # [9] "Ant_I_723" "Ant_J_724" "Ant_K_725" "Ant_L_725" "Ant_M_726" "Ant_N_726" "CC_A_727" "CC_B_724"
  648. # [17] "CC_C_723" "CC_D_727" "CC_F_718" "CC_G_720" "CC_H_720" "CC_I_723" "CC_J_724" "CC_K_725"
  649. # [25] "CC_L_725" "CC_M_726" "CC_N_726" "MB_A_727" "MB_C_723" "MB_D_727" "MB_E_718" "MB_F_718"
  650. # [33] "MB_G_720" "MB_H_720" "MB_I_723" "MB_J_724" "MB_K_725" "MB_L_725" "MB_M_726" "MB_N_726"
  651. # [41] "OL_A_727" "OL_B_724" "OL_C_723" "OL_E_718" "OL_F_718" "OL_G_720" "OL_H_720" "OL_I_723"
  652. # [49] "OL_J_724" "OL_K_725" "OL_L_725" "OL_M_726" "OL_N_726"
  653. design_rescue$Colony_ID <- as.factor(design_rescue$Colony_ID)
  654. design_rescue$Rescue <- as.factor(design_rescue$Rescue)
  655. design_rescue$Tissue <- as.factor(design_rescue$Tissue)
  657. dds_alltissues <- DESeqDataSetFromMatrix(countData = tissue_counts, colData = design_rescue , design = ~ Tissue + Rescue)
  658. # converting counts to integer mode
  659. keep <- rowSums(counts(dds_alltissues) >= 10) >= 6
  660. dds_alltissues <- dds_alltissues[keep,]
  662. # Visualization - let's first have a look at the data in a PCA plot. For this the data need to be transformed, for which several methods are described within the manual. Choose one, e.g. varianceStabilizingTransformation and carry out the command on the dds object. Then run the command plotPCA on the new transformed object, naming the grouping variables with intgroup
  663. dds_vst_allTissues <- varianceStabilizingTransformation(dds_alltissues)
  664. pcaData_allTissues <- plotPCA(dds_vst_allTissues, intgroup = c("Tissue", "Colony_ID"), returnData = T)
  665. # using ntop=500 top features by variance
  667. pca_allTissues <- ggplot(pcaData_allTissues, aes(x = PC1, y = PC2, label =name, color = Tissue,shape =Colony_ID )) +
  668. geom_point(size=3) +
  669. # geom_label_repel(label.size = 0)+
  670. # geom_text(aes(label=name), vjust=3, size=4,nudge_x = 3, nudge_y = 3)+
  671. xlab("PC1: 85% variance")+
  672. ylab("PC2 : 9% variance")+
  673. ggtitle("PCA all Tissues Rescue")+
  674. scale_shape_manual(values = c(0, 5, 10, 14, 1,8,17))+
  675. theme_minimal()+
  676. theme(panel.border = element_blank(),
  677. panel.grid.major = element_blank(),
  678. panel.grid.minor = element_blank(),
  679. axis.line = element_line(colour="black"),
  680. legend.position = "bottom",
  681. legend.text=element_text(size=12),
  682. axis.title = element_text(size = 12),
  683. axis.text = element_text(size = 12))
  684. pca_allTissues
  686. # Wald test in all tissues, but then test one to one each tissue, better all tissues that have DEGs with an LRT test
  687. # dds_alltissues <- DESeq(dds_alltissues)
  688. # Wald_rescue_results <- results(dds_alltissues, contrast = c("Rescue", "Yes", "No"), alpha = 0.05)
  689. # Wald_rescue_results_sig <- subset(Wald_rescue_results , padj < 0.05)
  690. # Wald_rescue_results_sig # DataFrame with
  691. # A reduced model without the rescue variable
  692. dds_allTissue_LTR_rescue <- DESeq(dds_alltissues, test="LRT", full= ~ Tissue + Rescue , reduced = ~Tissue)
  694. LRT_results_allTissue_rescue <- results(dds_allTissue_LTR_rescue)
  695. #write.csv(LRT_results_allTissue_rescue, "DEGs_alltissues.csv")
  696. LRT_results_allTissue_rescue_sig <- subset(LRT_results_allTissue_rescue , padj < 0.05)
  697. LRT_results_allTissue_rescue_sig # DataFrame with 9 rows and 6 columns
  698. as.data.frame(LRT_results_allTissue_rescue_sig)
  699. # baseMean log2FoldChange lfcSE stat pvalue padj
  700. # gene_930 278.50568 2.8281217 0.47015841 27.54054 1.538359e-07 0.001583740
  701. # gene_2555 117.04428 2.7623696 0.53417219 23.27896 1.401249e-06 0.004808619
  702. # gene_5245 2322.00176 1.0977149 0.22555155 22.50527 2.095677e-06 0.005393748
  703. #non-shared with MBs genes
  704. # gene_1997 34.03176 0.7625534 0.17878296 17.96929 2.244977e-05 0.029849548
  705. # gene_13318 789.81013 0.7336380 0.17270608 17.68300 2.609480e-05 0.029849548
  706. # gene_13632 631.05258 -0.2551454 0.05794553 19.29975 1.117209e-05 0.023003325
  707. # gene_14429 450.17155 0.3543307 0.08393271 17.79513 2.460122e-05 0.029849548
  708. # gene_14458 2677.47760 0.5531638 0.10852223 25.92235 3.554297e-07 0.001829574
  709. # gene_7603 95.96581 0.5167914 0.11947306 18.71925 1.514459e-05 0.025985597
  710. LRT_results_allTissue_rescue_sig@rownames
  711. alltissues_genes <- c("gene_1997", "gene_2555", "gene_930", "gene_5245", "gene_13318", "gene_13632", "gene_14429", "gene_14458", "gene_7603")
  713. LRT_results_allTissue_rescue_sig$gene <- rownames(LRT_results_allTissue_rescue_sig)
  714. df_allT <- lapply(LRT_results_allTissue_rescue_sig$gene, (x) {
  715. y <- plotCounts(dds_allTissue_LTR_rescue, returnData=TRUE, x, intgroup = "Rescue")
  716. y$feature <- x
  717. return(y)
  718. })
  720. df_allTs <- do.call(rbind, df_allT)
  721. df_allTs<- df_allTs %>%
  722. group_by(feature) %>%
  723. mutate(y.position=max(count))
  726. padj_allTs <- as.data.frame(LRT_results_allTissue_rescue_sig[[6]])
  727. padj_allTs$feature <- LRT_results_allTissue_rescue_sig[[7]]
  728. padj_allTs$`LRT_results_allTissue_rescue_sig[[6]]` <- lapply(padj_allTs$`LRT_results_allTissue_rescue_sig[[6]]`, signif, digits=3)
  730. large_padj_allTs <- merge(padj_allTs, df_allTs, by = "feature", all= FALSE)
  731. padj_allTs_small <- large_padj_allTs[!duplicated(large_padj_allTs$feature), ]
  733. plot_a_list <- function(master_list_with_plots, no_of_rows, no_of_cols) {
  735. patchwork::wrap_plots(master_list_with_plots,
  736. nrow = no_of_rows, ncol = no_of_cols)
  737. }
  740. # apply seperately for each gene showing the mean
  741. p_allT <- lapply(LRT_results_allTissue_rescue_sig$gene, function(gene) {
  742. ggplot(df_allTs[df_allTs[, "feature"]==gene,]) +
  743. geom_beeswarm(aes(x=Rescue, y=count, color= Rescue),priority = "none")+
  744. stat_summary(aes(x=Rescue, y=count), fun.y= mean, fun.ymin=mean, fun.ymax=mean, geom="crossbar", width=0.3, linewidth =0.3)+
  745. # annotate("text", label= "Padj=", x= 1,y=(padj_MBs_small[match( gene, padj_MBs_small$feature),5]))+
  746. geom_text(data = as.data.frame(padj_allTs_small),
  747. aes(label = (padj_allTs_small[match( gene,feature),2]),
  748. x = 1.45, y = (padj_allTs_small[match( gene, feature),5]))) +
  749. ylab("Normalized count")+
  750. scale_color_manual(values=c("blue3", "red3"))+
  751. theme_classic()+
  752. ggtitle(paste(gene))+
  753. theme(axis.title.y=element_blank(),
  754. legend.position="none")
  755. })
  758. # finally print your plots
  759. plot_a_list(p_allT, 3, 3)
  761. ##### 15 genes in MB, how are they in other tissues? ###
  762. MB_genes
  764. df_allT_15MBs <- lapply(MB_genes, (x) {
  765. y <- plotCounts(dds_allTissue_LTR_rescue, returnData=TRUE, x, intgroup =c("Rescue","Tissue"))
  766. y$feature <- x
  767. return(y)
  768. })
  770. p_allT_15 <- lapply(df_allT_15MBs, function(df) {
  771. ggplot(df) +
  772. geom_beeswarm(aes(x=Tissue, y=count, color = Rescue),priority = "none")+
  773. stat_summary(aes(x=Tissue, y=count, color = Rescue), fun.y= mean, fun.ymin=mean, fun.ymax=mean, geom="crossbar", width=0.3, linewidth =0.3)+
  774. ylab("Normalized count")+
  775. scale_color_manual(values=c("blue3", "red3"))+
  776. theme_classic()+
  777. ggtitle(paste(df$feature))+
  778. theme(axis.title.y=element_blank(),
  779. legend.position="none")
  780. })
  781. plot_a_list(p_allT_15, 5, 3)
  783. gene_14458 <- plotCounts(dds_allTissue_LTR_rescue, returnData=TRUE, "gene_14458", intgroup =c("Rescue","Tissue"))
  784. ggplot(data = gene_14458) +
  785. geom_beeswarm(aes(x=Tissue, y=count, color = Rescue),priority = "none")+
  786. stat_summary(aes(x=Tissue, y=count, color = Rescue), fun.y= mean, fun.ymin=mean, fun.ymax=mean, geom="crossbar", width=0.3, linewidth =0.3)+
  787. ylab("Normalized count")+
  788. scale_color_manual(values=c("blue3", "red3"))+
  789. theme_classic()+
  790. # ggtitle(paste(df$feature))+
  791. theme(legend.position="right")
  793. # A reduced model without the tissue variable
  794. dds_allTissue_LTR_tissue <- DESeq(dds_alltissues, test="LRT", full= ~ Tissue +Rescue, reduced = ~Rescue)
  796. LRT_results_allTissue_tissue <- results(dds_allTissue_LTR_tissue)
  797. LRT_results_allTissue_tissue_sig <- subset(LRT_results_allTissue_tissue , padj < 0.05)
  798. LRT_results_allTissue_tissue_sig # DataFrame with 9913 rows and 6 columns
  799. LRT_results_allTissue_tissue_sig$gene <- rownames(LRT_results_allTissue_tissue_sig)
  800. LRT_results_allTissue_frame <- as.data.frame(LRT_results_allTissue_tissue_sig)
  801. #write.xlsx(LRT_results_allTissue_frame, "DEGs_alltissues.xlsx")
  803. # DESEq per tissue
  805. ##### Optic lobes ####
  806. design_OL <- subset(design_rescue, Tissue =="OL")
  807. OL_counts <- tissue_counts[,colnames(tissue_counts) %in% design_OL$Sample_ID]
  808. dds_OL <- DESeqDataSetFromMatrix(countData = OL_counts, colData = design_OL , design = ~ Colony_ID + Rescue)
  809. # converting counts to integer mode
  810. keep_OL <- rowSums(counts(dds_OL) >= 10) >= 6
  811. dds_OL <- dds_OL[keep_OL,]
  813. # A reduced model without the rescue factor
  814. dds_OL_ltr_rescue <- DESeq(dds_OL, test="LRT", full= ~ Colony_ID + Rescue , reduced = ~ Colony_ID)
  816. LRT_results_rescue_OL <- results(dds_OL_ltr_rescue)
  817. #write.csv(LRT_results_rescue_OL, "DEGs_OL.csv")
  818. LRT_results_rescue_OL_sig <- subset(LRT_results_rescue_OL, padj < 0.05)
  819. LRT_results_rescue_OL_sig # 1 gene "affected"
  820. # log2 fold change (MLE): Rescue Yes vs No
  821. # LRT p-value: '~ Colony_ID + Rescue' vs '~ Colony_ID'
  822. # DataFrame with 1 row and 6 columns
  823. # baseMean log2FoldChange lfcSE stat pvalue padj
  824. # <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
  825. # gene_5301 399.255 1.3244 0.240348 28.644 8.6982e-08 0.000842073
  826. opticlobes_genes <- c("gene_5301")
  828. LRT_results_rescue_OL_sig$gene <- rownames(LRT_results_rescue_OL_sig)
  830. y <- plotCounts(dds_OL_ltr_rescue , returnData=TRUE,opticlobes_genes , intgroup = "Rescue")
  831. y$feature <- c("gene_5301")
  832. y.position=max(y$count)
  834. padj_OLs <- as.data.frame(LRT_results_rescue_OL_sig[[6]])
  835. padj_OLs$feature <- LRT_results_rescue_OL_sig[[7]]
  836. padj_OLs$`LRT_results_rescue_OL_sig[[6]]` <- lapply(padj_OLs$`LRT_results_rescue_OL_sig[[6]]`, signif, digits=3)
  838. large_padj_OLs <- merge(padj_OLs, y, by = "feature", all= FALSE)
  839. padj_OLs_small <- large_padj_OLs[!duplicated(large_padj_OLs$feature), ]
  841. plot_a_list <- function(master_list_with_plots, no_of_rows, no_of_cols) {
  843. patchwork::wrap_plots(master_list_with_plots,
  844. nrow = no_of_rows, ncol = no_of_cols)
  845. }
  847. # apply seperately for each gene showing the mean
  848. p_OL <- ggplot(y) +
  849. geom_beeswarm(aes(x=Rescue, y=count, color= Rescue),priority = "none")+
  850. stat_summary(aes(x=Rescue, y=count), fun.y= mean, fun.ymin=mean, fun.ymax=mean, geom="crossbar", width=0.3, linewidth =0.3)+
  851. # annotate("text", label= "Padj=", x= 1,y=(padj_MBs_small[match( gene, padj_MBs_small$feature),5]))+
  852. geom_text(data = as.data.frame(padj_OLs_small),
  853. aes(label = (padj_OLs_small[1,2]),
  854. x = 1.45, y = y.position)) +
  855. ylab("Normalized count")+
  856. scale_color_manual(values=c("blue3", "red3"))+
  857. theme_classic()+
  858. ggtitle(paste("gene_5301"))+
  859. theme(axis.title.y=element_blank(),
  860. legend.position="none")
  863. # finally print your plots
  864. plot_a_list(p_MB, 4, 4)
  866. ##### Mushroom Bodies ####
  867. design_MB <- subset(design_rescue, Tissue =="MB")
  868. MB_counts <- tissue_counts[,colnames(tissue_counts) %in% design_MB$Sample_ID]
  869. dds_MB <- DESeqDataSetFromMatrix(countData = MB_counts, colData = design_MB , design = ~ Colony_ID + Rescue)
  870. # converting counts to integer mode
  871. keep_MB <- rowSums(counts(dds_MB) >= 10) >= 6
  872. dds_MB <- dds_MB[keep_MB,]
  873. # A reduced model without the rescue factor
  874. dds_MB_ltr_rescue <- DESeq(dds_MB, test="LRT", full= ~ Colony_ID + Rescue , reduced = ~ Colony_ID)
  876. LRT_results_rescue_MB <- results(dds_MB_ltr_rescue)
  877. #write.csv(LRT_results_rescue_MB, "DEGs_MB.csv")
  878. LRT_results_rescue_MB_sig <- subset(LRT_results_rescue_MB, padj < 0.05)
  879. LRT_results_rescue_MB_sig
  880. # log2 fold change (MLE): Rescue Yes vs No
  881. # LRT p-value: '~ Colony_ID + Rescue' vs '~ Colony_ID'
  882. # DataFrame with 15 rows and 6 columns
  883. # baseMean log2FoldChange lfcSE stat pvalue padj
  884. # <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
  885. # gene_3661 156.4731 2.496080 0.426528 28.9972 7.24820e-08 9.90312e-05
  886. # gene_3662 98.3451 1.977841 0.435669 17.4028 3.02374e-05 2.10276e-02
  887. # gene_3426 250.8087 0.962072 0.225606 17.3690 3.07807e-05 2.10276e-02
  888. # gene_1607 167.2192 1.373933 0.242181 30.4875 3.36027e-08 6.42752e-05
  889. # gene_2555 56.0721 2.647648 0.446123 30.8862 2.73609e-08 6.42752e-05
  890. # ... ... ... ... ... ... ...
  891. # gene_13149 1019.5838 0.505851 0.114291 19.4021 1.05890e-05 9.43983e-03
  892. # gene_9895 64.7391 2.145827 0.369639 31.1385 2.40262e-08 6.42752e-05
  893. # gene_9900 169.3375 0.948923 0.234393 15.7326 7.29542e-05 4.65156e-02
  894. # gene_8499 46.5866 1.933660 0.399757 21.4861 3.56396e-06 4.26071e-03
  895. # gene_7997 137.1182 1.403238 0.301844 19.9517 7.94241e-06 8.44014e-03
  897. LRT_results_rescue_MB_sig@rownames
  898. MB_genes <-c( "gene_3661", "gene_3662", "gene_3426", "gene_1607", "gene_2555", "gene_1243" , "gene_930" , "gene_5298" , "gene_6089", "gene_5245", "gene_13149", "gene_9895", "gene_9900", "gene_8499", "gene_7997")
  899. library(ggVennDiagram)
  900. list_genes <- list(alltissues_genes, opticlobes_genes, MB_genes)
  901. Venn_genes <- ggVennDiagram(list_genes, label = c("count"), label_alpha = 0, category.names = c("All tissues", "OL", "MB"))+
  902. scale_fill_gradient(low = "#FFFFFF", high = "#FF0033")
  905. dds_vst_MB <- varianceStabilizingTransformation(dds_MB)
  906. pcaData_MB <- plotPCA(dds_vst_MB, intgroup = c("Rescue", "Colony_ID", "Pool"), returnData = T)
  908. # using ntop=500 top features by variance
  910. pca_MB <- ggplot(pcaData_MB, aes(x = PC1, y = PC2, color = Rescue, shape = Colony_ID)) +
  911. geom_point(size=3) +
  912. # geom_label_repel(label.size = 0)+
  913. xlab("PC1: 33% variance")+
  914. ylab("PC2 : 19% variance")+
  915. ggtitle("PCA MB Rescue")+
  916. # scale_color_manual(values = c("#C77CFF", "#00BFC4", "#7CAE00", "#F8766D"))+
  917. theme_minimal()+
  918. scale_shape_manual(values = c(0, 5, 10, 14, 1,8,17))+
  919. scale_color_manual(values=c("blue3", "red3"))+
  920. theme(panel.border = element_blank(),
  921. panel.grid.major = element_blank(),
  922. panel.grid.minor = element_blank(),
  923. axis.line = element_line(colour="black"),
  924. legend.position = "bottom",
  925. legend.text=element_text(size=12),
  926. axis.title = element_text(size = 12),
  927. axis.text = element_text(size = 12))
  928. pca_MB
  930. Venn_genes + pca_MB + plot_layout(widths = c(1,2))
  933. plotMA(LRT_results_rescue_MB, ylim=c(-2,2), alpha = 0.05) # logarithmic fold change log2(yes/no).
  934. # resLFC <- lfcShrink(dds_MB_ltr_rescue, coef="Rescue_Yes_vs_No", type="apeglm")
  935. # plotMA(resLFC, ylim=c(-2,2), alpha = 0.05) # shrunken log2 fold changes, which remove the noise associated with log2 fold changes from low count genes without requiring arbitrary filtering thresholds.
  937. #### Plot MB DEGs count ###
  939. #The counts are which normalizes counts by the estimated size factors (or normalization factors if these were used) and adds a pseudocount of 1/2 to allow for log scale plotting.
  941. #plotCounts(dds_MB_ltr_rescue, gene=which.min(LRT_results_rescue_MB$padj), intgroup="Rescue")
  942. LRT_results_rescue_MB_sig$gene <- rownames(LRT_results_rescue_MB_sig)
  943. df_MB <- lapply(LRT_results_rescue_MB_sig$gene, (x) {
  944. y <- plotCounts(dds_MB_ltr_rescue, returnData=TRUE, x, intgroup = "Rescue")
  945. y$feature <- x
  946. return(y)
  947. })
  951. df_MBs <- do.call(rbind, df_MB)
  952. df_MBs<- df_MBs %>%
  953. group_by(feature) %>%
  954. mutate(y.position=max(count))
  956. df_MBs$sample <- rep(sort_samples_data19_compl$shortID,length(df_3))
  957. df_threeruns$number_runs <- rep(sort_samples_data19_compl$number_runs,length(df_3))
  959. padj_MBs <- as.data.frame(LRT_results_rescue_MB_sig[[6]])
  960. padj_MBs$feature <- LRT_results_rescue_MB_sig[[7]]
  961. padj_MBs$`LRT_results_rescue_MB_sig[[6]]` <- lapply(padj_MBs$`LRT_results_rescue_MB_sig[[6]]`, signif, digits=3)
  963. large_padj_MBs <- merge(padj_MBs, df_MBs, by = "feature", all= FALSE)
  964. padj_MBs_small <- large_padj_MBs[!duplicated(large_padj_MBs$feature), ]
  966. plot_a_list <- function(master_list_with_plots, no_of_rows, no_of_cols) {
  968. patchwork::wrap_plots(master_list_with_plots,
  969. nrow = no_of_rows, ncol = no_of_cols)
  970. }
  972. # apply seperately for each gene showing the mean
  973. p_MB <- lapply(LRT_results_rescue_MB_sig$gene, function(gene) {
  974. ggplot(df_MBs[df_MBs[, "feature"]==gene,]) +
  975. geom_beeswarm(aes(x=Rescue, y=count, color= Rescue),priority = "none")+
  976. stat_summary(aes(x=Rescue, y=count), fun.y= mean, fun.ymin=mean, fun.ymax=mean, geom="crossbar", width=0.3, linewidth =0.3)+
  977. # annotate("text", label= "Padj=", x= 1,y=(padj_MBs_small[match( gene, padj_MBs_small$feature),5]))+
  978. geom_text(data = as.data.frame(padj_MBs_small),
  979. aes(label = (padj_MBs_small[match( gene,feature),2]),
  980. x = 1.45, y = (padj_MBs_small[match( gene, feature),5]))) +
  981. ylab("Normalized count")+
  982. scale_color_manual(values=c("blue3", "red3"))+
  983. theme_classic()+
  984. ggtitle(paste(gene))+
  985. theme(axis.title.y=element_blank(),
  986. legend.position="none")
  987. })
  990. # finally print your plots
  991. plot_a_list(p_MB, 4, 4)
  994. ##### Central Complex+ ####
  995. design_CC <- subset(design_rescue, Tissue =="CC")
  996. CC_counts <- tissue_counts[,colnames(tissue_counts) %in% design_CC$Sample_ID]
  997. dds_CC <- DESeqDataSetFromMatrix(countData = CC_counts, colData = design_CC , design = ~ Colony_ID + Rescue)
  998. # converting counts to integer mode
  999. keep_CC <- rowSums(counts(dds_CC) >= 10) >= 6
  1000. dds_CC <- dds_CC[keep_CC,]
  1001. # A reduced model without the rescue factor
  1002. dds_CC_ltr_rescue <- DESeq(dds_CC, test="LRT", full= ~ Colony_ID + Rescue , reduced = ~ Colony_ID)
  1004. LRT_results_rescue_CC <- results(dds_CC_ltr_rescue)
  1005. #write.csv(LRT_results_rescue_CC, "DEGs_CC.csv")
  1006. LRT_results_rescue_CC_sig <- subset(LRT_results_rescue_CC, padj < 0.05)
  1007. LRT_results_rescue_CC_sig
  1008. # DataFrame with 0 rows and 6 columns
  1010. ##### Antennae ####
  1011. design_Ant <- subset(design_rescue, Tissue =="Ant")
  1012. Ant_counts <- tissue_counts[,colnames(tissue_counts) %in% design_Ant$Sample_ID]
  1013. dds_Ant <- DESeqDataSetFromMatrix(countData = Ant_counts, colData = design_Ant , design = ~ Colony_ID + Rescue)
  1014. # converting counts to integer mode
  1015. keep_Ant <- rowSums(counts(dds_Ant) >= 10) >= 6
  1016. dds_Ant <- dds_Ant[keep_Ant,]
  1017. # A reduced model without the rescue factor
  1018. dds_Ant_ltr_rescue <- DESeq(dds_Ant, test="LRT", full= ~ Colony_ID + Rescue , reduced = ~ Colony_ID)
  1020. LRT_results_rescue_Ant <- results(dds_Ant_ltr_rescue)
  1021. #write.csv(LRT_results_rescue_Ant, "DEGs_antennae.csv")
  1022. LRT_results_rescue_Ant_sig <- subset(LRT_results_rescue_Ant, padj < 0.05)
  1023. LRT_results_rescue_Ant_sig
  1024. # 0
  1025. LRT_results_rescue_Ant_df<- as.data.frame(LRT_results_rescue_Ant)
  1029. #### GOterms upregulated in MB
  1031. ontology_type = "CC" # or CC or MF
  1032. GO_data_iup <- new("topGOdata", description="GO_Enrichment", ontology=ontology_type, allGenes=gene_list_iup, annot = annFUN.gene2GO, gene2GO=universe)
  1033. result_topGO_w01_iup <- runTest(GO_data_iup, algorithm = "weight01", statistic = "fisher")
  1034. result_table_w01_iup_CC <- GenTable(GO_data_iup, Fisher = result_topGO_w01_iup, orderBy = "Fisher", ranksOf = "Fisher", topNodes = numSigGenes(GO_data_iup))
  1035. result_table_w01_iup_CC
  1037. ontology_type = "MF" # or CC or MF
  1038. GO_data_iup <- new("topGOdata", description="GO_Enrichment", ontology=ontology_type, allGenes=gene_list_iup, annot = annFUN.gene2GO, gene2GO=universe)
  1039. result_topGO_w01_iup <- runTest(GO_data_iup, algorithm = "weight01", statistic = "fisher")
  1040. result_table_w01_iup_MF <- GenTable(GO_data_iup, Fisher = result_topGO_w01_iup, orderBy = "Fisher", ranksOf = "Fisher", topNodes = numSigGenes(GO_data_iup))
  1041. result_table_w01_iup_MF
  1044. ####### With the new annotation ####
  1045. setwd("~/Documents/Catagl_cognition/From_Eyal")
  1047. # Load InterProScan GO mappings
  1048. go_raw <- read.delim("go_raw.tsv", header = FALSE, stringsAsFactors = FALSE)
  1049. colnames(go_raw) <- c("transcript_id", "go_terms")
  1050. head(go_raw)
  1052. # Load gene-transcript mapping file
  1053. mapping <- read.csv("Cnig_gn3.1.fasta.transcripts.fasta.trinotate.csv", header = TRUE, stringsAsFactors = FALSE)
  1055. # Preview to confirm columns
  1056. head(mapping[, 1:2])
  1058. # Join GO annotations with gene mapping
  1059. go_merged <- merge(go_raw, mapping[, c("X.gene_id", "transcript_id")], by = "transcript_id")
  1060. head(go_merged)
  1062. library(dplyr)
  1063. library(tidyr)
  1065. # Clean and collapse
  1066. gene_go <- go_merged %>%
  1067. filter(go_terms != "") %>%
  1068. mutate(go_split = strsplit(go_terms, "[,|]")) %>%
  1069. unnest(go_split) %>%
  1070. group_by(X.gene_id) %>%
  1071. summarise(go_terms = paste(unique(go_split), collapse = ",")) %>%
  1072. ungroup()
  1074. # Remove any suffixes like (InterPro), (Pfam), etc.
  1075. gene_go$go_terms <- gsub("(.*?)", "", gene_go$go_terms)
  1077. # Build geneID2GO Object for topGO
  1078. geneID2GO <- setNames(strsplit(gene_go$go_terms, ","), gene_go$X.gene_id)
  1079. write.table(geneID2GO, file = "geneID2GO_universe.tsv", sep = "t", quote = FALSE, row.names = FALSE)
  1080. # generate list of the protein IDs
  1081. universe_genelist <- names(geneID2GO)
  1083. #### MBs rescuers vs non-rescuers ####
  1084. MB_up_genes <- as.data.frame(c("gene_3661", "gene_3662", "gene_3426", "gene_1607", "gene_2555", "gene_1243" , "gene_930" , "gene_5298" , "gene_6089", "gene_5245", "gene_13149", "gene_9895", "gene_9900", "gene_8499", "gene_7997"))
  1086. colnames(MB_up_genes) <- c("genes")
  1088. # generate a reference list of the protein IDs
  1089. #testset_genelist_idown <- readLines(LD_down_genes)
  1090. #determine matches between both datasets
  1091. gene_list_iup <- factor(as.integer(universe_genelist %in% MB_up_genes$genes))
  1092. names(gene_list_iup) <- universe_genelist
  1093. #set the ontology type you want to investigate, we first choose "Biological Processes"=BP, '93Cellular Components'94=CC, '93Molecular Function'94=MF
  1094. ontology_type = "BP" # or CC or MF
  1095. GO_data_iup <- new("topGOdata", description="GO_Enrichment", ontology=ontology_type, allGenes=gene_list_iup, annot = annFUN.gene2GO, gene2GO= geneID2GO)
  1097. result_topGO_w01_iup <- runTest(GO_data_iup, algorithm = "weight01", statistic = "fisher")
  1098. result_table_w01_iup_BP <- GenTable(GO_data_iup, Fisher = result_topGO_w01_iup, orderBy = "Fisher", ranksOf = "Fisher", topNodes = numSigGenes(GO_data_iup))
  1099. result_table_w01_iup_BP
  1100. #### Graphics ####
  1101. result_table_w01_iup_BP_sig <- result_table_w01_iup_BP[which(result_table_w01_iup_BP$Fisher<0.05),]
  1102. setwd("~/Documents/Catag_rescue/")
  1103. #write.csv(result_table_w01_idown_BP_sig, file=paste0("topGO_result_weight01_significant_OL_Corrprop_BP_down_cooks.csv")
  1104. colors <- colorRampPalette( brewer.pal( 12, "Paired" ) )( length(result_table_w01_iup_BP_sig$Term))
  1105. #Create the tag cloud
  1106. tagcloud(strmultline(result_table_w01_iup_BP_sig$Term), weights = - log(as.numeric(result_table_w01_iup_BP_sig$Fisher)),col=colors)
  1110. #### GOterms upregulated in MB
  1112. ontology_type = "CC" # or CC or MF
  1113. GO_data_iup <- new("topGOdata", description="GO_Enrichment", ontology=ontology_type, allGenes=gene_list_iup, annot = annFUN.gene2GO, gene2GO=geneID2GO)
  1114. result_topGO_w01_iup <- runTest(GO_data_iup, algorithm = "weight01", statistic = "fisher")
  1115. result_table_w01_iup_CC <- GenTable(GO_data_iup, Fisher = result_topGO_w01_iup, orderBy = "Fisher", ranksOf = "Fisher", topNodes = numSigGenes(GO_data_iup))
  1116. result_table_w01_iup_CC
  1118. ontology_type = "MF" # or CC or MF
  1119. GO_data_iup <- new("topGOdata", description="GO_Enrichment", ontology=ontology_type, allGenes=gene_list_iup, annot = annFUN.gene2GO, gene2GO=geneID2GO)
  1120. result_topGO_w01_iup <- runTest(GO_data_iup, algorithm = "weight01", statistic = "fisher")
  1121. result_table_w01_iup_MF <- GenTable(GO_data_iup, Fisher = result_topGO_w01_iup, orderBy = "Fisher", ranksOf = "Fisher", topNodes = numSigGenes(GO_data_iup))
  1122. result_table_w01_iup_MF
  1124. ############ reactome ##########
  1125. # Load Reactome mappings
  1126. reactome_raw <- read.delim("reactome_raw.tsv", header = FALSE, stringsAsFactors = FALSE)
  1127. colnames(reactome_raw) <- c("transcript_id", "reactome_id")
  1128. head(reactome_raw)
  1130. # Load gene-transcript mapping file
  1131. mapping <- read.csv("Cnig_gn3.1.fasta.transcripts.fasta.trinotate.csv", header = TRUE, stringsAsFactors = FALSE)
  1133. # Preview to confirm columns
  1134. head(mapping[, 1:2])
  1136. # Join GO annotations with gene mapping
  1137. reactome_merged <- merge(reactome_raw, mapping[, c("X.gene_id", "transcript_id")], by = "transcript_id")
  1138. head(reactome_merged)
  1140. # Clean and collapse
  1141. gene_reactome <- reactome_merged %>%
  1142. filter(reactome_id != "") %>%
  1143. mutate(react_split = strsplit(reactome_id, "[,|]")) %>%
  1144. unnest(react_split) %>%
  1145. group_by(X.gene_id) %>%
  1146. summarise(reactome_id = paste(unique(react_split), collapse = ",")) %>%
  1147. ungroup()
  1149. MB_up_reactome <- merge(MB_up_genes, gene_reactome, by.x = "genes", by.y="X.gene_id", all.x=T, all.y = F)
  1150. MB_up_reactome_list <- setNames(strsplit(MB_up_reactome$reactome_id, ","), MB_up_reactome$genes)
  1151. write.csv(MB_up_reactome, "MB_upDEGs_reactome_pathways.csv")
  1152. library(ReactomeContentService4R)
  1154. reactome_ids_gene9900 <- MB_up_reactome_list$gene_9900
  1155. # Remove "Reactome:" prefix
  1156. reactome_ids_gene9900clean <- sub("^Reactome:", "", reactome_ids_gene9900)
  1158. # Get pathway info
  1159. pathway_info <- lapply(reactome_ids_gene9900clean, function(id) {
  1160. tryCatch(
  1161. ReactomeContentService4R::getPathwaySummation(id),
  1162. error = function(e) NA
  1163. )
  1164. })
  1165. # Combine results
  1166. names(pathway_info) <- reactome_ids_gene9900clean
  1167. pathway_info

Script_Rescue.R at commit c2bd6a5, no license · at the source

Overview

Authors: Luisa Maria Jaimes-Nino1, Adi Bar2, Inon Scharf2, Susanne Foitzik1
ORCID iDs: Luisa Maria Jaimes-Nino, Inon Scharf, Susanne Foitzik
  1. Institute of Organismic and Molecular Evolution, Johannes Gutenberg University Mainz, 55128 Mainz, Germany
  2. School of Zoology, George S Wise Faculty of Life Sciences, Tel Aviv University, Tel Aviv 6997801, Israel
Institutions: Johannes Gutenberg University Mainz (Germany); Tel Aviv University (Israel)
Journal: The Journal of experimental biology, volume 229, issue 14, article jeb252086
Dates: received 19 December 2025; accepted 9 June 2026; published online 17 July 2026; in print July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1242/jeb.252086 · PMID 42332978 · PMCID PMC13405234 · OpenAlex W7165677441
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), other (organism), cellular / molecular (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions
Keywords: Altruistic behaviour, Brain activity, Gene expression, Helping behaviour, Juvenile hormone
MeSH: Ants*, Behavior, Animal*, Transcriptome*, Animals, Brain, Social Behavior (* major topic)
Topic: Insect and Arachnid Ecology and Behavior (Genetics, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: Deutsche Forschungsgemeinschaft (FO 298/31-1, JA 3614/1-1); Israel Science Foundation (699/24); Johannes Gutenberg University of Mainz
Citations: not cited yet (Europe PMC); 88 references in the paper

Abstract

Some social animals can recognize and respond to the distress of group members. While well documented in mammals, such behaviour has also independently evolved in ants. To uncover the molecular basis of this social trait, we compared gene expression in the central and peripheral nervous system of Cataglyphis nigra workers differing in rescue propensity. RNA-seq of mushroom bodies, optic lobes, central brain and antennae revealed nine genes upregulated across all tissues of rescuers, including those encoding cytochrome P450-9e2 (CYP9E2) and arginase, and immune-related genes such as the defensin (DEFA) and serine protease snake (snk) genes. The mushroom bodies showed the strongest signal, including genes involved in juvenile hormone signalling, polyamine synthesis and immune function. These results suggest that the immune pathways and polyamine metabolism modulate social responsiveness or alarm-cue detection in ants. In contrast, morphological and physiological traits did not differ between rescuers and non-rescuers, indicating that rescue behaviour propensity is governed by the activity of the nervous system.

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

ljaimesnino/RescueAnts

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: c2bd6a5378e072d2a806fbd6dfb7ea6737235af5, 24 May 2026
Languages: R (1)
Size: 2 files, 1 script
Software Heritage: not archived
Found in: the end of the paper
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: car (1 file), emmeans (1 file), ggplot2 (1 file), glmmTMB (1 file), patchwork (1 file), tidyverse (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
2 files
At the source: github.com/ljaimesnino/RescueAnts.git

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;
  • 10 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.

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, 4 authors, 5 keywords, 6 MeSH terms, 3 funders, 88 references.

Cite

This paper

Jaimes-Nino, L. M., Bar, A., Scharf, I., & Foitzik, S. (2026). Transcriptional predictors of rescue behaviour in ants. The Journal of experimental biology, 229(14), jeb252086. https://doi.org/10.1242/jeb.252086

BibTeX

@article{jaimesnino2026transcriptional,
author = {Jaimes-Nino, Luisa Maria and Bar, Adi and Scharf, Inon and Foitzik, Susanne},
title = {{Transcriptional predictors of rescue behaviour in ants}},
journal = {The Journal of experimental biology},
year = {2026},
month = jul,
volume = {229},
number = {14},
pages = {jeb252086},
publisher = {The Company of Biologists},
issn = {0022-0949},
doi = {10.1242/jeb.252086},
url = {https://doi.org/10.1242/jeb.252086},
pmid = {42332978},
pmcid = {PMC13405234}
}

RIS

TY - JOUR
AU - Jaimes-Nino, Luisa Maria
AU - Bar, Adi
AU - Scharf, Inon
AU - Foitzik, Susanne
TI - Transcriptional predictors of rescue behaviour in ants
T2 - The Journal of experimental biology
J2 - J Exp Biol
PY - 2026
DA - 2026/07/17
VL - 229
IS - 14
SP - jeb252086
SN - 0022-0949
PB - The Company of Biologists
DO - 10.1242/jeb.252086
UR - https://doi.org/10.1242/jeb.252086
LA - en
ER -

CSL-JSON

{
"id": "10.1242/jeb.252086",
"type": "article-journal",
"title": "Transcriptional predictors of rescue behaviour in ants",
"container-title": "The Journal of experimental biology",
"author": [
{
"family": "Jaimes-Nino",
"given": "Luisa Maria"
},
{
"family": "Bar",
"given": "Adi"
},
{
"family": "Scharf",
"given": "Inon"
},
{
"family": "Foitzik",
"given": "Susanne"
}
],
"container-title-short": "J Exp Biol",
"volume": "229",
"issue": "14",
"page": "jeb252086",
"DOI": "10.1242/jeb.252086",
"PMID": "42332978",
"PMCID": "PMC13405234",
"ISSN": "0022-0949",
"publisher": "The Company of Biologists",
"URL": "https://doi.org/10.1242/jeb.252086",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
17
]
]
}
}

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.1007/s12311-026-02006-1 [code]
Innovative 3D-Image Analysis of Cerebellar Vascularization Highlights Angiogenic Gene Dysregulations in a Murine Model of Apnea of Prematurity.
Journal: Cerebellum (London, England)
In common: glmmTMB, car, emmeans, 3 other tools, genetics / omics, 1 reference
[2] doi:10.3390/ijms27135713 [code]
Chronic Administration of Marinobufagenin in Mice Causes Hyperlocomotion and Decrease in Anxiety by Altering Monoamine Turnover Unaccompanied by Motor Deficits or Oxidative Stress.
Journal: International journal of molecular sciences
In common: glmmTMB, car, emmeans, 3 other tools, cellular / molecular
[3] 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: glmmTMB, emmeans, patchwork, 2 other tools, genetics / omics, 2 references
[4] doi:10.1016/j.celrep.2026.117505 [code]
Impaired spatial coding and neuronal hyperactivity in the medial entorhinal cortex of aged APP knock-in mice.
Journal: Cell reports
In common: glmmTMB, car, emmeans, 3 other tools
[5] doi:10.7554/elife.107088 [code]
Development of auditory and spontaneous movement responses to music over the first postnatal year.
Journal: eLife
In common: glmmTMB, car, emmeans, 2 other tools, 1 reference
[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: car, emmeans, patchwork, 2 other tools, genetics / omics, 2 references
[7] doi:10.1038/s41514-026-00391-9 [code]
Region-specific transcriptional signatures of brain aging in the absence of neuropathology at the single-cell level.
Journal: npj aging
In common: glmmTMB, patchwork, ggplot2, 1 other tool, genetics / omics, cellular / molecular, 1 reference
[8] doi:10.1038/s41586-026-10444-4 [code]
A brain reward circuit inhibited by next-generation weight-loss drugs in mice.
Journal: Nature
In common: glmmTMB, emmeans, patchwork, 2 other tools
[9] doi:10.1038/s41467-026-73858-8 [code]
A pegivirus associated with encephalitis in red-legged partridges shows neurotropism across avian species.
Journal: Nature communications
In common: car, emmeans, patchwork, 2 other tools, other
[10] doi:10.1126/sciadv.aeb8106 [code]
A thyroid hormone-mediated opsin switch initiates metamorphosis in a proto-vertebrate.
Journal: Science advances
In common: car, emmeans, patchwork, 2 other tools, other

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.