# This script was written for data analysis reported in the manuscript:
# Freshwater trematodes differ from marine trematodes in patterns connected with 
#   Division of Labor
# By Allison T Neal, Moira Stettner, Renytzabelle Ortega-Cotto, Daniel Dieringer,
#   and Lydia Reed
# Submitted to PeerJ in October 2023
# The script was written by Allison Neal and is provided to allow easy replication
#   and critical evaluation of the analysis performed.  The script author does not
#   claim to be an efficient or expert R coder, so please excuse any inefficiencies 
#   in the code.

# This script attempts to summarize what is different about the trematode 
#   colonies that have rediae with large pharynges 

### First, import the cleaned up data:

source('/Users/aneal1/OneDrive - Norwich University/2_Research/1_Trematode Research/Div of Labor/Analysis/Redial Sizes/DOLVTsizesImportCleanUp.R')

# Approximate volume of redia assuming a cylinder shape: pi*r^2*h
#   for all measurements (including damaged rediae, etc.)
data$vol<-pi*(data$width/2)^2*data$length
data$log.vol<-log(data$vol,10)
#   for reliable measurements only
data$volR<-pi*(data$widthR/2)^2*data$lengthR
data$log.volR<-log(data$volR,10)

# Calculate approximate volume of pharynx assuming a spherical pharynx: 4/3*pi*r^3 
#   for all measurements
data$pvol<-4/3*pi*(data$pharynx/2)^3  
data$log.pvol<-log(data$pvol,10)
  # for reliable measurements only
data$pvolR<-4/3*pi*(data$pharR/2)^3
data$log.pvolR<-log(data$pvolR,10)

# Calculate relative volume
data$rel.pvolR<-data$pvolR/data$volR

###############################################################################
# Identify outliers for pharynx size

# Get a list of all snail numbers for anlaysis
snails<-levels(as.factor(data$snail))

# Create an empty data frame to store information on which rediae have unusually large pharynges
large.pharynx <- data.frame(matrix(data = NA, nrow = length(snails), ncol = 6))
colnames(large.pharynx)<-c("snail","rel.range", "rel.IQR","n.total","n.large","list.large")

# Populate data frame with snail numbers
large.pharynx$snail <- snails

# For each snail number,identify the outliers as more than 1.5x the interquartile range
#   higher than the 3rd quartile
for(i in 1:length(snails)){
  temp <- data[data$snail==snails[i],]
  # Summarize pharynx sizes
  summary(temp$pharR)[1] -> minimum
  summary(temp$pharR)[6] -> maximum
  summary(temp$pharR)[3] -> med
  summary(temp$pharR)[2] -> Q1
  summary(temp$pharR)[5] -> Q3
  # Calculate interquartile range
  IQR <- Q3 - Q1
  # Calculate the range and IQR relative to the median
  #  (I thought this might be different for infections with outliers, but ended up
  #  not reporting these results)
  large.pharynx$rel.range[i] <- (maximum-minimum)/med
  large.pharynx$rel.IQR[i] <- IQR/med
  # Define cutoff for outliers
  outlier.min<-Q3 + 1.5*IQR
  # Identify and quantify outliers
  outliers <- temp$sample[which(temp$pharR > outlier.min)]
  large.pharynx$n.large[i] <- length(outliers)
  large.pharynx$n.total[i] <- sum(!is.na(temp$pharR))
  large.pharynx$list.large[i] <- paste(outliers,collapse = ', ')
}

######################################################################
### Look at the sizes of rediae with large pharynges relative to other rediae
#   in the infection

#  Create function to plot the outliers with the rest of the data
plot.outliers<-function(snail,dataset,summary){
  temp <- data[data$snail == snail,]
  # Get and clean up a list of outliers for the specified infection
  outliers <- summary$list.large[summary$snail == snail]
  outliers <- unlist(strsplit(outliers,","))
  outliers <- trimws(outliers, which = "both")
  # Create a column in the temporary dataframe for color; outliers are black, all others
  #   are gray
  temp$color<-rep("gray",length(temp$snail))
  temp$color[temp$sample %in% outliers]<-"black"
  # Plot length vs. width of rediae with the colors specifying outliers for pharynx size
  plot(x = temp$length, y = temp$width, col = temp$color, xlab='',ylab='')
  points(temp$length[temp$sample %in% outliers], temp$width[temp$sample %in% outliers]) #brings outliers to top
  mtext(expression(paste("length (", mu, "m)")),side = 1, line = 2.2)
  mtext(expression(paste("width (",mu,"m)")),side = 2, line = 2.2)
}

# Look at graphs for infections with 4 or more outliers
# Some of these were flagged as "interesting" (i.e. having some rediae with 
#   noticeably large pharynges) when the initial photos were taken
plot.outliers("P113",data,large.pharynx) # maybe "interesting", based on looking at photos
plot.outliers("P130",data,large.pharynx) # not "interesting", based on looking at photos
plot.outliers("P136",data,large.pharynx) # maybe "interesting based on looking at photos
plot.outliers("P177",data,large.pharynx) # interesting
plot.outliers("P178",data,large.pharynx) # interesting
plot.outliers("Ph96",data,large.pharynx) # not "interesting" based on looking at photos
plot.outliers("V168",data,large.pharynx) # interesting
plot.outliers("V182",data,large.pharynx) # interesting
plot.outliers("V183",data,large.pharynx) # interesting
plot.outliers("V184",data,large.pharynx) # interesting

# V172 and P153 were also identified as 'interesting' but didn't have many outliers
plot.outliers("V172",data,large.pharynx)
plot.outliers("P153",data,large.pharynx)

#####################################################################
# Calcualte what proportion of pharynx measurements were outliers
large.pharynx$p.outliers<-large.pharynx$n.large/large.pharynx$n.total
# Mark which snail numbers I thought looked interesting (i.e. trematode colonies 
#   in which some rediae have noticeably large pharynges) when looking at the photos
large.pharynx$interesting<-"N"
large.pharynx$interesting[large.pharynx$snail == "P153"] <- "Y"
large.pharynx$interesting[large.pharynx$snail == "P177"] <- "Y"
large.pharynx$interesting[large.pharynx$snail == "P178"] <- "Y"
large.pharynx$interesting[large.pharynx$snail == "V168"] <- "Y"
large.pharynx$interesting[large.pharynx$snail == "V172"] <- "Y"
large.pharynx$interesting[large.pharynx$snail == "V182"] <- "Y"
large.pharynx$interesting[large.pharynx$snail == "V183"] <- "Y"
large.pharynx$interesting[large.pharynx$snail == "V184"] <- "Y"
large.pharynx$interesting[large.pharynx$snail == "P136"] <- "M"
large.pharynx$interesting[large.pharynx$snail == "P113"] <- "M"
large.pharynx$interesting[large.pharynx$snail == "P120"] <- "M"
large.pharynx$interesting[large.pharynx$snail == "V141"] <- "M"

boxplot(large.pharynx$p.outliers~large.pharynx$interesting) # interesting infections have more outliers
summary(large.pharynx$p.outliers[large.pharynx$interesting == "Y"])  # interesting infections have 5-15% outliers;
#   most (75%) have at least 6% outliers (Q1 is 0.0617 for "interesting", 0 for "not interestiong")
summary(large.pharynx$p.outliers[large.pharynx$interesting == "N"])  # 'not interesting' infections
#  have 0-20% outliers, but most (75%) have <3% (Q3 is 0.03448 for "not interesting", 0.09607 for "interesting")


# Figure for publication
pdf(file = "/Users/aneal1/OneDrive - Norwich University/2_Research/1_Trematode Research/Div of Labor/DOL Drafts/Fig6LargePharynx",
    width=5,height=4.5)
layout(matrix(1:4,ncol=2))
par(mar = c(3.3,3.5,0.6,0.2))
# P153
plot.outliers("P153",data,large.pharynx)
title(main = "A", line = -1.2, adj = 0.02,cex.main = 1.5)
# P177
plot.outliers("P177",data,large.pharynx)
title(main = "C", line = -1.2, adj = 0.02,cex.main = 1.5)
# P178
plot.outliers("P178",data,large.pharynx)
title(main = "B", line = -1.2, adj = 0.02,cex.main = 1.5)
# V184
plot.outliers("V184",data,large.pharynx)
title(main = "D", line = -1.2, adj = 0.02,cex.main = 1.5)

dev.off()