## Housekeeping
cat("\014") # Clear console
rm(list=ls()) # Remove all variables

## Libraries
library(ggplot2) # general plotting
library(effsize) # cohen's d
library(pwr) # power analysis
library(lme4) # logistic and linear regression
library(stargazer) # auto latex outputs of lm, glm
library(tidyverse) # data transforming
library(ape) # read.dna command
library(pegas) # haplotype networks

## Set WD + load data
setwd("/Users/ba_whyte/Library/CloudStorage/GoogleDrive-ba.whyte@berkeley.edu/My Drive/0 - R/1 - Trematodes/ABC assay (code)")
data <- read.csv("ABC_ABX_Bites.csv")
data$PAIR <- paste(data$COLONY, data$ENEMY, sep = "")
data$ID <- paste(data$TRIAL, data$PAIR, sep = "-")
data$PAIR <- factor(data$PAIR, levels = c("AX","BX","AB","BA","AC","BC"))
data$TYPE <- factor(data$TYPE, levels = c("BL-Hirh","MB-Hirh","BL-Acha"))

## PAIR colors
pAB <- "#FFA630"
pAC <- "#0699EF"
pAX <- "#47E5BC"

## Figure 4: Haplotype network of H. rhigedana colonies in the conspecific assays ------------------------------------------------------------------------------

## Haplotype networks
dna_file <- "/Users/ba_whyte/Library/CloudStorage/GoogleDrive-ba.whyte@berkeley.edu/My Drive/0 - R/1 - Trematodes/ABC assay (code)/DNA/ManuallyTrimmedMAFFTNoOutgroup.fasta"
SEQS <- read.dna(dna_file, format="fasta")
region <- c("BL","BL","MB","BL","BL","MB","BL","BL","MB","BL","BL","MB")
h <- haplotype(SEQS, strict = TRUE)
d <- dist.dna(h, "TN93", as.matrix=TRUE) # Tamura & Nei (1993)
TCS <- haploNet(h, getProb = TRUE) # getProb TRUE = Templeton et al. (1992)
hf <- haploFreq(SEQS, fac = region, haplo = h)
setHaploNetOptions(pie.colors.function = c("#ffa630","#0699ef"), labels=FALSE)
plot(TCS, pie = hf, show.mutation = 0, fast = FALSE, size = 10, cex = 0.8)

## Table 1 - Statistical results of aggression assays --------------------------------------------------------------------------

## T-TEST by TYPE (BL-Hirh, MB-Hirh, BL-Acha)
t1 <- t.test(BITES.C ~ TYPE, data=data[data$ASSAY == "ABC",], type = "paired")
g1 <- as.numeric(data[data$ASSAY == "ABC" & data$TYPE == "BL-Hirh",]$BITES.C)
g2 <- as.numeric(data[data$ASSAY == "ABC" & data$TYPE == "MB-Hirh",]$BITES.C)
cohen.d(g1, g2, hedges.correction=TRUE, paired = TRUE) # effect size
pwr.t.test(n = 8, d = 1.4, sig.level = 0.05, power = NULL, type = "paired") # if n = 8, you still have strong power (0.92)
pwr.t.test(n = NULL, d = 1.4, sig.level = 0.05, power = 0.8, type = "paired") # for desired power, at least n > 6 is needed

## Table 2 - Regression model comparisons ----------------------------------------------------

## Upload data + z-score covariates
df <- read.csv("ABC_ABX_LMM_2.csv")
df <- df %>%
  mutate(zGERMS = scale(GERMS)) %>% # (z) Germinal balls
  mutate(zCOI.N = scale(COI.N)) %>% # (z) Hamming distance (just untransformed bp diffs)
  mutate(zCOI.K80 = scale(COI.K80)) %>% # (z) Kimura's 2-parameter distance 1980
  mutate(zCOI.TN93 = scale(COI.TN93)) # (z) Tamura & Nei 1993
## Logistic regression models
m1 <- glm(BITE.C ~ 1 + ENEMY + zGERMS + zCOI.TN93, data = df[df$ASSAY == "ABC",], family = "binomial") # ABC
m2 <- glmer(BITE.C ~ 1 + (1|COLONY) + ENEMY + zGERMS + zCOI.TN93, data = df[df$ASSAY == "ABC",], family = "binomial") # ABC + colony effect
m3 <- glm(BITE.C ~ 1 + ENEMY + zGERMS, data = df[df$ASSAY == "ABX",], family = "binomial") # ABX
m4 <- glmer(BITE.C ~ 1 + (1|COLONY) + ENEMY + zGERMS, data = df[df$ASSAY == "ABX",], family = "binomial") # ABX + colony effect
stargazer(m1,m2,m3,m4)