# 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 aims to compare morphology of small vs. large rediae
# The data analyzed were are notes taken from a visual examination of the rediae
#   listed.  This list of smallest and largest rediae was generated using the script
#   DOLVTlistSmallestLargest.R.

### Import data
morph<-read.csv('/Users/aneal1/OneDrive - Norwich University/2_Research/1_Trematode Research/Div of Labor/Analysis/Redial Sizes/DOLVTmorph.csv')
morph$size<-substr(morph$size,start = 2, stop = 3) # removes quotation marks from "sm" and "lg"

# The data lists whether the five largest and five smallest rediae in each infection
#  have 1) an anterior appendage or collar (app.ant), 2) a posterior appendage
#  (app.post), or 3) embryos/germ balls (germ).

# The possible values are Y (Yes, definitely clear), Y? (it looks like it but I'm not 100% 
#   sure), U (can't really tell at all from the picture), N? (probably not, but I'm not
#  100% sure) and N (no, clearly lacking structure).  I will convert these notes into
#  a numeric scale for how certain I am that the structure is present (4) or absent (0)
# or any level of uncertainty in between (1-3).

morph$app.ant[morph$app.ant=='N']<-0
morph$app.ant[morph$app.ant=='N?']<-1
morph$app.ant[morph$app.ant=='U']<-2
morph$app.ant[morph$app.ant=='Y?']<-3
morph$app.ant[morph$app.ant=='Y']<-4
morph$app.ant<-as.numeric(morph$app.ant) # for some reason this is storing as character data; convert to numeric

morph$app.post[morph$app.post=='N']<-0
morph$app.post[morph$app.post=='N?']<-1
morph$app.post[morph$app.post=='U']<-2
morph$app.post[morph$app.post=='Y?']<-3
morph$app.post[morph$app.post=='Y']<-4
morph$app.post<-as.numeric(morph$app.post)# for some reason this is storing as character data; convert to numeric

morph$germ[morph$germ=='N']<-0
morph$germ[morph$germ=='N?']<-1
morph$germ[morph$germ=='U']<-2
morph$germ[morph$germ=='Y?']<-3
morph$germ[morph$germ=='Y']<-4
morph$germ<-as.numeric(morph$germ)# for some reason this is storing as character data; convert to numeric


# To calculate an average score, I'll first have to get each redia labeled with
#  its infection number.

morph$snail<-rep(NA,length(morph$sample)) # create a new column in morph dataframe
# Extract the snail number from each sample number
for(i in 1:length(morph$sample)){
  if(substr(morph$sample[i],start=2,stop=3)=='Ph'){
    morph$snail[i]<-substr(morph$sample[i],start=2,stop=6) # Physidae samples have an extra charactr in their name
  }else{
    morph$snail[i]<-substr(morph$sample[i],start=2,stop=5) # All other snail numbers are 4 characters long
  }
}

### Summarize notes by calculating average scores for each category (small vs. large rediae,
#    different morphological structures/features)

# List the snails I want to analyze (it excludes some infections with missing/problematic data)
snails<-c('H102','H104','H107','H110','H127','H148','H156','H158','H070',
          'H096','P107','P113','P115','P116','P117','P119','P120','P121',
          'P122','P125','P126','P128','P130','P131','P135','P136','P137','P139',
          'P140','P141','P150','P152','P153','P154','P155','P157','P160','P162',
          'P163','P164','P165','P166','P170','P177','P178','Ph100','Ph107','Ph109',
          'Ph112','Ph114','Ph115','Ph116','Ph126','Ph159','Ph095','Ph096','V141',
          'V159','V168','V172','V182','V183','V184')

# Create empty data frame for storing average presence scores
results<-data.frame(matrix(data=NA, nrow = length(snails), ncol = 9))
colnames(results)<-c('snails','app.ant.sm','app.post.sm','germ.sm','app.ant.lg',
                     'app.post.lg','germ.lg','app.ant.p','app.post.p')
results$snails<-snails

# Calculates average presence scores and tests whether the scores differ b/w small and large rediae
for(i in 1:length(snails)){
  temp<-morph[morph$snail==snails[i],] # data set with only one snail's data
  # averages (median)
  results$app.ant.sm[i]<-median(temp$app.ant[temp$size == "sm"])
  results$app.ant.lg[i]<-median(temp$app.ant[temp$size == "lg"])
  results$app.post.sm[i]<-median(temp$app.post[temp$size == "sm"])
  results$app.post.lg[i]<-median(temp$app.post[temp$size == "lg"])
  results$germ.sm[i]<-median(temp$germ[temp$size == "sm"])
  results$germ.lg[i]<-median(temp$germ[temp$size == "lg"])
  # statistical test
  results$app.ant.p[i] <- wilcox.test(temp$app.ant~temp$size, alternative = "less")[3]
  results$app.post.p[i] <- wilcox.test(temp$app.post~temp$size, alternative = "less")[3]
}
# (This generates warnings because there are ties, which prevents the Wilcox
#   test from calculating an exact p-value)

##########################################################################
# I'm also interested in whether small rediae have evidence of germinal balls/
#   embryos at all.  Let's see what is the maximum score for small rediae for each 
#   infection

# Create empty vector for storing maximum presence scores
repro<-rep(NA,length(snails))

# Calculate and store maximum presence scores
for(i in 1:length(snails)){
  temp <- morph[morph$snail==snails[i],]
  repro[i] <- max(temp$germ[temp$size=='sm'], na.rm=T)
}

# List scores for each snail
cbind(snails,repro)