---
title: "MS_figure_plots"
author: "Carlos Coquinho"
date: "25/03/2020"
output:
  word_document: default
  pdf_document: default
  html_document: default
---

```{r set files & folder, results='hide'}
library(dplyr)
#Read settings file
set_info <- read.table(file="../new_file_list.txt", header=T, sep="\t")

## Excluding from analyses the files with bad data vib073, vib074, vib079, vib080, vib083, vib084, vib101, vib102, vib103, vib104,vib105, vib106
## The even and uneven file numbers correspond to pairs of feeding and pollinating stamens from the same flowers of which one was not useful.
## Remember to also remove the files from the folder to avoid error due different argument lengths.

set_info <- filter(set_info, !grepl('vib073|vib074|vib079|vib080|vib083|vib084|vib101|vib102|vib103|vib104|vib105|vib106', vib_file))
set_info

#Analysis of all files combined:
workdir<-c("../Nat_freq_data/vib_data")
filelist<-list.files(workdir)

#Number of files:
totnumfiles<-length(filelist)
paste("Total number of files:", totnumfiles)

```


```{r create_loop, results='hide'}
library(dplyr)

num.files.read<-totnumfiles
fns <- vector("list", length(filelist)) # vectors to store the multiple outputs
fminuses <- vector("list", length(filelist))
fpluses <- vector("list", length(filelist))
dratios <- vector("list", length(filelist))


#creating a loop to read list of files
for(i in 1:num.files.read){ 
  
  #Create path and name of file to be read:
  file.to.read<-paste(workdir,"/", filelist[i], sep="")
  
  
  #Read data:
  dat.temp <- read.table(file=file.to.read, sep="\t", skip = 5)
  #dat.temp <- read.table("../Nat_freq_data/vib_data/vib003.txt", sep="\t", skip = 5)
  
  #Rename the columns:
  names(dat.temp) <- c("frequency", "acceleration") ##Frequency (hz), Acceleration (m/s^2)
  
  #Filter data below 650Hz in frequency:
  dat.temp <- dat.temp[c(dat.temp$frequency<=400),]
  
  ## Obtaining the first mode or the first natural frequency (fn)
## fna1 <- srp1$`frequency (Hz)`[which.max(srp1$acceleration)] #returns the frequency value of maximum acceleration
fna.temp <- subset(dat.temp, acceleration == max(dat.temp$acceleration)) # returns the frequency and acceleration values of max acceleration
fn.temp <- fna.temp$frequency %>% round(digits = 3)
fns[[i]] <- fn.temp # assembling the fn`s into a vector

accdec.temp <- (fna.temp$acceleration)*0.7071068 # the decrease in acceleration we seek in the dataset for obtaining the f- and f+

# Creating a vector with simulated values of freq and acceleration.
# where the freq. is the same range as studied, acceleration is the created function function of Freq varying in a finer scale

freq.temp <- as.vector(seq(from = min(dat.temp$frequency), to = max(dat.temp$frequency), length.out = 50000))

xtemp <- dat.temp$frequency
ytemp <- dat.temp$acceleration

fun.temp <- splinefun(xtemp, ytemp, method = "fmm") #creates a function from the data

acc.temp <- as.vector(fun.temp(freq.temp))

simSpec.temp <- data.frame(freq.temp, acc.temp)
names(simSpec.temp) <- c("frequency", "acceleration")

# creating a interval comprising the acceleration values higher than the decreases, including the maximum value
interv.temp <- simSpec.temp%>%filter(acceleration > accdec.temp)

## calculating the f-`s from the fn`s
faminus.temp <- subset(interv.temp, frequency == min(interv.temp$frequency))
rAcminus.temp <- faminus.temp$acceleration %>% round(digits = 3)

fminus.temp <- faminus.temp$frequency %>% round(digits = 3)

faplus.temp <- subset(interv.temp, frequency == max(interv.temp$frequency))

rAcplus.temp <- faplus.temp$acceleration %>% round(digits = 3)
fplus.temp <- faplus.temp$frequency %>% round(digits = 3)

dratio.temp <- (fplus.temp - fminus.temp)/(2*fn.temp)
dratio.temp <- round(dratio.temp, digits = 3)

fns[[i]] <- fn.temp # assembling the fn`s into a vector
fminuses[[i]] <- fminus.temp # assembling the f-`s into a vector
fpluses[[i]] <- fplus.temp # assembling the f+`s into a vector
dratios[[i]] <- dratio.temp # assembling the damping rations into a vector
}

#creating vectors from the list of single values
fns <- as.vector(unlist(fns))
fminuses <-as.vector(unlist(fminuses))
fpluses <- as.vector(unlist(fpluses))
dratios <- as.vector(unlist(dratios))

#creating matrices from vectors and them dataframes from matrices
fns.matrix <- matrix(fns, ncol = 1)
fns.df <- as.data.frame(fns.matrix)

fminuses.matrix <- matrix(fminuses, ncol = 1)
fminuses.df <- as.data.frame(fminuses.matrix)

fpluses.matrix <- matrix(fpluses, ncol = 1)
fpluses.df <- as.data.frame(fpluses.matrix)

dratios.matrix <- matrix(dratios, ncol = 1)
dratios.df <- as.data.frame(dratios.matrix)

new.var.df <- bind_cols(fns.df, fminuses.df, fpluses.df, dratios.df)
colnames(new.var.df)<- c("fn", "fminus", "fplus", "dratio")
new.var.df


total.var.df <- bind_cols(set_info, new.var.df)
total.var.df <- select(total.var.df, -recovered_fn, -notes) # Removing unused columns from the dataset
total.var.df
```
```{r loading packages, results='hide'}
library(dplyr)
library(ggplot2)
library(Hmisc)
```

## Figure 1 - B and C

# Datasets
```{r datasets, results='hide'}

# S. rostratum pollination anthers (srp`s)
srp5 <- read.table("../Nat_freq_data/data/vib017.txt", sep="\t", skip = 5)
names(srp5) <- c("frequency (Hz)", "accelleration")
srp5 <- srp5[c(srp5$`frequency (Hz)`<=650),]

# S. rostratum feeding anthers (srf`s)
srf5 <- read.table("../Nat_freq_data/data/vib018.txt", sep="\t", skip = 5)
names(srf5) <- c("frequency (Hz)", "accelleration")
srf5 <- srf5[c(srf5$`frequency (Hz)`<=650),]

# S. fructo-tecto pollination anthers (sfp`s)
sfp5 <- read.table("../Nat_freq_data/data/vib039.txt", sep="\t", skip = 5)
names(sfp5) <- c("frequency (Hz)", "accelleration")
sfp5 <- sfp5[c(sfp5$`frequency (Hz)`<=650),]

# S. fructo-tecto feeding anthers (sff`s)
sff5 <- read.table("../Nat_freq_data/data/vib040.txt", sep="\t", skip = 5)
names(sff5) <- c("frequency (Hz)", "accelleration")
sff5 <- sff5[c(sff5$`frequency (Hz)`<=650),]
```

# Figure 1B
```{r Fig 1B, fig.height= 8, fig.width= 11}
library(ggplot2)
library(grid)

# Fig 1B S. fructo-tecto - pollinating and feeding stamen
fig1b.plot <- ggplot(NULL) +
  aes(x = `frequency (Hz)` , y = accelleration) +
  geom_line(data = sfp5, size = 1) +
  geom_line(data = sff5, size = 1, color = "darkgrey") +
  theme_bw()+
  xlab("Frequency (Hz)")+ 
  ylab(bquote('Acceleration ('*'m'~s^-2*')'))+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())+
  theme(axis.title.x = element_text(size=30, hjust = .44), 
                  axis.title.y = element_text(size=30))+
  theme(axis.text.y = element_text(color = "black", size = 24),
        axis.ticks.y = element_line(color = "black"))+
  theme(axis.text.x = element_text(color = "black", size = 24),
        axis.ticks.x = element_line(color = "black"))

pdf(file="Fig1bR.pdf", width=11, height=8)
fig1b.plot
dev.off()
fig1b.plot
```
Figure 1. B. Resonance pattern under white-noise vibration of feeding and pollinating anthers of the non-heterantherous Solanum fructo-tecto.


# Figure 1C
```{r Fig 1C, fig.height= 8, fig.width= 11}
library(ggplot2)
library(grid)

# Fig 1C - S. rostratum - pollinating and feeding stamen
fig1c.plot <- ggplot(NULL) +
  aes(x = `frequency (Hz)` , y = accelleration) +
  geom_line(data = srp5, size = 1) +
  geom_line(data = srf5, size = 1, color = "darkgrey") +
  theme_bw()+
  xlab("Frequency (Hz)")+ 
  ylab(bquote('Acceleration ('*'m'~s^-2*')'))+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())+
  theme(axis.title.x = element_text(size=30,hjust= .44), 
                  axis.title.y = element_text(size=30))+
  theme(axis.text.y = element_text(color = "black", size = 24),
        axis.ticks.y = element_line(color = "black"))+
  theme(axis.text.x = element_text(color = "black", size = 24),
        axis.ticks.x = element_line(color = "black"))

pdf(file="Fig1cR.pdf", width=11, height=8)
fig1c.plot
dev.off()
fig1c.plot
```
Figure 1. C.Resonance pattern under white-noise vibration (acceleration vs. frequency) of feeding and pollinating anthers of a typical heterantherous species, Solanum rostratum.


# Figure 2

First, create another a dataset selecting out the variables of interest from the full dataset "total.var.df".

```{r dataset Fig 2, results='hide'}
library(dplyr)

bdat<-read.csv(file="ColonyA1cit_to_het.csv")

df.plots <- select(total.var.df, flower_id, plant_sp, accession, anther_type, vib_file, fn, fminus, fplus, dratio)
df.plots <- as_tibble(df.plots)
df.plots <- rename(df.plots, Stamen = anther_type) #Renaming the column "anther_type"

df.plots$Stamen <-  recode_factor(df.plots$Stamen, 'pol' = "Pollinating", 'feed' = "Feeding")

df.plots$plant_sp <- recode_factor(df.plots$plant_sp, 'Solanum_rostratum' = "Solanum_rostratum", 'Solanum_fructo-tecto' = "S_fructo-tecto", 'Solanum_citrullifolium' = "S_citrullifolium", 'Solanum_heterodoxum' = "S_heterodoxum", 'Solanum_grayi_grandiflorum'= "S_grayi_grandiflorum", 'Solanum_grayi_grayi' = "S_grayi_grayi")

total.var.df
df.plots

mean_cl_boot(bdat$Frequency)
```

Then, we plot fns per categories (anther types within plant species).

```{r Fig 2, fig.height= 8, fig.width= 11}
library(grid)
library(ggplot2)

Fig2.plot <- ggplot(df.plots, 
                         aes(x = plant_sp, y = fn)) +
            stat_summary(fun.data = "mean_cl_boot",
                         position = position_dodge2(width = .2),
                         size = 1.5,
                         geom="errorbar", width=.15,
                         color = "darkgray",
                         aes(group= Stamen))+
            stat_summary(fun.y=mean,geom="point",
                         position = position_dodge2(width = .2),
                         shape = 21,
                         size=4.8, alpha=1,
                         stroke = 1.5,
                         aes(fill = Stamen))+
            scale_fill_manual(values = c("#440154FF", "#FDE725FF"))+
            xlab("Plant taxa")+ 
            ylab("Natural frequency (Hz)")+
            scale_x_discrete(labels=c("Solanum rostratum", 
                                      "S. fructo-tecto", 
                                      "S. citrullifolium",
                                      "S. heterodoxum", 
                                      "S. grayi grandiflorum", 
                                      "S. grayi grayi"))+
            theme_bw() +
            theme(panel.grid.major = element_blank(),
                  panel.grid.minor = element_blank())+
            theme(axis.title.x = element_text(size=20), 
                  axis.title.y = element_text(size=20)) +
            theme(axis.text.x = element_text(angle=0, hjust=0.5,
                                             vjust=0.3, size=16),
                  axis.text.y = element_text(size=16))+
            theme(axis.text.x = element_text(angle=15, hjust=.75,
                                             vjust=.9, size=16))+
                                             #aes(fontface = 3)))+
            theme(legend.position = c(0.175, 0.882))+
            #theme(legend.position = c(0.2, 0.8))+
            theme(legend.title = element_text(size=12))+
            theme(legend.text = element_text(size=12))+
            theme(legend.background = element_rect(size=0.5, 
                                                linetype="solid",
                                                colour = "black"))+
            theme(legend.title=element_blank())+
            geom_hline(yintercept = 400, 
                       linetype = "dotted", size = 1.1)+
            geom_hline(yintercept = 100, 
                       linetype = "dotted", size = 1.1)+
            geom_hline(yintercept = 344.87, size = 1.1,
                       linetype = "dashed")+
            theme(axis.text.y = element_text(color = "black"),
                  axis.ticks.y = element_line(color = "black"))+
            theme(axis.text.x = element_text(color = "black",
                                             face = "italic"),
                  axis.ticks.x = element_line(color = "black"))
Fig2.plot
pdf(file="Fig2R.pdf", width=11, height=8)
Fig2.plot
dev.off()


```
Figure 2. Resonant frequencies of different functional types of anther (feeding and pollinating) for the six Solanum taxa studied with the range of frequencies produced by bees while vibrating flowers to gather pollen indicated by area, illustrating the general morphology of flowers and stamens of each taxa.


# Figure S1 - Morphology and Natural Freq Correlation Plot
```{r plot 3, fig.height= 8, fig.width= 11}
library(dplyr)
library(corrplot)

corr.dat <- select(total.var.df, fn, w_corolla_diam_mm, h_corolla_diam_mm, flower_mass_g, stamen_length_mm, anther_length_mm, fil_length_mm, ldia_anther_mm, sdia_anther_mm, ldia_fil_mm, sdia_fil_mm, stamen_mass_g)

corr.morpho <- cor(corr.dat)
colnames(corr.morpho) <- c("Natural_frequency",
                           "Corolla_width_mm",
                           "Corolla_heigth_mm",
                           "Flower_mass_g",
                           "Stamen_length_mm",
                           "Anther_length_mm",
                           "Filament_length_mm",
                           "Anther_long_diameter_mm",
                           "Anther_short_diameter_mm",
                           "Filament_long_diameter_mm",
                           "Filament_short_diameter_mm",
                           "Stamen_mass_g")
rownames(corr.morpho) <- c("Natural_frequency",
                           "Corolla_width_mm",
                           "Corolla_heigth_mm",
                           "Flower_mass_g",
                           "Stamen_length_mm",
                           "Anther_length_mm",
                           "Filament_length_mm",
                           "Anther_long_diameter_mm",
                           "Anther_short_diameter_mm",
                           "Filament_long_diameter_mm",
                           "Filament_short_diameter_mm",
                           "Stamen_mass_g")


pdf(file="Fig3.pdf", width=11, height=8)
corrplot(corr.morpho, order="hclust")
dev.off()

corr_plot <- corrplot(corr.morpho, order="hclust")
corr_plot

```




