Analysis_Markdown.rmd

---
title: "Shell Thickness Analysis"
author: "Christopher Powers"
date: "2023-12-04"
output: html_document
---

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

# Thickness analysis

This file documents the analysis for the shell thickness data from both replicates of the OA experiment.
The first chunk loads required packages and defines the file locations.

```{r cars}
library(ggplot2)
library(plotly)
library(dplyr)
library(sjstats)
library(grid)
library(cowplot)
path <-("thicknessData")
file.names<-list.files(path = path, pattern = "csv$") #list all the file names in the folder to get only get 
```

# Water Chemistry Analysis
This chunk runs some basic analysis on the water chemistry data to generate table S2 and generates plot 5C.
This outputs a table in the working directory: Seawater_chemistry_table_Output_All.csv
```{r}

library(seacarb)

path2 = 'waterAnalysis/calibrations'
file.names<-list.files(path = path2, pattern = "csv$") #list all the file names in the folder to get only get the csv files
file.names
pH.cals <- data.frame(matrix(NA, nrow=length(file.names), ncol=3, dimnames=list(file.names,c("Date", "Intercept", "Slope")))) #generate a 3 column dataframe with specific column names

for(i in 1:length(file.names)) { # for every file in list start at the first and run this following function
  Calib.Data <-read.table(file.path(path2,file.names[i]), header=TRUE, sep="\t", na.string="NA", as.is=TRUE) #reads in the data files
  model <-lm(mVTris ~ TTris, data=Calib.Data) #runs a linear regression of mV as a function of temperature
  coe <- coef(model) #extracts the coeffecients
  summary(model)$r.squared
  plot(Calib.Data$mVTris, Calib.Data$TTris)
  pH.cals[i,2:3] <- coe #inserts them in the dataframe
  pH.cals[i,1] <- substr(file.names[i],1,8) #stores the file name in the Date column
}

colnames(pH.cals) <- c("Calib.Date",  "Intercept",  "Slope") #rename columns
pH.cals #view data
#constants for use in pH calculation 
R <- 8.31447215 #gas constant in J mol-1 K-1 
F <-96485.339924 #Faraday constant in coulombs mol-1
#read in probe measurements of pH, temperature, and salinity from tanks
daily <- read.csv("waterAnalysis/Daily.tsv", header=TRUE, sep="\t", na.strings="NA") #load data with a header, separated by commas, with NA as NA
SW.chem <- merge(pH.cals, daily, by="Calib.Date")
SW.chem$pH.total = 0
for (i in 1:nrow(SW.chem)) {
  SWS = pHnbs2sws(SW.chem[i,]$pH..meter., SW.chem[i,]$Salinity, SW.chem[i,]$Temperature)
  SW.chem[i,]$pH.total =  pHconv(flag=1, pH=SWS, T=SW.chem[i,]$Temperature, S=SW.chem[i,]$Salinity )
}


SW.chem = SW.chem[!is.na(SW.chem$TA),]
carb.output <- carb(flag=8, var1=SW.chem$pH.total, var2=SW.chem$TA/1000000, S= SW.chem$Salinity, T=SW.chem$Temperature, P=0, Pt=0, Sit=0, pHscale="T", kf="pf", k1k2="l", ks="d") #calculate seawater chemistry parameters using seacarb
carb.output$ALK <- carb.output$ALK*1000000 #convert to µmol kg-1
carb.output$CO2 <- carb.output$CO2*1000000 #convert to µmol kg-1
carb.output$HCO3 <- carb.output$HCO3*1000000 #convert to µmol kg-1
carb.output$CO3 <- carb.output$CO3*1000000 #convert to µmol kg-1
carb.output$DIC <- carb.output$DIC*1000000 #convert to µmol kg-1
carb.output <- carb.output[,-c(1,4,5,8,10:13)] #subset variables of interest
carb.output <- cbind(SW.chem$ID,  SW.chem$Day, SW.chem$Tank, SW.chem$Treatment, carb.output) #combine the sample information with the seacarb output
colnames(carb.output) <- c("ID",  "Day", "Tank","Treatment","Salinity",	"Temperature","pH",	"CO2",	"pCO2","HCO3",	"CO3",	"DIC", "TA",	"Aragonite.Sat" , "Calcite.Sat") #Rename columns to describe contents
View(carb.output)

write.table(carb.output, "waterAnalysis/Seawater_chemistry_table_Output_All.csv", sep=",", row.names = FALSE, quote=FALSE)
library(plyr)
# Total
carbo.melted <- melt(carb.output) #reshape the dataframe to more easily summarize all output parameters
mean.carb.output <-ddply(carbo.melted, .(Tank, variable), summarize, #For each subset of a data frame, apply function then combine results into a data frame.
                         N = length(na.omit(value)), #number of records
                         mean = (mean(value)),       #take the average of the parameters (variables) summarized by treatments
                         sem = (sd(value)/sqrt(N))) #calculate the SEM as the sd/sqrt of the count or data length
mean.carb.output # display mean and sem 

# Summer 2023
carbo.melted <- melt(carb.output[startsWith(carb.output$ID, "OA2"),]) #reshape the dataframe to more easily summarize all output parameters
mean.carb.output <-ddply(carbo.melted, .(Tank, variable), summarize, #For each subset of a data frame, apply function then combine results into a data frame.
                         N = length(na.omit(value)), #number of records
                         mean = (mean(value)),       #take the average of the parameters (variables) summarized by treatments
                         sem = (sd(value)/sqrt(N))) #calculate the SEM as the sd/sqrt of the count or data length
mean.carb.output # display mean and sem 

# Fall 2023
carbo.melted <- melt(carb.output[startsWith(carb.output$ID, "OA3"),]) #reshape the dataframe to more easily summarize all output parameters
mean.carb.output <-ddply(carbo.melted, .(Tank, variable), summarize, #For each subset of a data frame, apply function then combine results into a data frame.
                         N = length(na.omit(value)), #number of records
                         mean = (mean(value)),       #take the average of the parameters (variables) summarized by treatments
                         sem = (sd(value)/sqrt(N))) #calculate the SEM as the sd/sqrt of the count or data length
mean.carb.output # display mean and sem 


t1 = SW.chem[SW.chem$Tank == "T1",]
pH = c(8.1,7.9,7.8,7.7,7.6,7.5,7.4,7.3,7.2,7.1,7.0,6.9,6.8,6.7,6.6)

temp_df = as.data.frame(matrix(nrow=0, ncol=2))
colnames(temp_df) = c("pH", "CalcSat")
for (i in pH) {
  SW.chem_temp = t1
  SW.chem_temp$pH.total = i
  carb.output_temp <- carb(flag=8, var1=SW.chem_temp$pH.total, var2=SW.chem_temp$TA/1000000, S= SW.chem_temp$Salinity, T=SW.chem_temp$Temperature, P=0, Pt=0, Sit=0, pHscale="T", kf="pf", k1k2="l", ks="d") #calculate seawater chemistry parameters using seacarb
  carb.output_temp <- carb.output_temp[,-c(1,4,5,8,10:13)]
  carb.output_temp <- cbind(SW.chem_temp$ID,  SW.chem_temp$Day, SW.chem_temp$Tank, SW.chem_temp$Treatment, carb.output_temp) #combine the sample information with the seacarb output
  colnames(carb.output_temp) <- c("ID",  "Day", "Tank","Treatment","Salinity",	"Temperature","pH",	"CO2",	"pCO2","HCO3",	"CO3",	"DIC", "TA",	"Aragonite.Sat" , "Calcite.Sat") #Rename columns to describe contents
  temp_df = rbind(temp_df, c(mean(as.numeric(carb.output_temp$pH)), mean(as.numeric(carb.output_temp$Calcite.Sat))))
  colnames(temp_df) = c("pH", "CalcSat")
}
temp_df_t1 = temp_df

t2 = SW.chem[SW.chem$Tank == "T2",]
temp_df = as.data.frame(matrix(nrow=0, ncol=2))
colnames(temp_df) = c("pH", "CalcSat")
for (i in pH) {
  SW.chem_temp = t2
  SW.chem_temp$pH.total = i
  carb.output_temp <- carb(flag=8, var1=SW.chem_temp$pH.total, var2=SW.chem_temp$TA/1000000, S= SW.chem_temp$Salinity, T=SW.chem_temp$Temperature, P=0, Pt=0, Sit=0, pHscale="T", kf="pf", k1k2="l", ks="d") #calculate seawater chemistry parameters using seacarb
  carb.output_temp <- carb.output_temp[,-c(1,4,5,8,10:13)]
  carb.output_temp <- cbind(SW.chem_temp$ID,  SW.chem_temp$Day, SW.chem_temp$Tank, SW.chem_temp$Treatment, carb.output_temp) #combine the sample information with the seacarb output
  colnames(carb.output_temp) <- c("ID",  "Day", "Tank","Treatment","Salinity",	"Temperature","pH",	"CO2",	"pCO2","HCO3",	"CO3",	"DIC", "TA",	"Aragonite.Sat" , "Calcite.Sat") #Rename columns to describe contents
  temp_df = rbind(temp_df, c(mean(as.numeric(carb.output_temp$pH)), mean(as.numeric(carb.output_temp$Calcite.Sat))))
  colnames(temp_df) = c("pH", "CalcSat")
}
temp_df_t2 = temp_df

t3 = SW.chem[SW.chem$Tank == "T3",]
temp_df = as.data.frame(matrix(nrow=0, ncol=2))
colnames(temp_df) = c("pH", "CalcSat")
for (i in pH) {
  SW.chem_temp = t3
  SW.chem_temp$pH.total = i
  carb.output_temp <- carb(flag=8, var1=SW.chem_temp$pH.total, var2=SW.chem_temp$TA/1000000, S= SW.chem_temp$Salinity, T=SW.chem_temp$Temperature, P=0, Pt=0, Sit=0, pHscale="T", kf="pf", k1k2="l", ks="d") #calculate seawater chemistry parameters using seacarb
  carb.output_temp <- carb.output_temp[,-c(1,4,5,8,10:13)]
  carb.output_temp <- cbind(SW.chem_temp$ID,  SW.chem_temp$Day, SW.chem_temp$Tank, SW.chem_temp$Treatment, carb.output_temp) #combine the sample information with the seacarb output
  colnames(carb.output_temp) <- c("ID",  "Day", "Tank","Treatment","Salinity",	"Temperature","pH",	"CO2",	"pCO2","HCO3",	"CO3",	"DIC", "TA",	"Aragonite.Sat" , "Calcite.Sat") #Rename columns to describe contents
  temp_df = rbind(temp_df, c(mean(as.numeric(carb.output_temp$pH)), mean(as.numeric(carb.output_temp$Calcite.Sat))))
  colnames(temp_df) = c("pH", "CalcSat")
}
temp_df_t3 = temp_df

t1_subset = carb.output[carb.output$Tank == "T1",]
t2_subset = carb.output[carb.output$Tank == "T2",]
t3_subset = carb.output[carb.output$Tank == "T3",]
col_df = as.data.frame(matrix(c("High", "Moderate", "None")))
plot = ggplot() +  theme_classic() +
  geom_point(data=t1_subset, aes(x=pH, y=Calcite.Sat), colour="firebrick", size=3) + 
  geom_point(data=t2_subset, aes(x=pH, y=Calcite.Sat), colour="goldenrod", size=3) + 
  geom_point(data=t3_subset, aes(x=pH, y=Calcite.Sat), color="dodgerblue3", size=3) + 
  geom_segment(data=t1_subset,x = -Inf, xend=8.05, y = 1, yend=1, colour = "black", size = 1) +
  geom_segment(data=t1_subset,x = -Inf, xend=8.05, y = mean(t1_subset$Calcite.Sat), yend=mean(t1_subset$Calcite.Sat), colour = "firebrick", size = 1, linetype="dashed") +
  geom_segment(data=t1_subset,x = -Inf, xend=8.05, y = mean(t2_subset$Calcite.Sat), yend=mean(t2_subset$Calcite.Sat), colour = "goldenrod", size = 1, linetype="dashed") +
  geom_segment(data=t1_subset,x = -Inf, xend=8.05, y = mean(t3_subset$Calcite.Sat), yend=mean(t3_subset$Calcite.Sat), colour = "dodgerblue3", size = 1, linetype="dashed") +
  theme(text = element_text(size=28)) + ylab("Calcite Saturation") + 
  geom_segment(aes(x = 8.075 , y = 1.05, xend = 8.075, yend = 3),
                  arrow = arrow(length = unit(0.5, "cm")))+ 
  geom_segment(aes(x = 8.075 , y = 0.95, xend = 8.075, yend = -0.1),
                  arrow = arrow(length = unit(0.5, "cm"))) +
  annotate("text", x=8.115, y=0.5, label= "Dissolution", angle = 90, size=5) +
  annotate("text", x=8.115, y=2, label= "Precipitation", angle = 90, size=5)

legend_plot=ggplot(carb.output, aes(x=pH, y=Calcite.Sat, col=Tank)) + geom_line()+
  scale_color_manual(name="Treatment",labels=c(expression("High"~italic(p)*"CO"[2]),expression("Moderate"~italic(p)*"CO"[2]),expression("No"~italic(p)*"CO"[2]), "Saturation"),values=c("firebrick","goldenrod","dodgerblue3"))+ theme_classic() + theme(text = element_text(size=20))
legend=get_legend(legend_plot)
plot_grid(plot, legend, ncol=2, rel_widths=c(0.75,0.25))
```


## Compare all treatments, including the control treatment
This chunk will read in all of the thickness data and associate it with the metadata found in the metadata_supp.tsv file.
```{r pressure, echo=FALSE, results='hide'}
path <-("thicknessAnalysis")
file.names<-list.files(path = path, pattern = "csv$") #list all the file names in the folder to get only get 
# pH 7.1 and 8.1 live comparison
meta <- read.csv("metadata_supp.tsv", header=TRUE, strip.white = TRUE, sep="\t")
proloc = meta
shell <- data.frame(matrix(NA, nrow=0, ncol=6)) #generate a 3 column dataframe with specific column names
colnames(shell) <- c("cell", "pH", "Chamber", "State", "Diameter", "reproductive")

for(i in 1:length(file.names[1:length(file.names)])) { # for every file in list start at the first and run this following function
  names = strsplit(file.names[1:length(file.names)][i], "_", fixed = TRUE)[[1]]
    print(names)
    Data_full <-read.table(file.path(path,file.names[1:length(file.names)][i]), header=TRUE, sep=";", na.string="NA", as.is=TRUE) #reads in the data files
    Data = sample_n(Data_full, 10000)
  if (names[2] == "T1") {
    Data$pH = 7.2
  }
  if (names[2] == "T2") {
    Data$pH = 7.6
  }
  if (names[2] == "T3") {
    Data$pH = 8.1
  }
  if (names[2] == "C1") {
    Data$pH = "Control"
  }
  cell=paste(names[1],names[2],names[3], names[4])
  Data$cell=cell
  Data$Chamber = names[5]
  Data$State = names[4]
  if (cell != "OA2 T3 H5 Dead" ) { # This may be an elphidium shell that snuck into the analysis
    if (cell != "OA3 C1 H4 Live" ) { # 
  Data$Shell.Diameter = meta[meta$ID == cell,]$Shell.Diameter
  if (as.numeric(meta[meta$ID == cell,]$Proloculus.Diameter) > 55) { 
    Data$Repoductive = "Megalosphere"
  } else {
    Data$Repoductive = "Microsphere"
  }
  Data$Batch = names[1]
  
  shell2 <- as.data.frame(Data)
  shell <- rbind(shell, shell2)
    }
}
}

shell$ThickAdjust = (shell$Thickness..mm.*1000)/as.numeric(shell$Shell.Diameter)

shell[shell$pH == "Control",]$State = "Control"
```



# Figure S2
This chunk runs some analysis and generates Figure S2 and the associated stats
```{r}
pH7p2 = shell[shell$pH == 7.2 & shell$State == "Live",]

lm = aov(data = pH7p2, ThickAdjust ~ Repoductive)
summary(lm)
eta = effectsize::eta_squared(lm, partial=F)
eta
meg = pH7p2[pH7p2$Repoductive == "Megalosphere",]
meg$rel = meg$ThickAdjust/sum(meg[meg$Repoductive == "Megalosphere",]$ThickAdjust)
mic = pH7p2[pH7p2$Repoductive == "Microsphere",]
mic$rel = mic$ThickAdjust/sum(mic[mic$Repoductive == "Microsphere",]$ThickAdjust)
pH7p2 = rbind(mic, meg)
p1 = ggplot(pH7p2, aes(x=ThickAdjust, fill=Repoductive)) + 
  geom_histogram(data=(subset(pH7p2,Repoductive == "Megalosphere")), aes(y=..count../nrow(meg)), alpha=0.5, binwidth=0.002) +
  geom_histogram(data=(subset(pH7p2,Repoductive == "Microsphere")), aes(y=..count../nrow(mic)), alpha=0.5, binwidth=0.002) + theme_classic() + ggtitle(expression("A. High"~italic(p)*"CO"[2]~"Live")) + xlab("Normalized Thickness") + ylab("Relative Frequency")
anno1 = paste("Effect Size of Reproductive State:", round(eta[1,2], digits=2))
grob1 <- grobTree(textGrob(anno1, x=0.95,  y=0.9, hjust=1,
  gp=gpar(col="black", fontsize=8, fontface="italic")))
# Plot
p1 = p1 + annotation_custom(grob1)
p1 
ggsave(plot=p1, file="figures/FigS2A.svg", width=5,height=3)

# Dead cells
pH7p2 = shell[shell$pH == 7.2 & shell$State == "Dead",]

lm = aov(data = pH7p2, ThickAdjust ~ Repoductive)
summary(lm)
eta = effectsize::eta_squared(lm, partial=F)
eta
meg = pH7p2[pH7p2$Repoductive == "Megalosphere",]
meg$rel = meg$ThickAdjust/sum(meg[meg$Repoductive == "Megalosphere",]$ThickAdjust)
mic = pH7p2[pH7p2$Repoductive == "Microsphere",]
mic$rel = mic$ThickAdjust/sum(mic[mic$Repoductive == "Microsphere",]$ThickAdjust)
pH7p2 = rbind(mic, meg)
p2 = ggplot(pH7p2, aes(x=ThickAdjust, fill=Repoductive)) + 
  geom_histogram(data=(subset(pH7p2,Repoductive == "Megalosphere")), aes(y=..count../nrow(meg)), alpha=0.5, binwidth=0.002) +
  geom_histogram(data=(subset(pH7p2,Repoductive == "Microsphere")), aes(y=..count../nrow(mic)), alpha=0.5, binwidth=0.002) + theme_classic() + ggtitle(expression("B. High"~italic(p)*"CO"[2]~"Dead")) + xlab("Normalized Thickness") + ylab("Relative Frequency")

anno1 = paste("Effect Size of Reproductive State:", round(eta[1,2], digits=2))
grob1 <- grobTree(textGrob(anno1, x=0.95,  y=0.9, hjust=1,
  gp=gpar(col="black", fontsize=8, fontface="italic")))
# Plot
p2 = p2 + annotation_custom(grob1)
p2
ggsave(plot=p2, file="figures/FigS2B.svg", width=5,height=3)

# Dead cells
pH7p6 = shell[shell$pH == 7.6 & shell$State == "Dead",]

lm = aov(data = pH7p6, ThickAdjust ~ Repoductive)
summary(lm)
eta = effectsize::eta_squared(lm, partial=F)
eta
meg = pH7p6[pH7p6$Repoductive == "Megalosphere",]
meg$rel = meg$ThickAdjust/sum(meg[meg$Repoductive == "Megalosphere",]$ThickAdjust)
mic = pH7p6[pH7p6$Repoductive == "Microsphere",]
mic$rel = mic$ThickAdjust/sum(mic[mic$Repoductive == "Microsphere",]$ThickAdjust)
pH7p2 = rbind(mic, meg)
p3 = ggplot(pH7p6, aes(x=ThickAdjust, fill=Repoductive)) + 
  geom_histogram(data=(subset(pH7p6,Repoductive == "Megalosphere")), aes(y=..count../nrow(meg)), alpha=0.5, binwidth=0.002) +
  geom_histogram(data=(subset(pH7p6,Repoductive == "Microsphere")), aes(y=..count../nrow(mic)), alpha=0.5, binwidth=0.002) + theme_classic() + ggtitle(expression("C. Moderate"~italic(p)*"CO"[2]~"Dead"))+ xlab("Normalized Thickness") + ylab("Relative Frequency")

anno1 = paste("Effect Size of Reproductive State:", round(eta[1,2], digits=2))
grob1 <- grobTree(textGrob(anno1, x=0.95,  y=0.9, hjust=1,
  gp=gpar(col="black", fontsize=8, fontface="italic")))
# Plot
p3 = p3 + annotation_custom(grob1)
p3
ggsave(plot=p3, file="figures/FigS2C.svg", width=5,height=3)

pH8p1 = shell[shell$pH == 8.1 & shell$State == "Live",]

lm = aov(data = pH8p1, ThickAdjust ~ Repoductive)
summary(lm)
eta = effectsize::eta_squared(lm, partial=F)
eta

meg = pH8p1[pH8p1$Repoductive == "Megalosphere",]
meg$rel = meg$ThickAdjust/sum(meg[meg$Repoductive == "Megalosphere",]$ThickAdjust)
mic = pH8p1[pH8p1$Repoductive == "Microsphere",]
mic$rel = mic$ThickAdjust/sum(mic[mic$Repoductive == "Microsphere",]$ThickAdjust)
pH8p1 = rbind(mic, meg)
p4 = ggplot(pH8p1, aes(x=ThickAdjust, fill=Repoductive)) + 
  geom_histogram(data=(subset(pH8p1,Repoductive == "Megalosphere")), aes(y=..count../nrow(meg)), alpha=0.5, binwidth=0.002) +
  geom_histogram(data=(subset(pH8p1,Repoductive == "Microsphere")), aes(y=..count../nrow(mic)), alpha=0.5, binwidth=0.002) + theme_classic() + ggtitle(expression("D. No"~italic(p)*"CO"[2]~"Live")) + xlab("Normalized Thickness") + ylab("Relative Frequency")
anno1 = paste("Effect Size of Reproductive State:", round(eta[1,2], digits=2))
grob1 <- grobTree(textGrob(anno1, x=0.95,  y=0.9, hjust=1,
  gp=gpar(col="black", fontsize=8, fontface="italic")))
# Plot
p4 = p4 + annotation_custom(grob1)
p4
ggsave(plot=p4, file="figures/FigS2D.svg", width=5,height=3)
# Dead cells
pH8p1 = shell[shell$pH == 8.1 & shell$State == "Dead",]

lm = aov(data = pH8p1, ThickAdjust ~ Repoductive)
summary(lm)
eta = effectsize::eta_squared(lm, partial=F)
eta

meg = pH8p1[pH8p1$Repoductive == "Megalosphere",]
meg$rel = meg$ThickAdjust/sum(meg[meg$Repoductive == "Megalosphere",]$ThickAdjust)
mic = pH8p1[pH8p1$Repoductive == "Microsphere",]
mic$rel = mic$ThickAdjust/sum(mic[mic$Repoductive == "Microsphere",]$ThickAdjust)
pH7p2 = rbind(mic, meg)
p5 = ggplot(pH8p1, aes(x=ThickAdjust, fill=Repoductive)) + 
  geom_histogram(data=(subset(pH8p1,Repoductive == "Megalosphere")), aes(y=..count../nrow(meg)), alpha=0.5, binwidth=0.002) +
  geom_histogram(data=(subset(pH8p1,Repoductive == "Microsphere")), aes(y=..count../nrow(mic)), alpha=0.5, binwidth=0.002) + theme_classic() + ggtitle(expression("E. No"~italic(p)*"CO"[2]~"Dead"))+ xlab("Normalized Thickness") + ylab("Relative Frequency")

anno1 = paste("Effect Size of Reproductive State:", round(eta[1,2], digits=2))
grob1 <- grobTree(textGrob(anno1, x=0.95,  y=0.9, hjust=1,
  gp=gpar(col="black", fontsize=8, fontface="italic")))
# Plot
p5 = p5 + annotation_custom(grob1)
p5
ggsave(plot=p5, file="figures/FigS2E.svg", width=5,height=3)
grid.arrange(p1,p2,p3,p4,p5)
```

# Generate Figure S1A
This chunk generates a histogram that shows the two populations of cells. This is Figure S1A
```{r}
ggplot(proloc, aes(x=Proloculus.Diameter)) + 
  geom_histogram(aes(y= after_stat(density)), alpha=0.5, binwidth=5) + geom_density() +
  theme_classic() + xlab("proloculus size (uM)") + ylab("relative frequency") + geom_vline(xintercept=55, color = "red")
```

# Figure 2
Comparison of microspheric and megalospheric cells
Analysis of proloculus and chamber count
```{r}
proloc <- meta

proloc[as.numeric(proloc$Proloculus.Diameter) > 55,]$ReproductiveStage = "Megalosphere"
proloc[as.numeric(proloc$Proloculus.Diameter) <= 55,]$ReproductiveStage = "Microsphere"

proloc[proloc$tank == "Untreated",]$Control = "Control"
proloc$Treatment = "X"
proloc[proloc$Control == "Control",]$Treatment = "Untreated"
proloc[proloc$Control == "Live",]$Treatment = "Live"
proloc[proloc$Control == "Dead",]$Treatment = "Dead"
proloc$Growth = "X"
proloc[proloc$Control == "Control",]$Growth = "Not Live"
proloc[proloc$Control == "Live",]$Growth = "Live"
proloc[proloc$Control == "Dead",]$Growth = "Not Live"
proloc <- proloc[proloc$ID!="OA2 T3 H5 Dead",] # remove putative elphidium test
#proloc <- proloc[proloc$Processed...Y.N.=="Y" | proloc$Processed...Y.N.=="y",]
proloc <- proloc[proloc$ReproductiveStage!="",]
proloc %>% group_by(ReproductiveStage) %>% tally()
proloc$Nchamber = as.numeric(proloc$Nchamber)
proloc$Proloculus.Diameter = as.numeric(proloc$Proloculus.Diameter)
lm = aov(as.numeric(proloc$Nchamber)~proloc$Growth*proloc$ReproductiveStage)
TukeyHSD(lm)
summary(lm)
Treatment_p = summary(lm)[[1]][["Pr(>F)"]][1]
reproductive_p = summary(lm)[[1]][["Pr(>F)"]][2]

effectsize::eta_squared(lm, partial=F)
plot(lm)
g1=ggplot(proloc) + 
  geom_point(aes(x=Nchamber, y=Proloculus.Diameter, color=ReproductiveStage, shape=Treatment), size = 2.5) + 
  geom_density_2d(aes(x=Nchamber, y=Proloculus.Diameter, color=ReproductiveStage),alpha=.5,linemitre = 100) + 
  theme_classic() + 
  ylab(expression(paste("Proloculus Diameter (",mu,"m)"))) + xlab("Number of Chambers") +
  theme(text = element_text(size = 16)) + ggtitle("A.")+ labs(colour = "Life Stage") 
#g1 = g1 + theme(legend.position="none")
ggsave(plot = g1, file="figures/Fig2A.svg", width=4.5, height=3)

lm = aov(as.numeric(proloc$Nchamber)~proloc$ReproductiveStage)
summary(lm)
ncham = summary(lm)[[1]][["Pr(>F)"]][1]
lm = aov(as.numeric(proloc$Shell.Diameter)~proloc$ReproductiveStage)
summary(lm)
diameter = summary(lm)[[1]][["Pr(>F)"]][1]


g2=ggplot(proloc) + 
  geom_point(aes(x=Nchamber, y=as.numeric(Shell.Diameter), color=ReproductiveStage, shape=Treatment), size = 2.5) + 
  geom_smooth(aes(x=Nchamber, y=as.numeric(Shell.Diameter), color=ReproductiveStage),alpha=.5,linemitre = 100, method = "lm") + 
  theme_classic() + 
  ylab(expression(paste("Test Diameter (",mu,"m)"))) + xlab("Number of Chambers") +
  theme(text = element_text(size = 16)) + ggtitle("B.")+ labs(colour = "Life Stage")
anno1 = paste("p-values:", "\nChambers: < 0.001***\nDiameter:", "0.51")
grob1 <- grobTree(textGrob(anno1, x=0.95,  y=0.1, hjust=1,
  gp=gpar(col="black", fontsize=9, fontface="italic")))
g2 = g2 + annotation_custom(grob1)

legend1 = get_legend(g2 + theme(legend.box.margin = margin(0, 0, 0, 12)))
g2 = g2 #+  theme(legend.position="none")
ggsave(plot = g2, file="figures/Fig2B.svg", width=4.5, height=3)

print("dead microsphere")
d_micro = proloc[proloc$Treatment=="Not Growing" & proloc$ReproductiveStage == "Microsphere",]$Nchamber
d_micro = d_micro[!(is.na(d_micro))]
summary(d_micro)
sd(d_micro)/sqrt(length(d_micro))
print("live microsphere")
l_micro = proloc[proloc$Treatment=="Growing" & proloc$ReproductiveStage == "Microsphere",]$Nchamber
l_micro = l_micro[!(is.na(l_micro))]
summary(l_micro)
sd(l_micro)/sqrt(length(l_micro))
print("dead megalosphere")
d_megalo = proloc[proloc$Treatment=="Not Growing" & proloc$ReproductiveStage == "Megalosphere",]$Nchamber
d_megalo = d_megalo[!(is.na(d_megalo))]
summary(d_megalo)
sd(d_megalo)/sqrt(length(d_megalo))
print("live megalosphere")
l_megalo = proloc[proloc$Treatment=="Growing" & proloc$ReproductiveStage == "Megalosphere",]$Nchamber
l_megalo = l_megalo[!(is.na(l_megalo))]
summary(l_megalo)
sd(l_megalo)/sqrt(length(l_megalo))

print("proloculus microsphere")
summary(proloc[proloc$ReproductiveStage == "Microsphere",]$Proloculus.Diameter)
sd(proloc[proloc$ReproductiveStage == "Microsphere",]$Proloculus.Diameter)
print("proloculus megalosphere")
summary(proloc[proloc$ReproductiveStage == "Megalosphere",]$Proloculus.Diameter)
sd(proloc[proloc$ReproductiveStage == "Megalosphere",]$Proloculus.Diameter)

proloc[proloc$Control == "Control",]$Treatment = "Not Live"
proloc[proloc$Control == "Live",]$Treatment = "Live"
proloc[proloc$Control == "Dead",]$Treatment = "Not Live"
lm = aov(as.numeric(proloc$Nchamber)~proloc$Treatment*proloc$ReproductiveStage)
TukeyHSD(lm)
summary(lm)
Treatment_p = summary(lm)[[1]][["Pr(>F)"]][1]
reproductive_p = summary(lm)[[1]][["Pr(>F)"]][2]
g3 = ggplot(proloc) + 
  geom_boxplot(aes(x=ReproductiveStage, y=Nchamber, fill=Treatment), size = 1) + 
  theme_classic()  +
  ylab(expression(paste("Number of Chambers"))) + xlab("Life Stage") +
  theme(text = element_text(size = 16)) + ggtitle("C.") + scale_fill_manual(values = c( "green4", "gray25"))

anno1 = paste("p-values:", "\nTreatment: ",round(Treatment_p, digits=3),"**\nCell Type:", "<0.001***")
grob1 <- grobTree(textGrob(anno1, x=0.1,  y=0.9, hjust=0,
  gp=gpar(col="black", fontsize=9, fontface="italic")))
# Plot
g3 = g3 + annotation_custom(grob1)
g3

top = g1
bottom = plot_grid(g2, g3, ncol=2, rel_widths=c(1,1))
totalPlot = plot_grid(top, bottom, ncol=1)
totalPlot


lm = aov(as.numeric(proloc$Nchamber)~proloc$ReproductiveStage)
summary(lm)
lm = aov(as.numeric(proloc$Shell.Diameter)~proloc$ReproductiveStage)
summary(lm)

ggsave(plot = totalPlot, file="figures/Fig2.svg", width=10, height=7)
```

# Subset the data to only the microspheric data based on the microspheric/megalospheric analysis
This is necessary for the remaining analysis
```{r}
shell= shell[shell$Repoductive == "Microsphere",]
shell_total = shell
```

## Table 2
Comparison between the dead cells and the no-treatment control. 
This tests for the effects of dissolution

Generates a table of effect sizes for both pH and individual variation.
```{r}
shell2 = shell_total
print("Stats for the dead tests across pH levels")
dead = shell2[shell2$State == "Dead",]
dead8.1 = dead[dead$pH == 8.1,]
dead8.1$rel = dead8.1$ThickAdjust/sum(dead8.1[dead8.1$pH == 8.1,]$ThickAdjust)
dead7.6 = dead[dead$pH == 7.6,]
dead7.6$rel = dead7.6$ThickAdjust/sum(dead7.6[dead7.6$pH == 7.6,]$ThickAdjust)
dead7.2 = dead[dead$pH == 7.2,]
dead7.2$rel = dead7.2$ThickAdjust/sum(dead7.2[dead7.2$pH == 7.2,]$ThickAdjust)
control = shell2[shell2$State == "Control",]
control$rel = control$ThickAdjust/sum(control[control$pH == "Control",]$ThickAdjust)


live = shell2[shell2$State == "Live",]
live8.1 = live[live$pH == 8.1,]
live8.1$rel = live8.1$ThickAdjust/sum(live8.1[live8.1$pH == 8.1,]$ThickAdjust)
live7.6 = live[live$pH == 7.6,]
live7.6$rel = live7.6$ThickAdjust/sum(live7.6[live7.6$pH == 7.6,]$ThickAdjust)
live7.2 = live[live$pH == 7.2,]
live7.2$rel = live7.2$ThickAdjust/sum(live7.2[live7.2$pH == 7.2,]$ThickAdjust)

out_table = data.frame(matrix(NA, nrow=3, ncol=2))
out_table_cell = data.frame(matrix(NA, nrow=3, ncol=2))
colnames(out_table)=c("Live","Dead")
rownames(out_table)=c("none","moderate","high")

contrast = rbind(dead7.2, control)
lm = aov(contrast$ThickAdjust~contrast$pH+contrast$cell)
summary(lm)
effectsize::eta_squared(lm, partial=F)
out_table[3,2] = round(effectsize::eta_squared(lm, partial=F)[1,2], 3)
out_table_cell[3,2] = round(effectsize::eta_squared(lm, partial=F)[2,2], 3)

contrast = rbind(dead7.6, control)
lm = aov(contrast$ThickAdjust~contrast$pH+ contrast$cell)
summary(lm)
effectsize::eta_squared(lm, partial=F)
out_table[2,2] = round(effectsize::eta_squared(lm, partial=F)[1,2], 3)
out_table_cell[2,2] = round(effectsize::eta_squared(lm, partial=F)[2,2], 3)

contrast = rbind(dead8.1, control)
lm = aov(contrast$ThickAdjust~contrast$pH+ contrast$cell)
summary(lm)
effectsize::eta_squared(lm, partial=F)
out_table[1,2] = round(effectsize::eta_squared(lm, partial=F)[1,2], 3)
out_table_cell[1,2] = round(effectsize::eta_squared(lm, partial=F)[2,2], 3)

contrast = rbind(live7.2, control)
lm = aov(contrast$ThickAdjust~contrast$pH+ contrast$cell)
summary(lm)
effectsize::eta_squared(lm, partial=F)
out_table[3,1] = round(effectsize::eta_squared(lm, partial=F)[1,2], 3)
out_table_cell[3,1] = round(effectsize::eta_squared(lm, partial=F)[2,2], 3)

contrast = rbind(live7.6, control)
lm = aov(contrast$ThickAdjust~contrast$pH+ contrast$cell)
summary(lm)
effectsize::eta_squared(lm, partial=F)
out_table[2,1] = round(effectsize::eta_squared(lm, partial=F)[1,2], 3)
out_table_cell[2,1] = round(effectsize::eta_squared(lm, partial=F)[2,2], 3)

contrast = rbind(live8.1, control)
lm = aov(contrast$ThickAdjust~contrast$pH+ contrast$cell)
summary(lm)
effectsize::eta_squared(lm, partial=F)
out_table[1,1] = round(effectsize::eta_squared(lm, partial=F)[1,2], 3)
out_table_cell[1,1] = round(effectsize::eta_squared(lm, partial=F)[2,2], 3)
out_table
out_table_cell

```


# Figure 3
Comparison across all chambers analyzed
panels A and B account for comparisons across treatment
Panels C, D, and E account for differences across live vs dead within treatment
```{r}
shell = shell_total[shell_total$pH != "Control",]
semi.r = function(y, x, given){  # this function will compute the semi-partial r
  ryx  = cor(as.numeric(y), x)
  ryg  = cor(as.numeric(y), given)
  rxg  = cor(x, given)
  num  = ryx - (ryg*rxg)
  dnm  = sqrt( (1-rxg^2) )
  sp.r = num/dnm
  return(sp.r)
}
shell2 = shell
shell2$chamberCount=0
shell2[shell2$Chamber == "C1.csv",]$chamberCount = 1
shell2[shell2$Chamber == "C2.csv",]$chamberCount = 2
shell2[shell2$Chamber == "C3.csv",]$chamberCount = 3
shell2[shell2$Chamber == "C4.csv",]$chamberCount = 4
shell2[shell2$Chamber == "C5.csv",]$chamberCount = 5
shell2[shell2$Chamber == "C6.csv",]$chamberCount = 6
shell2[shell2$Chamber == "C7.csv",]$chamberCount = 7
shell2[shell2$Chamber == "C8.csv",]$chamberCount = 8

#shell2$batchCount=0
#shell2[shell2$Chamber == "OA2",]$batchCount = 1
#shell2[shell2$Chamber == "OA3",]$batchCount = 2

print("Stats for the live tests across pH levels")
live = shell2[shell2$State == "Live",]
live8.1 = live[live$pH == 8.1,]
live8.1$rel = live8.1$ThickAdjust/sum(live8.1[live8.1$pH == 8.1,]$ThickAdjust)
live7.6 = live[live$pH == 7.6,]
live7.6$rel = live7.6$ThickAdjust/sum(live7.6[live7.6$pH == 7.6,]$ThickAdjust)
live7.2 = live[live$pH == 7.2,]
live7.2$rel = live7.2$ThickAdjust/sum(live7.2[live7.2$pH == 7.2,]$ThickAdjust)
live = rbind(live8.1, live7.6, live7.2)
phPlot = ggplot(live, aes(x=ThickAdjust, fill=as.factor(pH))) + 
  geom_histogram(data=(subset(live,pH == 8.1)), aes(y=..count../nrow(live8.1)), alpha=0.5, binwidth=0.002) +
  geom_histogram(data=(subset(live,pH == 7.6)), aes(y=..count../nrow(live7.6)), alpha=0.5, binwidth=0.002) +
  geom_histogram(data=(subset(live,pH == 7.2)), aes(y=..count../nrow(live7.2)), alpha=0.5, binwidth=0.002) +
  scale_fill_manual(labels=c(expression("No"~italic(p)*"CO"[2]), expression("Moderate"~italic(p)*"CO"[2]), expression("High"~italic(p)*"CO"[2])), values = c("dodgerblue3", "goldenrod", "firebrick")) +
  theme_classic()+xlab("Normalized Thickness")+ylab("Relative Frequency") + xlim(min(shell2$ThickAdjust),max(shell2$ThickAdjust)) + ggtitle("A.") + 
  theme(legend.position="none", text = element_text(size = 10))
lm = aov(as.numeric(live$ThickAdjust) ~ live$pH + live$cell)
summary(lm)
effectsize::eta_squared(lm, partial=F)
eta = effectsize::eta_squared(lm, partial=F)[1,2]
anno1 = paste("Eta Squared\nTreatment:", round(eta, digits=3))
grob1 <- grobTree(textGrob(anno1, x=0.95,  y=0.9, hjust=1,
  gp=gpar(col="black", fontsize=8, fontface="italic")))
# Plot
p1 = phPlot + annotation_custom(grob1)


print("Stats for the dead tests across pH levels")
dead = shell2[shell2$State == "Dead",]
dead8.1 = dead[dead$pH == 8.1,]
dead8.1$rel = dead8.1$ThickAdjust/sum(dead8.1[dead8.1$pH == 8.1,]$ThickAdjust)
dead7.6 = dead[dead$pH == 7.6,]
dead7.6$rel = dead7.6$ThickAdjust/sum(dead7.6[dead7.6$pH == 7.6,]$ThickAdjust)
dead7.2 = dead[dead$pH == 7.2,]
dead7.2$rel = dead7.2$ThickAdjust/sum(dead7.2[dead7.2$pH == 7.2,]$ThickAdjust)
dead = rbind(dead8.1, dead7.6, dead7.2)
phPlot = ggplot(dead, aes(x=ThickAdjust, fill=as.factor(pH))) + 
  geom_histogram(data=(subset(dead,pH == 8.1)), aes(y=..count../nrow(dead8.1)), alpha=0.5, binwidth=0.002) +
  geom_histogram(data=(subset(dead,pH == 7.6)), aes(y=..count../nrow(dead7.6)), alpha=0.5, binwidth=0.002) +
  geom_histogram(data=(subset(dead,pH == 7.2)), aes(y=..count../nrow(dead7.2)), alpha=0.5, binwidth=0.002) +
  scale_fill_manual(labels=c(expression("No"~italic(p)*"CO"[2]), expression("Moderate"~italic(p)*"CO"[2]), expression("High"~italic(p)*"CO"[2])), values = c("dodgerblue3", "goldenrod", "firebrick")) +
  theme_classic()+xlab("Normalized Thickness")+ylab(NULL) + xlim(min(shell2$ThickAdjust),max(shell2$ThickAdjust)) + ggtitle("B.") + 
  theme(text = element_text(size = 10)) + labs(fill='')
lm = aov(as.numeric(dead$ThickAdjust) ~ dead$pH + dead$cell)
summary(lm)
effectsize::eta_squared(lm, partial=F)
eta = effectsize::eta_squared(lm, partial=F)[1,2]
anno1 = paste("Eta Squared\nTreatment:", round(eta, digits=3))
grob1 <- grobTree(textGrob(anno1, x=0.95,  y=0.9, hjust=1,
  gp=gpar(col="black", fontsize=8, fontface="italic")))
# Plot
p2 = phPlot + annotation_custom(grob1) + theme(legend.position="none")
legend1 = get_legend(phPlot + theme(legend.box.margin = margin(0, 0, 0, 12)))


print("Stats for the physiological effect at pH 7.2")
pH = rbind(live7.2, dead7.2)
lm = aov(pH$ThickAdjust ~ pH$State + pH$cell)
summary(lm)
effectsize::eta_squared(lm, partial=F)
phPlot = ggplot(pH, aes(x=ThickAdjust, fill=as.factor(State))) + 
  geom_histogram(data=subset(pH,State == "Live"), aes(y=..count../nrow(live7.2)), alpha=0.5, binwidth=0.002) +
  geom_histogram(data=subset(pH,State == "Dead"), aes(y=..count../nrow(dead7.2)),alpha=0.5, binwidth=0.002) +
  scale_fill_manual(values = c("green4","gray25")) +
  theme_classic()+xlab("Normalized Thickness")+ylab("Relative Frequency") + xlim(min(shell2$ThickAdjust),max(shell2$ThickAdjust)) + ggtitle(expression(C.)) + 
  theme(legend.position="none", text = element_text(size = 10))
effectsize::eta_squared(lm, partial=F)
eta = effectsize::eta_squared(lm, partial=F)[1,2]
anno1 = paste("Eta Squared\nLive Vs Dead:", round(eta, digits=3))
grob1 <- grobTree(textGrob(anno1, x=0.95,  y=0.9, hjust=1,
  gp=gpar(col="black", fontsize=8, fontface="italic")))
# Plot
p3 = phPlot + annotation_custom(grob1)


print("Stats for the physiological effect at pH 7.6")
pH = rbind(live7.6, dead7.6)
lm = aov(pH$ThickAdjust ~ pH$State + pH$cell)
summary(lm)
effectsize::eta_squared(lm, partial=F)
phPlot = ggplot(pH, aes(x=ThickAdjust, fill=as.factor(State))) + 
  geom_histogram(data=subset(pH,State == "Live"), aes(y=..count../nrow(live7.6)), alpha=0.5, binwidth=0.002) +
  geom_histogram(data=subset(pH,State == "Dead"), aes(y=..count../nrow(dead7.6)),alpha=0.5, binwidth=0.002) +
  scale_fill_manual(values = c("green4","gray25")) +
  theme_classic()+xlab("Normalized Thickness")+ylab(NULL) + xlim(min(shell2$ThickAdjust),max(shell2$ThickAdjust)) + ggtitle(expression(D.)) + 
  theme(legend.position="none", text = element_text(size = 10))
effectsize::eta_squared(lm, partial=F)
eta = effectsize::eta_squared(lm, partial=F)[1,2]
anno1 = paste("Eta Squared\nLive Vs Dead:", round(eta, digits=3))
grob1 <- grobTree(textGrob(anno1, x=0.95,  y=0.90, hjust=1,
  gp=gpar(col="black", fontsize=8, fontface="italic")))
# Plot
p4 = phPlot + annotation_custom(grob1)

print("Stats for the physiological effect at pH 8.1")
pH = rbind(live8.1, dead8.1)
lm = aov(pH$ThickAdjust ~ pH$State + pH$cell)
summary(lm)
effectsize::eta_squared(lm, partial=F)
phPlot = ggplot(pH, aes(x=ThickAdjust, fill=as.factor(State))) + 
  geom_histogram(data=subset(pH,State == "Live"), aes(y=..count../nrow(live8.1)), alpha=0.5, binwidth=0.002) +
  geom_histogram(data=subset(pH,State == "Dead"), aes(y=..count../nrow(dead8.1)),alpha=0.5, binwidth=0.002) +
  scale_fill_manual(values = c("green4", "gray25")) +
  theme_classic()+xlab("Normalized Thickness")+ylab(NULL) + xlim(min(shell2$ThickAdjust),max(shell2$ThickAdjust)) + ggtitle(expression(E.)) + 
  theme(text = element_text(size = 10)) + labs(fill='')
effectsize::eta_squared(lm, partial=F)
eta = effectsize::eta_squared(lm, partial=F)[1,2]
anno1 = paste("Eta Squared\nLive Vs Dead:", round(eta, digits=3))
grob1 <- grobTree(textGrob(anno1, x=0.95,  y=0.90, hjust=1,
  gp=gpar(col="black", fontsize=8, fontface="italic")))
# Plot
p5 = phPlot + annotation_custom(grob1) + theme(legend.position="none")
legend2 = get_legend(phPlot + theme(legend.box.margin = margin(0, 0, 0, 12)))
                     
top_row <- plot_grid(p1, p2, legend1, ncol=3, rel_widths=c(1,1,0.5))
bottom_row <- plot_grid(p3,p4,p5, legend2,ncol=4, rel_widths=c(1,1,1,0.35))

overall = plot_grid(top_row, bottom_row, ncol = 1)
overall
ggsave(plot = overall, file="figures/Fig3.svg", width=8, height=4.5)
```

# Figure 4: Summary of impact across chamber
This plot runs eta squared analysis within each chamber and treatment across the
live and dead to assess the impact of treatment at the chamber level.
```{r}
shell2 <- shell_total[shell_total$pH != "Control",]
thick_by_pH <- function(shell) {
  ph <- unique(shell$pH)
  chamb <- unique(shell$Chamber)
  state <- unique(shell$State)
  effect <- data.frame(matrix(NA, nrow=0, ncol=5))
  for (p in ph) {
    for (c in chamb) {
      C = shell[shell$Chamber == c,]
      C2 = C[C$pH == p,]
      lm = aov(C2$ThickAdjust~C2$State + C2$cell + C2$Batch)
      print(p)
      print(c)
      print(summary(lm))
      print(effectsize::eta_squared(lm, partial=F))
      se = effectsize::eta_squared(lm, partial=F)[1,2]
      subset1 = shell[shell$Chamber == c,]
      subset2 = subset1[subset1$pH == p,]
      subset3 = subset2[subset2$State == "Dead",]
      row = cbind(p, c, se)
      effect = rbind(effect, row)
    }
  }
  return(effect)
}

plot_thick <- function(e, plotname) {
  n = 25000
  
  C1_3d <- read.csv("3dHeatmap/OA2_T1_h1_Live_C1.csv", sep=";")
  C1_3d_sampled = sample(C1_3d$Vertex.Index, n, replace=F, prob=NULL)
  C1_3d = as.data.frame(C1_3d[C1_3d$Vertex.Index %in% C1_3d_sampled,])
  C1_3d$ef = e[3][1,]

  C2_3d <- read.csv("3dHeatmap/OA2_T1_h1_Live_C2.csv", sep=";")
  C2_3d_sampled = sample(C2_3d$Vertex.Index, n, replace=F, prob=NULL)
  C2_3d = C2_3d[C2_3d$Vertex.Index %in% C2_3d_sampled,]
  C2_3d$ef = e[3][2,]

  C3_3d <- read.csv("3dHeatmap/OA2_T1_h1_Live_C3.csv", sep=";")
  C3_3d_sampled = sample(C3_3d$Vertex.Index, n, replace=F, prob=NULL)
  C3_3d = C3_3d[C3_3d$Vertex.Index %in% C3_3d_sampled,]
  C3_3d$ef = e[3][3,]

  C4_3d <- read.csv("3dHeatmap/OA2_T1_h1_Live_C4.csv", sep=";")
  C4_3d_sampled = sample(C4_3d$Vertex.Index, n, replace=F, prob=NULL)
  C4_3d = C4_3d[C4_3d$Vertex.Index %in% C4_3d_sampled,]
  C4_3d$ef = e[3][4,]

  C5_3d <- read.csv("3dHeatmap/OA2_T1_h1_Live_C5.csv", sep=";")
  C5_3d_sampled = sample(C5_3d$Vertex.Index, n, replace=F, prob=NULL)
  C5_3d = C5_3d[C5_3d$Vertex.Index %in% C5_3d_sampled,]
  C5_3d$ef = e[3][5,]

  C6_3d <- read.csv("3dHeatmap/OA2_T1_h1_Live_C6.csv", sep=";")
  C6_3d_sampled = sample(C6_3d$Vertex.Index, n, replace=F, prob=NULL)
  C6_3d = C6_3d[C6_3d$Vertex.Index %in% C6_3d_sampled,]
  C6_3d$ef = e[3][6,]
  
  C7_3d <- read.csv("3dHeatmap/OA2_T1_h1_Live_C7.csv", sep=";")
  C7_3d_sampled = sample(C7_3d$Vertex.Index, n, replace=F, prob=NULL)
  C7_3d = C7_3d[C7_3d$Vertex.Index %in% C7_3d_sampled,]
  C7_3d$ef = e[3][7,]

  C8_3d <- read.csv("3dHeatmap/OA2_T1_h1_Live_C8.csv", sep=";")
  C8_3d_sampled = sample(C8_3d$Vertex.Index, n, replace=F, prob=NULL)
  C8_3d = C8_3d[C8_3d$Vertex.Index %in% C8_3d_sampled,]
  C8_3d$ef = e[3][8,]


  d3 = rbind(C1_3d, C2_3d, C3_3d, C4_3d, C5_3d, C6_3d, C7_3d, C8_3d)
  scene = list(camera = list(eye = list(x = 0.5, y = -2, z = 0.25)))
  fig <- plot_ly(d3, x = ~X.Coordinate..mm., y = ~Y.Coordinate..mm., z = ~Z.Coordinate..mm.,
               mode = "markers",
               marker = list(
    size = 0.5,
    color = ~ ef,
    colorbar = list(title = "Effect Size"),
    colorscale='Viridis',
    cmin = 0,
    cmax = 0.25
  ), type="scatter3d")%>%layout(title = plotname, 
                                scene = list(camera = list(eye = list(x = 0.5, y = -2, z = 0.25)),
                                             xaxis = list(title = '', showgrid = F, visible=F),
                                             yaxis = list(title = '', showgrid = F, visible=F),
                                             zaxis = list(title = '', showgrid = F, visible=F)))
  return(fig)
}

effect = thick_by_pH(shell2)


reticulate::use_miniconda('r-reticulate')
ph72 <- effect[effect$p == 7.2,]
fig <- plot_thick(ph72, "High pCO2")
fig
save_image(fig, "figures/Fig4A.svg")

ph76 <- effect[effect$p == 7.6,]
fig <- plot_thick(ph76, "Moderate pCO2")
fig
save_image(fig, "figures/Fig4B.svg")

ph81 <- effect[effect$p == 8.1,]
fig <- plot_thick(ph81, "No pCO2")
fig
save_image(fig, "figures/Fig4C.svg")

fig <- ggplot(effect, aes(x=p, y=as.numeric(se), fill=p)) + geom_boxplot(show.legend = TRUE) + theme_bw() + theme(text = element_text(size = 36)) + ylim(c(0,0.45)) + xlab(expression()) + ylab("Average Effect Size") + scale_x_discrete(labels = c('','','')) + scale_fill_manual(labels=c(expression("High"~italic(p)*"CO"[2]), expression("Moderate"~italic(p)*"CO"[2]), expression("No"~italic(p)*"CO"[2])), values = c("firebrick", "goldenrod", "dodgerblue3"))+ guides(fill=guide_legend(title="Treatment"))
fig
ggsave(plot=fig, file="figures/Fig4D.svg", width=10,height=6)
```