#### PROJECT: Mimulus cardinalis gene expression (Preston et al. manuscript from 2021)
#### PURPOSE: Make plots to visualize and compare RGR plasticity across genotypes and collection years
############# Data are from "RGR_families.csv" and "TPC data_cleaned.csv"
#### AUTHOR: Jill Preston and Seema Sheth
#### DATE LAST MODIFIED: 12-Nov-2021

#************************************************************************
# 1. PREPARING THE DATA
#************************************************************************

#set working directory
setwd("/")

Data <-read.table("RGR_families.csv", header=T, sep=",")
Data[1,]
attach(Data)

####Code for pairwise plots with lines linking the same gene.
####Need to first subset the data for N_2010 = Genotype (and then the same for N_2017 and S_2017)
####For each subset: x = Temp; y = Norm_expression; lines between = Gene; lines for average per Temp)

library(sm)
library(vioplot)

#create vector with unique gene names
family.names <- unique(Data$Family)
length(family.names)

#create data frame for N_2010
DataN_2010 <- Data[Data$Genotype=="N_2010",]
dim(DataN_2010)
class(DataN_2010)
head(DataN_2010)
tail(DataN_2010)
summary(DataN_2010)
unique(DataN_2010$Genotype)

#create data frame for N_2017
DataN_2017 <- Data[Data$Genotype=="N_2017",]
dim(DataN_2017)
class(DataN_2017)
head(DataN_2017)
tail(DataN_2017)
summary(DataN_2017)
unique(DataN_2017$Genotype)

#create data frame for S_2017
DataS_2017 <- Data[Data$Genotype=="S_2017",]
dim(DataS_2017)
class(DataS_2017)
head(DataS_2017)
tail(DataS_2017)
summary(DataS_2017)
unique(DataS_2017$Genotype)

#************************************************************************
# 2. PLOT EXPRESSION IN LINEAR SCALE AND MAKE VIOLIN PLOTS
#************************************************************************

#plot expression change in linear scale
color.N_2010 <- "cornflowerblue"
plot(DataN_2010$Temp=="T40", DataN_2010$RGR, col=color.N_2010, xlim=c(-0.8,1.8), ylim=c(0,0.7), xaxt="n", xlab="", ylab="RGR", pch=19, cex=0.5)
#axis(1, at=c(0,1), labels=c("T 20", "T 40"))
axis(1, at=c(0,1), labels=c(expression(paste(20*degree, "C")), expression(paste(40*degree, "C"))))
for(i in 1:length(family.names)){
	points(DataN_2010$Temp[DataN_2010$Family==family.names[i]]=="T40", DataN_2010$RGR[DataN_2010$Family==family.names[i]], type="l", col=color.N_2010)
}

#add violin plots
vioplot(DataN_2010$RGR[DataN_2010$Temp=="T20"], at=-0.2, add=T, col=color.N_2010, border=color.N_2010, rectCol="gray", lineCol="black", colMed="black", side="left", wex=1.3)
vioplot(DataN_2010$RGR[DataN_2010$Temp=="T20"], at=-0.2, add=T, col="transparent", border="transparent", rectCol="gray", lineCol="black", colMed="black", side="both", wex=1.3, lwd=1)
vioplot(DataN_2010$RGR[DataN_2010$Temp=="T40"], at=1.2, add=T, col=color.N_2010, border=color.N_2010, rectCol="gray", lineCol="black", colMed="black", side="right", wex=1.3)
vioplot(DataN_2010$RGR[DataN_2010$Temp=="T40"], at=1.2, add=T, col="transparent", border="transparent", rectCol="gray", lineCol="black", colMed="black", side="both", wex=1.3, lwd=1)

######################
#plot expression change in linear scale
color.N_2017 <- "purple"
plot(DataN_2017$Temp=="T40", DataN_2017$RGR, col=color.N_2017, xlim=c(-0.8,1.8), ylim=c(0,0.7), xaxt="n", xlab="", ylab="RGR", pch=19, cex=0.5)
#axis(1, at=c(0,1), labels=c("T 20", "T 40"))
axis(1, at=c(0,1), labels=c(expression(paste(20*degree, "C")), expression(paste(40*degree, "C"))))
for(i in 1:length(family.names)){
	points(DataN_2017$Temp[DataN_2017$Family==family.names[i]]=="T40", DataN_2017$RGR[DataN_2017$Family==family.names[i]], type="l", col=color.N_2017)
}

#add violin plots
vioplot(DataN_2017$RGR[DataN_2017$Temp=="T20"], at=-0.2, add=T, col=color.N_2017, border=color.N_2017, rectCol="gray", lineCol="black", colMed="black", side="left", wex=1.3)
vioplot(DataN_2017$RGR[DataN_2017$Temp=="T20"], at=-0.2, add=T, col="transparent", border="transparent", rectCol="gray", lineCol="black", colMed="black", side="both", wex=1.3, lwd=1)
vioplot(DataN_2017$RGR[DataN_2017$Temp=="T40"], at=1.2, add=T, col=color.N_2017, border=color.N_2017, rectCol="gray", lineCol="black", colMed="black", side="right", wex=1.3)
vioplot(DataN_2017$RGR[DataN_2017$Temp=="T40"], at=1.2, add=T, col="transparent", border="transparent", rectCol="gray", lineCol="black", colMed="black", side="both", wex=1.3, lwd=1)

######################
#plot expression change in linear scale
color.S_2017 <- "red"
plot(DataS_2017$Temp=="T40", DataS_2017$RGR, col=color.S_2017, xlim=c(-0.8,1.8), ylim=c(0,0.7), xaxt="n", xlab="", ylab="RGR", pch=19, cex=0.5)
#axis(1, at=c(0,1), labels=c("T 20", "T 40"))
axis(1, at=c(0,1), labels=c(expression(paste(20*degree, "C")), expression(paste(40*degree, "C"))))
for(i in 1:length(family.names)){
	points(DataS_2017$Temp[DataS_2017$Family==family.names[i]]=="T40", DataS_2017$RGR[DataS_2017$Family==family.names[i]], type="l", col=color.S_2017)
}

#add violin plots
vioplot(DataS_2017$RGR[DataS_2017$Temp=="T20"], at=-0.2, add=T, col=color.S_2017, border=color.S_2017, rectCol="gray", lineCol="black", colMed="black", side="left", wex=1.3)
vioplot(DataS_2017$RGR[DataS_2017$Temp=="T20"], at=-0.2, add=T, col="transparent", border="transparent", rectCol="gray", lineCol="black", colMed="black", side="both", wex=1.3, lwd=1)
vioplot(DataS_2017$RGR[DataS_2017$Temp=="T40"], at=1.2, add=T, col=color.S_2017, border=color.N_2017, rectCol="gray", lineCol="black", colMed="black", side="right", wex=1.3)
vioplot(DataS_2017$RGR[DataS_2017$Temp=="T40"], at=1.2, add=T, col="transparent", border="transparent", rectCol="gray", lineCol="black", colMed="black", side="both", wex=1.3, lwd=1)

#************************************************************************
# 3. ANOVA AND TUKEY TEST
#************************************************************************

# load packages
library(tidyverse)
library(multcomp)
library(cowplot)

# read in data, downloaded from Wooliver et al. 2020 Github repo: https://github.com/rwoolive/Cardinalis_TPC_evolution/blob/master/Processed%20data/TPC%20data_cleaned.csv
dat <- read.csv("TPC data_cleaned.csv") %>% filter((daytimeTemp==40|daytimeTemp==20) & (Pop=="N2"|(Pop=="S2" & Year==2017))&!is.na(RGR)) %>% dplyr::select(daytimeTemp,Pop,Year,Group,Group.ord,RGR) 

# make grouping_var column corresponding to population x year x temperature
dat = dat %>% mutate(grouping_var=as.factor(paste(Pop,Year,daytimeTemp,sep="-")),Pop=as.factor(Pop),Year=as.factor(Year),daytimeTemp=as.factor(daytimeTemp),Group=as.factor(Group),Group.ord=as.factor(Group.ord))

# inspect data
str(dat)

# perform ANOVA

# proper two-way ANOVA with interaction
analysis1=aov(RGR~daytimeTemp*Group,data=dat)
summary(analysis1)

# one-way ANOVA with interactions embedded into grouping variable
analysis2=aov(RGR~grouping_var,data=dat)
summary(analysis2)

# this website says to run the model as so; can't tell if analysis3 is any different from analysis2 https://cran.r-project.org/web/packages/multcomp/vignettes/multcomp-examples.pdf
analysis3=lm(RGR~grouping_var-1,data=dat)
summary(analysis3)

# perform Tukey test
summary(glht(analysis3, linfct=mcp(grouping_var="Tukey")))

# make violin plot
fig6=ggplot(dat, aes(x=grouping_var, y=RGR,group=grouping_var,fill=grouping_var)) + 
  geom_violin(trim=FALSE,show.legend = FALSE) +
  geom_boxplot(width=0.05,outlier.shape = NA,show.legend = FALSE) +
  ylab("Relative growth rate") +
  xlab("") +
  scale_fill_manual(values=c("cornflowerblue","cornflowerblue","purple","purple","red", "red")) +
  theme_minimal() +
  theme(axis.text=element_text(size=14),axis.title=element_text(size=18))

# save figure
ggsave(filename = "Fig6.tiff",
       plot = fig6,
       #bg = "transparent", 
       width = 11, height = 8.5, units = "in",dpi=600)