#### PROJECT: Mimulus cardinalis gene expression (Preston et al. manuscript from 2021)
#### PURPOSE: Visualize reaction norms and violin plots for thermal DEGs share across all three groups
############# Data are from "DEGs_HeatUP_Shared_July21.csv" and "DEGs_HeatDOWN_Shared_July21.csv"
#### AUTHOR: Jill Preston
#### DATE LAST MODIFIED: 12-Nov-2021

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

#set working directory
setwd("/")


### Load in required libraries

library(sm)
library(vioplot)

#************************************************************************
# 2. VIOLIN PLOTS AND REACTION NORMS - UPREGULATED GENES SHARED ACROSS GROUPS
#************************************************************************

Data4 <-read.table("DEGs_HeatUP_Shared_July21.csv", header=T, sep=",")
Data4[1,]

#create vector with unique gene names
gene.names <- unique(Data4$Gene)
length(gene.names)

#create data frame for N_2010
DataN_2010 <- Data4[Data4$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 <- Data4[Data4$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 <- Data4[Data4$Genotype=="S_2017",]
dim(DataS_2017)
class(DataS_2017)
head(DataS_2017)
tail(DataS_2017)
summary(DataS_2017)
unique(DataS_2017$Genotype)

#plot expression change in logarithmic scale for N_2010
color.N_2010 <- "cornflowerblue"
summary(DataN_2010$Norm_expression)
sum(DataN_2010$Norm_expression<=0) #how many zeros?
length(DataN_2010$Norm_expression) #total number of gene expression data points

plot(DataN_2010$Temp=="T40", log(DataN_2010$Norm_expression+1), col=color.N_2010, xlim=c(-0.8,1.8), ylim=c(0,10), xaxt="n", xlab="", ylab="Log( Expression +1 )", 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"))))
abline(h=0:10, col="gray90")
for(i in 1:length(gene.names)){
	points(DataN_2010$Temp[DataN_2010$Gene==gene.names[i]]=="T40", log(DataN_2010$Norm_expression[DataN_2010$Gene==gene.names[i]]+1), type="l", col=color.N_2010)
}
#add violin plots
vioplot(log(DataN_2010$Norm_expression+1)[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(log(DataN_2010$Norm_expression+1)[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(log(DataN_2010$Norm_expression+1)[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(log(DataN_2010$Norm_expression+1)[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 logarithmic scale for N_2017
color.N_2017 <- "purple"
summary(DataN_2017$Norm_expression)
sum(DataN_2017$Norm_expression<=0) #how many zeros?
length(DataN_2017$Norm_expression) #total number of gene expression data points

plot(DataN_2017$Temp=="T40", log(DataN_2017$Norm_expression+1), col=color.N_2017, xlim=c(-0.8,1.8), ylim=c(0,10), xaxt="n", xlab="", ylab="Log( Expression +1 )", 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"))))
abline(h=0:10, col="gray90")
for(i in 1:length(gene.names)){
	points(DataN_2017$Temp[DataN_2017$Gene==gene.names[i]]=="T40", log(DataN_2017$Norm_expression[DataN_2017$Gene==gene.names[i]]+1), type="l", col=color.N_2017)
}

#add violin plots
vioplot(log(DataN_2017$Norm_expression+1)[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(log(DataN_2017$Norm_expression+1)[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(log(DataN_2017$Norm_expression+1)[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(log(DataN_2017$Norm_expression+1)[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 logarithmic scale for S_2017
color.S_2017 <- "red"
summary(DataS_2017$Norm_expression)
sum(DataS_2017$Norm_expression<=0) #how many zeros?
length(DataS_2017$Norm_expression) #total number of gene expression data points

plot(DataS_2017$Temp=="T40", log(DataS_2017$Norm_expression+1), col=color.S_2017, xlim=c(-0.8,1.8), ylim=c(0,10), xaxt="n", xlab="", ylab="Log( Expression +1 )", 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"))))
abline(h=0:10, col="gray90")
for(i in 1:length(gene.names)){
	points(DataS_2017$Temp[DataS_2017$Gene==gene.names[i]]=="T40", log(DataS_2017$Norm_expression[DataS_2017$Gene==gene.names[i]]+1), type="l", col=color.S_2017)
}

#add violin plots
vioplot(log(DataS_2017$Norm_expression+1)[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(log(DataS_2017$Norm_expression+1)[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(log(DataS_2017$Norm_expression+1)[DataS_2017$Temp=="T40"], at=1.2, add=T, col=color.S_2017, border=color.S_2017, rectCol="gray", lineCol="black", colMed="black", side="right", wex=1.3)
vioplot(log(DataS_2017$Norm_expression+1)[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. VIOLIN PLOTS AND REACTION NORMS - DOWNREGULATED GENES SHARED ACROSS GROUPS
#************************************************************************

Data5 <-read.table("DEGs_HeatDOWN_Shared_July21.csv", header=T, sep=",")
Data5[1,]

#create vector with unique gene names
gene.names <- unique(Data5$Gene)
length(gene.names)

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

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

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

#plot expression change in logarithmic scale for N_2010
color.N_2010 <- "cornflowerblue"
summary(DataN_2010_down$Norm_expression)
sum(DataN_2010_down$Norm_expression<=0) #how many zeros?
length(DataN_2010_down$Norm_expression) #total number of gene expression data points

range(log(DataN_2010_down$Norm_expression))

plot(DataN_2010_down$Temp=="T40", log(DataN_2010_down$Norm_expression+1), col=color.N_2010, xlim=c(-0.8,1.8), ylim=c(0,12), xaxt="n", xlab="", ylab="Log( Expression +1 )", 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"))))
abline(h=0:12, col="gray90")
for(i in 1:length(gene.names)){
	points(DataN_2010_down$Temp[DataN_2010_down$Gene==gene.names[i]]=="T40", log(DataN_2010_down$Norm_expression[DataN_2010_down$Gene==gene.names[i]]+1), type="l", col=color.N_2010)
}

#add violin plots
vioplot(log(DataN_2010_down$Norm_expression+1)[DataN_2010_down$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(log(DataN_2010_down$Norm_expression+1)[DataN_2010_down$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(log(DataN_2010_down$Norm_expression+1)[DataN_2010_down$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(log(DataN_2010_down$Norm_expression+1)[DataN_2010_down$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 logarithmic scale for N_2017
color.N_2017 <- "purple"
summary(DataN_2017_down$Norm_expression)
sum(DataN_2017_down$Norm_expression<=0) #how many zeros?
length(DataN_2017_down$Norm_expression) #total number of gene expression data points

plot(DataN_2017_down$Temp=="T40", log(DataN_2017_down$Norm_expression+1), col=color.N_2017, xlim=c(-0.8,1.8), ylim=c(0,12), xaxt="n", xlab="", ylab="Log( Expression +1 )", 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"))))
abline(h=0:12, col="gray90")
for(i in 1:length(gene.names)){
	points(DataN_2017_down$Temp[DataN_2017_down$Gene==gene.names[i]]=="T40", log(DataN_2017_down$Norm_expression[DataN_2017_down$Gene==gene.names[i]]+1), type="l", col=color.N_2017)
}

#add violin plots
vioplot(log(DataN_2017_down$Norm_expression+1)[DataN_2017_down$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(log(DataN_2017_down$Norm_expression+1)[DataN_2017_down$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(log(DataN_2017_down$Norm_expression+1)[DataN_2017_down$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(log(DataN_2017_down$Norm_expression+1)[DataN_2017_down$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 logarithmic scale for S_2017
color.S_2017 <- "red"
summary(DataS_2017_down$Norm_expression)
sum(DataS_2017_down$Norm_expression<=0) #how many zeros?
length(DataS_2017_down$Norm_expression) #total number of gene expression data points

plot(DataS_2017_down$Temp=="T40", log(DataS_2017_down$Norm_expression+1), col=color.S_2017, xlim=c(-0.8,1.8), ylim=c(0,12), xaxt="n", xlab="", ylab="Log( Expression +1 )", 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"))))
abline(h=0:12, col="gray90")
for(i in 1:length(gene.names)){
	points(DataS_2017_down$Temp[DataS_2017_down$Gene==gene.names[i]]=="T40", log(DataS_2017_down$Norm_expression[DataS_2017_down$Gene==gene.names[i]]+1), type="l", col=color.S_2017)
}

#add violin plots
vioplot(log(DataS_2017_down$Norm_expression+1)[DataS_2017_down$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(log(DataS_2017_down$Norm_expression+1)[DataS_2017_down$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(log(DataS_2017_down$Norm_expression+1)[DataS_2017_down$Temp=="T40"], at=1.2, add=T, col=color.S_2017, border=color.S_2017, rectCol="gray", lineCol="black", colMed="black", side="right", wex=1.3)
vioplot(log(DataS_2017_down$Norm_expression+1)[DataS_2017_down$Temp=="T40"], at=1.2, add=T, col="transparent", border="transparent", rectCol="gray", lineCol="black", colMed="black", side="both", wex=1.3, lwd=1)