1_insect_mortality.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9982b90-c1ea-477f-9b06-c994d047e692",
   "metadata": {},
   "outputs": [],
   "source": [
    "#install.packages(c(\"tidyverse\", \"reshape2\", \"survminer\", \"survival\", \"splines\", \"coxme\", \"multcomp\"))\n",
    "library(tidyverse) # loads dplyr, ggplot2, stringr etc\n",
    "library(reshape2) # for melt and cast\n",
    "\n",
    "# The following packages are for analyzing mortality data and producing Kaplan-Meier plots:\n",
    "library(survminer)\n",
    "library(survival)\n",
    "library(splines)\n",
    "library(coxme)\n",
    "library(multcomp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76f02757-ca77-45aa-b0c7-0a355f902d6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "data <- read.csv(\"1_insect_mortality.csv\",\n",
    "                 header = TRUE)\n",
    "\n",
    "fit <- survfit(Surv(time, death) ~ treatment, data = data)\n",
    "print(fit, show.rmean=TRUE)\n",
    "plot(fit, conf.int=TRUE)\n",
    "summary(fit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0f23d2d-593b-40f1-85ce-ad1f07c2d18f",
   "metadata": {},
   "outputs": [],
   "source": [
    "survdiff(Surv(time,death)~treatment, data=data)\n",
    "\n",
    "comp<-pairwise_survdiff(Surv(time, death) ~ treatment,\n",
    "                        data = data, p.adjust.method = \"fdr\") #log rank compares the difference until an event (death) occurs\n",
    "comp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c003cd2-2f7c-4d23-9ea8-72735d4d3596",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Model with randomization effect by experiment\n",
    "set.seed(2)\n",
    "s1 = survfit(Surv(time,death)~treatment, data=data)\n",
    "print(s1, show.rmean=TRUE)\n",
    "plot(s1, conf.int=TRUE)\n",
    "summary(s1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0e51a81-e53b-49f0-a39e-658b14883ffe",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Statistics\n",
    "#Model random effect\n",
    "data$treatment <- as.factor(data$treatment)\n",
    "model <- coxme(Surv(time, death) ~ treatment+(1|biological_replicate), data=data)\n",
    "anova(model)\n",
    "\n",
    "summary(glht(model, linfct=mcp(treatment=\"Tukey\")), test = adjusted(type = \"bonferroni\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abf01b66-971d-4387-91b3-c98b3e4be439",
   "metadata": {},
   "outputs": [],
   "source": [
    "ggsurvplot(s1, data = data, \n",
    "           main = \"Survival Feeding (11 Days)\",\n",
    "           size = 1, ggtheme = theme_classic(), \n",
    "           linetype=c(\"solid\", \"dotted\", \"dotdash\", \"dashed\"), \n",
    "           xlab=\"Days\", ylab=\"% of survival\", palette=c(\"#887bb9\", \"#a35da5\", \"#69b19a\", \"#8ec372\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e06eef4b-b14a-49d7-ac1f-8d9765389f2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Midpoint and endpoint snapshot\n",
    "\n",
    "#Calculate survival at endpoint\n",
    "time_point = 11 # change for wanted timepoint\n",
    "mortality_at_timepoint = data.frame(Soil = c(\"C10\", \"C112\", \"S16\", \"S7\"), Mortality = (1 - summary(s1, times= time_point, extend = T)$surv)*100)\n",
    "\n",
    "\n",
    "p <- ggplot(mortality_at_timepoint, aes(x = Soil, y = Mortality)) +\n",
    "        geom_bar(position=\"stack\", stat = \"identity\") + \n",
    "        theme_classic() \n",
    "\n",
    "p\n",
    "\n",
    "\n",
    "#Mortality-benzoxazinoid correlation plots\n",
    "#change either use the data on Benzoxanzinoids or Phytohormones.\n",
    "benz = left_join(read.csv(\"4_benzoxazinoids.csv\", header = T), mortality_at_timepoint, by=\"Soil\")\n",
    "\n",
    "ggplot(benz, aes(x=Concentration_mg_g, y=Mortality))+\n",
    "geom_point(shape=16, aes(color=Soil))+\n",
    "stat_cor(method = \"spearman\")+\n",
    "geom_smooth(method = \"lm\", formula = y ~ x, se = T, color = \"black\") +\n",
    "scale_y_continuous(limits = c(0, 100))+\n",
    "ylab(\"Mortality (%)\") + xlab(\"Compound concentration (ng/g)\")+\n",
    "scale_color_manual(values=c(\"#887bb9\", \"#a35da5\", \"#69b19a\", \"#8ec372\"))+\n",
    "facet_wrap(~Treatment+Benzoxazinoid+Sample_Type, scales = \"free_x\", ncol = 5)\n",
    "\n",
    "#The plot needs to be saved in a large pdf to be able to fully visualize each facet."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.3.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}