---
title: "Population Analysis by Lineage - Non-Bottleneck Lineages related to Fig 1"
author: "Danielle Ireland"
date: "2025-06-24"
output:
  html_document: default
  pdf_document: default
---

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

## Load libraries

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

set.seed(228)
```

## Read in non-bottleneck lineage data


```{r }
data_NB<-as.data.frame(read_excel("Dataset_S1.xlsx", sheet="Dataset S1A"))
#select only the columns needed for further analysis
RWT<-data_NB[9:14]
RWT$FamilyName<-as.factor(RWT$FamilyName)

RWT_summary_stats<-RWT %>% 
  group_by(FamilyName) %>% 
  summarise(median=median(TtD), n=n())
```



#Plot Violin plots of RWT per lineage (Fig 1B) and perform Kruskal-Wallis omnibus test following by Dunn's posthoc
```{r RWT}
Fig1B<-ggplot(RWT, aes(x=reorder(FamilyName, TtD, FUN=median),y=TtD))+
  geom_violin(color='black')+
  ylim(0,300)+
  stat_summary(fun=median, geom="point", size=2, color='black')+
  scale_fill_grey()+
  theme(legend.position = "none")+
  xlab("Lineage")+
  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"))

Fig1B

## 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 1C)
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 frag}
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("H","P", "F", "S"))
Fig1C<-ggplot(fragRate_Freq, aes(fill=Div, x=FamilyName, y=Freq))+
  geom_bar(position="fill",stat="identity", color='black')+
  scale_fill_grey()+
  xlab("Lineage")+
  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"))
Fig1C

#statistical analysis: Create contigency table and calculate pairwise Fisher Exact tests
Frag_table <-xtabs(Freq~FamilyName+Div, data=fragRate_Freq)

library(RVAideMemoire)
Frag_ptable<-fisher.multcomp(Frag_table, p.method = "bonferroni")
Frag_ptable
```

## Plot # of middles per fragmentation event per lineage (Fig 1D). 
For visualization in the figure, 3 or more middles are grouped together as 3+. All raw counts are used for the statistical analysis.
```{r}
Trunks=RWT %>% filter(Middles>0)
trunks_Freq<-Trunks %>% group_by(FamilyName,Middles) %>% summarize(Freq=n())
trunks_Freq_combined<- trunks_Freq %>%
    group_by(FamilyName, Middles=replace(Middles, Middles %in% c('3','4','5','6'),'3+')) %>%
        summarise(Freq=sum(Freq))
trunks_Freq_combined$FamilyName<- factor(trunks_Freq_combined$FamilyName, levels=c("H","P", "F", "S"))
trunks_Freq_combined$Middles<- factor(trunks_Freq_combined$Middles)


Fig1D<-ggplot(trunks_Freq_combined, aes(fill=Middles, x=FamilyName, y=Freq))+
  geom_bar(position=position_fill(reverse=TRUE),stat="identity", color='black')+
  scale_fill_grey()+
  xlab("Lineage")+
  ylab('Proportion of fragmentations')+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))

Fig1D

#statistical analysis: Poisson GLM. Use emmeans to perform pairwise comparisons
library(lme4)
model_middles_pois<-glm(Middles~FamilyName, family=poisson(link="log"), data=Trunks)
summary(model_middles_pois)
library(emmeans)
middles_pwc<-pairs(emmeans(model_middles_pois, ~FamilyName), adjust="bonferroni")
middles_pwc
```

