---
title: "Population Analysis by Lineage - FS vs bottleneck Families. Related to Fig S7A and B"
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 division data for non-bottleneck lineages and only keep data for FS lineages.
```{r }
data_NB_original<-as.data.frame(read_excel("Dataset_S1.xlsx", sheet="Dataset S1A"))
data_NB<-data_NB_original %>%
  filter((Family==1) | (Family==4 ))
#create group variables (1=FS, 2=Bottleneck)
data_NB$Group<-1
RWT_NB<-data_NB[9:15]
```

# Read in bottleneck lineage data. Add Recode lineage and group names and isolate important variables
```{r }
data_B<-read_excel("Dataset_S1.xlsx", sheet="Dataset S1C")
data_B<-data_B %>% mutate(FamilyName=recode(Family,
                                            "0"="B1",
                                             "1"="B2"))
data_B$Group<-2
RWT_B<-data_B[9:15]
```


## Combine the data into one combined set
```{r }
RWT <-rbind(RWT_NB,RWT_B)
RWT$FamilyName<-as.factor(RWT$FamilyName)
RWT$Group<-as.factor(RWT$Group)
RWT_summary_stats<-RWT %>% 
  group_by(FamilyName) %>% 
  summarise(median=median(TtD), n=n())
```

##Plot Violin plots of RWT per Lineage (Fig S8A) and perform Kruskal-Wallis omnibus test following by Dunn's posthoc
```{r }
RWT_plot<-ggplot(RWT, aes(x=reorder(FamilyName, TtD, FUN=median),y=TtD,fill=Group)) +
  geom_violin(color='black') +
  stat_summary(fun=median, geom="point", size=2, color='black') +
  scale_fill_manual(values=c('blue','darkred')) +
  theme(legend.position = "none")+
  xlab("")+
  ylab("Reproductive waiting time (days)")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
RWT_plot

## Kruskal-wallis test to test significance
library(rstatix)
res.kruskal<-RWT %>% kruskal_test (TtD~FamilyName)
RWT_pwc<-RWT %>% 
  dunn_test(TtD~FamilyName, p.adjust.metho ="bonferroni")
RWT_pwc
```


## determining fission vs fragmentation rate (Fig S8B)
The RWT of the first division of a fragmentation event is always >5 days, but the subsequent divisions have RWT <= 5 days. Thus, we first remove the downstream fragmentation events to keep only the first division of the reproductive event. 

```{r }
RWT_fragRate<-RWT%>% filter(TtD>5)
RWT_fragRate$Div<-as.factor(ifelse(RWT_fragRate$Middles>0,"Fragment","Fission"))
fragRate_Freq<-RWT_fragRate %>% group_by(FamilyName,Div) %>% summarize(Freq=n())
fragRate_Freq$FamilyName<- factor(fragRate_Freq$FamilyName, levels=c("F", "S", "B1", "B2"))
DivType_plot<-ggplot(fragRate_Freq, aes(fill=Div, x=FamilyName, y=Freq))+
  geom_bar(position="fill",stat="identity", color='black')+
  scale_fill_grey()+
  xlab("Family")+
  ylab('Proportion of divisions')+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
DivType_plot

#statistical analysis
Frag_table <-xtabs(Freq~FamilyName+Div, data=fragRate_Freq)
library(RVAideMemoire)
Frag_ptable<-fisher.multcomp(Frag_table, p.method = "bonferroni")
Frag_ptable
```
