---
title: "Death Rate Analysis by Lineage - FS vs bottleneck Lineages. Related to Fig S8C"
author: "Danielle Ireland"
date: "2025-06-24"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


## Load libraries

```{r cars}
library(tidyverse)
library(dplyr)
library(readxl)
library(ggplot2)

set.seed(228)
```


## Read in non-bottleneck division data and only keep data for FS lineages

```{r }
div_NB_original<-as.data.frame(read_excel("Dataset_S1.xlsx", sheet="Dataset S1A"))
#only keep data from FS
div_NB_filtered<-div_NB_original %>%
  filter((Family==1) | (Family==4)) 

#create group variables (1=FS, 2=Bottleneck)
div_NB_filtered$Group<-1

div_NB<-div_NB_filtered[c(10, 11,12,14,15)]
#make a new column listing if this is a death (0=no, this is division). to be used to count the total number of worms in a generation
div_NB$Death=0
```

## Read in bottleneck lineage division data. Add Recode family and group names and isolate important variables
```{r }
div_B_original<-as.data.frame(read_excel("Dataset_S1.xlsx", sheet="Dataset S1C"))
div_B_original<-div_B_original %>% mutate(FamilyName=recode(Family,
                                            "0"="B1",
                                             "1"="B2"))
                                            
div_B_original$Group<-2
#subset only the columns we need
div_B<-div_B_original[c(10, 11,12,14,15)]
#make a new column listing if this is a death (0=no, this is division). to be used to count the total number of worms in a generation
div_B$Death=0
```

## Read in nonbottleneck DEATH data and only keep FS
```{r}
dth_NB_original<-as.data.frame(read_excel("Dataset_S1.xlsx", sheet="Dataset S1B"))
#only keep FS
dth_NB_filtered<-dth_NB_original %>%
  filter((Family==1) |(Family==4)) 

#create group variables (1=FS, 2=Bottleneck)
dth_NB_filtered$Group<-1
#subset only the columns we need
dth_NB<-dth_NB_filtered[c(10:14)]
#make a new column listing if this is a death (1=YES). to be used to count the total number of worms in a generation
dth_NB$Death=1
```

## Read in bottleneck lineage DEATH data. Add Recode family and group names and isolate important variables
```{r }
dth_B_original<-as.data.frame(read_excel("Dataset_S1.xlsx", sheet="Dataset S1D"))
dth_B_original<-dth_B_original %>% mutate(FamilyName=recode(Family,
                                            "0"="B1",
                                             "1"="B2"))
                                            
dth_B_original$Group<-3
#subset only the columns we need
dth_B<-dth_B_original[c(10:14)]
#make a new column listing if this is a death (1=YES). to be used to count the total number of worms in a generation
dth_B$Death=1
```

#combine data
```{r}
divdth <-rbind(div_NB,div_B, dth_NB, dth_B)
divdth$FamilyName<-as.factor(divdth$FamilyName)
divdth$Group<-as.factor(divdth$Group)
divdth$Type<-as.factor(ifelse(divdth$H.Gen!=0, 'H', 'T'))
divdth_grouped<-split(divdth,list(divdth$FamilyName,divdth$Type))
```

##Determine death rate per lineage
```{r}
dthNum<-function(x){
  colSums(x[6])}
AliveNum<-function(x){
  colSums(x[6]==0)}
NumDead<-sapply(divdth_grouped,dthNum)
NumAlive<-sapply(divdth_grouped,AliveNum)
DeathCount=as.data.frame(t(rbind(NumDead,NumAlive)))
DeathCount$FamilyName <- as.factor(c('B1', 'B2', 'F', 'S','B1', 'B2', 'F', 'S'))
DeathCount$Type <- as.factor(c('Head', 'Head', 'Head', 'Head', 'Tail', 'Tail', 'Tail', 'Tail'))
DeathCount$Percent<-100*DeathCount$NumDead/(DeathCount$NumDead+DeathCount$NumAlive)
DeathCount$Total<-DeathCount$NumDead+DeathCount$NumAlive
Death_plot<-ggplot(DeathCount, aes(x=(reorder(FamilyName, Percent)), y=Percent, fill=Type))+
  geom_bar(position='dodge', stat="identity", color='black')+
   scale_fill_grey()+
  xlab("Family")+
  ylab('% Death')+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black")
        )
Death_plot

#Perform pairwise Fisher's exact tests for head and tail offspring separately
library(RVAideMemoire)

#stats for heads
cont.tab_heads<-as.matrix(DeathCount%>% filter(Type=="Head") %>% select(NumDead,NumAlive))
Death_ptable_heads<-fisher.multcomp(cont.tab_heads, p.method = "bonferroni")
Death_ptable_heads

#stats for tails
cont.tab_tails<-as.matrix(DeathCount%>% filter(Type=="Tail") %>% select(NumDead,NumAlive))
Death_ptable_tails<-fisher.multcomp(cont.tab_tails, p.method = "bonferroni")
Death_ptable_tails

```
