##############################################################
####     Natural Growing Conditions Related to the Size   ####
####       and Nutritional Contents of Wild Rice Seeds    #### 
##############################################################
#### Figures and data analysis R script
#### Created by Jeremy Laurich
#### October 27, 2023
#### Data: Wild_Rice_Nutrition_Data



library(readr)
Consolidated_Data_R <- read_csv("C:/Users/jerem/OneDrive/Desktop/wild_Rice_sulfate/spreadsheets/Consolidated_Data_R_A.csv")
attach(Consolidated_Data_R)
# View(Consolidated_Data_R)

# setting colors and shapes for northeast and south/west ecoregions
colors <- c("#FF9900","#3399FF")
colors <- colors[as.factor(Ecoregion)]
shapes = c(19, 10) 
shapes <- shapes[as.factor(Ecoregion)]




##########################
######## FIGURE 2 ########
##########################

Geomass <- ggplot(data = Consolidated_Data_R,aes(x = Lon, y = Lat, color = Ave_Seed_Mass.mg,size = 1.5)) +
  geom_point()+theme(aspect.ratio=16/10, panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                     panel.background = element_blank())
Geomass + scale_color_gradient(low="light blue", high="black")




##########################
######## FIGURE 3 ########
##########################

dev.new(width = 20, height = 2, units = "cm")
par(mfrow = c(1,4))
boxplot(Ave_Seed_Mass.mg~Type,pch = 16,
        ,cex.names = 1.75, cex.axis = 1.75, cex.lab = 1.75,ylim=c(10,62), ylab="Average Seed Mass (mg)")

detach(Consolidated_Data_R)
library(readr)
L_R_Data <- read_csv("C:/Users/jerem/OneDrive/Desktop/wild_Rice_sulfate/spreadsheets/L_R_Data.csv")

attach(L_R_Data)
#View(L_R_Data)
boxplot(Mass_a~Type_a,ylab = NULL,pch = 16,
        ,cex.names = 1.75, cex.axis = 1.75, cex.lab = 1.75,ylim=c(10,62))
boxplot(Mass_b~Type_b,ylab = NULL,pch = 16,
        ,cex.names = 1.75, cex.axis = 1.75, cex.lab = 1.75,ylim=c(10,62))
boxplot(Mass_c~Type_c,ylab = NULL,pch = 16,
        ,cex.names = 1.75, cex.axis = 1.75, cex.lab = 1.75,ylim=c(10,62))

detach(L_R_Data)
attach(Consolidated_Data_R)




###########################
######## FIGURE 4A ########
###########################

library(ggplot2)
sp2<- ggplot(Consolidated_Data_R, aes(x = (Water_Sulfate.mgL), y = Ave_Seed_Mass.mg, color=Ecoregion, shape=Ecoregion)) + 
  geom_point(size = 4)+
  geom_smooth(method="lm", se = FALSE)+
  theme(panel.grid.major = element_blank(),text = element_text(size = 15),
        panel.grid.minor = element_blank(),panel.background = element_blank())+ scale_x_continuous(trans = 'log10')+
  xlab("Surface Water Sulfate (mg/L)") + ylab("Average Seed Mass (mg)")+ 
  scale_shape_manual(values=c(16, 17))+
  scale_color_manual(values=c('#CCCCCC','#333333')) + geom_text(x=100, y=0.06, label="A")+
  theme(panel.border  = element_rect(size = 0.5, colour = "black", linetype=1, fill = NA))
sp2




###########################
######## FIGURE 4B ########
###########################

attach(Consolidated_Data_R)

boxplot(Ave_Seed_Mass.mg~Sulfate_Quartiles)




##########################
######## FIGURE 5 ########
##########################

par(mfrow = c(1,1))
library(ggplot2)

library(ggpubr)
sp2<- ggplot(Consolidated_Data_R, aes(x = Ave_Seed_Mass.mg, y = Percent_Seed_Starch , color=Ecoregion,shape=Ecoregion)) + 
  geom_point(size = 4)+
  geom_smooth(method="lm", se = FALSE)+
  theme(panel.grid.major = element_blank(),text = element_text(size = 15),
        panel.grid.minor = element_blank(),panel.background = element_blank())+ scale_x_continuous(trans = 'log10')+
  xlab("Average Seed Mass (mg)") + ylab("Seed Starch %")+ 
  scale_shape_manual(values=c(16, 17))+
  scale_color_manual(values=c('#CCCCCC','#333333')) + geom_text(x=100, y=0.06, label="A")+
  theme(panel.border  = element_rect(size = 0.5, colour = "black", linetype=1, fill = NA))
sp2




##########################
######## FIGURE 6 ########
##########################

library(vioplot)
#Elements less than 250
vioplot(Seed_Ca, Seed_Fe, Seed_Mn, Seed_Na,Seed_Zn,names = c( "Ca", "Fe","Mn", "Na","Zn"),ylab = "Concentration(mg/Kg)", cex.names = 1.5, cex.axis = 1.5, cex.lab = 1.5)

#Elements greater than 250
vioplot(Seed_K,Seed_Mg,Seed_P,names = c("K","Mg", "P"),ylab = "Concentration(mg/Kg)",cex.names = 1.5, cex.axis = 1.5, cex.lab = 1.5)




##########################
######## FIGURE 7 ########
##########################

K <- ggplot(Consolidated_Data_R, aes(x = Ave_Seed_Mass.mg, y = Seed_K, color=Ecoregion, pch = Ecoregion)) + 
  geom_point(size = 2)+
  geom_smooth(method = 'lm', se = F)+
  theme(panel.grid.major = element_blank(),text = element_text(size = 12),
        panel.grid.minor = element_blank(),panel.background = element_blank(),axis.line = element_line(colour = "black"),panel.border = element_rect(colour = "black", fill=NA, size=1))+
  xlab(NULL) + ylab("Seed K (mg/kg)")+
  scale_color_manual(values=c('#CCCCCC','#333333'))+ggtitle("Group 1")
K


Mg <- ggplot(Consolidated_Data_R, aes(x = Ave_Seed_Mass.mg, y = Seed_Mg, color=Ecoregion,pch = Ecoregion)) + 
  geom_point(size = 2)+
  geom_smooth(method="lm", se = FALSE)+
  theme(panel.grid.major = element_blank(),text = element_text(size = 12),
        panel.grid.minor = element_blank(),panel.background = element_blank(),axis.line = element_line(colour = "black"),panel.border = element_rect(colour = "black", fill=NA, size=1))+
  xlab("Average Seed Mass (mg)") + ylab("Seed Mg (mg/kg)")+
  scale_color_manual(values=c('#CCCCCC','#333333'))
Mg
Fe <- ggplot(Consolidated_Data_R, aes(x = Ave_Seed_Mass.mg, y = Seed_Fe, color=Ecoregion,pch = Ecoregion)) + 
  geom_point(size = 2)+
  theme(panel.grid.major = element_blank(),text = element_text(size = 12),
        panel.grid.minor = element_blank(),panel.background = element_blank(),axis.line = element_line(colour = "black"),panel.border = element_rect(colour = "black", fill=NA, size=1))+
  xlab(NULL)+ylab("Seed Fe (mg/kg)")+
  scale_color_manual(values=c('#CCCCCC','#333333'))+ggtitle("Group 2")
Fe

Zn <- ggplot(Consolidated_Data_R, aes(x = Ave_Seed_Mass.mg, y = Seed_Zn, color=Ecoregion,pch = Ecoregion)) + 
  geom_point(size = 2)+
  theme(panel.grid.major = element_blank(),text = element_text(size = 12),
        panel.grid.minor = element_blank(),panel.background = element_blank(),axis.line = element_line(colour = "black"),panel.border = element_rect(colour = "black", fill=NA, size=1))+
  xlab("Average Seed Mass (mg)") + ylab("Seed Zn (mg/kg)")+
  scale_color_manual(values=c('#CCCCCC','#333333'))
Zn

Na<- ggplot(Consolidated_Data_R, aes(x = Ave_Seed_Mass.mg, y = Seed_Na, color=Ecoregion,pch = Ecoregion)) + 
  geom_point(size = 2)+
  theme(panel.grid.major = element_blank(),text = element_text(size = 12),
        panel.grid.minor = element_blank(),panel.background = element_blank(),axis.line = element_line(colour = "black"),panel.border = element_rect(colour = "black", fill=NA, size=1))+
  geom_smooth(method="lm", se = FALSE)+
  xlab(NULL) + ylab("Seed Na (mg/kg)")+
  scale_color_manual(values=c('#CCCCCC','#333333'))+ggtitle("Group 3")
Na

Si<- ggplot(Consolidated_Data_R, aes(x = Ave_Seed_Mass.mg, y = Seed_Si, color=Ecoregion,pch = Ecoregion)) + 
  geom_point(size = 2)+
  theme(panel.grid.major = element_blank(),text = element_text(size = 12),
        panel.grid.minor = element_blank(),panel.background = element_blank(),axis.line = element_line(colour = "black"),panel.border = element_rect(colour = "black", fill=NA, size=1))+
  geom_smooth(method="lm", se = FALSE)+
  xlab("Average Seed Mass (mg)") + ylab("Seed Si (mg/kg")+
  scale_color_manual(values=c('#CCCCCC','#333333'))
Si

library(ggpubr)
library(gridExtra)

library(cowplot)
plot_grid( K, Fe, Na, Mg, Zn, Si, 
           ncol = 3, nrow = 2)




###########################
######## FIGURE 8A ########
###########################

Ca <- ggplot(Consolidated_Data_R, aes(x = Ave_Seed_Mass.mg, y = Seed_Ca, color=Ecoregion,pch = Ecoregion)) + 
  geom_point(size = 2)+
  theme(panel.grid.major = element_blank(),text = element_text(size = 12),
        panel.grid.minor = element_blank(),panel.background = element_blank(),axis.line = element_line(colour = "black"),panel.border = element_rect(colour = "black", fill=NA, size=1))+
  xlab("Average Seed Mass (mg)") + ylab("Seed Ca (mg/kg)")+
  scale_color_manual(values=c('#CCCCCC','#333333'))
Ca




###########################
######## FIGURE 8B ########
###########################

Ca <- ggplot(Consolidated_Data_R, aes(x = Both_PW_Ca, y = Seed_Ca, color=Ecoregion,pch = Ecoregion)) + 
  geom_point(size = 2)+
  theme(panel.grid.major = element_blank(),text = element_text(size = 12),
        panel.grid.minor = element_blank(),panel.background = element_blank(),axis.line = element_line(colour = "black"),panel.border = element_rect(colour = "black", fill=NA, size=1))+
  xlab("Porewater Ca (mg/L)") + ylab("Seed Ca (mg/kg)")+
  scale_color_manual(values=c('#CCCCCC','#333333'))
Ca




###########################
######## FIGURE 8C ########
###########################

Ca <- ggplot(data = Consolidated_Data_R,aes(x = Lon, y = Lat, color = Both_PW_Ca,size = 1.5)) +
  geom_point()+theme(aspect.ratio=16/10, panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                     panel.background = element_blank())
Ca + scale_color_gradient(low="light blue", high="black")




#############################################################################################
######## APPENDIX FIGURES ###################################################################
#############################################################################################




##########################
######## FIGURE A2 #######
##########################

library(dplyr)

Seed_data_major <- Consolidated_Data_R %>% select(Seed_K, Seed_Ca, Seed_Mg, Seed_P)

library(GGally)
ggpairs(Seed_data_major)




###########################
######## FIGURE A3 ########
###########################

library(rpart) #for fitting decision trees
library(rpart.plot) #for plotting decision trees

tree <- rpart(Ave_Seed_Mass.mg ~ Type + Ecoregion + Water_Sulfate.mgL, data=Consolidated_Data_R, control=rpart.control(cp=.0001))

printcp(tree)

#identify best cp value to use
best <- tree$cptable[which.min(tree$cptable[,"xerror"]),"CP"]

#produce a pruned tree based on the best cp value
pruned_tree <- prune(tree, cp=best)

#plot the pruned tree
prp(pruned_tree,
    faclen=0, #use full names for factor labels
    extra=1, #display number of obs. for each terminal node
    roundint=F, #don't round to integers in output
    digits=5) #display 5 decimal places in output




###########################
######## Figure A4 ########
###########################

K <- ggplot(Consolidated_Data_R, aes(x = Ave_Seed_Mass.mg, y = Total_K, color=Ecoregion, pch = Ecoregion)) + 
  geom_point(size = 2)+
  geom_smooth(method = 'lm', se = F)+
  theme(panel.grid.major = element_blank(),text = element_text(size = 12),
        panel.grid.minor = element_blank(),panel.background = element_blank(),axis.line = element_line(colour = "black"),panel.border = element_rect(colour = "black", fill=NA, size=1))+
  xlab(NULL) + ylab("Seed K (µg)")+
  scale_color_manual(values=c('#CCCCCC','#333333'))+ggtitle("Group 1")
K


Mg <- ggplot(Consolidated_Data_R, aes(x = Ave_Seed_Mass.mg, y = Total_Mg, color=Ecoregion,pch = Ecoregion)) + 
  geom_point(size = 2)+
  geom_smooth(method="lm", se = FALSE)+
  theme(panel.grid.major = element_blank(),text = element_text(size = 12),
        panel.grid.minor = element_blank(),panel.background = element_blank(),axis.line = element_line(colour = "black"),panel.border = element_rect(colour = "black", fill=NA, size=1))+
  xlab("Average Seed Mass (mg)") + ylab("Seed Mg (µg)")+
  scale_color_manual(values=c('#CCCCCC','#333333'))
Mg
Fe <- ggplot(Consolidated_Data_R, aes(x = Ave_Seed_Mass.mg, y = Total_Fe, color=Ecoregion,pch = Ecoregion)) + 
  geom_point(size = 2)+
  theme(panel.grid.major = element_blank(),text = element_text(size = 12),
        panel.grid.minor = element_blank(),panel.background = element_blank(),axis.line = element_line(colour = "black"),panel.border = element_rect(colour = "black", fill=NA, size=1))+
  xlab(NULL)+ylab("Seed Fe (µg)")+geom_smooth(method="lm", se = FALSE)+
  scale_color_manual(values=c('#CCCCCC','#333333'))+ggtitle("Group 2")
Fe

Zn <- ggplot(Consolidated_Data_R, aes(x = Ave_Seed_Mass.mg, y = Total_Zn, color=Ecoregion,pch = Ecoregion)) + 
  geom_point(size = 2)+
  theme(panel.grid.major = element_blank(),text = element_text(size = 12),
        panel.grid.minor = element_blank(),panel.background = element_blank(),axis.line = element_line(colour = "black"),panel.border = element_rect(colour = "black", fill=NA, size=1))+
  xlab("Average Seed Mass (mg)") + ylab("Seed Zn (µg)")+geom_smooth(method="lm", se = FALSE)+
  scale_color_manual(values=c('#CCCCCC','#333333'))
Zn

Na<- ggplot(Consolidated_Data_R, aes(x = Ave_Seed_Mass.mg, y = Total_Na, color=Ecoregion,pch = Ecoregion)) + 
  geom_point(size = 2)+
  theme(panel.grid.major = element_blank(),text = element_text(size = 12),
        panel.grid.minor = element_blank(),panel.background = element_blank(),axis.line = element_line(colour = "black"),panel.border = element_rect(colour = "black", fill=NA, size=1))+
  geom_smooth(method="lm", se = FALSE)+
  xlab(NULL) + ylab("Seed Na (µg)")+
  scale_color_manual(values=c('#CCCCCC','#333333'))+ggtitle("Group 3")
Na

Si<- ggplot(Consolidated_Data_R, aes(x = Ave_Seed_Mass.mg, y = Total_Si, color=Ecoregion,pch = Ecoregion)) + 
  geom_point(size = 2)+
  theme(panel.grid.major = element_blank(),text = element_text(size = 12),
        panel.grid.minor = element_blank(),panel.background = element_blank(),axis.line = element_line(colour = "black"),panel.border = element_rect(colour = "black", fill=NA, size=1))+
  geom_smooth(method="lm", se = FALSE)+
  xlab("Average Seed Mass (mg)") + ylab("Seed Si (µg)")+
                                          scale_color_manual(values=c('#CCCCCC','#333333'))
Si
                                        
library(ggpubr)
library(gridExtra)
                                        
library(cowplot)
plot_grid( K, Fe, Na, Mg, Zn, Si, 
ncol = 3, nrow = 2)




#####################################################################################
#####################################################################################
############## DATA ANALYSIS ########################################################
#####################################################################################
#####################################################################################




# T test for average seed mass by ecoregion *** ####
t.test(log10(Ave_Seed_Mass.mg)~Ecoregion)

# T test Lakes vs. River for entire data set ###
t.test(log10(Ave_Seed_Mass.mg)~Type) 

# T test Lakes vs. River for Breda Lake vs. Petrel Creek ###
detach(Consolidated_Data_R)
attach(L_R_Data)

t.test(Mass_a~Type_a)

# T test Lakes vs. River for Shannon Lake vs. Shannon River ###
t.test(Mass_b~Type_b)

# T test Lakes vs. River for Lawrance Lake vs. Prarie River ###
t.test(Mass_c~Type_c)

detach(L_R_Data)
attach(Consolidated_Data_R)

# Linear model for Seed mass vs. Sediment bulk density ##
Bulkmod <- lm(Ave_Seed_Mass.mg~Sed_bulk_density)
summary(Bulkmod)

# Linear model for Seed mass vs. Loss on Ignition ##
LOImod <- lm(Ave_Seed_Mass.mg~L_O_I)
summary(LOImod)

# t test for Sediment bulk density vs. Waterbody type ##
t.test(Sed_bulk_density~Type)

# linear model for log 10 sulfate vs. seed mass##
smod <- lm(log10(Water_Sulfate.mgL)~ Ave_Seed_Mass.mg)
summary(smod)

# linear model for sulfate vs. Seed mass##
smod <- lm((Water_Sulfate.mgL)~ Ave_Seed_Mass.mg)
summary(smod)

#### Linear models for Table 3 ###
Mgmod <- lm(Seed_Mg~Ave_Seed_Mass.mg)
summary(Mgmod)

## linear mode for log10 sulfate vs. seed mass for northern ecoregions##

library(dplyr)
#Extract northeast data
Northeast_Data <- Consolidated_Data_R %>% filter(Ecoregion =="North")
head(Northeast_Data)
detach(Consolidated_Data_R)
attach(Northeast_Data)

Nsmod <- lm(log10(Water_Sulfate.mgL)~ Ave_Seed_Mass.mg)
summary(Nsmod)

## linear egression for log10 sulfate vs seed mass for southern/western ecoregions##
Southwest_Data <- Consolidated_Data_R %>% filter(Ecoregion == "South/West")

head(Southwest_Data)
detach(Northeast_Data)
attach(Southwest_Data)

Ssmod <- lm(log10(Water_Sulfate.mgL)~ Ave_Seed_Mass.mg)
summary(Ssmod)

# Sulfate kruskal willis test
detach(Southwest_Data)
attach(Consolidated_Data_R)
kruskal.test(Ave_Seed_Mass.mg~Sulfate_Quartiles)

# t test starch by ecoregions# 
t.test(Percent_Seed_Starch~Ecoregion)

## Starch vs seed mass linear regression ##
Starchmod <- lm(Percent_Seed_Starch~Ave_Seed_Mass.mg)
summary(Starchmod)

###Porewater calcium vs. Sulfate linear model##
Sulfate.Calcium.mod <- lm(log10(Water_Sulfate.mgL)~Both_PW_Ca)
summary(Sulfate.Calcium.mod)












Kmod <- lm(Seed_K~Ave_Seed_Mass.mg)
summary(Kmod)

Femod <- lm(Seed_Fe~Ave_Seed_Mass.mg)
summary(Femod)

Znmod <- lm(Seed_Zn~Ave_Seed_Mass.mg)
summary(Znmod)

Pmod <- lm(Seed_P~Ave_Seed_Mass.mg)
summary(Pmod)

Camod <- lm(Seed_Ca~Ave_Seed_Mass.mg)
summary(Camod)

Namod <- lm(Seed_Na~Ave_Seed_Mass.mg)
summary(Namod)

Simod <- lm(Seed_Si~Ave_Seed_Mass.mg)
summary(Simod)

Almod <- lm(Seed_Al~Ave_Seed_Mass.mg)
summary(Almod)

Mnmod <- lm(Seed_Mn~Ave_Seed_Mass.mg)
summary(Mnmod)

smod <- lm(Seed_S~Ave_Seed_Mass.mg)
summary(smod)

Bmod <- lm(Seed_B~Ave_Seed_Mass.mg)
summary(Bmod)

Crmod <- Ave_Seelm(Seed_Cr~Ave_Seed_Mass.mg)
summary(Crmod)

Cumod <- lm(Seed_Cu~Ave_Seed_Mass.mg)
summary(Cumod)


## regressions for appendix##

PWKmod <- lm(Seed_K~Pore_K)
summary(PWKmod)

PWFemod <- lm(Seed_Fe~Pore_Fe)
summary(PWFemod)

PWZnmod <- lm(Seed_Zn~Pore_Zn)
summary(PWZnmod)

PWPmod <- lm(Seed_P~Pore_P)
summary(PWPmod)

PWCamod <- lm(Seed_Ca~Pore_Ca)
summary(PWCamod)

PWNamod <- lm(Seed_Na~Pore_Na)
summary(PWNamod)

PWSimod <- lm(Seed_Si~ Pore_Si)
summary(PWSimod)

PWAlmod <- lm(Seed_Al~Pore_Al)
summary(PWAlmod)

PWMnmod <- lm(Seed_Mn~Pore_Mn)
summary(PWMnmod)

Pwsmod <- lm(Seed_S~Pore_S)
summary(PWsmod)

PWBmod <- lm(Seed_B~Pore_B)
summary(PWBmod)

PWCrmod <- lm(Seed_Cr~Pore_Cr)
summary(PWCrmod)

PWCumod <- lm(Seed_Cu~Pore_Cu)
summary(PWCumod)


# REGRESSIONS FOR SEDIMENT##



sedKmod <- lm(Seed_K~Sed_K)
summary(sedKmod)

sedFemod <- lm(Seed_Fe~Sed_Fe)
summary(sedFemod)

sedZnmod <- lm(Seed_Zn~Sed_Zn)
summary(sedZnmod)

sedPmod <- lm(Seed_P~Sed_P)
summary(sedPmod)

sedCamod <- lm(Seed_Ca~Sed_Ca)
summary(sedCamod)

sedNamod <- lm(Seed_Na~Sed_Na)
summary(sedNamod)

sedSimod <- lm(Seed_Si~Sed_Si)
summary(sedSimod)

sedMnmod <- lm(Seed_Mn~Sed_Mn)
summary(sedMnmod)

sedsmod <- lm(Seed_S~Sed_S)
summary(sedsmod)

sedBmod <- lm(Seed_B~Sed_B)
summary(sedBmod)

sedCumod <- lm(Seed_Cu~Sed_Cu)
summary(sedCumod)









#Average seed starch by region **
shapiro.test(log10(Percent_Seed_Starch))
hist(Percent_Seed_Starch)
t.test(Percent_Seed_Starch~Ecoregion)
wilcox.test(Percent_Seed_Starch~Ecoregion)
#seed starch by type
t.test(Percent_Seed_Starch~Type)





sedCamod <- lm(Ave_Seed_Mass.mg~Sed_Ca)
summary(sedCamod)

pwCamod <- lm(Ave_Seed_Mass.mg~PW_Ca)
summary(pwCamod)

t.test(Seed_Fe~Ecoregion)
t.test(Seed_Zn~Ecoregion)
t.test(Seed_Cu~Ecoregion)
t.test(Seed_P~Ecoregion)
t.test(Seed_Ca~Ecoregion)
t.test(Seed_Na~Ecoregion)
t.test(Seed_Si~Ecoregion)
t.test(Seed_K~Ecoregion)
t.test(Seed_Mg~Ecoregion)
t.test(Seed_Mn~Ecoregion)

t.test(Seed_Fe~Type)
t.test(Seed_Zn~Type)
t.test(Seed_Cu~Type)
t.test(Seed_P~Type)
t.test(Seed_Ca~Type)
t.test(Seed_Na~Type)
t.test(Seed_Si~Type)
t.test(Seed_K~Type)
t.test(Seed_Mg~Type)
t.test(Seed_Mn~Type)

t.test(Percent_Seed_Starch~Type)