PS117.Rmd

To do - 
Figure out SD and SE

2024-07-25 --  load packages 
```{r}
#plotting 
library (ggplot2) 
library (ggpattern)#for making cross hatches in the bars
library (ggplotify)
library (cowplot)
library (ggpubr) #for vertical legends 
library (gridExtra)
library (grid)
library (ggthemes)
library (patchwork) #allows for plot insets. e.g. map and the DFe and DMn in Figure 1
library (ggridges)
library (plotly)
library (svglite) #allows the export of .svg plots
#library (ggrepel) #makes point labels not touch each other in RA plots


#other
library(data.table)
library(dplyr)
library(magrittr)
library(reshape2)
library(stringr)
library(scales)
library(tidyr)
library(tidyverse)
library (patchwork) #for combining some plots


setwd("C:/Users/loayj/OneDrive - Woods Hole Oceanographic Institution/postdoc/other_projects/PS117/data/") #this is the directory where I keep all the data

#library(vegan) #for making NMDS analyses
```

2024-07-25 -- import data of interest 
```{r}
#import non-protein data
nonproteindata <- read.csv("LJ_ps117_allnonproteindata_20240424.csv", header = T)
nonproteindata_long <- melt (nonproteindata, id.vars = c(1:16)) #might need corrections here

#import protein data
proteindata_unnorm <- read.csv("all_ps117_taxon-function_peptides_non-normalized_injections-means_20230324.csv", header = T)

#arrange the 'uniq_cluster' annotations alphabetically in each row to reduce redundancies. 
proteindata_unnorm$uniq_cluster <- sapply(lapply(strsplit(proteindata_unnorm$uniq_cluster, split = ",\\s*"), sort), paste, collapse = ", ") 

### a quick test to show how to normalize to all of proteome, vs to group
# https://stackoverflow.com/questions/50332613/normalize-by-group-for-all-columns
#group <- c("A", "A", "B", "B", "C", "C", "C", "A")
#abundance <- c(1,3,4,5,6,7,8,0)
#all_norm <- c(1,3,4,5,6,7,8,0)
#grp_norm <- c(1,3,4,5,6,7,8,0)
#x <- data.frame(group, abundance, all_norm, grp_norm) %>% mutate_at (3, funs ((. / sum(.)))) # %>%   group_by(group) %>% #mutate_at (4, funs ((. / sum(.))


#counting number of peptides between samples etc
#proteindata <- proteindata_unnorm[c(40:69)]

#proteindata$zeroz <- rowSums(proteindata!=0)
```

2024-07-25 -- plotting functions and other useful things
```{r}
#https://r-charts.com/distribution/density-plot-group-ggplot2/ good colors 
#used in figure 1
x_vs_time <- function (yaxis, yaxislabel) {
  ggplot(nonproteindata, aes(x = daycount, 
                             y= {{yaxis}},
                             color = temperature_treatment, 
                             shape = iron_treatment)) + 
  stat_summary(fun = mean, geom = "point", show.legend = F, size = 2.5, stroke = 1, alpha = 0.9)+
  stat_summary(fun = mean, geom = "line", linewidth = 0.6, aes(linetype = iron_treatment))+
  stat_summary(fun.data =  mean_se, geom = "errorbar", size = 0.7, width = 0.2) +
  ylab (substitute(yaxislabel)) + 
  xlab (expression ("")) +
  scale_color_manual(name="Temperature",
                     breaks = c("LT","HT", "T0"),
                     labels = c("LT", "HT", "T0"),
                     values = c('black', "grey50", "black")) + #D55E00
  scale_linetype_manual(name = "Iron",
                        breaks = c("Fe","noFe"),   #this removes the T0 from the legend
                        labels = c("Fe", "noFe"),
                        values = c("solid", "dashed", "solid")) +
  scale_shape_manual (name = "Iron",
                       breaks = c("noFe","Fe", "T0"), 
                       labels = c("noFe", "Fe", "T0"),
                       values = c(1, 19, 5))+
  scale_x_continuous( labels = c("T0","T2","T4", "T6", "T8"))+
    facet_wrap(~bioassay)+
    theme_bw()+
        theme(axis.text=element_text(face = "bold", size = 10, color = "black"),
        axis.title.x=element_text(size = 10, color = "black"), 
        axis.title.y=element_text(size = 12, color = "black"),
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 10, face = "bold"),
        legend.position = "none", 
        panel.grid = element_blank()) 
}

#x_vs_time(NO3_uM, NO[3]^"-"~~(mu*M)) + ylim(5, 30) + theme (axis.text.x=element_blank())

#used in figure 1,2,3,4
T0_T8 <- function(dataset, yaxis, yaxislabel) {
   ggplot(filter ({{dataset}}, grepl ("T0|T8", timepoint)), 
         aes(x = tempfe_treatment, y = {{yaxis}})) +
    facet_wrap(~bioassay) +     
    geom_bar_pattern(stat = "summary", fun = "mean", size = 0.75, alpha = 0.9, 
                     aes(fill = tempfe_treatment, color = tempfe_treatment, 
                         pattern = tempfe_treatment, pattern_fill = tempfe_treatment, pattern_color = tempfe_treatment), 
                     pattern_density = 0.2, pattern_size = 0.2, pattern_spacing = 0.05) +
    stat_summary(fun.data = mean_se, geom = "errorbar", show.legend = F, linewidth = 0.7, width = 0.2, aes(color = tempfe_treatment)) +
     scale_x_discrete(limits = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"), 
                     labels = c("T0", "Ctrl", "+Fe", "HT", "HT +Fe")) +
    scale_fill_manual(breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      values = c('white', "white", "black", "white", "grey50")) +
    scale_color_manual(breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                        values = c('black', "black", "black", "grey30", "grey30")) +
    scale_pattern_manual(breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                         values = c('crosshatch', "none", "none", "none", "none")) +
    scale_pattern_fill_manual(breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                              values = c('white', NA, NA, NA, NA)) +
    scale_pattern_color_manual(breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                               values = c('grey30', NA, NA, NA, NA)) +
    xlab(expression("")) +
    ylab(substitute(yaxislabel)) + 
    theme_bw() +
    theme(axis.text = element_text(face = "bold", size = 9, color = "black"),
          axis.title.x = element_text(size = 10, color = "black"), 
          axis.title.y = element_text(size = 12, color = "black"),
          strip.background = element_rect(fill = "white"),
          strip.text.x = element_text(size = 10, face = "bold"),
          legend.position = "none", 
          panel.grid = element_blank()) 
}


```

2024-07-25 -- legend for figure 1 and maybe other figures
```{r}
#this legend is used in most of the plots, so it's good to have 
vertical_legend <- ggplot(filter (nonproteindata, grepl ("T0|8", full_treatment_name)),
       aes(x=iron_treatment, 
       y=fvfm,
       color = tempfe_treatment, 
       shape = tempfe_treatment))+
geom_point(size = 2, stroke = 1) + 
   
scale_color_manual(breaks = c("LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c(" Ctrl", " +Fe", " HT", " HT +Fe"),
                      values = c('black', 'black', 'grey50', "grey50", "black")) +
  scale_shape_manual(breaks = c("LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c(" Ctrl", " +Fe", " HT", " HT +Fe"),
                      values = c(1, 19,1, 19, 5)) +
  theme_bw() +
  
    theme(legend.title = element_blank(),
          legend.background = element_blank(),
        legend.text = element_text(size = 9, margin = margin(unit = 'cm')))+
    theme(legend.key.size = unit(0.1, "cm"))


vertical_legend

vertical_legend <- get_legend(vertical_legend)
vertical_legend <- as_ggplot(vertical_legend)
vertical_legend
```

#2024-08-27 -- Figure 1 
2024-05-24 -- map of Antarctica showing bioassay sites
```{r}
library(ggOceanMaps)
library (ggspatial)
#the numbers in the limits represent the range of the longs and lats 
#the first two numbers represent the longitude range to be shown, the secnd two numbers are the latitude range to be shown. 
bioassay_latlon <- data.frame(lon = c(-11.044,-0.011), lat =c(-70.25, -65), bioassay= c("Nearshore", "Offshore"))  #bioassay= c("BA2", "BA1")) 

#basemap(limits = c(-130, 100, -80,-50)) #for sqaure
#png("../figures/ps117_map_20241129.png", width= 2, height= 2, units = "in",  res = 600) 

map <- basemap(limits = c(-120, 40, -60,-50), #for sqaure #limits = c(-60), 
        bathymetry = TRUE, #this adds the bathymetry (depth info)
        glaciers = FALSE,
        grid.col = "NA")+ 
        labs(x = "", y = "") +
 # annotation_scale(location = "br") + 
        geom_spatial_point(data = bioassay_latlon, 
                           aes(x = lon, y = lat, shape = bioassay), 
                           color = "darkred", size = 3) +
  scale_shape_manual(values = c("Offshore" = 16, "Nearshore" = 17)) +
        geom_spatial_label(size = 3.5, fontface = "bold",  aes(label=bioassay_latlon$bioassay, x = bioassay_latlon$lon, #could also be "text"
                                                            y = bioassay_latlon$lat), hjust = 1.1, vjust = -0.2, color = "black") +
  theme(axis.text = element_blank(), 
              axis.ticks.y =element_blank(),
              axis.ticks.x =element_blank(),
              legend.position = "none") 
#dev.off()
map

##Additional information about maps: 
#https://mran.microsoft.com/snapshot/2021-03-14/web/packages/ggOceanMaps/vignettes/ggOceanMaps.html
#https://cran.r-project.org/web/packages/ggOceanMaps/ggOceanMaps.pdf
#bioassays <- data.frame(lon = c(165.4133), lat =c(-77.61895), bioassay= c("")) # mcmurdo ross sea TFG incubations
```

2024-11-11 -- figure 1, metals
```{r}
nonproteindata_long_modified <- nonproteindata_long %>%
  mutate(value = case_when(
    variable == "dMn_nM" ~ value * 10,
    variable == "dCo_nM" ~ value * 10,
    variable == "dFe56_nM" ~ value * 10,

    TRUE ~ value
  ))


dmetals <- ggplot(filter (nonproteindata_long_modified, grepl("dFe56|dMn|dCu|dNi|dCd|dCo|dZn", variable) & grepl ("0", timepoint)), 
       aes (x= variable,y = value, shape = bioassay))+
       stat_summary(fun = mean, geom = "point", size = 3)+
       stat_summary(fun.data =  mean_se, geom = "errorbar", show.legend = F, size = 0.7, width = 0.1) +
       labs(x = expression (""), y = expression ("Dissolved in situ concentration (nM)"))+
       theme_bw()+
  scale_x_discrete(limits = c("dFe56_nM", "dCo_nM", "dMn_nM", "dCd_nM", "dCu_nM", "dZn_nM", "dNi_nM"), 
                 labels = c("Fe (x10)", "Co (x10)", "Mn (x10)", "Cd", "Cu", "Zn", "Ni")) +
  #scale_y_log10()+
  
scale_shape_manual (name = "",
                       breaks = c("Nearshore","Offshore"), 
                       values = c(2, 1))+
  theme_bw()+
       theme(axis.text=element_text(face = "bold", size = 10, color = "black"),
        axis.title=element_text(size = 12,face = "bold", color = "black"), 
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 10, face = "bold"),
        legend.position = "none", #c(0.9,0.9), 
        panel.grid = element_blank())
dmetals






elementalratios <- read.csv("LJ_ps117_allnonproteindata_20240424.csv", header = T)


elementalratios$dFepFeratio <- elementalratios$dFe56_nM/elementalratios$pFe56_nM
elementalratios$dMnpMnratio <- elementalratios$dMn_nM/elementalratios$pMn_nM
elementalratios$dCupCuratio <- elementalratios$dCu_nM/elementalratios$pCu_nM
elementalratios$dZnpZnratio <- elementalratios$dZn_nM/elementalratios$pZn_nM
elementalratios$dNipNiratio <- elementalratios$dNi_nM/elementalratios$pNi_nM
elementalratios$dCdpCdratio <- elementalratios$dCd_nM/elementalratios$pCd_nM
elementalratios$dCopCoratio <- elementalratios$dCo_nM/elementalratios$pCo_nM
elementalratios$dNpNratio <- elementalratios$NOx_uM/elementalratios$PON_uM
elementalratios$dPpPratio <- elementalratios$PO4_uM/elementalratios$POP_nM*1000


elementalratios_long <- melt (elementalratios, id.vars = c(1:16)) #might need corrections here




elementalratio <- ggplot(filter (elementalratios_long, grepl("ratio|Si", variable) & grepl ("0", timepoint)), 
       aes (x= variable,y = value, shape = bioassay))+
       stat_summary(fun = mean, geom = "point", size = 3)+
       stat_summary(fun.data =  mean_se, geom = "errorbar", show.legend = F, size = 0.7, width = 0.1) +
       labs(x = expression (""), y = expression ("In situ dissolved:particulate concentration (M:M)"))+
       theme_bw()+
scale_x_discrete(limits = c("dNipNiratio", "dCupCuratio", "dPpPratio", "dNpNratio", "dCopCoratio", "dZnpZnratio", "dCdpCdratio", "dMnpMnratio", "dFepFeratio"),
                 labels = c("Ni", "Cu", "P", "N", "Co", "Zn", "Cd", "Mn", "Fe")) +
  scale_y_log10(breaks = 10^seq(-2, 10, by = 1),
                labels = trans_format("log10", math_format(10^.x))) +
    annotation_logticks(sides = "l") +
  geom_hline(yintercept=1, linetype = "dashed", size = 0.5) +

  
scale_shape_manual (name = "",
                       breaks = c("Nearshore","Offshore"), 
                       values = c(2, 1))+
  theme_bw()+
       theme(axis.text=element_text(face = "bold", size = 10, color = "black"),
        axis.title=element_text(size = 12,face = "bold", color = "black"), 
      #  legend.position = c(0.9,0.94), 
         legend.position = 'none', 

        legend.background = element_blank(),
        legend.box.background = element_blank(),
        panel.grid = element_blank())

elementalratio
```

2024-11-11 -- figure 1, nmds 
           -- The vegan package displays 'species' and 'sites'. The sites are the treatments, the 'species' are the variables (pigments, FC-groups, proteomics species etc) 
           -- First, make NMDS analyses from peptide, FC and pigment data. No need to re-run this! 
           -- Peptide NMDS -  2023-03-026 
```{r}
library(vegan)
library (goeveg)


c13_nmds <- c13_long_2 [c(1, 40:69)]

protein_nmds <- melt(protein_nmds, id.vars=c("peptide")) %>%
             rename(treatment = variable) %>%
             rename(unnorm_abundance = value )
sum(protein_nmds$unnorm_abundance) #this should equal 1.893E13 (total abundance of all peptides for all samples)

#normalize the peptide abundances for within each sample, irrespective of taxonomy
protein_nmds <- protein_nmds %>% #this creates a new column of normalized protein abundances grouped by treatment   
  group_by(treatment) %>% 
  mutate(norm_abundance = unnorm_abundance / sum(unnorm_abundance)) %>% 
  ungroup()

sum(protein_nmds$norm_abundance) #this should = 30 after normalization (number of replicates)

protein_nmds2 <- dcast(protein_nmds, treatment~peptide, fun.aggregate = sum,  value.var = "norm_abundance") #this creates MANY columns

protein_nmds2 <- protein_nmds2 %>% separate(treatment, into=c("bioassay", "timepoint" ,"temperature_treatment", "iron_treatment", "replicate"), sep = "_")

protein_nmds3 = protein_nmds2 [c(6:37477)] 

protein_nmds3 = as.matrix(protein_nmds3)

dimcheckMDS(protein_nmds3,
  distance = "bray",
  k = 4,
  trymax = 1000,
  autotransform = TRUE)

set.seed(123)

nmds_protein = metaMDS(protein_nmds3, try = 1000, trymax = 1000, k = 2, distance = "bray") # try = 500, trymax = 500
nmds_protein
scores(nmds_protein, display = "sites")


stressplot(nmds_protein, p.col = "black", l.col = "blue", title("Peptides")) #for the suplement

data.scores_protein <- as.data.frame(scores(nmds_protein, "sites")) #can select sites (treatments) or species (variable)

data.scores_protein$temperature_treatment <- protein_nmds2$temperature_treatment
data.scores_protein$iron_treatment <- protein_nmds2$iron_treatment
data.scores_protein$bioassay <- protein_nmds2$bioassay

data.scores_protein$bioassaytempfe <- paste(data.scores_protein$bioassay, data.scores_protein$temperature_treatment,data.scores_protein$iron_treatment, sep = "_")

data.scores_protein$bioassayfe <- paste(data.scores_protein$bioassay, data.scores_protein$iron_treatment, sep = "_")


# peptidenmds <- ggplot(data.scores_protein, aes(x = NMDS1, y = NMDS2))+
#           geom_point(size = 2, stroke = 1, 
#                     aes(fill = temperature_treatment,
#                     shape = bioassaytempfe,
#                     color = temperature_treatment))+
#     stat_ellipse(aes(fill=bioassay), alpha=.08, type='t', geom="polygon") +
# 
#           scale_shape_manual(values = c("Offshore_T0_T0" = 1, 
#                                         "Offshore_LT_noFe" = 1,
#                                         "Offshore_LT_Fe" = 19,
#                                         "Offshore_HT_Fe" = 19,
#                                         "Offshore_HT_noFe" = 1, 
#                                         "Nearshore_T0_T0" = 2, 
#                                         "Nearshore_LT_noFe" = 2,
#                                         "Nearshore_LT_Fe" = 17,
#                                         "Nearshore_HT_Fe" = 17,
#                                         "Nearshore_HT_noFe" =2 ))+
#   
#      annotate("text", x = max(data.scores_protein$NMDS1) * 0.9, 
#              y = 0.4, 
#              label = "Nearshore", 
#              size = 3, 
#              hjust = 0.8, fontface = "bold") + 
#    
#      annotate("text", x = max(data.scores_protein$NMDS1) * 0.9,
#              y = 0.8,
#              label = "Offshore",
#              size = 3, 
#              hjust = 6.9, fontface = "bold")+
#  theme_bw()+
# 
#           scale_color_manual(values = c("HT" = "grey50", "LT" = "black", "T0" = "darkred"))+
#           scale_fill_manual(values = c("HT" = "grey50", "LT" = "black", "T0" = "darkred"), guide = FALSE) +
#           ggtitle("Peptide-based NMDS")+
#           theme(axis.text = element_text(colour = "black", size = 10, face = "bold"), 
#                 axis.title = element_text(size = 10), 
#                # legend.position = "none", 
#                 panel.grid = element_blank(), 
#                 legend.key=element_blank(),
#                 plot.title = element_text( size = 8, vjust = -0.5))
# 
# peptidenmds





peptidenmds <- ggplot(data.scores_protein, aes(x = NMDS1, y = NMDS2))+
          geom_point(size = 2, stroke = 1, 
                    aes(fill = temperature_treatment,
                    shape = bioassaytempfe,
                    color = bioassaytempfe))+
    stat_ellipse(aes(fill=bioassay), alpha=.08, type='t', geom="polygon") +

          scale_shape_manual(values = c("Offshore_T0_T0" = 1, 
                                        "Offshore_LT_noFe" = 1,
                                        "Offshore_LT_Fe" = 19,
                                        "Offshore_HT_Fe" = 19,
                                        "Offshore_HT_noFe" = 1, 
                                        "Nearshore_T0_T0" = 2, 
                                        "Nearshore_LT_noFe" = 2,
                                        "Nearshore_LT_Fe" = 17,
                                        "Nearshore_HT_Fe" = 17,
                                        "Nearshore_HT_noFe" =2 ), 
                               name = "",
                             
 breaks = c("Offshore_HT_Fe", "Offshore_HT_noFe", "Offshore_LT_Fe", "Offshore_LT_noFe", "Offshore_T0_T0"), 
                     labels = c("HT +Fe", "HT", "+Fe", "Ctrl", "In situ")) +

  
  scale_color_manual(values = c("Offshore_T0_T0" = "darkred",
                                        "Offshore_LT_noFe" = "black",
                                        "Offshore_LT_Fe" = "black",
                                        "Offshore_HT_Fe" = "grey50",
                                        "Offshore_HT_noFe" = "grey50",
                                        "Nearshore_T0_T0" = "darkred", 
                                        "Nearshore_LT_noFe" = "black",
                                        "Nearshore_LT_Fe" = "black",
                                        "Nearshore_HT_Fe" = "grey50",
                                        "Nearshore_HT_noFe" = "grey50" ), 
                               name = "", 
                        breaks = c("Offshore_HT_Fe", "Offshore_HT_noFe", "Offshore_LT_Fe", "Offshore_LT_noFe", "Offshore_T0_T0"), 
                     labels = c("HT +Fe", "HT", "+Fe", "Ctrl", "In situ")) +
  
  
     annotate("text", x = max(data.scores_protein$NMDS1) * 0.9, 
             y = 0.4, 
             label = "Nearshore", 
             size = 3, 
             hjust = 0.8, fontface = "bold") + 
   
     annotate("text", x = max(data.scores_protein$NMDS1) * 0.9,
             y = 0.8,
             label = "Offshore",
             size = 3, 
             hjust = 7, fontface = "bold")+
 theme_bw()+

     
  
           scale_fill_manual(values = c("HT" = "grey50", "LT" = "black", "T0" = "darkred"), guide = FALSE) +
         # ggtitle("Peptide NMDS")+
          theme(axis.text = element_text(colour = "black", size = 10, face = "bold"), 
                axis.title = element_text(size = 10), 
                legend.position = c(0.92,0.14),
                panel.grid = element_blank(), 
                legend.key.spacing.y = unit(-1, "cm"),
                legend.key=element_blank(),
                legend.margin = margin(0, 0, 0, -10),
                legend.background = element_blank(),
                legend.box.background = element_blank(),
                plot.title = element_text( size = 8, vjust = -0.5))

peptidenmds

```

2024-11-11 -- figure 1, circlepack map, -- no need to re=run every time
           -- This figure divides all the protein abundances in the communities into major groups (taxon resolved), and presents them using circular packing (similar to a treemap)
           -- First we need to specificy the groups and combine the abundances of the groups 
```{r}
proteindata_unnorm$compartment <- ifelse (proteindata_unnorm$uniq_cluster  ==  "unassigned_cluster", "ribosome", 
                                  ifelse (proteindata_unnorm$uniq_grpnorm_compartment == "ribosomal", "ribosome",
                                  ifelse (grepl("ribosom|Ribosom", proteindata_unnorm$uniq_cluster_annotation) ,  "ribosome",
                                  
                                  ifelse (proteindata_unnorm$uniq_KO_pathway  ==  "Photosynthesis", "chloro", 
                                  ifelse (grepl("^clust_1820$|^clust_534$", proteindata_unnorm$uniq_cluster) ,  "chloro",
                                  ifelse (grepl("^clust_891$", proteindata_unnorm$uniq_cluster) ,  "other",

                                  ifelse (grepl("^clust_129$|^clust_2490$|^clust_1820$|^clust_534$",
                                                proteindata_unnorm$uniq_cluster) ,  "rubisco",

                                  ifelse (grepl("^clust_220$|^clust_115$|^clust_111$|^clust_2276$|^clust_2371$|^clust_277$|clust_555044|^clust_789$|^clust_2276, clust_2371$|^clust_100814, clust_115, clust_72123, clust_789$", proteindata_unnorm$uniq_cluster) ,  "histone",

                                  ifelse (grepl("Hsp|hsp", proteindata_unnorm$uniq_cluster_annotation) ,  "other",
                                  ifelse (grepl("ATP", proteindata_unnorm$uniq_cluster_annotation),  "other",
                                              proteindata_unnorm$uniq_grpnorm_compartment)))))))))) 


compartment <- proteindata_unnorm [c(4, 70, 5, 31:37, 40:43, 55:57)]   

compartment$taxon <- ifelse (compartment$F ==  "Polar-centric-Mediophyceae", "diatoms", 
                     ifelse (compartment$F ==  "Radial-centric-basal-Coscinodiscophyceae", "diatoms", 
                     ifelse (compartment$F ==  "Raphid-pennate", "diatoms", 
                     ifelse (compartment$F ==  "Araphid-pennate", "diatoms", 
                     ifelse (compartment$F ==  "Phaeocystaceae", "phaeocystis",
                     
                     ifelse (compartment$class_x ==  "Bacillariophyta_X", "diatoms", 
  
                     ifelse (compartment$C ==  "Ciliophora", "ciliates", 
                     ifelse (compartment$C ==  "Dinophyta", "dinophytes",
                     ifelse (compartment$C ==  "Haptophyta", "haptophytes",
                     
                     ifelse (compartment$B ==  "Eukaryota", "ambig_euk",
                   
                     ifelse (compartment$domain ==  "Bacteria", "bacteria",
                     ifelse (compartment$domain ==  "Viruses", "viruses_archea",
                     ifelse (compartment$domain ==  "Archaea", "viruses_archea",
                             
                     ifelse (compartment$domain ==  "Eukaryota", "other_euk",
                             
                     ifelse (compartment$domain ==  "ambiguous_domain", "ambig_domain",
                     ifelse (compartment$domain ==  "unassigned_tax", "unknown",
                     ifelse (compartment$domain ==  "not_in_annotation", "unknown",

                    "ALERT_somethingwentwrong!")))))))))))))))))

#check things to make sure it's all good
compartment2 <- compartment [c(18, 5:10)]   

#choose important columns for next steps
compartment2 <- compartment [c(2, 18, 11:17)]   


compartment2$compartment <- ifelse (compartment2$compartment ==  "rubisco", "rubisco", 
                     ifelse (compartment2$compartment ==  "ribosome", "ribosome", 
                     ifelse (compartment2$compartment ==  "chloro", "chloro", 
                     ifelse (compartment2$compartment ==  "histone", "histone",
                                        "other"))))

compartment2 <- melt(compartment2, id.vars=c("compartment", "taxon")) %>%
            rename(treatment = variable) %>%
            rename(unnnorm_abundance = value) 


#add the protein compartment abundances across samples and taxa
compartment_aggregated_group  <- aggregate(compartment2$unnnorm_abundance, 
                                 by=list(Category=compartment2$compartment, compartment2$treatment, compartment2$taxon), FUN=sum) %>%
                                 rename(unnorm_abundance = x) %>%
                        rename(compartment = Category) %>%
                        rename(treatment = Group.2) %>% 
                        rename(taxon = Group.3) 

#normalize by MS1 for each sample
compartment_aggregated_group  <- compartment_aggregated_group  %>%
                          group_by(treatment) %>% #normalize to treatment not by groups here
                          mutate(samplenorm_abundance = unnorm_abundance / sum(unnorm_abundance)) %>%
                          ungroup()

#get sample average
compartment_aggregated_group$treatment <- str_sub(compartment_aggregated_group$treatment, end=-3)


compartment_aggregated_group2  <- aggregate(compartment_aggregated_group$samplenorm_abundance, 
                                 by=list(Category=compartment_aggregated_group$compartment, compartment_aggregated_group$treatment,
                                         compartment_aggregated_group$taxon), FUN=mean) %>%
                                rename(mean_samplenormabundance = x) %>%
                                rename(compartment = Category) %>%
                                rename(treatment = Group.2) %>% 
                                rename(taxon = Group.3) 

#export these values
#this exported file will be used to make the indices and 'edges' file for circle pack plot
#write.csv(compartment_aggregated_group2, "ps117_circlepack_data_20241111.csv", row.names = F)
```

2024-11-15 -- circle pac, include all reps, not just T0, -- no need to re-run every time
```{r}
proteindata_unnorm$compartment <- ifelse (proteindata_unnorm$uniq_cluster  ==  "unassigned_cluster", "ribosome", 
                                  ifelse (proteindata_unnorm$uniq_grpnorm_compartment == "ribosomal", "ribosome",
                                  ifelse (grepl("ribosom|Ribosom", proteindata_unnorm$uniq_cluster_annotation) ,  "ribosome",
                                  
                                  ifelse (proteindata_unnorm$uniq_KO_pathway  ==  "Photosynthesis", "chloro", 
                                  ifelse (grepl("^clust_1820$|^clust_534$", proteindata_unnorm$uniq_cluster) ,  "chloro",
                                  ifelse (grepl("^clust_891$", proteindata_unnorm$uniq_cluster) ,  "other",

                                  ifelse (grepl("^clust_129$|^clust_2490$|^clust_1820$|^clust_534$",
                                                proteindata_unnorm$uniq_cluster) ,  "rubisco",

                                  ifelse (grepl("^clust_220$|^clust_115$|^clust_111$|^clust_2276$|^clust_2371$|^clust_277$|clust_555044|^clust_789$|^clust_2276, clust_2371$|^clust_100814, clust_115, clust_72123, clust_789$", proteindata_unnorm$uniq_cluster) ,  "histone",

                                  ifelse (grepl("Hsp|hsp", proteindata_unnorm$uniq_cluster_annotation) ,  "other",
                                  ifelse (grepl("ATP", proteindata_unnorm$uniq_cluster_annotation),  "other",
                                              proteindata_unnorm$uniq_grpnorm_compartment)))))))))) 


compartment <- proteindata_unnorm [c(4, 70, 5, 31:37, 40:69)]   

compartment$taxon <- ifelse (compartment$F ==  "Polar-centric-Mediophyceae", "diatoms", 
                     ifelse (compartment$F ==  "Radial-centric-basal-Coscinodiscophyceae", "diatoms", 
                     ifelse (compartment$F ==  "Raphid-pennate", "diatoms", 
                     ifelse (compartment$F ==  "Araphid-pennate", "diatoms", 
                     ifelse (compartment$F ==  "Phaeocystaceae", "phaeocystis",
                     
                     ifelse (compartment$class_x ==  "Bacillariophyta_X", "diatoms", 
  
                     ifelse (compartment$C ==  "Ciliophora", "ciliates", 
                     ifelse (compartment$C ==  "Dinophyta", "dinophytes",
                     ifelse (compartment$C ==  "Haptophyta", "haptophytes",
                     
                     ifelse (compartment$B ==  "Eukaryota", "ambig_euk",
                   
                     ifelse (compartment$domain ==  "Bacteria", "bacteria",
                     ifelse (compartment$domain ==  "Viruses", "viruses_archea",
                     ifelse (compartment$domain ==  "Archaea", "viruses_archea",
                             
                     ifelse (compartment$domain ==  "Eukaryota", "other_euk",
                             
                     ifelse (compartment$domain ==  "ambiguous_domain", "ambig_domain",
                     ifelse (compartment$domain ==  "unassigned_tax", "unknown",
                     ifelse (compartment$domain ==  "not_in_annotation", "unknown",

                    "ALERT_somethingwentwrong!")))))))))))))))))



#choose important columns for next steps
compartment2 <- compartment [c(2, 41, 11:40)]   


compartment2$compartment <- ifelse (compartment2$compartment ==  "rubisco", "rubisco", 
                     ifelse (compartment2$compartment ==  "ribosome", "ribosome", 
                     ifelse (compartment2$compartment ==  "chloro", "chloro", 
                     ifelse (compartment2$compartment ==  "histone", "histone",
                                        "other"))))

compartment2 <- melt(compartment2, id.vars=c("compartment", "taxon")) %>%
            rename(treatment = variable) %>%
            rename(unnnorm_abundance = value) 


compartment2$bioassay <- sub("_.*", "", compartment2$treatment)


#add the protein compartment abundances across samples and taxa
compartment_aggregated_group  <- aggregate(compartment2$unnnorm_abundance, 
                                 by=list(Category=compartment2$compartment, compartment2$bioassay, compartment2$taxon), FUN=sum) %>%
                                 rename(unnorm_abundance = x) %>%
                        rename(compartment = Category) %>%
                        rename(bioassay = Group.2) %>% 
                        rename(taxon = Group.3) 

#normalize by MS1 for each sample
compartment_aggregated_group  <- compartment_aggregated_group  %>%
                          group_by(bioassay) %>% #normalize to treatment not by groups here
                          mutate(samplenorm_abundance = unnorm_abundance / sum(unnorm_abundance)) %>%
                          ungroup()




compartment_aggregated_group2  <- aggregate(compartment_aggregated_group$samplenorm_abundance, 
                                 by=list(Category=compartment_aggregated_group$compartment, compartment_aggregated_group$bioassay,
                                         compartment_aggregated_group$taxon), FUN=mean) %>%
                                rename(mean_samplenormabundance = x) %>%
                                rename(compartment = Category) %>%
                                rename(bioassay = Group.2) %>% 
                                rename(taxon = Group.3) 

#export these values
#this exported file will be used to make the indices and 'edges' file for circle pack plot
#write.csv(compartment_aggregated_group2, "ps117_circlepack_data_alltreatments_20241115.csv", row.names = F)
```

2024-11-11 -- make the circle map 
           -- first we have to specify the nodes and leaves and edges etc. I manually did this on spreadsheets that I import here. 
```{r}
library(ggraph)
library(igraph)
library(data.tree)

edges <- read.csv ("ps117_circlepack_edges_20241111_nohardklor.csv")

vertices_offshore <- read.csv ("ps117_circlepack_offshoreallvertices_20241115_nohardklor.csv")
vertices_nearshore <- read.csv ("ps117_circlepack_nearshoreallvertices_20241115_nohardklor.csv") 

mygraph_offshore <- graph_from_data_frame(edges, vertices=vertices_offshore )
mygraph_nearshore <- graph_from_data_frame(edges, vertices=vertices_nearshore )



set.seed(2) #10, #1, 2
nearshorecircle <- ggraph(mygraph_nearshore, layout = 'circlepack', weight = size) +  #weight = size adjusts the circle size
      geom_node_circle(aes(fill = as.factor(colors), 
                         color = as.factor (circles),
                         group = depth))+
  #geom_node_circle()+ #this is the cicle around all the circles
  scale_fill_manual(values=c("allproteins" = "NA",
                            "unknown" = "black",
                            "ambig_domain" = "grey25",
                            "ambig_euk" = "#C0C0C0",
                            
                            "diatoms" = "NA",
                            "bacteria" = "#C0C0C0",
                            "other_haptophytes" = "#C0C0C0",
                            "ciliates" = "#C0C0C0",
                            "dinophytes" = "#C0C0C0",
                            
                            "other_euk" = "#C0C0C0",
                          
                            "phaeocystis" = "NA",
                           
                             "other" = "#DDCC77",
                             "histone" = "grey",
                             "ribosome" = "#56B4E9",
                             "photosyn" = "#009E73",
                             "rubisco" = "#E69F00"))+
  
         scale_color_manual(values=c("empty" = "NA",
                              "notempty" = "black")) +
         geom_node_text(aes(label = proteinlabel), size = 3) +
         geom_node_label(aes(label = taxonlabel), label.size = 0, fontface = "bold",
                         
            
                               size = ifelse(vertices_offshore$taxonlabel == "Unkn.", 1.8,
                               ifelse(vertices_offshore$taxonlabel == "Ambiguous domain", 2.3,
                               ifelse(vertices_offshore$taxonlabel == "Dinophytes", 2.3,
                               ifelse(vertices_offshore$taxonlabel == "Other eukaryotes", 2.3,
                               ifelse(vertices_offshore$taxonlabel == "Other hapto.", 2.3,
                               ifelse(vertices_offshore$taxonlabel == "Bacteria", 2,
                               ifelse(vertices_offshore$taxonlabel == "Ciliates", 2.3,
                                      3.3))))))),
                             
                        fill = ifelse(vertices_offshore$taxonlabel == "empty", "NA", "white"), 
                        color = ifelse(vertices_offshore$taxonlabel == "empty", "NA", "black"), 
                 
                        nudge_y = ifelse (vertices_offshore$taxonlabel == "Ambiguous eukaryotes", 0.0, 
                                  ifelse (vertices_offshore$taxonlabel == "Diatoms", -0.08, 0))) +

         theme_void() + 
         theme(legend.position = "none")+
            coord_fixed (ratio = 1)

nearshorecircle










set.seed(3) #10, #1, 2
offshorecircle <- ggraph(mygraph_offshore, layout = 'circlepack', weight = size) +  #weight = size adjusts the circle size
      geom_node_circle(aes(fill = as.factor(colors), 
                         color = as.factor (circles),
                         group = depth))+
  #geom_node_circle()+ #this is the cicle around all the circles
  scale_fill_manual(values=c("allproteins" = "NA",
                            "unknown" = "black",
                            "ambig_domain" = "grey25",
                            "ambig_euk" = "#C0C0C0",
                            
                            "diatoms" = "NA",
                            "bacteria" = "#C0C0C0",
                            "other_haptophytes" = "#C0C0C0",
                            "ciliates" = "#C0C0C0",
                            "dinophytes" = "#C0C0C0",
                            
                            "other_euk" = "#C0C0C0",
                          
                            "phaeocystis" = "NA",
                           
                             "other" = "#DDCC77",
                             "histone" = "grey",
                             "ribosome" = "#56B4E9",
                             "photosyn" = "#009E73",
                             "rubisco" = "#E69F00"))+
  
         scale_color_manual(values=c("empty" = "NA",
                              "notempty" = "black")) +
         geom_node_text(aes(label = proteinlabel), size = 3) +
         geom_node_label(aes(label = taxonlabel), label.size = 0, fontface = "bold",
                         
            
                               size = ifelse(vertices_offshore$taxonlabel == "Unkn.", 1.8,
                               ifelse(vertices_offshore$taxonlabel == "Ambiguous domain", 2.3,
                               ifelse(vertices_offshore$taxonlabel == "Dinophytes", 2.3,
                               ifelse(vertices_offshore$taxonlabel == "Other eukaryotes", 2.3,
                               ifelse(vertices_offshore$taxonlabel == "Other hapto.", 2.3,
                               ifelse(vertices_offshore$taxonlabel == "Bacteria", 2,
                               ifelse(vertices_offshore$taxonlabel == "Ciliates", 2.3,
                                      3.3))))))),
                             
                        fill = ifelse(vertices_offshore$taxonlabel == "empty", "NA", "white"), 
                        color = ifelse(vertices_offshore$taxonlabel == "empty", "NA", "black"), 
                 
                        nudge_y = ifelse (vertices_offshore$taxonlabel == "Ambiguous eukaryotes", 0.0, 
                                  ifelse (vertices_offshore$taxonlabel == "Diatoms", -0.1, 0.0))) +

         theme_void() + 
         theme(legend.position = "none")+
            coord_fixed (ratio = 1)

offshorecircle



```

2024-12-01  --  figure 1 construction
```{r}
#no3 <- x_vs_time(NO3_uM, NO[3]^"-"~~(mu*M)) + ylim(5, 30) + theme (axis.text.x=element_blank()) 
#po4 <- x_vs_time(PO4_uM, PO[4]^"3-"~~(mu*M))  + ylim(0, 2) +theme(strip.text.x = element_blank())
#pon <- T0_T8(nonproteindata, PON_uM, PON~~(mu*M)) + theme (axis.text.x=element_blank())
#poc <- T0_T8(nonproteindata, POC_uM, POC~~(mu*M)) +theme(strip.text.x = element_blank())

right <- plot_grid(NULL, 
          NULL,
          
          ncol = 1, align = "v", labels = c ("(C)", "(D)"))


left <- plot_grid (elementalratio, peptidenmds,   
           ncol = 1, 
           rel_heights = c(1, 1), rel_widths = c(1,1), align = "v", labels = c ("(A)", "(B)"))
left


png("../figures/ps117_figure1_20241201.png", width= 11, height= 9, units = "in",  res = 600) 
plot_grid (left, NULL, right, 
           ncol = 3, 
           rel_widths = c(0.75, 0.05,  0.7), align = "v") + 
 # draw_plot (map, x = 0.055, y = 0.52, width = 0.35, height = 0.35) + 
  draw_plot (nearshorecircle, x = 0.5, y = 0.53, width = 0.5, height = 0.5) + 
  draw_plot (offshorecircle, x = 0.48, y = -0.0, width = 0.5, height = 0.5) + 
  draw_label("Neashore", x = 0.67, y = 0.95, size = 14, hjust = 1, fontface = "bold") + 
  draw_label("Offshore", x = 0.67, y = 0.4, size = 14, hjust = 1, fontface = "bold") 
  
dev.off()
#svg("../figures/ps117_figure1new_20240830.svg", width= 13, height= 10) 



```

2024-12-01  --  figure 1 construction_v2
```{r}
left <- plot_grid (NULL, elementalratio,   
           ncol = 1, 
           rel_heights = c(0.7, 1), rel_widths = c(1,1), align = "v", labels = c ("(A)", "(B)"))
left


right <- plot_grid(peptidenmds, 
          NULL,
           rel_heights = c(0.8, 1),
          
          ncol = 1, align = "v", labels = c ("(C)", "(D)"))


right


png("../figures/ps117_figure1_20241201.png", width= 11, height= 9, units = "in",  res = 600) 
plot_grid (left, NULL, right, 
           ncol = 3, 
           rel_widths = c(0.6, 0.05,  0.7), align = "v") +
  draw_plot (map, x = 0.025, y = 0.57, width = 0.425, height = 0.425) + 
  draw_plot (offshorecircle, x = 0.39, y = 0.17, width = 0.42, height = 0.42) +
  draw_label("Offshore", x = 0.6, y = 0.18, size = 14, hjust = 1, fontface = "bold") +

  draw_plot (nearshorecircle, x = 0.64, y = -0.02, width = 0.42, height = 0.42)  +
  draw_label("Neashore", x = 0.96, y = 0.38, size = 14, hjust = 1, fontface = "bold") 
  
dev.off()
#svg("../figures/ps117_figure1new_20240830.svg", width= 13, height= 10) 



```









2024-05-25 -- normalized diatom and Phaeocystis protein abundance 
```{r}
allgroups <- proteindata_unnorm [c(32,34,35, 37, 40:69)]

allgroups <- melt(allgroups, id.vars=c("domain", "C","class", "F")) %>%
             rename(treatment = variable) %>%
             rename(unnorm_abundance = value )
sum(allgroups$unnorm_abundance) #this should equal 1.893E13 (total abundance of all peptides for all samples)

#normalize the peptide abundances for within each sample, irrespective of taxonomy
allgroups <- allgroups %>% #this creates a new column of normalized protein abundances grouped by treatment   
  group_by_at(c("treatment")) %>% #this can be normalized to treatment and/or class
  mutate(norm_abundance = unnorm_abundance / sum(unnorm_abundance)) %>% 
  ungroup()
sum(allgroups$norm_abundance) #this should = 30 after normalization (number of replicates)

allgroups$taxon <- ifelse (grepl("pennate", allgroups$F),  "Pennate diatom", 
                    ifelse (grepl("centric", allgroups$F), "Centric diatom", 
                    ifelse (grepl("Phaeocystaceae", allgroups$F), "Phaeocystis",
                    allgroups$F))) 

allgroups_aggregated <- aggregate(allgroups$norm_abundance, by=list( allgroups$treatment, allgroups$taxon), FUN=sum) %>%
                  rename(norm_abundance = x) %>%
                  rename(treatment = Group.1) %>%
                  rename(taxon = Group.2) %>% 
                  separate(treatment, into =  c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"))  

allgroups_aggregated <- allgroups_aggregated %>%
  mutate(tempfe_treatment = paste(temperature_treatment, iron_treatment, sep = "_")) %>%
  mutate(tempfe_treatment = str_replace(tempfe_treatment, "T0_T0", "T0"))

T0_T8(filter(allgroups_aggregated, grepl ("Pennate", taxon)), norm_abundance, "Pennate diatom fraction") + theme(axis.text.x = element_blank())

pennategroup <- T0_T8(filter(allgroups_aggregated, grepl ("Pennate", taxon)), norm_abundance, "Pennate diatom fraction") + theme(axis.text.x = element_blank())

centricgroup <- T0_T8(filter(allgroups_aggregated, grepl ("Centric", taxon)), norm_abundance, "Centric diatom fraction") + theme(axis.text.x = element_blank()) + theme(strip.text.x = element_blank())

phaeogroup <- T0_T8(filter(allgroups_aggregated, grepl ("Phaeocystis", taxon)), norm_abundance, "Phaeocystis fraction") +theme(strip.text.x = element_blank()) 
```







#2024-08-27 -- Figure 1/2 leftoevers 
            -- (this will need to get re-shuffled, NMDS is in sup. now). is this pennate diatoms of total protein? or normalized to total diaotm?
          
2024-05-25 -- protein heatmap
```{r}
allgroups <- proteindata_unnorm [c(32,34,35, 40:69)]

allgroups <- melt(allgroups, id.vars=c("domain", "C","class")) %>%
             rename(treatment = variable) %>%
             rename(unnorm_abundance = value )
sum(allgroups$unnorm_abundance) #this should equal 1.893E13 (total abundance of all peptides for all samples)

#normalize the peptide abundances for within each sample, irrespective of taxonomy
allgroups <- allgroups %>% #this creates a new column of normalized protein abundances grouped by treatment   
  group_by(treatment) %>% 
  mutate(norm_abundance = unnorm_abundance / sum(unnorm_abundance)) %>% 
  ungroup()
sum(allgroups$norm_abundance) #this should = 30 after normalization (number of replicates)

allgroups$class_a <- ifelse (allgroups$class ==  "Bacillariophyta", "Bacillariophyta", 
           
                    ifelse (allgroups$C ==  "Dinophyta", "Dinophyta",
                                    
                    ifelse (allgroups$class ==  "Haptophyta", "Haptophyta",
                    ifelse (allgroups$class ==  "Prymnesiophyceae", "Haptophyta",
                    ifelse (allgroups$class ==  "Pavlovophyceae", "Haptophyta",
                                    
                    ifelse (allgroups$class ==  "Cryptophyceae", "Cryptophyta",
                    ifelse (allgroups$class ==  "Cryptophyta", "Cryptophyta",
                                    
                    ifelse (allgroups$class ==  "Pelagophyceae", "Pelagophyceae",
                                    
                    ifelse (allgroups$class ==  "Chlorophyceae", "Chlorophyta",
                    ifelse (allgroups$class ==  "Chlorophyta", "Chlorophyta",
                    ifelse (allgroups$class ==  "Pyramimonadales", "Chlorophyta",
                    ifelse (allgroups$class ==  "Prasinophyceae", "Chlorophyta",
                            
                    ifelse (allgroups$class ==  "Spirotrichea", "Ciliophora",
                    ifelse (allgroups$class ==  "Prostomatea", "Ciliophora",  
                    ifelse (allgroups$class ==  "Oligohymenophorea", "Ciliophora",
                    ifelse (allgroups$class ==  "Litostomatea", "Ciliophora",   
                    ifelse (allgroups$class ==  "Heterotrichea", "Ciliophora", 
                    ifelse (allgroups$class ==  "Colpodea", "Ciliophora", 
                    ifelse (allgroups$class ==  "Ciliophora", "Ciliophora", 
                            
                    ifelse (allgroups$class == "Eukaryota", "Other Eukaryota",
                    ifelse (allgroups$domain == "Bacteria", "Bacteria",
                    ifelse (allgroups$domain == "Viruses", "Viruses",
                    "Other")))))))))))))))))))))) #this takes a second

allgroups_aggregated <- aggregate(allgroups$norm_abundance, by=list( allgroups$treatment, allgroups$class_a), FUN=sum) %>%
                  rename(norm_abundance = x) %>%
                  rename(treatment = Group.1) %>%
                  rename(class = Group.2) %>% 
                  separate(treatment, into =  c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"))  

allgroups_aggregated$class <- reorder (allgroups_aggregated$class , allgroups_aggregated$norm_abundance)

allgroups_aggregated$tempfe_treatment<- paste (allgroups_aggregated$temperature_treatment, allgroups_aggregated$iron_treatment, sep = "_")

allgroups_aggregated_means <- aggregate(allgroups_aggregated$norm_abundance, by=list (allgroups_aggregated$tempfe_treatment, allgroups_aggregated$bioassay, allgroups_aggregated$class), FUN=mean) %>%
  rename (treatment = Group.1, 
          bioassay = Group.2, 
          class = Group.3, 
          norm_abundance = x)

allgroups_aggregated_means$treatment2 <- ifelse (allgroups_aggregated_means$treatment ==  "T0_T0", "T0",
                                                 allgroups_aggregated_means$treatment)

fig2_proteintopmeans <- ggplot(filter(allgroups_aggregated_means, grepl ("Bacillar|Hapto|Bacteria", class)), 
                               aes (x = treatment2,  y = class, fill= norm_abundance))+ 

  geom_tile(color = "black", lwd = 0.7) +
  scale_fill_viridis_c(values = rescale(c(0.00001, 0.2, 0.55)), option= "E")+ #H D, E 
   labs(x = "", y = "") +
   facet_grid(.~ bioassay)+
   scale_x_discrete(limits=c("T0", "LT_noFe","LT_Fe", "HT_noFe", "HT_Fe"), expand = expansion(mult = 0.01))+
   theme(axis.text.x=element_blank(),
         axis.text.y=element_text(face = "bold", size = 10, color = "black"), 
         axis.ticks  = element_blank(), 
         strip.background =element_rect(fill="white", color = "black"),
         strip.text.x = element_text(size = 10,  face = "bold"),   
         legend.box.margin=margin(l=-10), 
         legend.title = element_blank())+ 
     scale_y_discrete(expand = expansion(mult = 0.2))
fig2_proteintopmeans

fig2_proteinbottommeans <- ggplot(filter(allgroups_aggregated_means, grepl ("Dinophy|Ciliop|Cryptophy|Pelagophy|Chlorophy|Viruses", class)), 
       aes (x = treatment2,  y = class, fill= norm_abundance))+ 
  geom_tile(color = "black", lwd = 0.7) +
  scale_fill_viridis_c(values = rescale(c(0.00001, 0.004, 0.034)), option= "E")+ #H D, E 
  # scale_fill_gradientn(colours = c("white","#999933","darkgreen"), 
  #                      values = rescale(c(0.0,0.002,0.034)),
  #                      guide = "colorbar", limits=c(0.000001,0.034),
  #                      name = "")+
   labs(x = "", y = "")+
   scale_x_discrete(limits=c("T0", "LT_noFe","LT_Fe", "HT_noFe", "HT_Fe"), 
                      labels = c("T0", "Ctrl", "+Fe", "HT", "HT +Fe"))+
   theme(axis.text =element_text(face = "bold", size = 9, color = "black"),
        axis.text.y=element_text(face = "bold", size = 10, color = "black"),
        panel.border = element_blank(),
        panel.background = element_blank())+ 
      #  axis.ticks  = element_blank()) +
  facet_grid(.~ bioassay)+
  theme(strip.background =element_rect(fill="white"), 
        strip.text = element_blank(), 
        plot.margin=margin(b=0, unit="cm"), 
        legend.box.margin=margin(l=-10), 
        legend.title = element_blank())  + 
  scale_y_discrete(expand = c(0.002, 0.002))
   #theme(plot.margin = margin(t = -5, r = 0, l = 0, unit = "mm"))
fig2_proteinbottommeans


```


Figure 2 misc. and oldies 
2023-05-07 -- after making all the NMDS analysis and combining the results, 
              make nmds plots 2023-04-21 (see chunks fig2nmdsA,fig2nmdsB,fig2nmdsC... below for calculations)
 2024-03-26 -- Voronoi tree map. #maybe for later
           -- one thing that would be good is to export the treemap coordinates and plot in ggplot. 
```{r}
#https://cran.r-project.org/web/packages/WeightedTreemaps/WeightedTreemaps.pdf
library(WeightedTreemaps)

# load example data
data(mtcars)
mtcars$car_name = gsub(" ", "\n", row.names(mtcars))

mtcars$newname <- ifelse(mtcars$gear == 4, "nitrogen metabolism", 
                  ifelse(mtcars$gear == 3, "information processing/translation", 
                         "photosynthesis"))

#Generate the treemap. It will return a list of polygons and metadata. The columns of the data frame that are used as levels of the treemap need to be specified. Different parameters like the initial shape, or the maximum number of iterations are optional.

treempadata <- filter (allgroups_aggregated_means, grepl ("T0", treatment), grepl ("Offshore", bioassay))

# generate treemap; set seed to obtain same pattern every time
tm <- voronoiTreemap(
  data = treempadata,
  levels = c("class"),
  cell_size = "norm_abundance",
  shape = "circle",
  seed = 123
  
)
#Draw the treemap
drawTreemap(tm,label_size = 3, 
            label_color = "black", 
            label_level = 1,
            border_size = 2, 
           # border_level = 2,
            border_color = "black")
```


- Pigments heatmap
```{r}
pigment_fract <- nonproteindata [c(15,62:67)] #use the 'fract' columns of all groups 
pigment_fract <- na.omit(pigment_fract)
pigment_fract2 <- melt(pigment_fract, id.vars=c("full_treatment_name")) %>%
            rename(taxon_group = variable) %>%
            rename(fract = value )

pigment_fract2 <- pigment_fract2 %>% separate(full_treatment_name, c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate")) 

pigment_fract2$fulltreatment <- paste(pigment_fract2$temperature_treatment,pigment_fract2$iron_treatment, sep = "_")

pigment_fract2 <- pigment_fract2 %>% separate(taxon_group, c("bla", "taxon_group")) 

pigment_fract2$taxon_group <- reorder(pigment_fract2$taxon_group, pigment_fract2$fract) #this reorders the data so when it's plotted, the most abundant fraction is on the bottom #only wroks when there are no NA's

#dev.off()
#pigmentheatmap
pigment_means <- aggregate(pigment_fract2$fract, by=list( pigment_fract2$fulltreatment, pigment_fract2$bioassay, pigment_fract2$taxon_group), FUN=mean) %>%
  rename (fulltreatment = Group.1, 
          bioassay = Group.2, 
          taxon_group = Group.3, 
          fract = x)
pigment_means$taxon_group2 <-   ifelse (pigment_means$taxon_group ==  "diatoms", "Bacillariophyta", 
                   ifelse (pigment_means$taxon_group ==  "haptophytes", "Haptophyta", 
                   ifelse (pigment_means$taxon_group ==  "pelagophytes", "Pelagophyceae", 
                   ifelse (pigment_means$taxon_group ==  "cryptophytes", "Cryptophyta",
                   ifelse (pigment_means$taxon_group ==  "chlorophytes", "Chlorophyta", 
                   ifelse (pigment_means$taxon_group ==  "dinoflagellates", "Dinophyta",
                           "else" 
                    ))))))

pigment_means$taxon_group2 <- reorder(pigment_means$taxon_group2, pigment_means$fract) #this reorders the data so when it's plotted, the most abundant fraction is on the bottom #only wroks when there are no NA's

pigment_means$treatment2 <- ifelse (pigment_means$fulltreatment ==  "LT_noFe", "LT_noFe", 
                            ifelse (pigment_means$fulltreatment ==  "LT_Fe", "LT_Fe",
                            ifelse (pigment_means$fulltreatment ==  "HT_noFe", "HT_noFe", 
                            ifelse (pigment_means$fulltreatment ==  "HT_Fe", "HT_Fe",
                            "T0"))))

#x = factor(treatment2, level = tempFe_order)

fig2_pigmentmeans <- ggplot(pigment_means, aes(x = treatment2,
                                              y = taxon_group2, 
                                              fill = fract)) +  
    geom_tile(color = "black", lwd = 0.7)+
    #coord_fixed()+ #this makes it all squares 
    scale_fill_gradientn(colours = c("deepskyblue2","white","darkred"),
                         values = rescale(c(0.0000001,0.35,.9)),
                         guide = "colorbar", limits=c(0.0000001,0.9),
                       name = "")+
# guides (fill = guide_legend (title.position = "right"))+
  
    labs(x = "", y =  "")+
  ggtitle("Chl-a fraction") +
  theme (plot.title = element_text(size = 12, vjust = -8, hjust = 0.01))+
    scale_x_discrete(limits=c("T0", "LT_noFe","LT_Fe", "HT_noFe", "HT_Fe"))+
    facet_grid(.~ bioassay, scales = "free_x") + 
    theme(axis.text.x= element_blank(),
        axis.text.y=element_text(face = "bold", size = 11, color = "black"),
        panel.border = element_blank(),
        panel.background = element_blank(), 
        axis.ticks  = element_blank()) +
    theme(strip.background =element_rect(fill="white"), 
        strip.text = element_text(size = 11, face = "bold"),
        legend.box.margin=margin(l=-10))+
        theme(plot.margin=grid::unit(c(-2.5,0,0,0), "mm")) 

fig2_pigmentmeans



#stacked bar chart
# ggplot (pigment_fract2, aes(
#   x = fulltreatment,
#   y = fract, 
#   fill = taxon_group)) +  
#   geom_bar(position = "fill", stat = "identity") + #position = "fill" does percentage, stack does raw values
#   facet_grid(.~ bioassay)+
# theme_bw()+
#    theme (axis.text.x = element_text(angle = 90)) +
#    ylab ("Contribution to total pigment (fraction)")+
#     xlab ("")+
#   scale_x_discrete(limits = c('T0_T0', 'LT_noFe', 'LT_Fe', 'HT_noFe', 'HT_Fe'))+
#    scale_fill_manual(name="",
#                      values = c("diatoms" = "#117733",
#                                 "dinoflagellates" = "#0072B2",
#                                 "haptophytes" = "#CC79A7",
#                                 "prasinophytes" = "#88CCEE",
#                                 "cryptophytes" = "grey", #000000",
#                                 "pelagophytes" = "#E69F00"))+
#   theme(axis.text.x=element_text(face = "bold", size = 18, color = "black", vjust = 1, hjust = 1, angle = 45),
#         axis.title.y=element_text(size=20,face = "bold", color = "black"), 
#         axis.text.y=element_text(face = "bold", size = 15, color = "black"),
#         strip.background =element_rect(fill = "white"),
#         legend.text = element_text(size = 15), 
#         strip.text.x = element_text(size = 18, face = "bold"))
```

frag pseudo nitz etc. 2023-06-24
```{r}
allgroups <- proteindata_unnorm [c(37, 38, 40:69)]

allgroups <- melt(allgroups, id.vars=c("F", "G")) %>%
            rename(treatment = variable) %>%
            rename(unnnorm_abundance = value) 
            #separate(treatment, into=c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"), sep = "_")

allgroups$F <- ifelse (allgroups$F ==  "Polar-centric-Mediophyceae", "Centric diatom", 
                 ifelse (allgroups$F ==  "Radial-centric-basal-Coscinodiscophyceae", "Centric diatom", 
               #  ifelse (compartment$F ==  "Bacillariophyta_X", "Bacillariophyta", 
                 ifelse (allgroups$F ==  "Raphid-pennate", "Pennate diatom", 
                 ifelse (allgroups$F ==  "Araphid-pennate", "Pennate diatom", 
                 ifelse (allgroups$F ==  "Phaeocystaceae", "Phaeocystaceae",
                                          allgroups$F)))))

#normalize to group 
aggregated_group  <- aggregate(allgroups$unnnorm_abundance, 
                                  by=list(Category=allgroups$G, allgroups$treatment, allgroups$F), FUN=sum) %>%
                        rename(unnorm_abundance = x) %>%
                        rename(taxon_G = Category) %>%
                        rename(treatment = Group.2) %>% 
                        rename(taxon_F = Group.3) 

aggregated_group  <- aggregated_group  %>% 
                          group_by(treatment,taxon_F) %>% 
                          mutate(groupnorm_abundance = unnorm_abundance / sum(unnorm_abundance)) %>% 
                          ungroup()

aggregated_group <- aggregated_group %>% 
  separate(treatment, into=c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"), sep = "_") %>%
  mutate(tempfe_treatment = paste(temperature_treatment, iron_treatment, sep = "_")) %>%
  mutate(tempfe_treatment = str_replace(tempfe_treatment, "T0_T0", "T0"))



protein_taxgroups <- function (taxonomicgroup, yaxislabel) { 
ggplot(filter (allgroups_aggregated, grepl (taxonomicgroup, taxon)), 
   aes(x= tempfe_treatment,
       y= norm_abundance, 
       color = temperature_treatment, 
       shape = iron_treatment))+
  facet_wrap(~bioassay)+           
  stat_summary(fun = mean, geom = "point", size = 3.5, stroke = 1.2, alpha = 0.5)+
  stat_summary(fun.data =  mean_se, geom = "errorbar", show.legend = F, linewidth = 1, width = 0.1)+
  scale_x_discrete(limits=c("T0", "LT_noFe","LT_Fe", "HT_noFe", "HT_Fe"), 
                   labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"))+
      scale_color_manual(name="Temperature",
                     breaks = c("LT","HT"),
                     labels = c("LT", "HT"),
                     values = c('black', 'darkorange2','black')) +
  scale_shape_manual(name = "Iron",
                       breaks = c("Fe","noFe", "T0"),
                       labels = c("Fe", "noFe", "T0"),
                       values = c(19,1,5)) +
  xlab(expression ("")) +
  ylab (substitute(yaxislabel))+ 
  theme_bw()+
  theme(axis.text.y =element_text(face = "bold", size = 11, color = "black"),
        axis.text.x=element_text(face = "bold", size = 11, color = "black", vjust = 1.2, hjust = 1,  angle = 45),
        axis.title.y =element_text(size = 12,face = "bold", color = "black"), 
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 12, face = "bold"),
        legend.position = "none") +
    ylim (0,NA)
}

pennategroup <- protein_taxgroups("Pennate", "Pennate diatoms") + theme(axis.text.x = element_blank()) + ggtitle("Protein fraction")
centricgroup <- protein_taxgroups("Centric", "Centric diatoms") + theme(axis.text.x = element_blank()) 
```

groups heatmap
```{r}

###
#try to do this in a heatmap 

compartment_aggregated_group_means <- aggregate(compartment_aggregated_group$groupnorm_abundance, by=list (compartment_aggregated_group$bioassay, compartment_aggregated_group$compartment, compartment_aggregated_group$timepoint, compartment_aggregated_group$temperature_treatment, compartment_aggregated_group$iron_treatment, compartment_aggregated_group$taxon_F), FUN=mean)
  
  rename (treatment = Group.1, 
          bioassay = Group.2, 
          class = Group.3, 
          norm_abundance = x)


compartment_aggregated_group_means$treatment2 <- paste (compartment_aggregated_group_means$Group.4, compartment_aggregated_group_means$Group.5, sep = "_")


  
ggplot(filter(compartment_aggregated_group_means, grepl ("Bacilla|Phaeocy", Group.6) & grepl ("^chloro$", Group.2)), 
       aes (x = treatment2,  y = Group.6, fill= x))+ 
  geom_tile(color = "black", lwd = 0.7) +
  scale_fill_gradientn(colours = c("deepskyblue2","white","darkred"),
                       values = rescale(c(0.05,0.4,.72)),
                       guide = "colorbar", limits=c(0.05,0.72),
                       name = "Fraction")+
   labs(x = "", y = "") +
   scale_x_discrete(limits=c("T0_T0", "LT_noFe","LT_Fe", "HT_noFe", "HT_Fe"))+
   theme(axis.text.x=element_text(face = "bold", size = 11, color = "black", vjust = 1.2, hjust = 1,  angle = 45),
        axis.text.y=element_text(face = "bold", size = 11, color = "black"),
        panel.border = element_blank(),
        panel.background = element_blank(), 
        axis.ticks  = element_blank()) +
  facet_grid(.~ Group.1)+
  theme(strip.background =element_rect(fill="white"), 
        strip.text = element_blank(), 
        plot.margin=margin(b=0, unit="cm"), 
        legend.box.margin=margin(l=-10)) + 
  theme(plot.margin=grid::unit(c(0,0,0,0), "mm")) 

```

Fig 2 legend  2023-03-13.
```{r}
nmdsplotdata$tempfe_treatment <- paste (nmdsplotdata$temperature_treatment, nmdsplotdata$iron_treatment, sep = "_")

legendA <- ggplot(nmdsplotdata, 
       aes(x=NMDS1, 
       y=NMDS2,
       color = tempfe_treatment, 
       shape = tempfe_treatment))+
  geom_point(size = 4, alpha =0.65, stroke =1.5)+
  scale_color_manual(breaks = c("T0_T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c('purple', 'black', 'black','darkorange2', "darkorange2")) +
  scale_shape_manual(breaks = c("T0_T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c(1, 1, 19,1, 19)) +
  theme_bw()+
  theme(legend.title = element_blank(),
      #  legend.position = "top",
        #legend.justification = "left",
        legend.text = element_text(size = 10), 
        legend.margin  = unit(c(0, 0, 0, 0), "cm"))
legendA

legendA <- get_legend(legendA)
fig2_legendA <- as_ggplot(legendA)
fig2_legendA

###

legendB <- ggplot(nmdsplotdata, 
       aes(x=NMDS1, 
       y=NMDS2,
       shape = bioassay))+
  geom_point(size = 4, alpha =0.65, stroke =1.5)+
  scale_shape_manual(breaks = c("BA1", "BA2"),
                      labels = c("BA1", "BA2"),
                      values = c(1, 2)) +
  theme_bw()+
  theme(legend.title = element_blank(),
        legend.text = element_text(size = 10), 
        #legend.position = "top",
        #legend.justification = "left",
        legend.margin  = unit(c(0, 0, 0, 0), "cm"))
 
legendB <- get_legend(legendB)
fig2_legendB <- as_ggplot(legendB)
fig2_legendB

###
#fig2legend <- ggdraw (legendB) + draw_plot (legendA, x = 0.32 , y = 0) 
#fig2legend <- ggdraw (legendA) + draw_plot (legendB, x = -0.1, y = 0.13) 
#fig2legend
```

Figure 2 Plot : 2023-03-13
```{r}
nmds <- plot_grid (pigment_nmds, fc_nmds, peptide_nmds,   NULL,
                  ncol=1, 
                  rel_heights = c(1,1,1,0.12),
                  labels = c("A", "B","C", ""))

heatmaps <- plot_grid(fig2_pigmentmeans, fig2_proteintopmeans, fig2_proteinbottommeans, 
          ncol=1,
          labels = c("D","E", ""), 
          align = "v", 
          rel_heights = c(1.25,0.65,1.1))

# photogroups <- plot_grid (pennatephotofraction, centricphotofraction, phaeophotofraction, 
#                           align = 'v', 
#                           ncol = 1)
# 
# y.grobphotogroups <- textGrob("Photosynthetic protein fraction", 
#                    gp=gpar( fontsize=16, fontface="bold"), rot=90)
# 
# allphotogroups <- grid.arrange(arrangeGrob(photogroups, left = y.grobphotogroups))
# 
# labelz <- plot_grid (NULL, NULL, NULL, ncol=1, 
#                     labels = c("F", "G", "H"))
# allphotogroups_labeled <- plot_grid (labelz, allphotogroups,
#            ncol = 2, rel_widths =c(0.05, 1))
# 
# allphotogroups_labeled
# 
# 
# groups <-  plot_grid(pennategroup, centricgroup, phaeogroup,
#             align='v',
#             ncol =1, 
#             rel_heights = c(1,1,1.25))                    

#y.grob <- textGrob("Fraction of total protein abundance", 
                   #gp=gpar( fontsize=16, fontface="bold"), rot=90, hjust = 0.5) #fontface="bold",
#allgroups <- grid.arrange(arrangeGrob(groups, left = y.grob))
#labelzz <- plot_grid (NULL, NULL, NULL, ncol=1, 
                    #labels = c("I", "J", "K"))

#allgroup_labeled <- plot_grid (labelzz, allgroups,
          # ncol = 2, rel_widths =c(0.05, 1))
#allgroup_labeled

#plot_grid (allphotogroups_labeled, allgroup_labeled, ncol = 1, 
                         #align = "v",  axis = "l")
           
rightpanel <- plot_grid (NULL, NULL, 
           #NULL, pennatephotofraction, 
           #NULL, centricphotofraction, 
           #NULL, phaeophotofraction, 
           NULL, pennategroup, 
           NULL, centricgroup, 
           NULL,  phaeogroup,
           ncol =2, #0.05
           rel_widths = c(0.05,1), rel_heights = c(0.05, 1.1,1,1.2), 
          labels = c("", "",
                     "F", "", 
                     "G", "",
                     "H", ""),
                     #"I", "",
                    # "J", "",  "K", "")
                     align = "v")

#png("ps117_figure2_highres_20230421.png", width=20, height=10, units="in", res=900) 
fig2 <- plot_grid (nmds, NULL, heatmaps, NULL, rightpanel, 
            ncol = 5,
            rel_widths = c(0.7, 0.05, 1.2, 0.05, 1),
            rel_heights =  c(1,1,1,1,1),
            axis = "bt",
            align = "vh")
#dev.off()

#png("ps117_figure2_highres_20230508.png", width=20, height=10, units="in", res=900) 
ggdraw(fig2) + draw_plot (fig2_legendB, x = 0.245, y = 0.45, width = 0.0001, height = 0.001) +
               draw_plot (fig2_legendA, x = 0.25, y = 0.35, width = 0.0001, height = 0.001) 
                
#dev.off()
```

Normalized diatom and Phaeocystis protein abundance - 2023-07-06
```{r}
allgroups <- proteindata_unnorm [c(32,34,35, 37, 40:69)]

allgroups <- melt(allgroups, id.vars=c("domain", "C","class", "F")) %>%
             rename(treatment = variable) %>%
             rename(unnorm_abundance = value )
sum(allgroups$unnorm_abundance) #this should equal 1.893E13 (total abundance of all peptides for all samples)

#normalize the peptide abundances for within each sample, irrespective of taxonomy
allgroups <- allgroups %>% #this creates a new column of normalized protein abundances grouped by treatment   
  group_by_at(c("treatment")) %>% #this can be normalized to treatment and/or class
  mutate(norm_abundance = unnorm_abundance / sum(unnorm_abundance)) %>% 
  ungroup()
sum(allgroups$norm_abundance) #this should = 30 after normalization (number of replicates)

allgroups$taxon <- ifelse (grepl("pennate", allgroups$F),  "Pennate diatom", 
                    ifelse (grepl("centric", allgroups$F), "Centric diatom", 
                    ifelse (grepl("Phaeocystaceae", allgroups$F), "Phaeocystis",
                    allgroups$F))) 

allgroups_aggregated <- aggregate(allgroups$norm_abundance, by=list( allgroups$treatment, allgroups$taxon), FUN=sum) %>%
                  rename(norm_abundance = x) %>%
                  rename(treatment = Group.1) %>%
                  rename(taxon = Group.2) %>% 
                  separate(treatment, into =  c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"))  

allgroups_aggregated <- allgroups_aggregated %>%
  mutate(tempfe_treatment = paste(temperature_treatment, iron_treatment, sep = "_")) %>%
  mutate(tempfe_treatment = str_replace(tempfe_treatment, "T0_T0", "T0"))


#allgroupsstats <- filter(allgroups_aggregated, grepl("Pennate|Centric|Haptophyta|Bacillariophyta|Phaeocys", taxon))
#write.csv(allgroupsstats, "pennate-centric-hapt_stat.csv", row.names = F)


protein_taxgroups <- function (taxonomicgroup, yaxislabel) { 
ggplot(filter (allgroups_aggregated, grepl (taxonomicgroup, taxon)), 
   aes(x= tempfe_treatment,
       y= norm_abundance, 
       color = temperature_treatment, 
       shape = iron_treatment))+
  facet_wrap(~bioassay)+           
  stat_summary(fun = mean, geom = "point", size = 3.5, stroke = 1.2, alpha = 0.5)+
  stat_summary(fun.data =  mean_se, geom = "errorbar", show.legend = F, linewidth = 1, width = 0.1)+
  scale_x_discrete(limits=c("T0", "LT_noFe","LT_Fe", "HT_noFe", "HT_Fe"), 
                   labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"))+
      scale_color_manual(name="Temperature",
                     breaks = c("LT","HT"),
                     labels = c("LT", "HT"),
                     values = c('black', 'darkorange2','black')) +
  scale_shape_manual(name = "Iron",
                       breaks = c("Fe","noFe", "T0"),
                       labels = c("Fe", "noFe", "T0"),
                       values = c(19,1,5)) +
  xlab(expression ("")) +
  ylab (substitute(yaxislabel))+ 
  theme_bw()+
  theme(axis.text.y =element_text(face = "bold", size = 11, color = "black"),
        axis.text.x=element_text(face = "bold", size = 11, color = "black", vjust = 1.2, hjust = 1,  angle = 45),
        axis.title.y =element_text(size = 12,face = "bold", color = "black"), 
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 12, face = "bold"),
        legend.position = "none") +
    ylim (0,NA)
}

pennategroup <- protein_taxgroups("Pennate", "Pennate diatom fraction") 
centricgroup <- protein_taxgroups("Centric", "Centric diatom fraction") 
phaeogroup <- protein_taxgroups("Phaeocystis", "Phaeocystis fraction") 
```

Figure 2 Plot long : 2023-07-06
```{r}
nmds <- plot_grid (pigment_nmds, fc_nmds, peptide_nmds,   NULL,
                  ncol=1, 
                  rel_heights = c(1,1,1,0.12),
                  labels = c("A", "B","C", ""))

heatmaps <- plot_grid(fig2_pigmentmeans, fig2_proteintopmeans, fig2_proteinbottommeans, 
          ncol=1,
          labels = c("D","E", ""), 
          align = "v", 
          rel_heights = c(1.25,0.65,1.1))


rightpanel <- plot_grid (pennategroup, NULL, centricgroup, NULL,  phaeogroup,
           ncol =5, 
           rel_widths   = c(1, 0.05, 1, 0.05,1),
          labels = c("F", "", 
                     "G", "",
                     "H", ""),
                    
                     align = "vh")

rightpanel


a <- plot_grid (nmds, NULL, heatmaps, 
           ncol = 3, 
           rel_widths = c(0.7, 0.05, 1.2))


b <- plot_grid (a, rightpanel, 
           nrow = 2, 
           rel_heights = c(1,0.4))



#png("ps117_figure2_20230706.png", width=12, height=11, units="in", res=300) 
ggdraw(b) +               draw_plot (fig2_legendA, x = 0.387, y = 0.53, width = 0.0001, height = 0.001) +
  draw_plot (fig2_legendB, x = 0.38, y = 0.63, width = 0.0001, height = 0.001) 
                
dev.off()
```

fig2nmdsA - The vegan package displays 'species' and 'sites'. The sites are the treatments, the 'species' are the variables (pigments, FC-groups, proteomics species etc) 
First, make NMDS analyses from peptide, FC and pigment data. No need to re-run this! 
Peptide NMDS -  2023-03-026 
```{r}
protein_nmds <- proteindata_unnorm [c(1, 40:69)]

protein_nmds <- melt(protein_nmds, id.vars=c("peptide")) %>%
             rename(treatment = variable) %>%
             rename(unnorm_abundance = value )
sum(protein_nmds$unnorm_abundance) #this should equal 1.893E13 (total abundance of all peptides for all samples)

#normalize the peptide abundances for within each sample, irrespective of taxonomy
protein_nmds <- protein_nmds %>% #this creates a new column of normalized protein abundances grouped by treatment   
  group_by(treatment) %>% 
  mutate(norm_abundance = unnorm_abundance / sum(unnorm_abundance)) %>% 
  ungroup()

sum(protein_nmds$norm_abundance) #this should = 30 after normalization (number of replicates)

protein_nmds2 <- dcast(protein_nmds, treatment~peptide, fun.aggregate = sum,  value.var = "norm_abundance") #this creates MANY columns

protein_nmds2 <- protein_nmds2 %>% separate(treatment, into=c("bioassay", "timepoint" ,"temperature_treatment", "iron_treatment", "replicate"), sep = "_")

protein_nmds3 = protein_nmds2 [c(6:37477)] 

protein_nmds3 = as.matrix(protein_nmds3)

dimcheckMDS(protein_nmds3,
  distance = "bray",
  k = 4,
  trymax = 1000,
  autotransform = TRUE)

set.seed(123)

nmds_protein = metaMDS(protein_nmds3, try = 1000, trymax = 1000, k = 2, distance = "bray") # try = 500, trymax = 500
nmds_protein
scores(nmds_protein, display = "sites")


stressplot(nmds_protein, p.col = "black", l.col = "blue", title("Peptides")) #for the suplement

data.scores_protein <- as.data.frame(scores(nmds_protein, "sites")) #can select sites (treatments) or species (variable)

data.scores_protein$temperature_treatment <- protein_nmds2$temperature_treatment
data.scores_protein$iron_treatment <- protein_nmds2$iron_treatment
data.scores_protein$bioassay <- protein_nmds2$bioassay

data.scores_protein$bioassaytempfe <- paste(data.scores_protein$bioassay, data.scores_protein$temperature_treatment,data.scores_protein$iron_treatment, sep = "_")
```

fig2nmdsB - Flow cytometry NMDS - 2023-03-26
```{r}
FC <- nonproteindata [c(15,69,71,73,75,77,79,81,83,85,87,89,91,93,95,97,99,101,103,105, 107,109,111,113)] %>% na.omit() %>% separate(full_treatment_name, c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate")) 

FC_nmds <- FC [c(6:28)] 

FC_nmds = as.matrix(FC_nmds)

set.seed(123)
nmds_FC = metaMDS(FC_nmds, try = 10000, distance = "bray") # try = 500, trymax = 500
nmds_FC

data.scores_FC <- as.data.frame(scores(nmds_FC, "sites"))
stressplot(nmds_FC, p.col = "black", l.col = "blue", title("Flow Cytometry")) #for the suplement

data.scores_FC$temperature_treatment <- FC$temperature_treatment
data.scores_FC$iron_treatment <- FC$iron_treatment
data.scores_FC$bioassay <- FC$bioassay

data.scores_FC$bioassaytempfe <- paste(data.scores_FC$bioassay, data.scores_FC$temperature_treatment,data.scores_FC$iron_treatment, sep = "_")
```

fig2nmdsC - Pigment NMDS - 2023-03-26 
```{r}
pigment_nmds <- nonproteindata [c(15,49:54)] %>% na.omit() %>% separate(full_treatment_name, c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate")) 

pigment_nmds2 <- pigment_nmds [c(6:11)] 

pigment_nmds2 = as.matrix(pigment_nmds2)

set.seed(123)
nmds_pigment = metaMDS(pigment_nmds2, try = 1000, distance = "bray") # try = 500, trymax = 500
nmds_pigment

data.scores_pigment <- as.data.frame(scores(nmds_pigment, "sites"))
stressplot(nmds_pigment, p.col = "black", l.col = "blue", title("Pigments")) #for the suplement

data.scores_pigment$temperature_treatment <- pigment_nmds$temperature_treatment
data.scores_pigment$iron_treatment <- pigment_nmds$iron_treatment
data.scores_pigment$bioassay <- pigment_nmds$bioassay

data.scores_pigment$bioassaytempfe <- paste(data.scores_pigment$bioassay, data.scores_pigment$temperature_treatment,data.scores_pigment$iron_treatment, sep = "_")
```

fig2nmdsD - Then combine all three analyses outputs into one file for plotting - 20220421
```{r}
#data.scores_FC$analysis <- "FC"
#data.scores_pigment$analysis <- "pigment"
#data.scores_protein$analysis <- "peptide"

#ps117_nmds_20230421 <- rbind (data.scores_FC, data.scores_pigment, data.scores_protein)
#write.csv(ps117_nmds_20230421, "ps117_nmds_20230421.csv", row.names = F)
```

Figure 2 - Individual replicates - supplement ? 
```{r}
allgroupshm$treatmentrep <- paste (allgroupshm$temperature_treatment, allgroupshm$iron_treatment, allgroupshm$replicate, sep = "_")


treatmentrep_order <- c('T0_T0_A', 'T0_T0_B', 'T0_T0_C', 'T0_T0_D',
                     'LT_noFe_A', 'LT_noFe_B','LT_noFe_C',
                     'LT_Fe_A', 'LT_Fe_B', 'LT_Fe_C', 
                     'HT_noFe_A', 'HT_noFe_B','HT_noFe_C',
                     'HT_Fe_A', 'HT_Fe_B', 'HT_Fe_C')

  allgroupstop <- filter (allgroupshm, grepl("Bacillar|Other|Hapto|Bacteria", class))

ggplot(allgroupstop, aes( x = factor(treatmentrep, level = treatmentrep_order), y = class, fill= abundance))+ 
  geom_tile(aes(fill = norm_abundance),color = "black", lwd = 0.7 ) + 
  #coord_fixed()+ #this makes it all squares 
  scale_fill_gradient2(low = "deepskyblue2", 
                       high = "darkred",
                       mid = "white", 
                       midpoint = 0.25,
                       name = "normalized abundance") +
   ylab ("")+
   xlab ("")+
   theme(axis.text.x=element_text(face = "bold", size = 8, color = "black", vjust = 1,            hjust = 1, angle = 45),
        axis.text.y=element_text(face = "bold", size = 15, color = "black"),
        panel.border = element_blank(),
        panel.background = element_blank(), 
        axis.ticks.y  = element_blank()) +
   guides(fill = guide_colourbar(barwidth = 1, barheight = 15)) +
  facet_grid(.~ bioassay, scales = "free_x") #free_x removes replicates that aren't there


allgroupsbottom <- filter (allgroupshm, !grepl("Bacillar|Other|Hapto|Bacteria", class))
ggplot(allgroupsbottom, aes( x = factor(treatmentrep, level = treatmentrep_order), y = class, fill= abundance))+   geom_tile(aes(fill = norm_abundance),color = "black", lwd = 0.7 ) + 
  #coord_fixed()+ #this makes it all squares 
  scale_fill_gradient2(low = "deepskyblue2", 
                       high = "darkred",
                       mid = "white", 
                       midpoint = 0.003,
                       name = "normalized abundance") +
   ylab ("")+
   xlab ("")+
   theme(axis.text.x=element_text(face = "bold", size = 8, color = "black", vjust = 1,            hjust = 1, angle = 45),
        axis.text.y=element_text(face = "bold", size = 15, color = "black"),
        panel.border = element_blank(),
        panel.background = element_blank(), 
        axis.ticks.y  = element_blank()) +
   guides(fill = guide_colourbar(barwidth = 1, barheight = 15)) +
  facet_grid(.~ bioassay, scales = "free_x") #free_x removes replicates that aren't there

```

2024-05-24 -- nitrate, phosphate, silicate, POC, Si:No3 drawdown. 
```{r}
no3 <- x_vs_time(NO3_uM, NO[3]^"-"~~(mu*M)) + ylim(5, 30) + theme (axis.text.x=element_blank())
po4 <- x_vs_time(PO4_uM, PO[4]^"3-"~~(mu*M))  + ylim(0, 2) +theme(strip.text.x = element_blank())


pon <- T0_T8(nonproteindata, PON_uM, PON~~(mu*M)) + theme (axis.text.x=element_blank())
poc <- T0_T8(nonproteindata, POC_uM, POC~~(mu*M)) +theme(strip.text.x = element_blank())
chla <- T0_T8(nonproteindata, total_chla_ug_L, chla~~(mu*M)) +theme(strip.text.x = element_blank())

# pon <- x_vs_time(PON_uM, PON~~(mu*M))  + theme (axis.text.x=element_blank())
# poc <- x_vs_time(POC_uM, POC~~(mu*M)) +theme(strip.text.x = element_blank())

```

2024-05-25 -- figure 1 construction
```{r}
#png("../figures/ps117_figure1_20240808.png", width= 13, height= 6, units="in", res=600) 

a <- ggdraw (plot_grid(NULL, NULL, no3, NULL, pon, 
                  inset, NULL, po4, NULL, poc,
                 
        labels = c("(A)","", "(C)","", "(E)",
                   "(B)","", "(D)", "", "(F)"),
          label_size = 10,
          ncol = 5,
          rel_widths = c(0.5, 0.07, 1 ,0.07,1),
          align = "v")) +
  draw_plot (map, x = -0.015, y = 0.53, width = 0.25, height = 0.45)  +
  draw_plot (vertical_legend , x = 0.29, y= 0.1, width = 0.01, height = 0.01)  #y= 0.58

#dev.off ()

#if we want to add a map legend 
# ggdraw (x) + draw_plot (map, x = 0.15, y = 0.77, width = 0.22, height = 0.22)  +
#           draw_plot (vertical_legend , x = 0.068, y= 0.789, width = 0.05, height = 0.05) +  
#           draw_plot (maplegendvert  , x = 0.33, y= 0.79, width = 0.15, height = 0.15)   
#     # draw_plot (maplegendvert , x = 0.09, y= 0.84, width = 0.05, height = 0.05) +  
#           #draw_plot (vertical_legend , x = 0.36, y= 0.74, width = 0.15, height = 0.15)   
```

supplement ? 
```{r}
#things that go in the supplement 

#fvfm <- T0_T8_FenoFe (fvfm, F[v]/F[m]) + ylim(0.2, NA)
#chla <- T0_T8_FenoFe (total_chla_ug_L, Total~chl-a~~(mu*g~L^-1))
#si <- x_vs_time (Si_uM, Si~mu*M) + ylim(35,75) + theme(strip.text.x = element_blank())


si_no3_drawdown <- read.csv("LJ_ps117_si_consumption.csv", header = T)

sino3drawdown <- ggplot(filter (si_no3_drawdown, grepl ("8|0", daycount)), 
       aes(x=iron_treatment, 
       y=SiNO3ratio,
       color = temperature_treatment, 
       shape = iron_treatment))+ 
  stat_summary(fun = mean, geom = "point", size = 3.5, stroke = 1.3, alpha = 0.8)+
  stat_summary(fun.data =  mean_sd, geom = "errorbar", show.legend = F, linewidth = 1, width = 0.1)+
  scale_color_manual(name = "Temperature",
                     breaks = c("LT","HT"),
                     labels = c("LT", "HT"),
                    values = c('#0072B2', "#D55E00", "black")) +
  
  
  scale_shape_manual(name = "Iron",
                       breaks = c("Fe","noFe", "T0"), #this removes the T0 from the legend
                       labels = c("+Fe", "noFe", "T0"),
                       values = c(19,1,5)) +
  scale_x_discrete(limits=c("T0", "noFe","Fe"))+
  ylab (expression (Si:NO[3]^"-"~drawdown~~(M:M))) + 
 # ylab (expression (Delta*Si:Delta*NO[3]~(M:M))) + 
  xlab (expression ("")) +
  theme_bw()+
  theme(axis.text=element_text(face = "bold", size = 10, color = "black"),
        axis.title.x=element_text(size = 12,face = "bold", color = "black"), 
        axis.title.y=element_text(size=12, color = "black"), 
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 10, face = "bold"),
        legend.position = "none", 
        panel.grid = element_blank()) +
  ylim(0,2.5)+
  facet_wrap(~bioassay)
sino3drawdown
```












#2024-07-26 -- Figure 3 
            -- NP drawdown and PON:POP ratio, ribosomes etc. 
```{r}
#the drawdown is calculated between day 4 and day 8 based on the NO3 and PO4 plots in Fig 1, showing when the incubations are acclimated. 
npdrawdown <- read.csv("LJ_ps117_NO3_PO4_consumption.csv")
npdrawdown$bioassayFe <- paste (npdrawdown$bioassay, npdrawdown$iron_treatment, sep = "_")

npdrawdown$bioassay <- ifelse (npdrawdown$bioassay == "BA1" , "Offshore", 
                      ifelse (npdrawdown$bioassay == "BA2", "Nearshore", 
                          "somethingwentwrong"))



NP_drawdown <- T0_T8(filter (npdrawdown, grepl("8", timepoint)), (total_Ndrawdown/total_Pdrawdown), NO[3]^"-":PO[4]^"3-"~drawdown~~(M:M)) + 
  geom_hline (aes(yintercept = 16), col = "black", lty = 2, lwd=0.6) + 
  scale_x_discrete(limits = c("LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"), 
                     labels = c("Ctrl", "+Fe", "HT", "HT +Fe"))

particulatenp <- T0_T8(filter (nonproteindata, grepl("8|0", timepoint)), ((1000*PON_uM)/POP_nM), Particulate~N:P~~(M:M)) +   
  geom_hline (aes(yintercept = 16), col = "black", lty = 2, lwd=0.6) + 
  theme(strip.text.x = element_blank()) 

# #old script
# ggplot((filter(npdrawdown, grepl("8", timepoint))) , aes(
#     x = tempfe_treatment,
#     y = (total_Ndrawdown/total_Pdrawdown),
#     shape = iron_treatment, 
#     color = temperature_treatment))+
#    stat_summary(fun = mean, geom = "point", show.legend = T, size = 4, stroke = 1.5, alpha = 0.65)+
#    stat_summary(fun.data =  mean_se, geom = "errorbar", size = 1, width = 0.1)+
#    ylab (expression (NO[3]^"-":PO[4]^"3-"~drawdown~~(M:M))) + 
#   xlab("") +
#    facet_wrap (~bioassay) +
#    
#    scale_color_manual(name="Temperature",
#                      breaks = c("LT","HT"),
#                      labels = c("LT", "HT"),
#                      values = c('black', 'darkorange2','black')) +
#   scale_shape_manual(name = "Iron",
#                       breaks = c("Fe","noFe", "T0"), 
#                       labels = c("+Fe", "-Fe", "T0"),
#                       values = c(19,1,5)) 
#   scale_x_discrete(limits=c("HT_noFe","HT_Fe", "LT_noFe", "LT_Fe"))

```

2024-07-26 -- NP from literature (prelim data)
```{r}
alldata <- read.csv("LJ_prelim_data_20221206.csv")
#alldata$limitation  <- paste(alldata$NO3_lim, alldata$PO4_lim)
#alldata2 <- filter(alldata, !grepl("y",limitation))

stoik_temp <- function (ratio, hline, yaxislabel) {
    ggplot(filter(alldata, grepl("polar",provenance)), aes(
    x = temperature,
    y = {{ratio}}))+
  #  stat_summary(fun = mean, geom = "point", show.legend = T, size = 4, stroke = 1.5, alpha = 0.65)+
  # stat_summary(fun.data =  mean_se, geom = "errorbar", size = 1, width = 0.1)+
  geom_point(size  = 1.5, alpha = 0.6, stroke = 1.1, shape = 4) +
    theme_bw() +
    ylab (substitute(yaxislabel)) +
    xlab (expression("Temperature (°C)")) +
    theme(axis.text=element_text(face = "bold", size = 9, color = "black"),
        axis.title.x=element_text(size = 12,face = "bold", color = "black"), 
        axis.title.y=element_text(size=12, color = "black"), 
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 10, face = "bold"),
        legend.position = "none", 
        panel.grid = element_blank())+
 geom_hline (aes(yintercept = hline), col = "black", lty = 2, lwd=0.6 )}

#stoik_temp(NP_MM, 16, N:P~~~(M:M)) 
 metaanaNP <- stoik_temp(NP_MM, 16, Particulate~N:P~~~(M:M)) + facet_wrap(.~ type) + theme ( legend.position = "none") 
 
```

2024-07-26 -- sample-normalized, taxon-agnostic ribosome abundance
```{r}
#First make sure that all the ribosomal peptides are annotated as such, and that the groups are properly annotated 
proteindata_unnorm$compartment <- ifelse (proteindata_unnorm$uniq_cluster  ==  "unassigned_cluster", "ribosome", 
                                  ifelse (proteindata_unnorm$uniq_grpnorm_compartment == "ribosomal", "ribosome",
                                  ifelse (grepl("ribosom|Ribosom", proteindata_unnorm$uniq_cluster_annotation) ,  "ribosome",
                                   
                                  ifelse (proteindata_unnorm$uniq_KO_pathway  ==  "Photosynthesis", "chloro", 
                                  ifelse (grepl("^clust_1820$|^clust_534$", proteindata_unnorm$uniq_cluster) ,  "chloro",
                                  ifelse (grepl("^clust_891$", proteindata_unnorm$uniq_cluster) ,  "nitrogen",

                                          
                                  ifelse (grepl("^clust_129$|^clust_2490$|^clust_1820$|^clust_534$", proteindata_unnorm$uniq_cluster) ,  "rubisco",
                                  ifelse (grepl("Hsp|hsp", proteindata_unnorm$uniq_cluster_annotation) ,  "HSP",
                                  ifelse (grepl("ATP", proteindata_unnorm$uniq_cluster_annotation),  "ATP",
                                              proteindata_unnorm$uniq_grpnorm_compartment))))))))) #

compartment <- proteindata_unnorm [c(70, 37, 40:69  )]     

compartment <- melt(compartment, id.vars=c("compartment", "F")) %>%
            rename(treatment = variable) %>%
            rename(unnnorm_abundance = value) 
            #separate(treatment, into=c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"), sep = "_")


#notmalize to total protein
compartment_aggregated_total <- aggregate(compartment$unnnorm_abundance, 
                                  by=list(Category=compartment$compartment, compartment$treatment), FUN=sum) %>%
                        rename(unnorm_abundance = x) %>%
                        rename(compartment = Category) %>%
                        rename(treatment = Group.2) 
                       
compartment_aggregated_total <- compartment_aggregated_total %>% 
                          group_by(treatment) %>% 
                          mutate(totalproteinnorm_abundance = unnorm_abundance / sum(unnorm_abundance)) %>% 
                          ungroup()

compartment_aggregated_total <- compartment_aggregated_total %>% 
  separate(treatment, into=c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"), sep = "_") %>%
  mutate(tempfe_treatment = paste(temperature_treatment, iron_treatment, sep = "_")) %>%
  mutate(tempfe_treatment = str_replace(tempfe_treatment, "T0_T0", "T0"))



#quick plot of the major compartments under different ocnditions
ggplot(filter (compartment_aggregated_total, grepl ("^chloro$|^ribosome$|^rubisco$", compartment)),
       aes(x= tempfe_treatment,
       y= totalproteinnorm_abundance*100,
       fill = compartment))+
  facet_wrap(~bioassay) +
  stat_summary(fun = mean, geom = "bar", width=0.75,position=position_dodge(width=0.75), color = "black", alpha = 1)+
  stat_summary(fun.data =  mean_se, geom = "errorbar", show.legend = F, linewidth = 1, position=position_dodge(width=0.75), width = 0.2, color = "black") +
  scale_fill_manual(name="",
                     breaks = c("chloro", "ribosome", "rubisco"),
                    labels = c("Photosynthesis", "Ribosome", "RuBisCO"),
                     values = c('#009E73', '#56B4E9', '#E69F00')) +
  
  
  scale_x_discrete(limits=c("T0", "LT_noFe","LT_Fe", "HT_noFe", "HT_Fe"),
                   labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe")) +
  xlab (expression(""))+
  ylab("Percent of all matched proteins") +
  theme_bw()+
  theme(axis.text=element_text(face = "bold", size = 12, color = "black"),
        axis.title.x=element_text(size = 14,face = "bold", color = "black"), 
        axis.title.y=element_text(size=16, color = "black", face ="bold"), 
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 12, face = "bold"))

ribosome <- T0_T8 (filter (compartment_aggregated_total, grepl ("^ribosome$", compartment)), totalproteinnorm_abundance, "Ribosomal fraction of total protein") +
  theme(axis.text.x = element_blank())
  

#old plot with orange dots etc. 
# samplenormcompartment <- function (compartments, yaxislabel) { 
# ggplot(filter (compartment_aggregated_total, grepl (compartments, compartment)),
#        aes(x= iron_treatment,
#        y= totalproteinnorm_abundance,
#        color = temperature_treatment, 
#        shape = iron_treatment)) + 
#   facet_wrap(~bioassay)+ 
#   stat_summary(fun = mean, geom = "point", size = 4, stroke = 1.5, alpha = 0.5)+
#   stat_summary(fun.data =  mean_se, geom = "errorbar", show.legend = F, linewidth = 1, width = 0.1)+
#   scale_color_manual(name="Temperature",
#                      breaks = c("LT","HT"),
#                      labels = c("LT", "HT"),
#                      values = c('black', 'darkorange2','black')) +
#   scale_shape_manual(name = "Iron",
#                        breaks = c("Fe","noFe", "T0"), 
#                        labels = c("Fe", "noFe", "T0"),
#                        values = c(19,1,5)) +
#   scale_x_discrete(limits=c("T0","noFe","Fe"), 
#                    labels = c("T0", "-Fe", "+Fe"))+
#                  # labels = expression(bold("T0"), bold("T8")~"(-Fe)", bold("T8")~"(+Fe)"))+
#   xlab(expression ("")) +
#   ylab (substitute(yaxislabel))+ 
# 
#   theme_bw()+
#   theme(axis.text=element_text(face = "bold", size = 12, color = "black"),
#         axis.title =element_text(size = 14, color = "black"), 
#         strip.background =element_rect(fill = "white"),
#         strip.text.x = element_text(size = 12, face = "bold"),
#         legend.position = "none") 
# }

#ribosome <- samplenormcompartment ("^ribosome$", "Ribosomal fraction of total protein") +   ylim (0, 0.22)
#ribosome
```

2024-07-26 -- figure 3 construction 
```{r}
#png("../figures/ps117_figure3_20240726.png", width=11, height=7.5, units="in", res=600) 


plot_grid (ribosome, NULL, NP_drawdown, 
           particulatenp, NULL, metaanaNP,
           rel_widths = c(1,0.04,1), 
           ncol = 3, 
           align = "vh", 
            axis ="b", 
           labels = c ("(A)", "", "(C)", 
                       "(B)", "", "(D)"))
         
#dev.off()
```


#2024-07-27 -- Figure 4
            -- pMetal:pN and DE protein data 
2024-07-29  -- metal plots - use only metal:pN for main paper, maybe metal:pP  and metal:pC in supplement 
```{r}
pFepN <- T0_T8 (nonproteindata, (pFe57_nM)/(PON_uM),p ^{57}*"Fe":pN~~(mM:M)) + theme(axis.text.x = element_blank()) + 
  theme(axis.title.y=element_text(colour="#c34007"))


pZnpN <- T0_T8(nonproteindata, pZn_nM/(PON_uM), pZn:pN~~(mM:M)) + theme(strip.text.x = element_blank(), axis.text.x = element_blank()) + theme(axis.title.y=element_text(colour="#0072B2"))

pCupN <- T0_T8(nonproteindata, pCu_nM/(PON_uM), pCu:pN~~(mM:M)) + theme(strip.text.x = element_blank(), axis.text.x = element_blank()) + theme(axis.title.y=element_text(colour="orange"))

pCdpN <- T0_T8(nonproteindata, pCd_nM/(PON_uM), pCd:pN~~(mM:M)) + theme(strip.text.x = element_blank(), axis.text.x = element_blank()) + theme(axis.title.y=element_text(colour="grey50"))

pMnpN <- T0_T8(nonproteindata, pMn_nM/(PON_uM), pMn:pN~~(mM:M)) + theme(strip.text.x = element_blank(), axis.text.x = element_blank()) + theme(axis.title.y=element_text(colour="pink"))

pNipN <- T0_T8(nonproteindata, pNi_nM/(PON_uM), pNi:pN~~(mM:M)) + theme(strip.text.x = element_blank()) + theme(axis.title.y=element_text(colour="#009E73"))
#totalFe57 <- metalplot (PFe57_nM+DFe57_nM, Total ^{57}*Fe~(nM))
```

2024-08-2 -- massage DE and fold change analysis outputs 
```{r}
#Read in the DE analyses for BA1 and B2, and combine
# Offshore <- read.csv("ps117_BA1_clusterDE_20230328.csv", header = T)
# colnames(Offshore)[3:14] <- paste("Offshore", colnames(Offshore[,c(3:14)]), sep = "_")
# 
# Offshorenotax <- read.csv ("ps117_Offshorenotaxon_clusterDE_20240802.csv")
# colnames(Offshorenotax)[2:13] <- paste("Offshore", colnames(Offshorenotax[,c(2:13)]), sep = "_")
# Offshorenotax$taxon_F <- "All taxa"
# 
# OffshoreDE <- rbind (Offshore, Offshorenotax, all = TRUE)
# 
# 
# Nearshore <- read.csv("ps117_BA2_clusterDE_20230328.csv", header = T)
# colnames(Nearshore)[3:14] <- paste("Nearshore", colnames(Nearshore[,c(3:14)]), sep = "_")
# 
# Nearshorenotax <- read.csv ("ps117_Nearshorenotaxon_clusterDE_20240804.csv")
# colnames(Nearshorenotax)[2:13] <- paste("Nearshore", colnames(Nearshorenotax[,c(2:13)]), sep = "_")
# Nearshorenotax$taxon_F <- "All taxa"
# NearshoreDE <- rbind (Nearshore, Nearshorenotax, all = TRUE)


##



#do this for the diatoms grouped as one
Offshore <- read.csv("ps117_offshore_diatomphaeo_clusterDE_20241120.csv", header = T)
colnames(Offshore)[3:14] <- paste("Offshore", colnames(Offshore[,c(3:14)]), sep = "_")

Offshorenotax <- read.csv ("ps117_Offshorenotaxon_clusterDE_20240802.csv")
colnames(Offshorenotax)[2:13] <- paste("Offshore", colnames(Offshorenotax[,c(2:13)]), sep = "_")
Offshorenotax$taxon_F <- "All taxa"

OffshoreDE <- rbind (Offshore, Offshorenotax, all = TRUE)


Nearshore <- read.csv("ps117_nearshore_diatomphaeo_clusterDE_20241120.csv", header = T)
colnames(Nearshore)[3:14] <- paste("Nearshore", colnames(Nearshore[,c(3:14)]), sep = "_")

Nearshorenotax <- read.csv ("ps117_Nearshorenotaxon_clusterDE_20240804.csv")
colnames(Nearshorenotax)[2:13] <- paste("Nearshore", colnames(Nearshorenotax[,c(2:13)]), sep = "_")
Nearshorenotax$taxon_F <- "All taxa"
NearshoreDE <- rbind (Nearshore, Nearshorenotax, all = TRUE)


##
allDE <- merge (OffshoreDE, NearshoreDE, by = c("cluster", "taxon_F"), all = TRUE)



#pick the three groups I'm interested in, and change the names of clusters to include the names I'm interested in.                     
#BA_DE <- filter (allDE, grepl ("centric|pennate|Phaeocyst|All taxa", taxon_F))
BA_DE <- filter (allDE, grepl ("diatom|Phaeocyst|All taxa", taxon_F))


#to annotate clusters 
#annotations <- fread("datannotations/annotation_allTFG.grpnorm_mmetsp_fc_pn_reclassified.edgeR.csv") 

#annotations2 <- annotations [,c(2, 4:9, 12:16, 27:31, 44)]


#annotations3 <- annotations2 %>% distinct(across(everything()))

#give taxa sensible names
BA_DE$taxon_F <- ifelse (BA_DE$taxon_F  ==  "centric", "Centric diatoms", 
                  ifelse (BA_DE$taxon_F  ==  "pennate", "Pennate diatoms", 
                  ifelse (BA_DE$taxon_F  ==  "diatom", "Diatoms", 
                  ifelse(BA_DE$taxon_F  ==  "Phaeocystaceae", "Phaeocystis",
                  BA_DE$taxon_F))))

#do some organization 
BA_DE$Offshore_ironDE <-  ifelse (BA_DE$Offshore_Fe_FDR < 0.05, "*", "")
BA_DE$Nearshore_ironDE <-  ifelse (BA_DE$Nearshore_Fe_FDR < 0.05, "*", "")

#kog_KOG0985_Vesicle coat protein clathrin, heavy chain

#this will be the end all be all cluster annotations. 
#This is based on consensus annotations as well as manual blasting, checking and double checking

#https://www.genome.jp/entry/pathway+ko00195

cluster_combinedannotations <- c(
  "clust_332" = "LHC",
      "clust_3986" = "LHC", # double check
      "clust_5925" = "LHC",
      "clust_1044" = "LHC",
      "clust_195" = "LHC",
      "clust_1236" = "LHC",
      "clust_3440" = "LHC",
      "clust_3658" = "LHC",
      "clust_3703" = "LHC",
      "clust_3551" = "LHC",
  
   "clust_42" = "PSII",
     "clust_1978" = "PSII",
      "clust_498" = "PSII",
      "clust_211" = "PSII",
      "clust_1784" = "PSII",
      "clust_4256" = "PSII",
      "clust_5714" = "PSII",
      "clust_1550" = "PSII",

  "clust_1906" = "OEC",
      "clust_2763" = "OEC",
  
  "clust_613" = "Cyt-b6f",
      "clust_1497" = "Cyt-b6f",
      "clust_1668" = "Cyt-b6f",
      "clust_1822" = "Cyt-b6f",
  
    "clust_55" = "PSI",
      "clust_2815" = "PSI",
      "clust_1564" = "PSI",
      "clust_2439" = "PSI",
      "clust_2460" = "PSI",
      "clust_3396" = "PSI",
  
  "clust_1820" = "Pcyn",
  "clust_534" = "Fld",
  "clust_231" = "Rho",
  
    "clust_2460" = "Fd",
    "clust_808" = "Fd",

    "clust_216" = "Pyrophosphatase",  "clust_522" = "Pyrophosphatase",  "clust_868" = "Pyrophosphatase",  "clust_522" = "Pyrophosphatase",

  "clust_2490" = "Rubisco",
  "clust_655" = "CA_clust_655",
  "clust_6573" = "Cd-CA",
  "clust_1525" = "SLC-4",
  "clust_1886" = "SLC",
  "clust_411" = "NR",
  "clust_3563" = "NiR",

  "clust_417" = "AMT",
  "clust_891" = "Glutamate Synthase",
  "clust_1372" = "UT",
  "clust_320" = "Urea carboxylase",
  "clust_6019" = "Ornithine cyclodeaminase",
  "clust_955" = "Mg chltase",
      "clust_416" = "Mg chltase",
      "clust_2922" = "Mg chltase",

  "clust_5562" = "Cu/Zn SOD",
  "clust_35741" = "Cu/Zn SOD",

  "clust_2501" = "Fe/Mn SOD",
  "clust_1289" = "DoxX", 
  "clust_520" = "Alkyl hydroperoxide redct-peroxiredoxin", 
      "clust_9714" = "Peroxiredoxin",
  
  "clust_2244" = "PEPC",
  
  "clust_2088" = "FAS", 
      "clust_6071" = "FAS", 
      "clust_11151" = "FAS",
      "clust_149204, clust_3492, clust_7032, clust_87" = "FAS", 

  "clust_155" = "RNR", #this is interesting
  
  "clust_2934" = "Multicopper Oxidase",
  "clust_753" = "Fructose bisphosphatase",
  "clust_997" = "APRT",
  "clust_10365" = "PET8",

  "clust_1322" = "Porphobilinogen deaminase",
  
  "clust_1814" = "ISIP",
      "clust_343" = "ISIP",
      "clust_1154" = "ISIP")

cluster_combinedannotationslongname <- c(
  "clust_332" = "Light harvesting complex",
      "clust_3986" = "Light harvesting complex", # double check
      "clust_5925" = "Light harvesting complex",
      "clust_1044" = "Light harvesting complex",
      "clust_195" = "Light harvesting complex",
      "clust_1236" = "Light harvesting complex",
      "clust_3440" = "Light harvesting complex",
      "clust_3658" = "Light harvesting complex",
      "clust_3703" = "Light harvesting complex",
      "clust_3551" = "Light harvesting complex",
  
   "clust_42" = "Photosystem II",
     "clust_1978" = "Photosystem II",
      "clust_498" = "Photosystem II",
      "clust_211" = "Photosystem II",
      "clust_1784" = "Photosystem II",
      "clust_4256" = "Photosystem II",
      "clust_5714" = "Photosystem II",
      "clust_1550" = "Photosystem II",

  "clust_1906" = "Oxyen evolving complex",
      "clust_2763" = "Oxyen evolving complex",
  
  "clust_613" = "Cytochrome b6f",
      "clust_1497" = "Cytochrome b6f",
      "clust_1668" = "Cytochrome b6f",
      "clust_1822" = "Cytochrome b6f",
  
    "clust_55" = "Photosystem I",
      "clust_2815" = "Photosystem I",
      "clust_1564" = "Photosystem I",
      "clust_2439" = "Photosystem I",
      "clust_2460" = "Photosystem I",
      "clust_3396" = "Photosystem I",
  
  "clust_1820" = "Plastocyanin",
  "clust_534" = "Flavodoxin",
  "clust_231" = "Rhodopsin",
  
    "clust_2460" = "Ferredoxin",
    "clust_808" = "Ferredoxin",

    "clust_216" = "Pyrophosphatase",  "clust_522" = "Pyrophosphatase",  "clust_868" = "Pyrophosphatase",  "clust_522" = "Pyrophosphatase",

  "clust_2490" = "Rubisco",
  "clust_655" = "CA_clust_655",
  "clust_6573" = "Cd carbonic anhydrase",
  "clust_1525" = "SLC-4",
  "clust_1886" = "SLC",
  "clust_411" = "Nitrate reductase",
  "clust_3563" = "Nitrite/sulphite reductase",

  "clust_417" = "Ammonium transporter",
  "clust_891" = "Glutamate Synthase",
  "clust_1372" = "Urea transporter",
  "clust_320" = "Urea carboxylase",
  "clust_6019" = "Ornithine cyclodeaminase",
  "clust_955" = "Mg chelatase",
      "clust_416" = "Mg chelatase",
      "clust_2922" = "Mg chelatase",

  "clust_5562" = "Cu/Zn superoxide dismutase",
  "clust_35741" = "Cu/Zn superoxide dismutase",
  
  "clust_2501" = "Fe/Mn superoxide dismutase",
  "clust_1289" = "DoxX", 
  "clust_520" = "Alkyl hydroperoxide redct-peroxiredoxin", 
      "clust_9714" = "Peroxiredoxin",
  
  "clust_2244" = "PEPC",
  
  "clust_2088" = "Fasciclin", 
      "clust_6071" = "Fasciclin", 
      "clust_11151" = "Fasciclin",
      "clust_149204, clust_3492, clust_7032, clust_87" = "Fasciclin", 

  "clust_155" = "Ribonucleotide reductase", #this is interesting
  
  "clust_2934" = "Multicopper Oxidase",
  "clust_753" = "Fructose bisphosphatase",
  "clust_997" = "APRT",
  "clust_10365" = "PET8",

  "clust_1322" = "Porphobilinogen deaminase",
  
  "clust_1814" = "Iron stress induced proteins",
      "clust_343" = "Iron stress induced proteins",
      "clust_1154" = "Iron stress induced proteins"
)

#'Porphobilinogen deaminase', 
              
#https://www.frontiersin.org/journals/cellular-and-infection-microbiology/articles/10.3389/fcimb.2014.00052/full#B40
#https://pubmed.ncbi.nlm.nih.gov/21338418/
#https://pubmed.ncbi.nlm.nih.gov/21371140/
#https://onlinelibrary.wiley.com/doi/10.1111/j.1365-2958.2011.07593.x

#clust_155
BA_DE$cluster_annotation_combined <- cluster_combinedannotations[BA_DE$cluster]

BA_DE$cluster_annotation_combined_long <- cluster_combinedannotationslongname[BA_DE$cluster]

#metals these proteins are known to contain 
cluster_metal <- c(
  "PSI" = "Fe",
  "PSII" = "Fe", 
  "Fd" = "Fe",
  "OEC"= "Mn",
  "Cyt-b6f" = "Fe", 
  "Pcyn" = "Cu",
  "Cu/Zn SOD" = "Cu", 
  "Fe/Mn SOD" = "Fe", 
  "NiR" =  "Fe", 
  "Cd-CA" = "Cd")

cluster_metal2 <- c(
  "Cu/Zn SOD" = "Zn", 
  "Fe/Mn SOD" = "Mn")

BA_DE$cluster_metal <- cluster_metal[BA_DE$cluster_annotation_combined]
BA_DE$cluster_metal2 <- cluster_metal2[BA_DE$cluster_annotation_combined]
```

2024-08-08 -- sliding plots for groups of clusters (e.g. all the proteins in PSI as one row)
```{r}
#https://pmc.ncbi.nlm.nih.gov/articles/PMC6332867/ ubiquitin SOD
cluster_combined_order <- c('LHC', 
                            'PSII', 
                            'OEC', 
                            'Cyt-b6f', 
                            'Pcyn',  
                            'PSI', 
                            'Fld', 
                            'Fd', 
                            'Rho',
                            'Mg chltase',
                            'ISIP',
                            'RNR',
                            'FAS',
                            'Cu/Zn SOD',
                            'Fe/Mn SOD',
                            'DoxX',
                            'Cd-CA', 
                            'NR',
                            'NiR',
                            'UT', 
                            'AMT')

cluster_combined_order_long <- c('Light harvesting complex', 
                            'Photosystem II', 
                            'Oxyen evolving complex', 
                            'Cytochrome b6f', 
                            'Plastocyanin',  
                            'Photosystem I', 
                            'Flavodoxin', 
                            'Ferredoxin', 
                            'Rhodopsin',
                            'Mg chelatase',
                            'Iron stress induced proteins',
                            'Ribonucleotide reductase',
                            'Fasciclin',
                            'Cu/Zn superoxide dismutase',
                            'Fe/Mn superoxide dismutase',
                            'DoxX',
                            'Cd carbonic anhydrase', 
                            'Nitrate reductase',
                            'Nitrite/sulphite reductase',
                            'Urea transporter', 
                            'Ammonium transporter')

clustersofinteres_combined <- paste(cluster_combined_order, collapse = "|")

clustersofinteres_combined_long <- paste(cluster_combined_order_long, collapse = "|")


BA_DE_metal <- BA_DE [c("cluster_annotation_combined", "cluster_metal", "cluster_metal2")]

BA_DE_metal <- melt(BA_DE_metal,  id.vars =  "cluster_annotation_combined")

pannel1 <- ggplot(filter(BA_DE_metal, grepl(clustersofinteres_combined, cluster_annotation_combined)),
                  aes(x = value, 
                      y = factor(cluster_annotation_combined, levels = rev(cluster_combined_order)),
                      fill = value)) + 
    scale_x_discrete(limits = c("Fe", "Zn", "Cu", "Cd", "Mn", "Ni")) +
    geom_tile(width = 0.7, height = 0.8, lwd = 0.7, color = "black") + 
    scale_fill_manual(values = c("Fe" = "#c34007", 
                                 "Zn" = "#0072B2", 
                                 "Cu" = "orange", 
                                 "Cd" = "grey50", 
                                 "Mn" = "pink", 
                                 "Ni" = "blue")) +
    xlab(expression(paste(""))) +
    coord_equal() +
    theme_bw() +
    theme(axis.text.x = element_text(size = 9, face = "bold", color = "black"),
          axis.title.y = element_blank(),
          axis.text.y = element_text(size = 12, face = "bold", color = "black"), 
          legend.position = "none") + 
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))

pannel1
       
####
####
#



BA_DE <- BA_DE %>%
  mutate(taxon_combined = ifelse(taxon_F %in% c("Centric diatoms", "Pennate diatoms"), "Diatoms", taxon_F))

BA_DE <- BA_DE %>%
  mutate(taxon_combined = factor(taxon_combined, levels = c("Diatoms", "Phaeocystis", "All taxa")))

DE_slide <- function(effect, significance, title) {
ggplot(filter(BA_DE, grepl(clustersofinteres_combined, cluster_annotation_combined)), 
      aes(x = {{effect}}, 
             y = factor(cluster_annotation_combined, levels = rev(cluster_combined_order)), 
             #y = reorder (cluster_combinedannotations, {{effect}}), #this is for reordering from positive to negative DE
             color = taxon_F, 
             shape = interaction(taxon_F, {{significance}}))) + 
    geom_point(alpha = 0.5, size = 4, stroke = 1) +
    
    # scale_color_manual(name = "",
    #                    values = c("Centric diatoms" = 'black',
    #                               "Pennate diatoms" = "chocolate4", ##999933
    #                               "All taxa" = 'black',
    # 
    #                               "Phaeocystis" = 'black'),
    #                    breaks = c("Centric diatoms", "Pennate diatoms"),
    #                   labels = c("Centric", "Pennate")) +
    # scale_shape_manual(name = "",
    #                    values = c("Centric diatoms.*" = 19, "Centric diatoms." = 1,
    #                               "Pennate diatoms.*" = 19, "Pennate diatoms." = 1,
    #                               "Phaeocystis.*" = 19, "Phaeocystis." = 1,
    #                               "All taxa.*" = 19, "All taxa." = 1), #15,0
    #                    guide = "none") +
    # 
  
  
   scale_color_manual(name = "",
                        values = c("All taxa" = 'black',
                                   "Diatoms" = "black",
                                 "Phaeocystis" = 'black'), guide = "none")+

    scale_shape_manual(name = "",
                       values = c("Diatoms.*" = 19, "Diatoms." = 1,
                                  "Phaeocystis.*" = 19, "Phaeocystis." = 1,
                                  "All taxa.*" = 19, "All taxa." = 1), #15,0
                       guide = "none") +

  
  ylab ("") + 
    xlab(expression(paste("Log"[2]~fold~change))) +
    ggtitle(title) +
    theme_bw() +
    theme(plot.title = element_text(hjust = 0, size = 10, color = "black", face = "bold"),
          axis.text.x = element_text(size = 10, face = "bold", color = "black"),
          axis.title.x = element_text(size = 12, face = "bold", color = "black"), 
          axis.text.y = element_text(size = 12, face = "bold", color = "black"), 
          legend.text = element_text(size = 9, color = "black", face = "bold"), 
          strip.background =element_rect(fill = "white"),
          strip.text.x = element_text(size = 10, face = "bold")) +
    geom_vline(xintercept = c(0), col = "black", lty = 2, linewidth = 0.7) + 
    facet_wrap(~taxon_combined, ncol = 3)
  }


FeDE_Nearshore <- DE_slide (Nearshore_Fe_logFC, Nearshore_ironDE, "Nearshore Fe-addition effect" ) + 
theme(legend.position = "none")+
  guides(colour = guide_legend(override.aes = list(size=3))) +
   theme(legend.margin=margin(t=-0.7,  r=0, b=0, l=0, unit="cm"), axis.text.y = element_blank()) +
   xlim(-7,8.2) +   
  theme(plot.margin = margin(0, 0, 0, 0))

  FeDE_Nearshore

  
  
DE_slide <- function(effect, significance, title) {
ggplot(filter(BA_DE, grepl(clustersofinteres_combined_long, cluster_annotation_combined_long)), 
      aes(x = {{effect}}, 
             y = factor(cluster_annotation_combined_long, levels = rev(cluster_combined_order_long)), 
             #y = reorder (cluster_combinedannotations, {{effect}}), #this is for reordering from positive to negative DE
             color = taxon_F, 
             shape = interaction(taxon_F, {{significance}}))) + 
    geom_point(alpha = 0.5, size = 4, stroke = 1) +
    
    
   scale_color_manual(name = "",
                        values = c("All taxa" = 'black',
                                   "Diatoms" = "black",
                                 "Phaeocystis" = 'black'), guide = "none")+

    scale_shape_manual(name = "",
                       values = c("Diatoms.*" = 19, "Diatoms." = 1,
                                  "Phaeocystis.*" = 19, "Phaeocystis." = 1,
                                  "All taxa.*" = 19, "All taxa." = 1), #15,0
                       guide = "none") +
    
    ylab ("") + 
    xlab(expression(paste("Log"[2]~fold~change))) +
    ggtitle(title) +
    theme_bw() +
    theme(plot.title = element_text(hjust = 0, size = 10, color = "black", face = "bold"),
          axis.text.x = element_text(size = 10, face = "bold", color = "black"),
          axis.title.x = element_text(size = 12, face = "bold", color = "black"), 
          axis.text.y = element_text(size = 12, face = "bold", color = "black"), 
          legend.text = element_text(size = 9, color = "black", face = "bold"), 
          strip.background =element_rect(fill = "white"),
          strip.text.x = element_text(size = 10, face = "bold")) +
    geom_vline(xintercept = c(0), col = "black", lty = 2, linewidth = 0.7) + 
    facet_wrap(~taxon_combined, ncol = 3)
  }
  
  
  
  

FeDE_Offshore <- 
  DE_slide (Offshore_Fe_logFC, Offshore_ironDE, "Offshore Fe-addition effect" ) + 
   theme (legend.position = "none") +
    guides(colour = guide_legend(override.aes = list(size=3))) + 
 theme(plot.margin = margin(0, 0, 0, 0))+
 xlim(-7,8.2)

 FeDE_Offshore

```

2024-08-08 -- figure 4 construction
```{r}
plot_grid(FeDE_Nearshore,
          FeDE_Offshore,
          ncol = 1, 
          align = "hv")

a <- plot_grid (FeDE_Nearshore, 
                FeDE_Offshore,
                labels = c("(G)", "(H)"),
                ncol = 1, 
                align = "hv") 

b <- plot_grid (pFepN, 
                pZnpN,
                pCupN,
                pCdpN, 
                pMnpN, 
                pNipN, 
                ncol = 1,
                labels = c ("(A)", "(B)" ,"(C)", 
                       "(D)", "(E)", "(F)"), 
                align = "v")
png("../figures/ps117_figure4_20241129.png", width=19, height=14, units="in", res=600) 
ggdraw (plot_grid (b, NULL,  a,
        rel_widths = c(0.6,0.06, 1),
        ncol = 3, 
        align = "vh") + 
        draw_plot (pannel1, x = 0.325, y = 0.512, width = 0.25, height = 0.455)) 
 # draw_grob(grid::rectGrob(gp = grid::gpar(col = "#c34007", fill = "#c34007")), x = 0.04, y = 0.955, width = 0.01, height = 0.015) + 
   # draw_grob(grid::rectGrob(gp = grid::gpar(col = "#0072B2", fill = "#0072B2")), x = 0.04, y = 0.806, width = 0.01, height = 0.015) + 
    #  draw_grob(grid::rectGrob(gp = grid::gpar(col = "orange", fill = "orange")), x = 0.04, y = 0.64, width = 0.01, height = 0.015)+

     # draw_grob(grid::rectGrob(gp = grid::gpar(col = "grey50", fill = "grey50")), x = 0.04, y = 0.474, width = 0.01, height = 0.015) + 
       # draw_grob(grid::rectGrob(gp = grid::gpar(col = "pink", fill = "pink")), x = 0.04, y = 0.308, width = 0.01, height = 0.015) + 
         # draw_grob(grid::rectGrob(gp = grid::gpar(col = "#009E73", fill = "#009E73")), x = 0.04, y = 0.142, width = 0.01, height = 0.015)






dev.off()
```


2024-08-08 -- sliding plots for each cluster individually
```{r}
cluster_annotations <- c(
  "clust_1814" = "ISIP-1",
  "clust_343" = "ISIP-2A",
  "clust_1154" = "ISIP-3",
  "clust_1820" = "Plastocyanin",
  "clust_534" = "Flavodoxin",
  "clust_231" = "Rhodopsin",

  "clust_332" = "LHC",
      "clust_3986" = "LHC", # double check
      "clust_5925" = "LHC",
      "clust_1044" = "LHC",
      "clust_195" = "LHC",
      "clust_1236" = "LHC",
      "clust_3440" = "LHC",
      "clust_3658" = "LHC",
      "clust_3703" = "LHC",
      "clust_3551" = "LHC",
  
    "clust_216" = "Pyrophosphatase",  "clust_522" = "Pyrophosphatase",  "clust_868" = "Pyrophosphatase",  "clust_522" = "Pyrophosphatase",


  
  "clust_55" = "PSI PsaA/PsaB",
  "clust_2815" = "PSI RC",
  "clust_1906" = "Mn-cluster",
      "clust_2763" = "Mn-cluster",
  "clust_613" = "Cyt-b6",
  "clust_1497" = "Cyt-f",
  "clust_42" = "PSII RC",
  "clust_1978" = "PSII PsbU",
  "clust_498" = "PSII CP47",
  "clust_211" = "PSII PsbC",
  "clust_2490" = "Rubisco",
  "clust_655" = "CA_clust_655",
  "clust_6573" = "Cd-CA",
  "clust_1525" = "SLC-4",
  "clust_1886" = "SLC",
  "clust_411" = "Nitrate reductase",
  "clust_417" = "AMT",
  "clust_891" = "Glutamate Synthase",
  "clust_1372" = "Urea transporter",
  "clust_320" = "Urea carboxylase",
  "clust_6019" = "Ornithine cyclodeaminase",
  "clust_955" = "Mg chelatase",
      "clust_416" = "Mg chelatase",
      "clust_2922" = "Mg chelatase",

  "clust_5562" = "Cu-Zn SOD", "clust_2501" = "Fe-Mn SOD", "clust_1289" = "DoxX", 
  "clust_520" = "Alkyl hydroperoxide redct-peroxiredoxin", 
      "clust_9714" = "Peroxiredoxin",
  
  "clust_2244" = "PEPC",
  "clust_2088" = "Fasciclin", "clust_6071" = "Fasciclin2",
  "clust_2460" = "Ferredoxin",
  "clust_2934" = "Multicopper Oxidase",
  "clust_753" = "Fructose bisphosphatase",
  "clust_3563" = "Nitrite/Sulphite reductase",
  "clust_997" = "APRT",
  "clust_10365" = "PET8",
  "clust_1550" = "Cyt-b559",
  "clust_1322" = "Porphobilinogen deaminase"
)






cluster_metal <- c(
  "PSI PsaA/PsaB" = "Fe",
  "PSII PsbU" = "Fe", 
  "PSII CP47" = "Fe", 
  "PSII RC" = "Fe",
  "PSI RC" = "Fe",
  "Ferredoxin" = "Fe",
  "Mn-cluster"= "Mn",
  "Cyt-b6" = "Fe", 
  "Cyt-f" = "Fe", 
  "Plastocyanin" = "Cu",
  "Cu-Zn SOD" = "Cu", 
  "Fe-Mn SOD" = "Fe", 
  "Nitrite/Sulphite reductase" =  "Fe")

cluster_metal2 <- c(
  "Cu-Zn SOD" = "Zn", 
  "Fe-Mn SOD" = "Mn")


BA_DE$cluster_metal <- cluster_metal[BA_DE$cluster_annotation]
BA_DE$cluster_metal2 <- cluster_metal2[BA_DE$cluster_annotation]




cluster_order <- c('LHC', #split into PSI and PSII LHCs?
                   
                   "PSII PsbU", 
                   'PSII CP47',
                   'PSII RC',
                   'Cyt-b559',

                   'Mn-cluster', 
                   'Plastocyanin',  
                   'Rhodopsin',
                   'Flavodoxin', 
                   'Ferredoxin', 
                   "Cyt-b6", 
                   'Cyt-f',
                   'Mg chelatase',

                   'PSI RC',
                   'PSI PsaA/PsaB',
                 
                   'ISIP-1',
                   'ISIP-2A', 
                   'ISIP-3',
                   
                   'Fasciclin', 'Fasciclin2',

                   'Cu-Zn SOD',
                   'Fe-Mn SOD',
                    "DoxX",
                   'Cd-CA', 
                   
                   'Nitrate reductase',
                   'Nitrite/Sulphite reductase', 
                   'Urea transporter', "AMT")


clustersofinterest <- paste(cluster_order, collapse = "|")


```

```{r}
#pfam_PF00006_ATP synthase alpha/beta family, nucleotide-binding domain ?
#proteasome 
# porins - clustclust_3629

protein <- filter (proteindata_unnorm, grepl ("clust_110|clust_11709", uniq_cluster))


clusters <- proteindata_unnorm [c(4, 37, 40:69)]

clusters2 <- melt(clusters, id.vars=c("uniq_cluster", "F")) %>%
            rename (treatment = variable) %>%
            rename (unnorm_abundance = value ) %>% 
            separate(treatment, into=c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"), sep = "_")

clusters2$F <- ifelse (clusters2$F  ==  "Polar-centric-Mediophyceae", "Centric diatom", 
                  ifelse (clusters2$F  ==  "Radial-centric-basal-Coscinodiscophyceae", "Centric diatom", 
                  ifelse( clusters2$F  ==  "Raphid-pennate", "Pennate diatom", 
                  ifelse(clusters2$F  ==  "Araphid-pennate", "Pennate diatom", 
                  ifelse(clusters2$F  ==  "Phaeocystaceae", "Phaeocystis",
                  clusters2$F)))))

clusters2_aggregated <- aggregate(clusters2$unnorm_abundance, 
                                  by=list(Category=clusters2$uniq_cluster, clusters2$F, clusters2$bioassay), FUN=sum) %>%
                        rename(unnorm_abundance = x) %>%
                        rename(cluster = Category) %>%
                        rename(taxon_F = Group.2) %>%
                        rename(treatment = Group.3)

clusters3 <- dcast(clusters2_aggregated, taxon_F+cluster~treatment) #%>% mutate_at (3, funs ((. / sum(.))))

clusters3 <- clusters3 %>% #this creates a new column of normalized protein abundances, either normalized to total protein or group normalized. 
  group_by(taxon_F) %>% 
  mutate(Offshore_grpnorm = Offshore / sum(Offshore)) %>% 
  ungroup()

clusters3 <- clusters3 %>% #this creates a new column of normalized protein abundances, either normalized to total protein or group normalized. 
  group_by(taxon_F) %>% 
  mutate(Nearshore_grpnorm = Nearshore / sum(Nearshore)) %>% 
  ungroup()
#because some clusters are not found in any taxa, we have NA values in the grpnorm column. (dividing by 0 )

#Step 2: get the DE and fold change for taxon-normalized clusters
BA_DE$taxon_F <- ifelse (BA_DE$taxon_F  ==  "centric", "Centric diatom", 
                  ifelse (BA_DE$taxon_F  ==  "pennate", "Pennate diatom", 
                  ifelse(BA_DE$taxon_F  ==  "Phaeocystaceae", "Phaeocystis",
                  BA_DE$taxon_F)))

alldata <- merge(BA_DE, clusters3, by = c("cluster", "taxon_F"))
alldata <- alldata [c(1,27,2, 30, 31, 3:26)]
```

#2024-11-05 -- PSI:PSII ratios
2024-11-24 -- total PSI:PSII ratios
```{r}
#https://www.pnas.org/doi/full/10.1073/pnas.1810886116
#https://www.pnas.org/doi/full/10.1073/pnas.0711370105
#https://onlinelibrary.wiley.com/doi/full/10.1111/j.1529-8817.2011.01098.x#b55
#https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2022GB007382   #Nick might be a good reviewer, or Ben twining 
#https://www.nature.com/articles/nature02954  PSII:PSI 10:1 in coastal vs offshore

norm_proteomicsdata <- proteindata_unnorm  %>% mutate_at (40:69, funs ((. / sum(.)))) #this normalizes the data to total peptide abundance, 


PSII <- filter(norm_proteomicsdata, grepl("^clust_1784$", uniq_cluster))

PSII2 <- PSII [c(4, 40:69)]

PSII3 <- melt(PSII2, id.vars=c("uniq_cluster")) %>%
            rename(treatment = variable) %>%
            rename(norm_abundance = value ) %>% 
            separate(treatment, into =  c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"))  
 
PSII3$tempfe_treatment <- paste (PSII3$temperature_treatment, PSII3$iron_treatment, sep = "_")
  
PSII4 <- aggregate(PSII3$norm_abundance, by=list(PSII3$bioassay, PSII3$tempfe_treatment, PSII3$replicate, PSII3$uniq_cluster), FUN = sum) %>%
          rename (bioassay = Group.1, 
          tempfe_treatment = Group.2, 
          replicate = Group.3, 
          cluster = Group.4, 
          PSII_norm_abundance = x)

PSII4$normabunancadjusted <- ifelse (PSII4$cluster ==  "clust_1784", PSII4$PSII_norm_abundance*5.598e-11, 
                                    ifelse(PSII4$cluster == "clust_2763", PSII4$PSII_norm_abundance*6.6667e-11, 
                                           ifelse(PSII4$cluster == "clust_5856", PSII4$PSII_norm_abundance*4.347e-11, 
                                         "else")))
                                    
PSII4$tempfe_treatment <- ifelse (PSII4$tempfe_treatment ==  "T0_T0", "T0",
                                                 PSII4$tempfe_treatment)

alldata_PSII <- merge (nonproteindata, PSII4, by = c("bioassay", "tempfe_treatment", "replicate"), all=TRUE)





#  "clust_55" = "PSI",
#       "clust_2815" = "PSI",
#       "clust_1564" = "PSI",
#       "clust_2439" = "PSI",
#       "clust_2460" = "PSI", #21,342da
#       "clust_3396" = "PSI"



PSI <- filter(norm_proteomicsdata, grepl("^clust_2460$", uniq_cluster))

PSI2 <- PSI [c(4, 40:69)]

PSI3 <- melt(PSI2, id.vars=c("uniq_cluster")) %>%
            rename(treatment = variable) %>%
            rename(norm_abundance = value ) %>% 
            separate(treatment, into =  c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"))  
 
PSI3$tempfe_treatment <- paste (PSI3$temperature_treatment, PSI3$iron_treatment, sep = "_")
  
PSI4 <- aggregate(PSI3$norm_abundance, by=list(PSI3$bioassay, PSI3$tempfe_treatment, PSI3$replicate, PSI3$uniq_cluster), FUN = sum) %>%
          rename (bioassay = Group.1, 
          tempfe_treatment = Group.2, 
          replicate = Group.3, 
          cluster = Group.4, 
          PSI_norm_abundance = x)

#clust_55 mw = 83,710 Da 

# "clust_2460" = "PSI", #21,342da
#https://www.bioline.com/media/calculator/01_04.html

#1 ug of psbv  = 1.195 e-8 moles



PSI4$normabunancadjusted <- ifelse (PSI4$cluster ==  "clust_55", PSI4$PSI_norm_abundance*1.195e-8, 
                                    ifelse(PSI4$cluster == "clust_2460", PSI4$PSI_norm_abundance*4.686e-11, 
                                           ifelse(PSI4$cluster == "clust_5856", PSI4$PSI_norm_abundance*4.347e-11, 
                                         "else")))
                                    
PSI4$tempfe_treatment <- ifelse (PSI4$tempfe_treatment ==  "T0_T0", "T0",
                                                 PSI4$tempfe_treatment)

alldata_PSI <- merge (nonproteindata, PSI4, by = c("bioassay", "tempfe_treatment", "replicate"), all=TRUE)



alldata_PSII$psIIpsIratio <- alldata_PSI$normabunancadjusted/alldata_PSII$normabunancadjusted


psratio_boxplot <- ggplot(filter (alldata_PSII, grepl ("Fe|", full_treatment_name) & grepl ("BA2|", bioassay)), 
       aes( x = iron_treatment, 
          y = psIIpsIratio))+
    geom_boxplot() + 
    geom_point(size = 3, alpha =0.65, stroke =1.5, aes (shape = tempfe_treatment))+
        scale_shape_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c(5, 1, 19,1, 19)) +
      labs (x = "", y = "PSI:PSII" )+
   scale_x_discrete(limits = c("noFe", "Fe" ), 
                          labels=c('no added Fe', '+Fe'))+
   theme_bw()+
  facet_wrap(~bioassay) +
 theme(axis.text=element_text(face = "bold", size = 12, color = "black"),
        axis.title.y=element_text(size=12, color = "black", face = "bold")) + 
  ylim (0,1.2)+
      theme_bw()+
        theme(axis.text=element_text(face = "bold", size = 10, color = "black"),
        axis.title.y=element_text(size = 14, color = "black"),
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 10, face = "bold"),
        legend.position = "none", 
        panel.grid = element_blank()) 


psratio_Fe56 <- ggplot(filter (alldata_PSII, grepl ("_T0_", full_treatment_name) & grepl ("BA2|", bioassay)), 
       aes(x = (pFe56_nM)/PON_uM, 
          y = psIIpsIratio))+
      geom_point(size = 3, alpha =0.65, stroke =1.5, aes (shape = tempfe_treatment))+

          scale_shape_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c(5, 1, 19,1, 19)) +
  
      xlab (expression ("in situ"~ p ^{56}*"Fe":pN~~(mM:M))) +
     ylab (expression("PSI:PSII"))+ 
   ylim (0,1.2)+

     theme_bw()+
     theme(axis.text=element_text(face = "bold", size = 12, color = "black"),
        axis.title.y=element_text(size=12, color = "black", face = "bold")) +
      theme_bw()+
        theme(axis.text=element_text(face = "bold", size = 10, color = "black"),
        axis.title.x=element_text(size = 12, color = "black"),
        axis.title.y=element_text(size = 12, color = "black"),
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 10, face = "bold"),
        legend.position = "none",
        panel.grid = element_blank())

psratio_Fe57 <- ggplot(filter (alldata_PSII, !grepl ("_T0_", full_treatment_name) & grepl ("BA2|", bioassay)), 
       aes(x = (pFe57_nM)/PON_uM, 
          y = psIIpsIratio))+
      geom_point(size = 3, alpha =0.65, stroke =1.5, aes (shape = tempfe_treatment))+

          scale_shape_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c(5, 1, 19,1, 19)) +
  
      xlab (expression (p ^{57}*"Fe":pN~~(mM:M))) +
     ylab (expression("PSI:PSII"))+ 
    ylim (0,1.2)+

     theme_bw()+
     theme(axis.text=element_text(face = "bold", size = 12, color = "black"),
        axis.title.y=element_text(size=12, color = "black", face = "bold")) +
      theme_bw()+
        theme(axis.text=element_text(face = "bold", size = 10, color = "black"),
        axis.title.x=element_text(size = 12, color = "black"),
        axis.title.y=element_text(size = 12, color = "black"),
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 10, face = "bold"),
        legend.position = "none",
        panel.grid = element_blank())





```



PSI : PSII figure construction 
```{r}
ratiometalconc <- plot_grid(NULL, psratio_Fe56,
                            NULL, psratio_Fe57,
          ncol = 2, 
          align = "hv", 
          rel_widths = c(0.05,1),
          labels = c ("(B)", "",   "(C)"))


#png("../figures/ps117_PSIPSII_20241105.png", width=14, height=7, units="in", res=600) 

plot_grid (psratio_boxplot, NULL, ratiometalconc,
                labels = c("(A)", ""),
                ncol = 3, 
                rel_widths = c(1,0.03,1.5),
                align = "v") 
#dev.off()
```

       

#2024-08-20 -- Figure 5 Metal Calculations 
            -- Calculate theoretical maximum possible amount of Cd-CA that can be made
```{r}
norm_proteomicsdata <- proteindata_unnorm  %>% mutate_at (40:69, funs ((. / sum(.)))) #this normalizes the data to total peptide abundance, 
```
     
2024-11-05 -- setting up the data
```{r}
#This is the full calculation
#beta <- ((gamma * Pi_ms_signalratio) * (bulkprotein) * (1/MPi) * (sigma * nEjPi)) / particulateEj
#have to convert the pEj_nM to PEj_M (because we calculated the metal from protein in moles, we need to convert that to nM)

# beta = percentage of metalJ found in proteini; for example, what percent of pFe is in PSII? 
# gamma = correction factor for protein quant. do we normalize to hardklor? TIC, match abundance etc..; 0-1
# Pi_ms_signalratio = (summed intensities of peptides matching to proteini (i.e. Pi))/ (summed intensities of all matched peptides)
# bulkprotein = protein concentration per L of seawater (BCA, or from PON etc). right now it's from PON
# MPi = moles of Pi per ug protein
# sigma = metallation coifficient; 0-1, 1 = 100% metallation, 0 = no metallation. I do the caclulations with a range of this. 
# nEjPi = moles of metal (Ej) per moles of metalloprotein (Pi)  (e.g. 1 = one metal per protein)
# particulateEj =  molar concentration of particulate element of interest, mol L-1 


#1) gamma - the data is already normalized to total matched peptide abundance. can renormalize to gamma_hardklor or gamma_tic
gamma_df <- read.csv("../data/ps117_all_norm_factors_20240818.csv") %>% 
  group_by(bioassay, replicate, tempfe_treatment, timepoint) %>%
  summarize(hardklor = mean(all_possible_peptides, na.rm = TRUE), 
            TIC = mean(TIC, na.rm = TRUE), 
            matchedpeptides = mean (matched_peptides, na.rm = T))

gamma_df$gamma_hardklor <- gamma_df$matchedpeptides/gamma_df$hardklor
gamma_df$gamma_tic <- gamma_df$matchedpeptides/gamma_df$TIC

metalcalcs_df <- merge ((filter(nonproteindata, grepl("T0|T8", timepoint))), gamma_df, by = c("bioassay", "tempfe_treatment", "replicate", 
 "timepoint"), all = T) %>% 
                select(c(1:3, 17, 23, 36:43, 105, 106))

metalcalcs_df$pFe_nM <- metalcalcs_df$pFe56_nM+metalcalcs_df$pFe57_nM



#2) Pi_ms_signalratio - this will depend on the protein of interest (proteini)
#First make a list of metalloproteins of interest: 
cluster <- c ("clust_1906",
              "clust_2763",
              "clust_43119",#, #testing with one protein for now
              "clust_42",
              "clust_211",
              "clust_498",
              "clust_55",
              "clust_1820",
              "clust_6573",
              "clust_1784",
              "clust_1497",
              "clust_2460" ,
              "clust_1978")

#plastocyanin
#cluster <- c ("clust_1820")
              
proteini <- read.csv("../data/ps117_proteini_20241207.csv")


proteini <- data.frame(cluster)


proteini$protein_name <-  ifelse (cluster == 'clust_1906', "PSBO_OEC_PSII", 
                              ifelse (cluster == "clust_2763", "OEC_psbq_PSII", 
                              ifelse (cluster == "clust_43119", "OEC_psbq_PSII", 
                              ifelse(cluster == 'clust_42', "D1_PSII",
                              ifelse(cluster == 'clust_211', "D2_PSII",
                              ifelse(cluster == 'clust_6573', "CdCA",
                              ifelse(cluster == 'clust_1784', "PsbV_PSII",
                              ifelse(cluster == 'clust_1497', "cytf",
                              ifelse(cluster == 'clust_2460', "PsaD_PSI",
                              ifelse(cluster == 'clust_55', "PsaA/B_PSI",
                              ifelse(cluster == 'clust_2934', "multi_Cu_oxidase",
                              ifelse(cluster == 'clust_1820', "Pcyn",
                              ifelse(cluster == 'clust_1978', "PsbU_PSII",
                              ifelse(cluster == 'clust_498', "CP47_PSII",
                                      'somethingwrong'))))))))))))))

#3) bulkprotein - calculate using PON. this is the total amount of particulate protein
metalcalcs_df$bulkprotein <- metalcalcs_df$PON_ugL *4.78

#4) MPi - moles of proteini per ug protein
#https://www.bioline.com/media/calculator/01_04.html
proteini$MPi <- ifelse (cluster == 'clust_1906', 3.0303e-11, 
                         ifelse (cluster == "clust_2763", 5.882e-11, 
                              ifelse (cluster == "clust_43119", 5.882e-11, 

                              ifelse(cluster == 'clust_42', 2.632e-11,
                              ifelse(cluster == 'clust_211', 2.532e-11, ##this needs to be fixed, clust_211  contains another protien that is not D2, so look for D2 proteins here! look for "psbD in uniq_kegg_desc
                              ifelse(cluster == 'clust_6573', 1.454e-11,
                              ifelse(cluster == 'clust_1784', 5.598e-11,
                              ifelse(cluster == 'clust_1497', 3.195e-11,
                              ifelse(cluster == 'clust_2460', 4.686e-11,
                              ifelse(cluster == 'clust_1820', 1.0204e-10,
                              ifelse(cluster == 'clust_1978', "notsure",
                              ifelse(cluster == 'clust_498', 1.786e-11,
                              ifelse(cluster == 'clust_2934', 1.768e-11,

                              ifelse(cluster == 'clust_55', 1.22e-11,
                                      'somethingwrong'))))))))))))))

#5) #this accounts for number of protein copies in each group. for example, there are 3 Fe atoms in PSII, but if there are three PsaDs in PSII, then the nEjPi_Cd for PsaD will be 1
proteini$nEjPi_Cd <- ifelse (cluster == 'clust_1906', 0,
                              ifelse (cluster == "clust_43119", 0,
                              ifelse(cluster == 'clust_42', 0,
                              ifelse(cluster == 'clust_211', 0,
                              ifelse(cluster == 'clust_6573', 1,
                              ifelse(cluster == 'clust_1784', 0,
                              ifelse(cluster == 'clust_55', 0,
                              ifelse(cluster == 'clust_2763', 0,


                              ifelse(cluster == 'clust_1497', 0,
                              ifelse(cluster == 'clust_2460', 0,
                              ifelse(cluster == 'clust_1820', 0,
                              ifelse(cluster == 'clust_1978', 0,
                              ifelse(cluster == 'clust_498', 0,
                                      'somethingwrong')))))))))))))

#https://pmc.ncbi.nlm.nih.gov/articles/PMC9502191/ copper paper. 50-70% Cu in plastocyanin
proteini$nEjPi_Cu <- ifelse (cluster == 'clust_1906', 0,
                              ifelse (cluster == "clust_43119", 0,
                              ifelse(cluster == 'clust_42', 0,
                              ifelse(cluster == 'clust_211', 0,
                              ifelse(cluster == 'clust_6573', 0,
                              ifelse(cluster == 'clust_1784', 0,
                              ifelse(cluster == 'clust_1497', 0,
                              ifelse(cluster == 'clust_2460', 0,
                              ifelse(cluster == 'clust_1820', 1,
                              ifelse(cluster == 'clust_2934', 4,
                              ifelse(cluster == 'clust_55', 0,
                              ifelse(cluster == 'clust_2763', 0,

                              ifelse(cluster == 'clust_1978', 0,
                              ifelse(cluster == 'clust_498', 0,
                                      'somethingwrong'))))))))))))))

proteini$nEjPi_Mn <- ifelse (cluster == 'clust_1906', 2, #there are 4Mn total in Mn-cluster, but there are two psbO  proteins in OEC
                              ifelse (cluster == "clust_2763", 4,
                              ifelse (cluster == "clust_43119",4,  
                              ifelse(cluster == 'clust_42', 4,#there is one D1 protein in PSII
                              ifelse(cluster == 'clust_211', 4,
                              ifelse(cluster == 'clust_6573', 0,
                              ifelse(cluster == 'clust_1784', 4,
                              ifelse(cluster == 'clust_1497', 0,
                              ifelse(cluster == 'clust_2460', 0,
                              ifelse(cluster == 'clust_1820', 0,
                              ifelse(cluster == 'clust_1978', 0,
                              ifelse(cluster == 'clust_55', 0,
                              ifelse(cluster == 'clust_498', 0,
                                      'somethingwrong')))))))))))))

proteini$nEjPi_Fe <- ifelse (cluster == 'clust_1906', 1.5, #assuming 3Fe atoms in PSII 
                              ifelse (cluster == "clust_2763", 3,
                              ifelse(cluster == 'clust_42', 3,
                              ifelse(cluster == 'clust_211', 3,
                              ifelse(cluster == 'clust_6573', 0,
                              ifelse(cluster == 'clust_1784', 3,
                              ifelse(cluster == 'clust_1497', 6,
                              ifelse(cluster == 'clust_2460', 12,
                              ifelse(cluster == 'clust_55', 24,
                              ifelse(cluster == 'clust_1820', 0,
                              ifelse(cluster == 'clust_1978', 0,
                              ifelse(cluster == 'clust_498', 3,
                              ifelse(cluster == 'clust_43119', 0,
                                      'somethingwrong')))))))))))))

# proteini$nEjPi_chla <- ifelse (cluster == 'clust_1906', 40, #assuming 3Fe atoms in PSII 
#                               ifelse (cluster == "clust_2763", 40,
#                               ifelse(cluster == 'clust_42', 40,
#                               ifelse(cluster == 'clust_211', 40,
#                               ifelse(cluster == 'clust_6573', 0,
#                               ifelse(cluster == 'clust_1784', 40,
#                               ifelse(cluster == 'clust_1497', 0,
#                               ifelse(cluster == 'clust_2460', 100,
#                               ifelse(cluster == 'clust_55', 100,
#                               ifelse(cluster == 'clust_1820', 0,
#                               ifelse(cluster == 'clust_1978', 40,
#                               ifelse(cluster == 'clust_498', 40,
#                                       'somethingwrong'))))))))))))


#https://www.bioline.com/media/calculator/01_04.html
#clust_1906, psbO, OEC = 3.0303e-11, two proteins per OEC(4Mn)
#clust_42, D1 protein = 2.632e-11
#clust_211, D2 protein = 2.532e-11 #this needs to be fixed, the 211 cluster contains another protien that is not D2, so look for D2 proteins here! look for "psbD in uniq_kegg_desc
#clust_6573, cdCa = 1.454e-11, 1Cd per protein
#clust_1784, PsbV; cytochrome c550, PSII = 5.598e-11, 3 Fe per protein
#clust_1497, ETC, apocytochrome f = 3.195e-11, 5-7Fe per protein
#clust_2460, PSI PsaD = 4.686e-11, 12 Fe per protein
#clust_1820, plastocyanin = 1.0204e-10, 1 Cu per protein
#clust_1978, psbU PSII
#clust_498, CP47, ~6 per PSII
```
     
2024-11-05 -- do the calculation                
```{r}
#first, choose the clusters that represent the metalloproteins of interest, and all the intensities of all the peptides in those clusters
ms_signalratio <- filter(norm_proteomicsdata, uniq_cluster %in% proteini$uniq_cluster) %>% 
                  filter(!grepl ( "PSBC, photosystem II 44" , uniq_kegg_desc)) %>%
                  filter(!grepl("clust_43119", uniq_cluster) | !grepl("unassigned_PFams_desc", uniq_PFams_desc)) %>%

                  select(c(4, 40:69)) %>% #then select relavant columns 
                  melt(id.vars=c("uniq_cluster")) %>%
                  rename(treatment = variable, norm_abundance = value) %>%
                  separate(treatment, into =  c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate")) %>% 
                  mutate(tempfe_treatment = paste(temperature_treatment, iron_treatment, sep = "_")) %>%
                  
             
                  group_by(bioassay, tempfe_treatment, replicate, uniq_cluster) %>% #group things so we can add the abundance of each peptide within each cluster and within each treatment
                  summarize(Pi_ms_signalratio = sum(norm_abundance, na.rm = TRUE), .groups = 'drop') %>% #add things up
                  #rename(cluster = uniq_cluster) %>% 
                  mutate(tempfe_treatment = ifelse(tempfe_treatment == "T0_T0", "T0", tempfe_treatment))

     # mutate (protein_name = ifelse (uniq_cluster == 'clust_1906', "PSBO_OEC_PSII", 
     #                          ifelse (uniq_cluster == "clust_2763", "OEC_psbq_PSII", 
     #                          ifelse (uniq_cluster == "clust_43119", "OEC_psbq_PSII", 
     #                          ifelse(uniq_cluster == 'clust_42', "D1_PSII",
     #                          ifelse(uniq_cluster == 'clust_211', "D2_PSII",
     #                          ifelse(uniq_cluster == 'clust_6573', "CdCA",
     #                          ifelse(uniq_cluster == 'clust_1784', "PsbV_PSII",
     #                          ifelse(uniq_cluster == 'clust_1497', "cytf",
     #                          ifelse(uniq_cluster == 'clust_2460', "PsaD_PSI",
     #                          ifelse(uniq_cluster == 'clust_55', "PsaA/B_PSI",
     #                          ifelse(uniq_cluster == 'clust_498', "CP47_PSII",
     #                          ifelse(cluster == 'clust_2934', "multi_Cu_oxidase",
     #                          ifelse(uniq_cluster == 'clust_1820', "Pcyn",
     #                          ifelse(uniq_cluster == 'clust_1978', "PsbU_PSII",
     #                                  'somethingwrong'))))))))))))))) %>%



#Now we have a dataframe with values that represent Pi_ms_signalratio, so we merge it with the rest of the information to do the calculation 
metalcalcs_df2 <-  merge (metalcalcs_df, ms_signalratio,  by = c("bioassay", "tempfe_treatment", "replicate"), all = T)

metalcalcs_df3 <-  merge (metalcalcs_df2, proteini,  by = c("uniq_cluster"), all = T)

metalcalcs_df3$MPi <- as.numeric(metalcalcs_df3$MPi) #need to revisit this to make sure the numerics are good 
metalcalcs_df3$nEjPi_Mn <- as.numeric(metalcalcs_df3$nEjPi_Mn)
metalcalcs_df3$nEjPi_Fe <- as.numeric(metalcalcs_df3$nEjPi_Fe)
metalcalcs_df3$nEjPi_Cu <- as.numeric(metalcalcs_df3$nEjPi_Cu)
#metalcalcs_df3$nEjPi_chla <- as.numeric(metalcalcs_df3$nEjPi_chla)
#metalcalcs_df3$total_chla_ug_L <- as.numeric(metalcalcs_df3$total_chla_ug_L)


### now that the data is in order, we can perform the calculation 

#a reminder of the paramters
#beta <- ((gamma * Pi_ms_signalratio) * (bulkprotein) * (1/MPi) * (sigma * nEjPi)) / particulateEj
#have to convert the pEj_nM to PEj_M (because we calculated the metal from protein in moles, we need to convert that to nM)

# beta = percentage of metalJ found in proteini; for example, what percent of pFe is in PSII? 
# gamma = correction factor for protein quant. do we normalize to hardklor? TIC, match abundance etc..; 0-1
# Pi_ms_signalratio = (summed intensities of peptides matching to proteini (i.e. Pi))/ (summed intensities of all matched peptides)
# bulkprotein = protein concentration per L of seawater (BCA, or from PON etc). right now it's from PON
# MPi = moles of Pi per ug protein
# sigma = metallation coifficient; 0-1, 1 = 100% metallation, 0 = no metallation. I do the caclulations with a range of this. 
# nEjPi = moles of metal (Ej) per moles of metalloprotein (Pi)  (e.g. 1 = one metal per protein)
# particulateEj =  molar concentration of particulate element of interest, mol L-1 

# Function to perform the calculation
nEjPi <- "nEjPi_Fe" 
particulateEj <- "pFe_nM" 
gamma <- "gamma_hardklor"

perform_calculation <- function(sigma_value, df) {
  # Randomly sample one Pi_ms_signalratio from the replicates
  sampled_ratio <- sample(df$Pi_ms_signalratio, 1, replace = TRUE) #instead of picking one ratio, randomly select a replicate with every calculation 

#100* is for % 
#the 1e9 is to go to nM 
  beta <- 100*(1e9*df[[gamma]] * (sampled_ratio * df$MPi) * (df$bulkprotein) * (sigma_value * df[[nEjPi]])) / (df[[particulateEj]])
  return(mean(beta, na.rm = TRUE))  # Return mean for simplicity
}

# Calculate and store results for each treatment separately
set.seed(123)  # For reproducibility

results <- metalcalcs_df3 %>%
  group_by(tempfe_treatment, bioassay, uniq_cluster) %>%
  do({
    treatment_df <- .
    sigma_results <- replicate(2000, {
      sigma_random <- runif(1, 1, 1)  # Random sigma between 0.6 and 1
      perform_calculation(sigma_random, treatment_df)
    })
    data.frame(tempfe_treatment = unique(treatment_df$tempfe_treatment), 
               bioassay = unique(treatment_df$bioassay), 
               results = sigma_results)
  }) %>%
filter(!is.na(results))  # Remove NA values



results <- results %>%
  separate(tempfe_treatment, into = c("temp_treatment", "iron_treatment"), sep = "_")

#png("../figures/ps117_Cu-pcyn_20241117.png", width= 8, height= 8, units="in", res=600) 
ggplot(filter (results, !grepl ("T0", temp_treatment) & !grepl ("CdCA", protein_name)), aes(x = protein_name, y = results, shape = iron_treatment, color =temp_treatment)) +
  #geom_violin(aes( fill = iron_treatment), alpha = 0.5, size = 1, trim = F, scale = 'width') + 
 #geom_flat_violin (scale = 'width', trim = F, aes( fill = iron_treatment), alpha = 0.5, size = 1)  +
# geom_boxplot(width = 0.2, position = position_nudge(x = 0), alpha = 0.7, size = 0.5, show.legend = F ) +
    #geom_point(size = 3, stroke = 1 )+

  facet_wrap(~ bioassay) +
  theme_bw() + 
  scale_shape_manual (name = "",
                       breaks = c("noFe","Fe", "T0"), 
                       labels = c("noFe", "+Fe", "T0"),
                       values = c(1, 19, 5))+
  ylab (expression ("Fe allocation to PSII and PSII (%)")) + 
  xlab (expression ("")) +
 #  scale_x_discrete (limits=c("Pcyn", "CP47_PSII"))+
 #                   labels = c("T0" ))+
  
#  scale_x_discrete (limits=c("CP47_PSII", "D1_PSII", "D2_PSII", "PSBO_OEC_PSII", "PsaA/B_PSI", "PsaD_PSI"))+
 #                   labels = c("T0" ))+
  
  
  
  stat_summary(fun = mean, geom = "point", show.legend = F, size = 3, stroke = 1, alpha = 0.9)+
  stat_summary(fun.data =  mean_sd, geom = "errorbar", size = 0.7, width = 0.2) +
  scale_color_manual(name="Temperature",
                     breaks = c("LT","HT", "T0"),
                     labels = c("LT", "HT", "T0"),
                     values = c('black', "grey50", "black")) + #D55E00

  
  
  #scale_x_discrete(labels = c("D1 (PSII)",  "PsbQ (OEC)", "PsbO (OEC)")) +  
 # scale_x_discrete(labels = c("CP47",  "D1", "D2", "PsbO")) +  
  # scale_x_discrete(labels = c("PsaA/PsaB",  "PsaD")) +  


    theme_bw()+
# scale_y_continuous(limits = c(0, 150), breaks = seq(0, 150, by = 25)) +
  theme(axis.text=element_text(face = "bold", size = 10, color = "black"),
        axis.title.x=element_text(size = 10, color = "black"), 
        axis.title.y=element_text(size = 13, color = "black", face = "bold"),
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 10, face = "bold"),
       # legend.position = "none", 
        panel.grid = element_blank()) 

#dev.off()

```

2024-12-07 -- new way of doing it 
```{r}
norm_proteomicsdata <- proteindata_unnorm  %>% mutate_at (40:69, funs ((. / sum(.)))) #this normalizes the data to total peptide abundance, 

#This is the full calculation
#beta <- ((gamma * Pi_ms_signalratio) * (bulkprotein) * (1/MPi) * (sigma * nEjPi)) / particulateEj
#have to convert the pEj_nM to PEj_M (because we calculated the metal from protein in moles, we need to convert that to nM)

# beta = percentage of metalJ found in proteini; for example, what percent of pFe is in PSII? 
# gamma = correction factor for protein quant. do we normalize to hardklor? TIC, match abundance etc..; 0-1
# Pi_ms_signalratio = (summed intensities of peptides matching to proteini (i.e. Pi))/ (summed intensities of all matched peptides)
# bulkprotein = protein concentration per L of seawater (BCA, or from PON etc). right now it's from PON
# MPi = moles of Pi per ug protein
# sigma = metallation coifficient; 0-1, 1 = 100% metallation, 0 = no metallation. I do the caclulations with a range of this. 
# nEjPi = moles of metal (Ej) per moles of metalloprotein (Pi)  (e.g. 1 = one metal per protein)
# particulateEj =  molar concentration of particulate element of interest, mol L-1 


#1) gamma - the data is already normalized to total matched peptide abundance. can renormalize to gamma_hardklor or gamma_tic
gamma_df <- read.csv("../data/ps117_all_norm_factors_20240818.csv") %>% 
  group_by(bioassay, replicate, tempfe_treatment, timepoint) %>%
  summarize(hardklor = mean(all_possible_peptides, na.rm = TRUE), 
            TIC = mean(TIC, na.rm = TRUE), 
            matchedpeptides = mean (matched_peptides, na.rm = T))

gamma_df$gamma_hardklor <- gamma_df$matchedpeptides/gamma_df$hardklor
gamma_df$gamma_tic <- gamma_df$matchedpeptides/gamma_df$TIC

metalcalcs_df <- merge ((filter(nonproteindata, grepl("T0|T8", timepoint))), gamma_df, by = c("bioassay", "tempfe_treatment", "replicate", 
 "timepoint"), all = T) %>% 
                select(c(1:3, 17, 23, 36:43, 105, 106))

metalcalcs_df$pFe_nM <- metalcalcs_df$pFe56_nM+metalcalcs_df$pFe57_nM

metalcalcs_df$bulkprotein <- metalcalcs_df$PON_ugL *4.78



ms_signalratio <- filter(norm_proteomicsdata, uniq_cluster %in% proteini$uniq_cluster) %>% 
                  filter(!grepl ( "PSBC, photosystem II 44" , uniq_kegg_desc)) %>%
                  filter(!grepl("clust_43119", uniq_cluster) | !grepl("unassigned_PFams_desc", uniq_PFams_desc)) %>%

                  select(c(4, 40:69)) %>% #then select relavant columns 
                  melt(id.vars=c("uniq_cluster")) %>%
                  rename(treatment = variable, norm_abundance = value) %>%
                  separate(treatment, into =  c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate")) %>% 
                  mutate(tempfe_treatment = paste(temperature_treatment, iron_treatment, sep = "_")) %>%
                  
             
                  group_by(bioassay, tempfe_treatment, replicate, uniq_cluster) %>% #group things so we can add the abundance of each peptide within each cluster and within each treatment
                  summarize(Pi_ms_signalratio = sum(norm_abundance, na.rm = TRUE), .groups = 'drop') %>% #add things up
                  #rename(cluster = uniq_cluster) %>% 
                  mutate(tempfe_treatment = ifelse(tempfe_treatment == "T0_T0", "T0", tempfe_treatment))


#Now we have a dataframe with values that represent Pi_ms_signalratio, so we merge it with the rest of the information to do the calculation 
metalcalcs_df2 <-  merge (metalcalcs_df, ms_signalratio,  by = c("bioassay", "tempfe_treatment", "replicate"), all = T)

metalcalcs_df3 <-  merge (metalcalcs_df2, proteini,  by = c("uniq_cluster"), all = T)

metalcalcs_df3$MPi <- as.numeric(metalcalcs_df3$MPi) #need to revisit this to make sure the numerics are good 
metalcalcs_df3$nEjPi_Mn <- as.numeric(metalcalcs_df3$nEjPi_Mn)
metalcalcs_df3$nEjPi_Fe <- as.numeric(metalcalcs_df3$nEjPi_Fe)
metalcalcs_df3$nEjPi_Cu <- as.numeric(metalcalcs_df3$nEjPi_Cu)
#metalcalcs_df3$nEjPi_chla <- as.numeric(metalcalcs_df3$nEjPi_chla)
#metalcalcs_df3$total_chla_ug_L <- as.numeric(metalcalcs_df3$total_chla_ug_L)


### now that the data is in order, we can perform the calculation 

#a reminder of the paramters
#beta <- ((gamma * Pi_ms_signalratio) * (bulkprotein) * (1/MPi) * (sigma * nEjPi)) / particulateEj
#have to convert the pEj_nM to PEj_M (because we calculated the metal from protein in moles, we need to convert that to nM)

# beta = percentage of metalJ found in proteini; for example, what percent of pFe is in PSII? 
# gamma = correction factor for protein quant. do we normalize to hardklor? TIC, match abundance etc..; 0-1
# Pi_ms_signalratio = (summed intensities of peptides matching to proteini (i.e. Pi))/ (summed intensities of all matched peptides)
# bulkprotein = protein concentration per L of seawater (BCA, or from PON etc). right now it's from PON
# MPi = moles of Pi per ug protein
# sigma = metallation coifficient; 0-1, 1 = 100% metallation, 0 = no metallation. I do the caclulations with a range of this. 
# nEjPi = moles of metal (Ej) per moles of metalloprotein (Pi)  (e.g. 1 = one metal per protein)
# particulateEj =  molar concentration of particulate element of interest, mol L-1 

# Function to perform the calculation
nEjPi <- "nEjPi_Fe" 
particulateEj <- "pFe_nM" 
gamma <- "gamma_hardklor"

perform_calculation <- function(sigma_value, df) {
  # Randomly sample one Pi_ms_signalratio from the replicates
  sampled_ratio <- sample(df$Pi_ms_signalratio, 1, replace = TRUE) #instead of picking one ratio, randomly select a replicate with every calculation 

#100* is for % 
#the 1e9 is to go to nM 
  beta <- 100*(1e9*df[[gamma]] * (sampled_ratio * df$MPi) * (df$bulkprotein) * (sigma_value * df[[nEjPi]])) / (df[[particulateEj]])
  return(mean(beta, na.rm = TRUE))  # Return mean for simplicity
}

# Calculate and store results for each treatment separately
set.seed(123)  # For reproducibility

results <- metalcalcs_df3 %>%
  group_by(tempfe_treatment, bioassay, uniq_cluster) %>%
  do({
    treatment_df <- .
    sigma_results <- replicate(500, {
      sigma_random <- runif(1, 0.8, 1)  # Random sigma between 0.6 and 1
      perform_calculation(sigma_random, treatment_df)
    })
    data.frame(tempfe_treatment = unique(treatment_df$tempfe_treatment), 
               bioassay = unique(treatment_df$bioassay), 
               results = sigma_results)
  }) %>%
filter(!is.na(results))  # Remove NA values

results <- results %>%
                    separate(tempfe_treatment, into = c("temp_treatment", "iron_treatment"), sep = "_")


results <- merge(results, proteini[, c("uniq_cluster", "protein_name")], by = "uniq_cluster", all.x = TRUE)
  


#png("../figures/ps117_Cu-pcyn_20241117.png", width= 8, height= 8, units="in", res=600) 
ggplot(filter (results, !grepl ("T0", temp_treatment)), aes(x = protein_name, y = results, shape = iron_treatment, color =temp_treatment)) +
  facet_wrap(~ bioassay) +
  theme_bw() + 
  scale_shape_manual (name = "",
                       breaks = c("noFe","Fe", "T0"), 
                       labels = c("noFe", "+Fe", "T0"),
                       values = c(1, 19, 5))+
  ylab (expression ("Fe allocation to PSII and PSII (%)")) + 
  xlab (expression ("")) +
  stat_summary(fun = mean, geom = "point", show.legend = F, size = 3, stroke = 1, alpha = 0.9)+
  stat_summary(fun.data =  mean_sd, geom = "errorbar", size = 0.7, width = 0.2) +
  scale_color_manual(name="Temperature",
                     breaks = c("LT","HT", "T0"),
                     labels = c("LT", "HT", "T0"),
                     values = c('black', "grey50", "black")) + #D55E00
  #scale_x_discrete(limits = c("clust_498",  "clust_42", "clust_211")) +  
   scale_x_discrete (limits=c("CP47_PSII", "D1_PSII", "D2_PSII", "PSBO_OEC_PSII", "PsaA/B_PSI", "PsaD_PSI"))+
  theme_bw()+
  theme(axis.text=element_text(face = "bold", size = 10, color = "black"),
        axis.title.x=element_text(size = 10, color = "black"), 
        axis.title.y=element_text(size = 13, color = "black", face = "bold"),
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 10, face = "bold"),
       # legend.position = "none", 
        panel.grid = element_blank()) 

#dev.off()



#info about proteins and whatnot
#https://www.bioline.com/media/calculator/01_04.html
#clust_1906, psbO, OEC = 3.0303e-11, two proteins per OEC(4Mn)
#clust_42, D1 protein = 2.632e-11
#clust_211, D2 protein = 2.532e-11 #this needs to be fixed, the 211 cluster contains another protien that is not D2, so look for D2 proteins here! look for "psbD in uniq_kegg_desc
#clust_6573, cdCa = 1.454e-11, 1Cd per protein
#clust_1784, PsbV; cytochrome c550, PSII = 5.598e-11, 3 Fe per protein
#clust_1497, ETC, apocytochrome f = 3.195e-11, 5-7Fe per protein
#clust_2460, PSI PsaD = 4.686e-11, 12 Fe per protein
#clust_1820, plastocyanin = 1.0204e-10, 1 Cu per protein
#clust_1978, psbU PSII
#clust_498, CP47, ~6 per PSII
```




pick a random protein abundance from the triplicates every time a calcualtion is made
```{r}
# Define the perform_calculation function
perform_calculation <- function(sigma_value, df) {
  # Randomly sample one Pi_ms_signalratio
  sampled_ratio <- sample(df$Pi_ms_signalratio, 1, replace = TRUE)
  
  # Perform the calculation using the sampled value
  beta_Mn_OEC <- (1e9 * (sampled_ratio * MPi) * (df$bulkprotein) * (sigma_value * nEjPi)) / (df$pMn_nM)
  
  # Return the mean of beta_Mn_OEC
  return(mean(beta_Mn_OEC, na.rm = TRUE))  # Return mean for simplicity
}

# Calculate and store results for each treatment separately
set.seed(123)  # For reproducibility

results <- alldata_PiEj %>%
  group_by(iron_treatment, bioassay, temperature_treatment) %>%
  do({
    treatment_df <- .
    sigma_results <- replicate(2000, {
      sigma_random <- runif(1, 0.6, 1)  # Random sigma between 0.6 and 1
      
      # Call the perform_calculation function with the random sigma and treatment dataframe
      perform_calculation(sigma_random, treatment_df)
    })
    
    # Return a dataframe with the results
    data.frame(iron_treatment = unique(treatment_df$iron_treatment), 
               bioassay = unique(treatment_df$bioassay), 
               results = sigma_results)
  }) %>%
  filter(!is.na(results))  # Remove NA values



  
ggplot(results, aes(x = iron_treatment, y = results)) +
  geom_violin(scale = 'width', trim = F) + 
  stat_summary(fun = "mean",
               geom = "crossbar", 
               width = 0.5,
               colour = "red")+
  
    #geom_jitter(alpha = 0.1, height = 0, width = 0.1)+

#  geom_point(alpaha = 0.2)+
 # geom_histogram(alpha = 0.7, bins = 100) +
  #geom_density(alpha = 0.5)+
 # geom_vline(data = mean_results, aes(xintercept = mean_result, color = iron_treatment),
            # linetype = "dashed", size = 1) +
 # geom_errorbarh(data = mean_results, aes(xmin = mean_result - sem, xmax = mean_result + sem, y = 0),
                # height = 0.2, color = "black") +
  facet_wrap(~ bioassay) +
  theme_bw()+ 
  theme()



# Plot histograms with mean lines and error bars
ggplot(results, aes(y = iron_treatment, x = results, fill = iron_treatment)) +
  geom_density_ridges()+
  geom_point(alpha = 0.2)+
  theme_ridges()+
 # geom_histogram(alpha = 0.7, bins = 100) +
  #geom_density(alpha = 0.5)+
 # geom_vline(data = mean_results, aes(xintercept = mean_result, color = iron_treatment),
            # linetype = "dashed", size = 1) +
 # geom_errorbarh(data = mean_results, aes(xmin = mean_result - sem, xmax = mean_result + sem, y = 0),
                # height = 0.2, color = "black") +
  facet_wrap(~ bioassay) +
  labs(title = 'Histogram of beta_Mn_OEC Values by Treatment', x = 'beta_Mn_OEC', y = 'Frequency') +
  theme_bw()+ 

```
           
```{r}
#total moles of Cd-CA, assuming 1mole Cd  = 1mol Cd-CA
ggplot(filter (nonproteindata, grepl ("T0|T8", full_treatment_name) & grepl ("BA2|", bioassay)), 
       aes(x = pCd_nM, 
        y = pCd_nM *1 ,
        color = tempfe_treatment, 
        shape = tempfe_treatment))+
    geom_point(size = 4, alpha =0.65, stroke =1.5)+
    scale_color_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c('purple', 'black', 'black','darkorange2', "darkorange2")) +
   scale_shape_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c(5, 1, 19,1, 19)) +
   labs (x = "PCd (nM)", y = "Theoretical maximum Cd-CA (nM) " )+
   facet_wrap(~bioassay) +
   theme_bw()+
   theme(axis.title = element_text (size = 14, face = "bold"), 
        axis.text = element_text(size = 12, color = 'black'), 
        legend.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 10, margin = margin(r = 0.5, unit = 'cm')), 
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 14, face = "bold"))

#Convert moles of Cd-CA to ug. Assuming One mole Cd-CA = 68.76 kg
#1 mole Cd-CA = 68760 g  = 6.876 e+10 ug
 #– This is calculated based on Cd-CA AA sequence from Lane et al. 2005 (Cd-CA = 68.76 kDa = 68.76 kg mole -1)
#1 nmole Cd-CA = 6.876 e+1 ug


nonproteindata$theor_CdCA_ug <- nonproteindata$pCd_nM * 68.76

ggplot(filter (nonproteindata, grepl ("T0|2|4|6|8", full_treatment_name) & grepl ("BA2|", bioassay)), 
       aes(x = pCd_nM, 
        y = theor_CdCA_ug ,
        color = tempfe_treatment, 
        shape = tempfe_treatment))+
    geom_point(size = 4, alpha =0.65, stroke =1.5)+
    scale_color_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c('purple', 'black', 'black','darkorange2', "darkorange2")) +
   scale_shape_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c(5, 1, 19,1, 19)) +
   labs (x = "pCd (nM)", y = "Theoretical maximum Cd-CA (ug/L) " )+
   facet_wrap(~bioassay) +
   theme_bw()+
   theme(axis.title = element_text (size = 14, face = "bold"), 
        axis.text = element_text(size = 12, color = 'black'), 
        legend.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 10, margin = margin(r = 0.5, unit = 'cm')), 
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 14, face = "bold"))
```

Is the theoretical maximum reasonable? 
```{r}
ggplot(filter (nonproteindata, grepl ("T0|T8|", full_treatment_name) & grepl ("BA2|", bioassay)), 
       aes(#x = protein_ug_L*3, #this is to adjust for the wrong protein calculation - need to go back and double check all of this
          x = PON_ugL*4.78,
          y = theor_CdCA_ug, 
          color = tempfe_treatment, 
         shape = tempfe_treatment))+
      scale_color_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c('purple', 'black', 'black','darkorange2', "darkorange2")) +
   scale_shape_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c(5, 1, 19,1, 19)) +
       
    geom_point(size = 4, alpha =0.65, stroke =1.5)+
   labs (x = "Total protein estimate (ug/L) PON*4.78", y = "Theoretical maximum Cd-CA (ug/L) " )+
   facet_wrap(~bioassay) +
   theme_bw()+
   theme(axis.title = element_text (size = 14, face = "bold"), 
        axis.text = element_text(size = 12, color = 'black'), 
        legend.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 10, margin = margin(r = 0.5, unit = 'cm')), 
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 14, face = "bold"))  
  # xlim (0,1200) + 
  # ylim (0,1200)


#Fraction of total protein that could be Cd-Ca
ggplot(filter (nonproteindata, grepl ("T0|T8|", full_treatment_name) & grepl ("BA2|", bioassay)), 
       aes(#x = protein_ug_L*3, #this is to adjust for the wrong protein calculation - need to go back and double check all of this
          y = 100*(theor_CdCA_ug/(PON_ugL*4.78)),
           #y=(theor_CdCA_ug/(protein_ug_L)),
          x = bioassay))+
    geom_boxplot() + 
   geom_point(size = 3, alpha =0.65, stroke =1.5, aes (color = tempfe_treatment, 
         shape = tempfe_treatment))+

      scale_color_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c('purple', 'black', 'black','darkorange2', "darkorange2")) +
      scale_shape_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c(5, 1, 19,1, 19)) +
       
   labs (x = "", y = "Maximum contribution of Cd-CA to total protein (%)" )+
  # facet_wrap(~bioassay) +
   theme_bw()+
   theme(axis.title = element_text (size = 14, face = "bold"), 
        axis.text = element_text(size = 12, color = 'black'), 
        legend.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 10, margin = margin(r = 0.5, unit = 'cm')), 
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 14, face = "bold"))  

```

Look at Cd-CA 
```{r}

cdca <- filter(norm_proteomicsdata, grepl("clust_6573", uniq_cluster))
cdca2 <- cdca [c(4, 40:69)]

cdca3 <- melt(cdca2, id.vars=c("uniq_cluster")) %>%
            rename(treatment = variable) %>%
            rename(norm_abundance = value ) %>% 
            separate(treatment, into =  c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"))  
 
cdca3$tempfe_treatment <- paste (cdca3$temperature_treatment, cdca3$iron_treatment, sep = "_")
  
cdca4 <- aggregate(cdca3$norm_abundance, by=list(cdca3$bioassay, cdca3$tempfe_treatment, cdca3$replicate), FUN = sum) %>%
          rename (bioassay = Group.1, 
          tempfe_treatment = Group.2, 
          replicate = Group.3, 
          cdca_norm_abundance = x)

ggplot(cdca4, aes(x = bioassay, 
                  y = 100*(cdca_norm_abundance)))+
    geom_boxplot() + 
    geom_point(size = 3, alpha =0.65, stroke =1.5, aes (color = tempfe_treatment, 
         shape = tempfe_treatment))+

      scale_color_manual(name="",
                      breaks = c("T0_T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c('black', 'black', 'black','darkorange2', "darkorange2")) +
      scale_shape_manual(name="",
                      breaks = c("T0_T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c(5, 1, 19,1, 19)) +
       
   labs (x = "", y = "Contribution of Cd-CA to total matched protein abundance (%)" )+
  # facet_wrap(~bioassay) +
   theme_bw()+
   theme(axis.title = element_text (size = 14, face = "bold"), 
        axis.text = element_text(size = 12, color = 'black'), 
        legend.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 10, margin = margin(r = 0.5, unit = 'cm')), 
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 14, face = "bold"))  

cdca4$tempfe_treatment <- ifelse (cdca4$tempfe_treatment ==  "T0_T0", "T0",
                                                 cdca4$tempfe_treatment)

alldata2 <- merge (nonproteindata, cdca4, by = c("bioassay", "tempfe_treatment", "replicate"), all=TRUE)


#I was thinking just to calculate the proportion of the metal that can be accounted for by protein.

#step 1: Calculate moles of pCd from mass of Cd-CA
#each kg of Cd-CA = 0.0145 moles

#1 ug of cdca = 1.454 e-11 moles

#calculate ug/L of cdca: cdca_norm_abundance*PON_ug_L*4.78)
#calculate moles/L of cdca from ug/L : 1.454e-11 *(cdca_norm_abundance*PON_ug_L*4.78)

#calculate nM pCd from moles/L cd-ca 1e9 *(1.454e-11 *(cdca_norm_abundance*PON_ug_L*4.78))
cdcacalc <-  ggplot(filter (alldata2, grepl ("Fe", full_treatment_name) & grepl ("BA2|", bioassay)), 
       aes( x = bioassay, 
                       

           #y = (cdca_norm_abundance*PON_ugL*4.78)))+
           y =0.6*(1e9 *(1.454e-11 *(cdca_norm_abundance*PON_ugL*4.78)))/pCd_nM*100))+
            #y = (1e10 *(1.454e-11*(cdca_norm_abundance*PON_ug_L*4.78)))))+

   geom_boxplot() + 
    geom_point(size = 3, alpha =0.65, stroke =1.5, aes (color = tempfe_treatment, 
         shape = tempfe_treatment)) +

      scale_color_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c('black', 'black', 'black','black', "black")) +
      scale_shape_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c(5, 1, 19,1, 19)) +
      labs (x = "", y = "PCd accounted for by Cd-CA (%)" )+
  # facet_wrap(~bioassay) +
   theme_bw()+
    theme(axis.text=element_text(face = "bold", size = 12, color = "black"),
          axis.title.y= element_text(size=12, color = "black", face = "bold"), 
          legend.position = "left")
#png("cbiomesCadmiumcluster_20240624.png", width= 6, height= 7, units="in", res=700) 

cdcacalc
#dev.off()
```

Look at Mn cluster
```{r}
#clust_5856
#clust_1906 - psbO this is what we use 
#clust_2763 - psbQ
oec <- filter(norm_proteomicsdata, grepl("clust_1906", uniq_cluster))

oec2 <- oec [c(4, 40:69)]

oec3 <- melt(oec2, id.vars=c("uniq_cluster")) %>%
            rename(treatment = variable) %>%
            rename(norm_abundance = value ) %>% 
            separate(treatment, into =  c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"))  
 
oec3$tempfe_treatment <- paste (oec3$temperature_treatment, oec3$iron_treatment, sep = "_")
  
oec4 <- aggregate(oec3$norm_abundance, by=list(oec3$bioassay, oec3$tempfe_treatment, oec3$replicate, oec3$uniq_cluster), FUN = sum) %>%
          rename (bioassay = Group.1, 
          tempfe_treatment = Group.2, 
          replicate = Group.3, 
          cluster = Group.4, 
          oec_norm_abundance = x)

oec4$normabunancadjusted <- ifelse (oec4$cluster ==  "clust_1906", oec4$oec_norm_abundance*3.0303e-11, 
                                    ifelse(oec4$cluster == "clust_2763", oec4$oec_norm_abundance*6.6667e-11, 
                                           ifelse(oec4$cluster == "clust_5856", oec4$oec_norm_abundance*4.347e-11, 
                                         "else")))
                                    
oec4$tempfe_treatment <- ifelse (oec4$tempfe_treatment ==  "T0_T0", "T0",
                                                 oec4$tempfe_treatment)

alldata_oec <- merge (nonproteindata, oec4, by = c("bioassay", "tempfe_treatment", "replicate"), all=TRUE)


#I was thinking just to calculate the proportion of the metal that can be accounted for by protein.
#Calculate moles of pMn from mass of Mn-cluster

#each kg of OEC-cluster = 0.0303moles

#1 ug of oec = 3.0303 e-11 moles

#calculate ug/L of oec: oec_norm_abundance*PON_ug_L*4.78)
#calculate moles/L of oec from ug/L : 3.0303e-11 *(oec_norm_abundance*PON_ug_L*4.78)

#calculate nM pMn from moles/L Mnclust 1e9 *(3.0303e-11 *(oec_norm_abundance*PON_ug_L*4.78))

ggplot(filter (alldata_oec, grepl ("Fe|", full_treatment_name) & grepl ("BA2|", bioassay)), 
       aes( x = iron_treatment, 
          y = (0.5*0.4*(1e9*4* as.numeric(normabunancadjusted)*PON_ugL*4.78 )/pMn_nM*100)))+
            #y = (1e10 *(1.454e-11*(cdca_norm_abundance*PON_ug_L*4.78)))))+

   geom_boxplot() + 
 stat_summary(fun = "mean",
               geom = "crossbar", 
               width = 0.5,
               colour = "red")+
  geom_point(size = 3, alpha =0.65, stroke =1.5, aes (color = tempfe_treatment, 
         shape = tempfe_treatment))+

      scale_color_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c('black', 'black', 'black','red', "red")) +
      scale_shape_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c(5, 1, 19,1, 19)) +
      labs (x = "", y = "PMn accounted for by Mn-cluster (%)" )+
  facet_wrap(~bioassay) +
   theme_bw()+
 theme(axis.text=element_text(face = "bold", size = 12, color = "black"),
        axis.title.y=element_text(size=12, color = "black", face = "bold"))
 
#png("cbiomesMncluster_20240624.png", width= 6, height= 7, units="in", res=700) 

mnclustcalc    
#dev.off()


ggplot(filter (alldata_oec, grepl ("Fe|", full_treatment_name) & grepl ("BA2|", bioassay)), 
       aes( x = iron_treatment, 
           #y = (cdca_norm_abundance*PON_ug_L*4.78))+
          y = 0.5*(3*0.4*(1e9* as.numeric(normabunancadjusted)*PON_ugL*4.78 ))/(pFe56_nM+pFe57_nM)*100)) +
            #y = (1e10 *(1.454e-11*(cdca_norm_abundance*PON_ug_L*4.78)))))+

   geom_boxplot() + 
    geom_point(size = 3, alpha =0.65, stroke =1.5, aes (color = tempfe_treatment, 
         shape = tempfe_treatment))+

      scale_color_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c('black', 'black', 'black','red', "red")) +
      scale_shape_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c(5, 1, 19,1, 19)) +
      labs (x = "", y = "PFe accounted for by PSII from Mn cluster calcs (%)" )+
   facet_wrap(~bioassay) +
   theme_bw()+
 theme(axis.text=element_text(face = "bold", size = 12, color = "black"),
        axis.title.y=element_text(size=12, color = "black", face = "bold"))

```

PSII metal (2-3 Fe)
```{r}
#https://nph.onlinelibrary.wiley.com/doi/epdf/10.1111/j.1469-8137.1990.tb00505.x
PSII <- filter(norm_proteomicsdata, grepl("^clust_1784$", uniq_cluster))

PSII2 <- PSII [c(4, 40:69)]

PSII3 <- melt(PSII2, id.vars=c("uniq_cluster")) %>%
            rename(treatment = variable) %>%
            rename(norm_abundance = value ) %>% 
            separate(treatment, into =  c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"))  
 
PSII3$tempfe_treatment <- paste (PSII3$temperature_treatment, PSII3$iron_treatment, sep = "_")
  
PSII4 <- aggregate(PSII3$norm_abundance, by=list(PSII3$bioassay, PSII3$tempfe_treatment, PSII3$replicate, PSII3$uniq_cluster), FUN = sum) %>%
          rename (bioassay = Group.1, 
          tempfe_treatment = Group.2, 
          replicate = Group.3, 
          cluster = Group.4, 
          PSII_norm_abundance = x)

#psbv mw = 17,863 Da = 2.9662207e-20g , https://www.uniprot.org/uniprotkb/P48263/entry
#https://www.bioline.com/media/calculator/01_04.html

#1 ug of psbv  = 5.598 e-11 moles



PSII4$normabunancadjusted <- ifelse (PSII4$cluster ==  "clust_1784", PSII4$PSII_norm_abundance*5.598e-11, 
                                    ifelse(PSII4$cluster == "clust_2763", PSII4$PSII_norm_abundance*6.6667e-11, 
                                           ifelse(PSII4$cluster == "clust_5856", PSII4$PSII_norm_abundance*4.347e-11, 
                                         "else")))
                                    
PSII4$tempfe_treatment <- ifelse (PSII4$tempfe_treatment ==  "T0_T0", "T0",
                                                 PSII4$tempfe_treatment)

alldata_PSII <- merge (nonproteindata, PSII4, by = c("bioassay", "tempfe_treatment", "replicate"), all=TRUE)


ggplot(filter (alldata_PSII, grepl ("Fe|", full_treatment_name) & grepl ("BA2|", bioassay)), 
       aes( x = iron_treatment, 
          y = (3*0.4*(1e9* as.numeric(normabunancadjusted)*PON_ugL*4.78 ))/(pFe56_nM+pFe57_nM)*100))+
    geom_boxplot() + 
    geom_point(size = 3, alpha =0.65, stroke =1.5, aes (color = tempfe_treatment, 
         shape = tempfe_treatment))+

      scale_color_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c('black', 'black', 'black','red', "red")) +
      scale_shape_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c(5, 1, 19,1, 19)) +
      labs (x = "", y = "PFe accounted for by psbv (%)" )+
   theme_bw()+
  facet_wrap(~bioassay) +
 theme(axis.text=element_text(face = "bold", size = 12, color = "black"),
        axis.title.y=element_text(size=12, color = "black", face = "bold"))

#https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4989028/
#https://website.whoi.edu/gfd/wp-content/uploads/sites/14/2018/10/Twining__Baines_ANN_REV_MAR_SCI_2013_264044.pdf

#https://www.pnas.org/doi/full/10.1073/pnas.1810886116#t01
#12 Fe atoms per PSI protein complex compared to 3 Fe atoms per PSII complex (Raven et al., 1999)
```

ETC metal (5-7 Fe)
```{r}
  # "clust_613" = "Cyt b6f",
  #     "clust_1497" = "Cyt b6f", 31.3 kd
  #     "clust_1668" = "Cyt b6f",
  #     "clust_1822" = "Cyt b6f",


#https://nph.onlinelibrary.wiley.com/doi/epdf/10.1111/j.1469-8137.1990.tb00505.x
etc <- filter(norm_proteomicsdata, grepl("^clust_1497$", uniq_cluster))

etc2 <- etc [c(4, 40:69)]

etc3 <- melt(etc2, id.vars=c("uniq_cluster")) %>%
            rename(treatment = variable) %>%
            rename(norm_abundance = value ) %>% 
            separate(treatment, into =  c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"))  
 
etc3$tempfe_treatment <- paste (etc3$temperature_treatment, etc3$iron_treatment, sep = "_")
  
etc4 <- aggregate(etc3$norm_abundance, by=list(etc3$bioassay, etc3$tempfe_treatment, etc3$replicate, etc3$uniq_cluster), FUN = sum) %>%
          rename (bioassay = Group.1, 
          tempfe_treatment = Group.2, 
          replicate = Group.3, 
          cluster = Group.4, 
          etc_norm_abundance = x)

#psbv mw = 17,863 Da = 2.9662207e-20g , https://www.uniprot.org/uniprotkb/P48263/entry
#https://www.bioline.com/media/calculator/01_04.html

#1 ug of psbv  = 5.598 e-11 moles



etc4$normabunancadjusted <- ifelse (etc4$cluster ==  "clust_1497", etc4$etc_norm_abundance*3.195e-11, 
                                    ifelse(etc4$cluster == "clust_2763", etc4$etc_norm_abundance*6.6667e-11, 
                                           ifelse(etc4$cluster == "clust_5856", etc4$etc_norm_abundance*4.347e-11, 
                                         "else")))
                                    
etc4$tempfe_treatment <- ifelse (etc4$tempfe_treatment ==  "T0_T0", "T0",
                                                 etc4$tempfe_treatment)

alldata_etc <- merge (nonproteindata, etc4, by = c("bioassay", "tempfe_treatment", "replicate"), all=TRUE)


ggplot(filter (alldata_etc, grepl ("Fe|", full_treatment_name) & grepl ("BA2|", bioassay)), 
       aes( x = iron_treatment, 
          y = (6*0.4*(1e9* as.numeric(normabunancadjusted)*PON_ugL*4.78 ))/(pFe56_nM+pFe57_nM)*100))+
    geom_boxplot() + 
    geom_point(size = 3, alpha =0.65, stroke =1.5, aes (color = tempfe_treatment, 
         shape = tempfe_treatment))+

      scale_color_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c('black', 'black', 'black','red', "red")) +
      scale_shape_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c(5, 1, 19,1, 19)) +
      labs (x = "", y = "PFe accounted for by cytf_ETC (%)" )+
   theme_bw()+
  facet_wrap(~bioassay) +
 theme(axis.text=element_text(face = "bold", size = 12, color = "black"),
        axis.title.y=element_text(size=12, color = "black", face = "bold"))


```

PSI metal (12 Fe)
```{r}
#  "clust_55" = "PSI",
#       "clust_2815" = "PSI",
#       "clust_1564" = "PSI",
#       "clust_2439" = "PSI",
#       "clust_2460" = "PSI", #21,342da
#       "clust_3396" = "PSI"



PSI <- filter(norm_proteomicsdata, grepl("^clust_2460$", uniq_cluster))

PSI2 <- PSI [c(4, 40:69)]

PSI3 <- melt(PSI2, id.vars=c("uniq_cluster")) %>%
            rename(treatment = variable) %>%
            rename(norm_abundance = value ) %>% 
            separate(treatment, into =  c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"))  
 
PSI3$tempfe_treatment <- paste (PSI3$temperature_treatment, PSI3$iron_treatment, sep = "_")
  
PSI4 <- aggregate(PSI3$norm_abundance, by=list(PSI3$bioassay, PSI3$tempfe_treatment, PSI3$replicate, PSI3$uniq_cluster), FUN = sum) %>%
          rename (bioassay = Group.1, 
          tempfe_treatment = Group.2, 
          replicate = Group.3, 
          cluster = Group.4, 
          PSI_norm_abundance = x)

#clust_55 mw = 83,710 Da 

# "clust_2460" = "PSI", #21,342da
#https://www.bioline.com/media/calculator/01_04.html

#1 ug of psbv  = 1.195 e-8 moles



PSI4$normabunancadjusted <- ifelse (PSI4$cluster ==  "clust_55", PSI4$PSI_norm_abundance*1.195e-8, 
                                    ifelse(PSI4$cluster == "clust_2460", PSI4$PSI_norm_abundance*4.686e-11, 
                                           ifelse(PSI4$cluster == "clust_5856", PSI4$PSI_norm_abundance*4.347e-11, 
                                         "else")))
                                    
PSI4$tempfe_treatment <- ifelse (PSI4$tempfe_treatment ==  "T0_T0", "T0",
                                                 PSI4$tempfe_treatment)

alldata_PSI <- merge (nonproteindata, PSI4, by = c("bioassay", "tempfe_treatment", "replicate"), all=TRUE)


ggplot(filter (alldata_PSI, grepl ("Fe|", full_treatment_name) & grepl ("BA2|", bioassay)), 
       aes( x = iron_treatment, 
          y = (12*0.4*(1e9*as.numeric(normabunancadjusted)*PON_ugL*4.78 ))/(pFe56_nM+pFe57_nM)*100))+
    geom_boxplot() + 
    geom_point(size = 3, alpha =0.65, stroke =1.5, aes (color = tempfe_treatment, 
         shape = tempfe_treatment))+

      scale_color_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c('black', 'black', 'black','red', "red")) +
      scale_shape_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c(5, 1, 19,1, 19)) +
      labs (x = "", y = "PFe accounted for by clust_2460 (PSI) (%)" )+
   theme_bw()+
  facet_wrap(~bioassay) +
 theme(axis.text=element_text(face = "bold", size = 12, color = "black"),
        axis.title.y=element_text(size=12, color = "black", face = "bold"))

#https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4989028/
#https://website.whoi.edu/gfd/wp-content/uploads/sites/14/2018/10/Twining__Baines_ANN_REV_MAR_SCI_2013_264044.pdf

#https://www.pnas.org/doi/full/10.1073/pnas.1810886116#t01
#12 Fe atoms per PSI protein complex compared to 3 Fe atoms per PSII complex (Raven et al., 1999)
```


Plastocyanin and copper
```{r}
#https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9502191/#:~:text=It%20is%20estimated%20that%20about,affinity%20Fe%20uptake%20%5B65%5D.
# is estimated that about 50–70% of intracellular metabolic Cu is allocated to plastocyanin in diatoms and green algae [46,63,64]. 

#https://www.frontiersin.org/journals/marine-science/articles/10.3389/fmars.2022.975184/full
# The chlorophyll-derived physiological parameter, rETRmax, provides an estimate of the maximum photosynthetic electron transport rate and is positively correlated with growth rate in the diatom (Kim and Price, 2017). Pooled data from the LCuA1nM, CCuA and ancestral populations showed a significant, positive relationship between growth rate and rETRmax (F(1,4) = 47.1, p=0.002, Figure 4A), suggesting growth rate depended on photosynthetic electron flow and photosynthetic capacity. In T. oceanica, plastocyanin, the most abundant Cu-containing protein, transfers electrons from cytochrome b6f to photosystem I (Peers and Price, 2006) and is critical to maintaining efficient photosynthetic electron flow. Its abundance decreases to about 1 μmol L-1 CV when Cu is limiting, but it makes up about 70% of the total QCu under these conditions, much higher than when Cu is sufficient. Preferential retention of Cu in plastocyanin relative to other Cu-containing proteins underscores its overall importance to cellular metabolism. Kim and Price (2017) showed rETRmax was positively correlated with QCu in Cu-limited T. oceanica, consistent with the function of Cu-containing plastocyanin in photosynthesis. In the present study, rETRmax was negatively correlated with QCu implying that faster electron transport rates were surprisingly associated with lower cellular Cu (Figure 4B). Note that QCu and rETRmax data in Figure 4 were obtained from three different populations (CCuA, LCuA1nM and ancestral) that were adapted to different growth conditions, while Kim and Price (2017) assayed a single population under different environmental conditions. Faster rates of rETRmax and growth were sustained by less Cu in LCuA1nM cells, raising the question how this diatom restored plastocyanin-dependent electron flow in the absence of additional Cu.


#https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0280827
# Plastocyanin (PETE), a replacement protein for cytochrome c6, is mainly found in low-iron adapted diatoms. PETE possesses copper as a cofactor and is therefore regulated by copper availability [54]. Based on the low expression of cytochrome c6 in T. oceanica cultures in high- and low-iron conditions, it has been suggested that PETE replaces cytochrome c6 permanently [14]. This theory was further supported by studies in the North Pacific, where a permanent replacement of cytochrome c6 in low-iron adapted diatoms and a temporary replacement of cytochrome c6 in coastal diatoms was suggested. In contrast to previous findings [14], our study did not show differences in PETE transcript levels between the three treatments. These differences could stem from a lower iron concentration in the previous study [14]. However, the long half-life supports the hypothesis of a permanent replacement of cytochrome c6 by PETE in T. oceanica (Figs 5B and 6D) [13–15].

# https://link.springer.com/article/10.1186/GB-2012-13-7-R66
# While the cytochrome c6 genes CYTC6A and CYTC6B are found to be weakly expressed, the plastocyanin gene PETE shows high expression under high-iron conditions with a characteristic decrease in low-iron conditions as seen from many constitutively expressed chloroplast genes in the course of the chlorosis response. This suggests a constitutive use of plastocyanin (PETE, p175) instead of cytochrome c6 for photosynthetic electron transport and is consistent with prior findings [11, 12].




# 
# 
# https://aslopubs.onlinelibrary.wiley.com/doi/pdf/10.4319/lo.2013.58.2.0613
# 
# One interpretation is that our species does not
# use a multicopper oxidase as part of a high-affinity Fe
# uptake system (i.e., does not need to produce DA to
# acquire Cu under Fe limitation). 





#  "clust_55" = "PSI",
#       "clust_2815" = "PSI",
#       "clust_1564" = "PSI",
#       "clust_2439" = "PSI",
#       "clust_2460" = "PSI", #21,342da
#       "clust_3396" = "PSI"



PSI <- filter(norm_proteomicsdata, grepl("^clust_1820$", uniq_cluster))

PSI2 <- PSI [c(4, 40:69)]

PSI3 <- melt(PSI2, id.vars=c("uniq_cluster")) %>%
            rename(treatment = variable) %>%
            rename(norm_abundance = value ) %>% 
            separate(treatment, into =  c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"))  
 
PSI3$tempfe_treatment <- paste (PSI3$temperature_treatment, PSI3$iron_treatment, sep = "_")
  
PSI4 <- aggregate(PSI3$norm_abundance, by=list(PSI3$bioassay, PSI3$tempfe_treatment, PSI3$replicate, PSI3$uniq_cluster), FUN = sum) %>%
          rename (bioassay = Group.1, 
          tempfe_treatment = Group.2, 
          replicate = Group.3, 
          cluster = Group.4, 
          PSI_norm_abundance = x)

#clust_55 mw = 83,710 Da 

# "clust_2460" = "PSI", #21,342da
#https://www.bioline.com/media/calculator/01_04.html

#1 ug of psbv  = 1.195 e-8 moles



PSI4$normabunancadjusted <- ifelse (PSI4$cluster ==  "clust_1820", PSI4$PSI_norm_abundance*1.0204e-10, 
                                    ifelse(PSI4$cluster == "clust_2460", PSI4$PSI_norm_abundance*4.686e-11, 
                                           ifelse(PSI4$cluster == "clust_5856", PSI4$PSI_norm_abundance*4.347e-11, 
                                         "else")))
                                    
PSI4$tempfe_treatment <- ifelse (PSI4$tempfe_treatment ==  "T0_T0", "T0",
                                                 PSI4$tempfe_treatment)

alldata_PSI <- merge (nonproteindata, PSI4, by = c("bioassay", "tempfe_treatment", "replicate"), all = T) %>% 
                filter(grepl ("0|8", daycount))
               


ggplot(filter (alldata_PSI, grepl ("Fe|", full_treatment_name) & grepl ("BA2|", bioassay)), 
       aes( x = iron_treatment, 
          y = (1*0.4*(1e9*as.numeric(normabunancadjusted)*PON_ugL*4.78 ))/(pCu_nM)*100))+
    geom_boxplot() + 
    geom_point(size = 3, alpha =0.65, stroke =1.5, aes (color = tempfe_treatment, 
         shape = tempfe_treatment))+

      scale_color_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c('black', 'black', 'black','red', "red")) +
      scale_shape_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c(5, 1, 19,1, 19)) +
      labs (x = "", y = "PFe accounted for by clust_2460 (PSI) (%)" )+
   theme_bw()+
  facet_wrap(~bioassay) +
 theme(axis.text=element_text(face = "bold", size = 12, color = "black"),
        axis.title.y=element_text(size=12, color = "black", face = "bold"))



```


3D LCMS figure
```{r}
#https://lgatto.github.io/2020-02-17-RProt-Berlin/raw-ms-data-mzr-and-msnbase.html
#https://www.kaggle.com/code/desertman/plotting-total-ion-chromatogram-of-a-mzml-file/notebook
#BiocManager::install("mzR")

library(mzR)

xfile <- openMSfile("C:/Users/loayj/OneDrive - Woods Hole Oceanographic Institution/Desktop/211216_1037_097_LJ_S27_99.mzML") #this is the directory where I keep all the data
peakTable <- header(xfile) # this gets table of details for each spectra
spectraList <- spectra(xfile) # this gets the spectra values.
# Conver the RT from sec to min
peakTable$retentionTime <- round(peakTable$retentionTime/60,digits = 2)
head(peakTable)







peakTable2 <- peakTable[peakTable$basePeakIntensity >=100000, ]


ggplot(peakTable2, 
       aes( y = retentionTime, 
          x = basePeakMZ, 
          color = basePeakIntensity))+
    scale_color_viridis_c(option= "D")+ #H

  scale_x_continuous(limits = c(600,610))+
   #scale_x_continuous(limits = c(50,60))+

    geom_point()+
  theme_minimal()+
  theme(legend.position = "none")



# Convert ggplot to plotly object and add 3D histogram-like shape
p_plotly <- plotly::plot_ly(peakTable2, 
                            x = ~retentionTime, 
                            y = ~basePeakMZ, 
                            z = ~basePeakIntensity,
                            color = ~basePeakIntensity,
                            colors = viridis::viridis(256),
                            type = "scatter3d", #mesh3d
                            mode = "markers") %>%

plotly::layout(scene = list(xaxis = list(title = 'Retention Time',
                                           range = c(0, 75)),
                              yaxis = list(title = 'Base Peak MZ',
                                           range = c(600, 610)),
                              zaxis = list(title = 'Base Peak Intensity')))


p_plotly




#%>%
  plotly::add_trace(
    type = "histogram2dcontour", 
    x = ~retentionTime, 
    y = ~basePeakMZ, 
    z = ~basePeakIntensity,
    colorscale = "Viridis",
    opacity = 0.6
  ) %>%
  plotly::layout(scene = list(xaxis = list(title = 'Retention Time'),
                              yaxis = list(title = 'Base Peak MZ'),
                              zaxis = list(title = 'Base Peak Intensity')))

# Display the plot
p_plotly








```


#2024-05-29 -- Figure 4 - Sample-normalized major protein groups and N:P  

Plot the rubisco, photosynthesis, ribosomes 
```{r}
samplenormcompartment <- function (compartments, yaxislabel) { 
ggplot(filter (compartment_aggregated_total, grepl (compartments, compartment)),
       aes(x= iron_treatment,
       y= totalproteinnorm_abundance,
       color = temperature_treatment, 
       shape = iron_treatment)) + 
  facet_wrap(~bioassay)+ 
  stat_summary(fun = mean, geom = "point", size = 4, stroke = 1.5, alpha = 0.5)+
  stat_summary(fun.data =  mean_se, geom = "errorbar", show.legend = F, linewidth = 1, width = 0.1)+
  scale_color_manual(name="Temperature",
                     breaks = c("LT","HT"),
                     labels = c("LT", "HT"),
                     values = c('black', 'darkorange2','black')) +
  scale_shape_manual(name = "Iron",
                       breaks = c("Fe","noFe", "T0"), 
                       labels = c("Fe", "noFe", "T0"),
                       values = c(19,1,5)) +
  scale_x_discrete(limits=c("T0","noFe","Fe"), 
                   labels = c("T0", "-Fe", "+Fe"))+
                 # labels = expression(bold("T0"), bold("T8")~"(-Fe)", bold("T8")~"(+Fe)"))+
  xlab(expression ("")) +
  ylab (substitute(yaxislabel))+ 

  theme_bw()+
  theme(axis.text=element_text(face = "bold", size = 12, color = "black"),
        axis.title =element_text(size = 14, color = "black"), 
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 12, face = "bold"),
        legend.position = "none") 
}


compartment_aggregated_total$bioassay <- ifelse (compartment_aggregated_total$bioassay == "BA1", "Offshore", 
                                          ifelse (compartment_aggregated_total$bioassay == "BA2", "Nearshore", 
                                                   "somethingwentwrong"))


photosyn <- samplenormcompartment ("^chloro$", "Fraction of total matched proteins") +  
   ylim (0, 0.22) + 
  ggtitle ("Photosynthetic proteins")

rubisco <- samplenormcompartment ("^rubisco$", "Fraction of total matched proteins") +  
  ggtitle ("RuBisCO") +   ylim (0, 0.22)


ribosome <- samplenormcompartment ("^ribosome$", "Fraction of total matched proteins") +  
  ggtitle ("Ribosomal proteins") +   ylim (0, 0.22)

    #theme (axis.text.x = element_blank())  

ribosome

#png("ps117_figureS4majorproteingrops.png", width=7, height=11, units="in", res=400) 

ggdraw (plot_grid(photosyn, rubisco, ribosome, 
          ncol = 1, 
          labels = c("A", "B", "C"))) + 
  draw_plot (vertical_legend, x = 0.8, y = 0.0, width = 0.25, height = 0.18)  

#dev.off()
    

ribosome
```

#Figure 3 - Taxon-normalized major protein groups 2023-05-26
```{r}
#Ribosomes, photosynthesis, rubisco

#chlorotest <- filter (proteindata_unnorm, grepl ("^ribosome$", compartment))

proteindata_unnorm$compartment <- ifelse (proteindata_unnorm$uniq_cluster  ==  "unassigned_cluster", "ribosome", 
                                  ifelse (proteindata_unnorm$uniq_grpnorm_compartment == "ribosomal", "ribosome",
                                  ifelse (grepl("ribosom|Ribosom", proteindata_unnorm$uniq_cluster_annotation) ,  "ribosome",
                                   
                                  ifelse (proteindata_unnorm$uniq_KO_pathway  ==  "Photosynthesis", "chloro", 
                                  ifelse (grepl("^clust_1820$|^clust_534$", proteindata_unnorm$uniq_cluster) ,  "chloro",
                                  ifelse (grepl("^clust_891$", proteindata_unnorm$uniq_cluster) ,  "nitrogen",

                                          
                                  ifelse (grepl("^clust_129$|^clust_2490$|^clust_1820$|^clust_534$", proteindata_unnorm$uniq_cluster) ,  "rubisco",
                                  ifelse (grepl("Hsp|hsp", proteindata_unnorm$uniq_cluster_annotation) ,  "HSP",
                                  ifelse (grepl("ATP", proteindata_unnorm$uniq_cluster_annotation),  "ATP",
                                              proteindata_unnorm$uniq_grpnorm_compartment))))))))) # proteindata_unnorm$uniq_cluster))) #


compartment <- proteindata_unnorm [c(70, 37, 40:69  )]                           

compartment <- melt(compartment, id.vars=c("compartment", "F")) %>%
            rename(treatment = variable) %>%
            rename(unnnorm_abundance = value) 
            #separate(treatment, into=c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"), sep = "_")

compartment$F <- ifelse (compartment$F ==  "Polar-centric-Mediophyceae", "Centric diatom", 
                 ifelse (compartment$F ==  "Radial-centric-basal-Coscinodiscophyceae", "Centric diatom", 
               #  ifelse (compartment$F ==  "Bacillariophyta_X", "Bacillariophyta", 
                 ifelse (compartment$F ==  "Raphid-pennate", "Pennate diatom", 
                 ifelse (compartment$F ==  "Araphid-pennate", "Pennate diatom", 
                 ifelse (compartment$F ==  "Phaeocystaceae", "Phaeocystaceae",
                                          compartment$F)))))

#normalize to group 
compartment_aggregated_group  <- aggregate(compartment$unnnorm_abundance, 
                                  by=list(Category=compartment$compartment, compartment$treatment, compartment$F), FUN=sum) %>%
                        rename(unnorm_abundance = x) %>%
                        rename(compartment = Category) %>%
                        rename(treatment = Group.2) %>% 
                        rename(taxon_F = Group.3) 

compartment_aggregated_group  <- compartment_aggregated_group  %>% 
                          group_by(treatment,taxon_F) %>% 
                          mutate(groupnorm_abundance = unnorm_abundance / sum(unnorm_abundance)) %>% 
                          ungroup()

compartment_aggregated_group <- compartment_aggregated_group %>% 
  separate(treatment, into=c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"), sep = "_") %>%
  mutate(tempfe_treatment = paste(temperature_treatment, iron_treatment, sep = "_")) %>%
  mutate(tempfe_treatment = str_replace(tempfe_treatment, "T0_T0", "T0"))


grpnorm_compartment <- function (compartments, taxon, yaxislabel) { 
ggplot(filter (compartment_aggregated_group, grepl (compartments, compartment) & grepl (taxon, taxon_F)),
       aes(x= iron_treatment,
       y= groupnorm_abundance,
       color = temperature_treatment, 
       shape = iron_treatment)) + 
  facet_wrap(~bioassay)+
  stat_summary(fun = mean, geom = "point", size = 4, stroke = 1.5, alpha = 0.5)+
  stat_summary(fun.data =  mean_se, geom = "errorbar", show.legend = F, linewidth = 1, width = 0.1)+
  scale_color_manual(name="Temperature",
                     breaks = c("LT","HT"),
                     labels = c("LT", "HT"),
                     values = c('black', 'darkorange2','black')) +
  scale_shape_manual(name = "Iron",
                       breaks = c("Fe","noFe", "T0"), 
                       labels = c("Fe", "noFe", "T0"),
                       values = c(19,1,5)) +
  scale_x_discrete(limits=c("T0","noFe","Fe"), 
                   labels = c("T0", "-Fe", "+Fe"))+
                 # labels = expression(bold("T0"), bold("T8")~"(-Fe)", bold("T8")~"(+Fe)"))+
  xlab(expression ("")) +
  ylab (substitute(yaxislabel))+ 

  theme_bw()+
  theme(axis.text=element_text(face = "bold", size = 13, color = "black"),
        axis.title.x=element_text(size = 14,face = "bold", color = "black"), 
        axis.title.y=element_text(size=14.5, color = "black", face = "bold"), 
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 12, face = "bold"),
        legend.position = "none") 
}




#pennates
pennatephoto <- grpnorm_compartment ("^chloro$", "Pennate di", atop ("Taxon-normalized", paste ("Photosynthetic protein fraction"))) +
  ylim(0,0.33) + 
  theme (axis.text.x = element_blank())+
  ggtitle("Pennate diatoms")

pennaterubisco <- grpnorm_compartment ("^rubisco$", "Pennate di", atop ("Taxon-normalized", paste ("RuBisCO protein fraction"))) +
  ylim(0,0.07) + 
  theme (axis.text.x = element_blank()) 
 
pennateribosome <- grpnorm_compartment ("^ribosome$", "Pennate di", atop ("Taxon-normalized", paste ("Ribosomal protein fraction"))) +
  ylim(0,0.23) 


#centric
centricphoto <- grpnorm_compartment ("^chloro$", "Centric di", " ") + 
  ylim(0,0.33) + 
  theme (axis.text.x = element_blank()) + 
   ggtitle("Centric diatoms")

centricrubisco <- grpnorm_compartment ("^rubisco$", "Centric di", "") + 
  ylim(0,0.35) + 
  theme (axis.text.x = element_blank()) 
   
centricribosome <- grpnorm_compartment ("^ribosome$", "Centric di", "") +
  ylim(0,0.23)  
  
  

#phaeocystis
phaeochloro <- grpnorm_compartment ("^chloro$", "Phaeocys", "") +
  ylim(0, 0.33) + 
  theme (axis.text.x = element_blank()) + 
  ggtitle("Phaeocystis")


phaeorubisco <- grpnorm_compartment ("^rubisco$", "Phaeocys", "")  +
  theme (axis.text.x = element_blank()) + 
  ylim(0,NA)
    
phaeoribosome <- grpnorm_compartment ("^ribosome$", "Phaeocys", "") +
  ylim(0,0.23) 

```


```{r}
ribo <- filter (compartment_aggregated_group, grepl ("ribosome", compartment))
ribo <- ribo[c(2:7, 9:10)]


riboprotein <- merge (allgroups_aggregated, ribo, by = c("tempfe_treatment", "replicate", "taxon_F"))


ggplot ((filter(riboprotein, grepl ("Centric|Pennate|Phaeocystis", taxon_F))), aes (
         x = groupnorm_abundance, 
         y = unnorm_abundance, 
         shape = iron_treatment.x, 
         color = taxon_F))+
  geom_point() + 
  scale_y_log10() + 
    facet_wrap(~iron_treatment.x+temperature_treatment.x)

```


testing out this figure of normalizing chloro to ribosomes, rubisco to chloro etc. 
```{r}
test1 <- filter (compartment_aggregated_group, grepl ("Centric|Pennate|Phaeocystis", taxon_F)) 

test2 <- dcast(test1, taxon_F~compartment)
```



cbiomes3 
```{r}
pennatephoto <- grpnorm_compartment ("^chloro$", "Pennate di", atop ("Taxon-normalized", paste ("Photosynthetic fraction"))) +
  ylim(0,0.33) + 
  theme (axis.text.x = element_blank())+
  ggtitle("Pennate diatoms")

pennaterubisco <- grpnorm_compartment ("^rubisco$", "Pennate di", atop ("Taxon-normalized", paste ("RuBisCO fraction"))) +
  ylim(0,0.07) + 
  theme (axis.text.x = element_blank()) 
 
pennateribosome <- grpnorm_compartment ("^ribosome$", "Pennate di", atop ("Taxon-normalized", paste ("Ribosomal fraction"))) +
  ylim(0,0.23) 



#png("cbiomes3.png", width=13, height=6.3, units="in", res=900) 

plot_grid (pennatephoto, NULL, centricphoto , NULL, phaeochloro,
           pennaterubisco, NULL, centricrubisco, NULL, phaeorubisco, 
           pennateribosome, NULL, centricribosome, NULL, phaeoribosome,
           rel_widths = c(1,0.02,1,0.02,1),
           ncol=5, 
           align = "vh", axis ="b")



###
pennateribo <- grpnorm_compartment ("^ribosome$", "Pennate di", atop ("Taxon-normalized ribosomal fraction")) +
  ylim(0,0.33) + 
  ggtitle("Pennate diatoms")

centricribo <- grpnorm_compartment ("^ribosome$", "Centric di", atop ("Taxon-normalized ribosomal fraction")) +
  ylim(0,0.33) + 
  ggtitle("Centric diatoms")


phaeoribo <- grpnorm_compartment ("^ribosome$", "Phaeocy", atop ("Taxon-normalized ribosomal fraction")) +
  ylim(0,0.33) + 
  ggtitle("Phaeocystis")


#png("cbiomestalk.png", width=14, height=4, units="in", res=900) 

ggdraw (plot_grid (pennateribo, NULL, centricribo , NULL, phaeoribo,
           rel_widths = c(1,0.02,1,0.02,1),
           ncol=5, 
           align = "vh") + 
          
draw_plot (vertical_legend, x = 0.838, y = 0.13, width = 0.25, height = 0.18)) 


#dev.off()



```

diatom and Phaeocystis protein abundance - 2023-06-19. unnoramlized
```{r}
allgroups <- proteindata_unnorm [c(32,34,35, 37, 40:69)]

allgroups <- melt(allgroups, id.vars=c("domain", "C","class", "F")) %>%
             rename(treatment = variable) %>%
             rename(unnorm_abundance = value )
sum(allgroups$unnorm_abundance) #this should equal 1.893E13 (total abundance of all peptides for all samples)

allgroups$taxon <- ifelse (grepl("pennate", allgroups$F),  "Pennate diatom", 
                    ifelse (grepl("centric", allgroups$F), "Centric diatom", 
                    ifelse (grepl("Phaeocystaceae", allgroups$F), "Phaeocystis",
                    allgroups$F))) 

allgroups_aggregated <- aggregate(allgroups$unnorm_abundance, by=list(allgroups$treatment, allgroups$taxon), FUN=sum) %>%
                  rename(unnorm_abundance = x) %>%
                  rename(treatment = Group.1) %>%
                  rename(taxon_F = Group.2) %>% 
                  separate(treatment, into =  c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"))  

allgroups_aggregated <- allgroups_aggregated %>%
  mutate(tempfe_treatment = paste(temperature_treatment, iron_treatment, sep = "_")) %>%
  mutate(tempfe_treatment = str_replace(tempfe_treatment, "T0_T0", "T0"))


#allgroupsstats <- filter(allgroups_aggregated, grepl("Pennate|Centric|Haptophyta|Bacillariophyta|Phaeocys", taxon))
#write.csv(allgroupsstats, "pennate-centric-hapt_stat.csv", row.names = F)


protein_taxgroups <- function (taxonomicgroup, yaxislabel) { 
ggplot(filter (allgroups_aggregated, grepl (taxonomicgroup, taxon_F)), 
   aes(x= iron_treatment,
       y= unnorm_abundance, 
       color = temperature_treatment, 
       shape = iron_treatment))+
  facet_wrap(~bioassay)+           
  stat_summary(fun = mean, geom = "point", size = 3.5, stroke = 1.2, alpha = 0.5)+
  stat_summary(fun.data =  mean_se, geom = "errorbar", show.legend = F, linewidth = 1, width = 0.1)+
  scale_x_discrete(limits=c("T0", "noFe","Fe"), 
                   labels = c("T0", "-Fe", "+Fe"))+
      scale_color_manual(name="Temperature",
                     breaks = c("LT","HT"),
                     labels = c("LT", "HT"),
                     values = c('black', 'darkorange2','black')) +
  scale_shape_manual(name = "Iron",
                       breaks = c("Fe","noFe", "T0"),
                       labels = c("Fe", "noFe", "T0"),
                       values = c(19,1,5)) +
  xlab(expression ("")) +
  ylab (substitute(yaxislabel))+ 
  theme_bw()+
   theme(axis.text=element_text(face = "bold", size = 12, color = "black"),
        axis.title.x=element_text(size = 14,face = "bold", color = "black"), 
        axis.title.y=element_text(size=16, color = "black"), 
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 12, face = "bold"),
        legend.position = "none")  +
    ylim (0,NA)
}




pennate <- protein_taxgroups("Pennate", "Total pennate diatom protein")
centric <- protein_taxgroups("Centric", "Total centric diatom protein") 
phaeo <- protein_taxgroups("Phaeocystis", "Total Phaeocystis protein") 


unnormprotein <- plot_grid(pennate, NULL, centric, NULL,  phaeo, 
         rel_widths = c(1,0.02,1,0.02,1), 
         nrow = 1, 
         labels = c ("A","", "B","", "C"))






###
pennateribo <- grpnorm_compartment ("^ribosome$", "Pennate di", atop ("Pennate diatom ribosomal fraction")) +
  ylim(0,0.33) 

centricribo <- grpnorm_compartment ("^ribosome$", "Centric di", atop ("Centric diatom ribosomal fraction")) +
  ylim(0,0.33) 


phaeoribo <- grpnorm_compartment ("^ribosome$", "Phaeocy", atop ("Phaeocystis ribosomal fraction")) +
  ylim(0,0.33) 


#png("cbiomestalk.png", width=14, height=4, units="in", res=900) 

x <- ggdraw (plot_grid (pennateribo, NULL, centricribo , NULL, phaeoribo,
           rel_widths = c(1,0.02,1,0.02,1),
           ncol=5, 
           align = "vh", 
           labels = c("D", "", "E", "", "F")) + 
          
draw_plot (vertical_legend, x = 0.838, y = 0.15, width = 0.25, height = 0.18)) 


#dev.off()
png("ps117_proteinribosome_20230625.png", width=14, height=8.5, units="in", res=900) 

plot_grid (unnormprotein, NULL, x, 
           align = "vh", 
           rel_heights = c(1,0.2,1),
           
           nrow = 3)
dev.off()
```

Figure 
```{r}
a <- plot_grid (pennatephoto, NULL, centricphoto , NULL, phaeochloro,
           pennaterubisco, NULL, centricrubisco, NULL, phaeorubisco, 
           pennateribosome, NULL, centricribosome, NULL, phaeoribosome,
           rel_widths = c(1,0.02,1,0.02,1),
           ncol=5, 
           
           labels = c("A", "", "D", "", "G", 
                      "B", "",  "E", "", "H",
                      "C", "",  "F", "", "I"),
          
           align = "vh", axis ="b")


#png("ps117_taxonphotoriborubisco_20230530.png", width=15, height=10, units="in", res=900) 

ggdraw (a) + draw_plot (vertical_legend, x = 0.838, y = 0, width = 0.25, height = 0.18)  

# plot_grid (a, vertical_legend, ncol =  2, 
#            rel_widths = c(4, 0.3), 
#            rel_heights = c(1,0.5))

#dev.off()
```








```{r}
cbiomes4
#png("cbiomes4.png", width= 5, height=13, units="in", res=900) 
plot_grid(Fe57PON, 
                    MnPON, 
                    CdPON, 
                    ZnPON, ncol = 1)
#dev.off()

#png("cbiomes4b.png", width= 4, height=13, units="in", res=900) 
plot_grid (mnclustcalc, cdcacalc, 
           ncol = 1, 
           align = "hv", axis = "lr",
           labels = c("", "") )



png("cbiomes4c.png", width= 8, height=13, units="in", res=700) 
plot_grid(FeDE_BA1, FeDE_BA2,
          ncol = 2,
          rel_widths = c(1.2,0.75))

dev.off()




```


Figure 5 _ 20230611
```{r}
#png("ps117_fig5_DE_20230611.png", width=10, height=6, units="in", res=900) 

#dev.off()
metals <- plot_grid(NULL, Fe57PON, 
                    NULL, MnPON, 
                    NULL, CdPON, 
                    NULL, CuPON,
                    NULL, ZnPON,
  
         labels = c("A", "" , 
                    "B", "", 
                    "C", "", 
                    "D", "",
                    "E", ""), 
          ncol = 2, 
         rel_widths = c(0.1,1),
          align = "hv")

de <- plot_grid(FeDE_BA1, FeDE_BA2,
          ncol = 2,
          labels = c("F", "G"),
          rel_widths = c(1,0.75))

calcs <- plot_grid (mnclustcalc, cdcacalc, 
           ncol = 1, 
           align = "hv", axis = "lr",
           labels = c("H", "I") )

#png("ps117_figure5_20230713.png", width=16, height=10, units="in", res=400)

ggdraw(plot_grid (metals, de, NULL, calcs,  ncol =4, rel_widths = c(1,3, 0.1,1), align = "hv") + 
            draw_plot (vertical_legend , x = 0.2, y= 0.5, width = 0.05, height = 0.05)) 

#dev.off()


```


Figure 5
```{r}
#png("ps117_fig5_DE_20230605.png", width=10, height=6, units="in", res=900) 
plot_grid(FeDE_BA1, FeDE_BA2,
          ncol = 2,
          rel_widths = c(1,0.75))

dev.off()


#ForErin
metals <- plot_grid(ZnPON,
          CdPON,
          CuPON,
          MnPON, 
          labels = c("A", "B", "C", "D"), 
          ncol = 1, 
          align = "hv")

calcs <- plot_grid (mnclustcalc, cdcacalc, 
           ncol = 1, 
           align = "hv", axis = "lr",
           labels = c("E", "F") )


#png("ps117_metals_20230526.png", width=12, height=10, units="in", res=900) 

plot_grid (metals, NULL, calcs,
           ncol = 3, 
           align = "v", axis = "b",
           rel_widths = c(1,0.1,1))

```

Figure 6 
```{r}
metals <- plot_grid(Fe57PON, NULL, totalFe57,
          CdPON, NULL,   MnPON, 
          CuPON, NULL, ZnPON,
          ncol = 3, 
          rel_widths = c(1,0.1,1),
          labels = c("A", "",  "B",
                     "C", "", "D", 
                     "E", "" ,  "F"), 
       
          align = "hv")


metalcalcs <- plot_grid (cdcacalc, mnclustcalc,
           ncol = 1, 
          labels = c("G", "H") )


#png("ps117_fig6_metals_20230605.png", width=13, height=8, units="in", res=900) 

plot_grid (metals, NULL, metalcalcs, 
           ncol = 3, 
           rel_widths = c(1,0.1,0.5))
 
#dev.off()


```


Misc. 
RA plots - I don't think these will make it into the paper 
```{r}
# make list of labels of interest for RA plot  
label_list <- c("ISIP-1", "ISIP-2A", "ISIP-3", "Fld", "Pcyn", "LHC", "Ukn", "DoxX", "Rhodopsin", "Fe-Mn SOD", "Cu-Zn SOD",  "PSI PsaA/PsaB", "PSII RC, D1", "PSII CP47", "Cyt f", "Cyt b6", "Fasciclin", "Rhodopsin", "OEC_PSBO", "OEC_PSBQ", "CA_clust_655", "CA_clust_6573", "Fructose bisphosphatase", "NO2/SO3 reductase" )


#Step 3: combine FC and abundance 
taxon_F <- filter(alldata, grepl ("Phaeocys|diatom", taxon_F))
taxon_F$newabund <- abs(log2(taxon_F$BA1_grpnorm))

foldchangevsabundance <- function(title, xaxis, yaxis, FDR) {
       ggplot(taxon_F, aes (x = {{xaxis}},y = {{yaxis}}, color = (taxon_F), 
       label = ifelse({{FDR}} <0.05 & cluster_annotation %in% label_list, cluster_annotation, "")))+
    geom_point(size = 2, aes(alpha = {{FDR}} < 0.05)) +
    geom_label_repel (size = 3, max.overlaps = Inf, fontface = "bold", show.legend = F, box.padding   = 1, 
                 point.padding = 0.01, force = 10) +
    scale_alpha_manual(values = c(0.1,0.7), guide = 'none') +
    scale_colour_manual(values = c('black','orange4', "seagreen"), name = "") +
   # scale_x_log10()+
    xlab (expression ("Group-normalized protein abundance")) +
    ylab (expression(Fe-Log[2]~fold~change))+ 
    labs(title = title)+
    theme_bw()+
    theme(axis.title=element_text(size = 14,face = "bold", color = "black"), #axes formatting
          axis.text=element_text(face = "bold", size = 12, color = "black", vjust = 0.5, hjust = 1)) + 
    theme(legend.text = element_text(size = 11, face = "bold"))  
}

foldchangevsabundance ("BA1", BA1_grpnorm, BA1_Fe_logFC, BA1_Fe_FDR) + theme (legend.position = "none",   legend.box.margin=margin(-10,-10,-10,-10), legend.background = element_blank()) + xlim(0, 0.03)

foldchangevsabundance ("BA1", BA1_grpnorm, BA1_Fe_logFC, BA1_Fe_FDR) + theme (legend.position = c(0.01,1), legend.justification = c(0.01,1),   legend.box.margin=margin(-10,-10,-10,-10), legend.background = element_blank())


BA2 <- foldchangevsabundance ("BA2", BA2_grpnorm, BA2_Fe_logFC, BA2_Fe_FDR ) + theme (legend.position = "none")

#png("ps117_figure3_highres_20230328.png", width=10, height=11, units="in", res=900) 

plot_grid(BA1, BA2,
          ncol = 1)

#dev.off()

```

Cd-CA tree and whatnot
```{r}
#Consider not splitting the data into bioassays, the Cd and CA trends become strong. 
T0_T8_FenoFe (PCd_nM/(PON_uM),PCd:PON~~(mM:M))
T0_T8_FenoFe (PCd_nM/(POP_nM),PCd:POP~~(M:M))

cluster <- proteindata_unnorm [c(4, 37, 40:69)] #cluster, taxon_F, all reps. 

cluster$F <- ifelse (cluster$F ==  "Polar-centric-Mediophyceae", "Centric diatom", 
                 ifelse (cluster$F ==  "Radial-centric-basal-Coscinodiscophyceae", "Centric diatom", 
                 ifelse (cluster$F ==  "Bacillariophyta_X" & cluster$uniq_cluster == "clust_6573", "Centric diatom", 
                 ifelse (cluster$F ==  "Raphid-pennate", "Pennate diatom", 
                 ifelse (cluster$F ==  "Araphid-pennate", "Pennate diatom", 
                 ifelse (cluster$F ==  "Phaeocystaceae", "Phaeocystaceae",
                                          cluster$F))))))

cluster <- melt(cluster, id.vars=c("uniq_cluster", "F")) %>%
            rename(treatment = variable) %>%
            rename(unnnorm_abundance = value) 


#normalize to group 
cluster_aggregated_group  <- aggregate(cluster$unnnorm_abundance, 
                                  by=list(Category=cluster$uniq_cluster, cluster$treatment, cluster$F), FUN=sum) %>%
                        rename(unnorm_abundance = x) %>%
                        rename(cluster = Category) %>%
                        rename(treatment = Group.2) %>% 
                        rename(taxon_F = Group.3) 

cluster_aggregated_group  <- cluster_aggregated_group  %>% 
                          group_by(treatment,taxon_F) %>% 
                          mutate(groupnorm_abundance = unnorm_abundance / sum(unnorm_abundance)) %>% 
                          ungroup() %>% 
  separate(treatment, into=c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"), sep = "_") %>%
  mutate(tempfe_treatment = paste(temperature_treatment, iron_treatment, sep = "_")) %>%
  mutate(tempfe_treatment = str_replace(tempfe_treatment, "T0_T0", "T0"))


#normalize to total 
cluster_aggregated_total  <- aggregate(cluster$unnnorm_abundance, 
                                  by=list(Category=cluster$uniq_cluster, cluster$treatment), FUN=sum) %>%
                        rename(unnorm_abundance = x) %>%
                        rename(cluster = Category) %>%
                        rename(treatment = Group.2) 

cluster_aggregated_total  <- cluster_aggregated_total  %>% 
                          group_by(treatment) %>% 
                          mutate(groupnorm_abundance = unnorm_abundance / sum(unnorm_abundance)) %>% 
                          ungroup() %>% 
  separate(treatment, into=c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate"), sep = "_") %>%
  mutate(tempfe_treatment = paste(temperature_treatment, iron_treatment, sep = "_")) %>%
  mutate(tempfe_treatment = str_replace(tempfe_treatment, "T0_T0", "T0"))


ggplot(filter (cluster_aggregated_total, grepl ("clust_6573", cluster)), 
                 aes(x = iron_treatment,
                 y = groupnorm_abundance,
                 color = temperature_treatment, 
                 shape = iron_treatment)) + 
  geom_point()+
  stat_summary(fun = mean, geom = "point", size = 4, stroke = 1.5, alpha = 0.5)+
  stat_summary(fun.data =  mean_se, geom = "errorbar", show.legend = F, linewidth = 1, width = 0.1)+
  scale_color_manual(name = "Temperature",
                     breaks = c("LT","HT"),
                     labels = c("LT", "HT"),
                     values = c('black', 'darkorange2','black')) +
  scale_shape_manual(name = "Iron",
                       breaks = c("Fe","noFe", "T0"), 
                       labels = c("Fe", "noFe", "T0"),
                       values = c(19,1,5)) +
  scale_x_discrete(limits = c("T0","noFe","Fe"), 
                   labels = c("T0", "-Fe", "+Fe")) + 
  xlab(expression ("")) +
  ylab (substitute(yaxislabel))+ 
  facet_wrap(~bioassay)+
  theme_bw()+
  theme(axis.text=element_text(face = "bold", size = 12, color = "black"),
        axis.title.x=element_text(size = 14,face = "bold", color = "black"), 
        axis.title.y=element_text(size=16, color = "black"), 
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 12, face = "bold"),
        legend.position = "none") 

CA_centric <- protein_clusters ("centric|Bacill", "CA_clust_6573", protein = c("CA_clust_6573"), shape = c(19,15)) + ggtitle ("Centric diatoms - clust_6573") + theme(legend.position = "none") 

CA_pennate <- protein_clusters ("pennate", "CA_clust_655", protein = c("CA_clust_655"), shape = c(19,15)) + ggtitle ("Pennate diatoms - clust_655") + theme(legend.position = "none") 

CA_phaeo <- protein_clusters ("Phaeo", "CA_clust_655", protein = c("CA_clust_655"), shape = c(19,15)) + ggtitle ("Phaeocystis - clust_655") + theme(legend.position = "none")

#protein_clusters ("pennate", "Cation Efflux protein", protein = c("Cation Efflux protein"), shape = c(19,15)) look into efflux proteins
#protein_clusters_hm ("centric|pennate", "Carbonic|CA", "BA1")
#protein_clusters_hm ("centric|pennate", "Carbonic|CA", "BA2")
#Maybe look atr SLC4 ? 
```

plot
```{r}
plot_grid( metals_Cd,  CA_centric,  
           CA_pennate, CA_phaeo, 
           nrow = 4, 
           align = 'hv')
#and heatmap
```

Stats
```{r}
#Fe effect is low temp Fe vs no Fe
#temp effect is low temp no Fe, vs high temp no Fe
#tempxFe effect is is Fe at high temp vs Fe at low temp


#total_cells_mL
#total_chla_ug_L
#POC_uM, NO3, PO4

#fract_diatoms
#fract_haptophytes
BA1 <- filter (nonproteindata, grepl ("BA1", bioassay) & grepl ("6", daycount))
x <- aov (as.numeric(Si_uM)~as.factor(iron_treatment)*as.factor(temperature_treatment), data = BA1)
TukeyHSD(x)

summary (x)


# centric
# pennate
# phaeo

pencentpha <- read.csv ("pennate-centric-hapt_stat.csv") #this data can also be used for the main diatom and hapto trends. 
pencentpha2 <- filter (pencentpha, grepl ("BA2", bioassay) & grepl ("T8", timepoint))
z <- aov (as.numeric(centric)~as.factor(iron_treatment)*as.factor(temperature_treatment), data = pencentpha2)
summary (z)
TukeyHSD(z)


#diatomfraction
#haptophytefraction
# diatomhaptoprotein <- read.csv ("diatomhapto_protein_stats.csv")
# diatomhaptoprotein2 <- filter (diatomhaptoprotein, grepl ("BA2", bioassay) & grepl ("T8", timepoint))
# z <- aov (as.numeric(haptophytefraction)~as.factor(iron_treatment)*as.factor(temperature_treatment), data = diatomhaptoprotein2)
# summary (z)
# TukeyHSD(z)


#si_no3_drawdown <- read.csv("LJ_ps117_si_consumption.csv", header = T)
#si_no3_drawdown2 <- filter (si_no3_drawdown, grepl ("BA2", bioassay) & grepl ("8", daycount))
#y <- aov (as.numeric(SiNO3ratio)~as.factor(iron_treatment)*as.factor(temperature_treatment), data = si_no3_drawdown2)
#summary (y)
#TukeyHSD(y)
```




#Supplemental figures
Dissolved nutrients
```{r}
dsi <- x_vs_time (Si_uM, Si(OH)[4]~~(mu*M)) + ylim(0,NA)

dmn<- x_vs_time (DMn_nM, DMn~~(n*M)) + ylim(0,NA)
dfe57 <- x_vs_time (DFe57_nM,D^{57}*"Fe" ~~ (nM)) + ylim(0,NA)
dco <- x_vs_time (DCo_nM, DCo~~(n*M)) + ylim(0,NA)
dni <- x_vs_time (DNi_nM, DNi~~(n*M)) + ylim(0,NA)
dcu <- x_vs_time (DCu_nM, DCu~~(n*M)) + ylim(0,NA)
dzn <- x_vs_time (DZn_nM, DZn~~(n*M)) + ylim(0,NA)
dcd <- x_vs_time (DCd_nM, DCd~~(n*M)) + ylim(0,NA)

dfe56 <- ggplot(filter (nonproteindata_long, !grepl ("0.35|0.14", value) & grepl("DFe56", variable) & grepl ("T0|", full_treatment_name)), 
       aes (x= daycount,
       y = value, 
       color = temperature_treatment, 
       shape = iron_treatment))+
  
      stat_summary(fun = mean, geom = "point", show.legend = F, size = 4, stroke = 1.5, alpha = 0.5)+
      stat_summary(fun = mean, geom = "line", linewidth = 0.9, aes(linetype = iron_treatment))+
      stat_summary(fun.data =  mean_sd, geom = "errorbar", size = 1, width = 0.1) +
      
  scale_color_manual(name="Temperature",
                     breaks = c("LT","HT"),
                     labels = c("LT", "HT"),
                     values = c('black', "darkorange", "blue")) +
  scale_linetype_manual(name = "Iron",
                        breaks = c("Fe","noFe"),   #this removes the T0 from the legend
                        labels = c("Fe", "noFe"),
                        values = c("solid", "dotted", "solid")) +
  scale_shape_manual (name = "Iron",
                       breaks = c("noFe","Fe", "T0"), 
                       labels = c("noFe", "Fe", "T0"),
                       values = c(1, 19, 5))+
  scale_x_continuous( labels = c("T0","T2","T4", "T6", "T8"))+
  
   labs(x = expression (""), y = D^{56}*"Fe" ~~ (nM)) + 
       facet_wrap(~bioassay) + 
       theme_bw()+
  theme_bw()+
        theme(axis.text=element_text(face = "bold", size = 12, color = "black"),
        axis.title.x=element_text(size = 14,face = "bold", color = "black"), 
        axis.title.y=element_text(size=16, color = "black"), 
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 12, face = "bold"),
        legend.position = "none") +ylim(0,NA)


#png("ps117_figureS_dnutrients_highres_20230617.png", width=18, height=9, units="in", res=900) 
ggdraw (plot_grid(dsi, NULL,  dni, NULL,  dzn, 
          dcu, NULL,  dfe56,NULL,  dcd, 
          dfe57,NULL,  dmn, NULL,  dco, 
          rel_widths = c(1,0.04,1,0.04),
          ncol = 5, align = "vh", 
          
       labels = c("A", "", "B", "", "C", 
                  "D", "", "E", "", "F", 
                  "G", "", "H", "", "I")) + 
   draw_plot (vertical_legend , x = 0.94, y= 0.08, width = 0.05, height = 0.05))
#dev.off()
```

fvfm_mn
```{r}
fvfmmn <- read.csv ("fvfmmn.csv")


#png("ps117_fvfmmMn_20230617.png", width= 8, height=5, units="in", res=900) 

ggplot(fvfmmn, aes (x = treatment, y = fvfm))+
  geom_point (size = 5)+
                
  xlab("") +
  ylab (expression (F[v]/F[m])) +
  theme_bw()+
    
  theme(axis.title =element_text(size = 18,face = "bold", color = "black"), 
        axis.text = element_text(size = 13, color = "black"))
#dev.off()
```

ctd data 
```{r}
library (oce)

ba1 <- read.oce("station 17 Cast 1 BioAssay #1_rerun.cnv")
ba1_v2 <- subset (ba1, depth > 4)
summary (ba1)

ba1_v3 = ba1_v2%>%ctdTrim(method = "downcast")%>%ctdDecimate(p = 1)


ba1_v3 = ba1_v3@data%>%
  as_data_frame()

colnames(ba1_v3) [13] <- "par"

ba1temp <- ggplot(data = ba1_v3%>%na.omit(), 
       aes(x = temperature, y = depth))+
  geom_path( col = "black", size = 0.8)+
  scale_y_reverse()+
  scale_x_continuous(position = "top")+
  theme_bw()+
  ggtitle("Bioassay 1")  +
  theme(axis.text = element_text(size = 12, colour = 1),
        axis.title = element_text(size = 14, colour = 1), 
        plot.title = element_text(size = 14, face = "bold"))+
  labs(x = expression(~Temperature~(degree*C)), y = "Depth (m)") + 
  geom_hline (aes(yintercept = 20.3), col = "blue", lty = 2) 


ba1par <- ggplot(data = ba1_v3%>%na.omit(), 
       aes(x = par, y = depth))+
  geom_path( col = "black", size = 0.8)+
  scale_y_reverse()+
  scale_x_continuous(position = "top")+
  theme_bw()+
  ggtitle("")  +
  theme(axis.text = element_text(size = 12, colour = 1),
        axis.title = element_text(size = 14, colour = 1), 
        plot.title = element_text(size = 14, face = "bold"))+
  labs(x = expression(PAR~(mu*M~m^-2~S^-1)), y = "Depth (m)") + 
  geom_hline (aes(yintercept = 20.3), col = "blue", lty = 2)

ba1oxygen <- ggplot(data = ba1_v3%>%na.omit(), 
       aes(x = oxygen, y = depth))+
  geom_path( col = "black", size = 0.8)+
  scale_y_reverse()+
  scale_x_continuous(position = "top")+
  theme_bw()+
  ggtitle("")  +
  theme(axis.text = element_text(size = 12, colour = 1),
        axis.title = element_text(size = 14, colour = 1), 
        plot.title = element_text(size = 14, face = "bold"))+
  labs(x = expression(Oxygen~(mu*M)), y = "Depth (m)") + 
  geom_hline (aes(yintercept = 20.3), col = "blue", lty = 2)


##BA2
ba2 <- read.oce("Station 36 Cast 1 BioAssay #2-rerun.cnv")
ba2_v2 <- subset (ba2, depth > 4)
summary (ba2)

ba2_v3 = ba2_v2%>%ctdTrim(method = "downcast")%>%ctdDecimate(p = 1)


ba2_v3 = ba2_v3@data%>%
  as_data_frame()

colnames(ba2_v3) [13] <- "par"

ba2temp <- ggplot(data = ba2_v3%>%na.omit(), 
       aes(x = temperature, y = depth))+
  geom_path( col = "black", size = 0.8)+
  scale_y_reverse()+
  scale_x_continuous(position = "top")+
  theme_bw()+
  ggtitle("Bioassay 2")  +
  theme(axis.text = element_text(size = 12, colour = 1),
        axis.title = element_text(size = 14, colour = 1), 
        plot.title = element_text(size = 14, face = "bold"))+
  labs(x = expression(~Temperature~(degree*C)), y = "Depth (m)") + 
  geom_hline (aes(yintercept = 20.3), col = "blue", lty = 2) + 
  ggtitle("Bioassay 2")

ba2par <- ggplot(data = ba2_v3%>%na.omit(), 
       aes(x = par, y = depth))+
  geom_path( col = "black", size = 0.8)+
  scale_y_reverse()+
  scale_x_continuous(position = "top")+
  theme_bw()+
  ggtitle("")  +
  theme(axis.text = element_text(size = 12, colour = 1),
        axis.title = element_text(size = 14, colour = 1), 
        plot.title = element_text(size = 14, face = "bold"))+
  labs(x = expression(PAR~(mu*M~m^-2~S^-1)), y = "Depth (m)") + 
  geom_hline (aes(yintercept = 20.3), col = "blue", lty = 2)

ba2oxygen <- ggplot(data = ba2_v3%>%na.omit(), 
       aes(x = oxygen, y = depth))+
  geom_path( col = "black", size = 0.8)+
  scale_y_reverse()+
  scale_x_continuous(position = "top")+
  theme_bw()+
   ggtitle("")  +
  theme(axis.text = element_text(size = 12, colour = 1),
        axis.title = element_text(size = 14, colour = 1), 
        plot.title = element_text(size = 14, face = "bold"))+
  labs(x = expression(Oxygen~(mu*M)), y = "Depth (m)") + 
  geom_hline (aes(yintercept = 20.3), col = "blue", lty = 2)


#png("ps117_figureS_ctd_20230603.png", width=10, height=13, units="in", res=900) 
#pdf("ps117_figure1_highres_20230519.pdf", width=18.2, height=8.4) 

plot_grid(ba1temp, ba1par, ba1oxygen, 
          NULL, NULL, NULL, 
          ba2temp, ba2par, ba2oxygen, 
          ncol = 3,
          rel_heights = c(1,0.1,1))
        
#dev.off()
```


Figure S_57Fe (show particulate vs dissolved Fe57)
```{r}

x <- ggplot(filter (nonproteindata, !grepl ("T0", iron_treatment)), aes(x= DFe57_nM, 
              
              y=  PFe57_nM, 
                color = temperature_treatment, 
                shape = iron_treatment)) + 
  facet_wrap(~bioassay)+
  geom_point(size =4, alpha =0.6, stroke = 1.5)+
 
  scale_color_manual(name = "Temperature",
                     breaks = c("LT","HT"),
                     labels = c("LT", "HT"),
                     values = c('black', 'darkorange2','black')) +
  scale_shape_manual(name = "Iron",
                       breaks = c("Fe","noFe", "T0"), 
                       labels = c("Fe", "noFe", "T0"),
                       values = c(19,1,5)) +
  xlab (expression (D^{57}*"Fe" ~~ (nM))) + 
  ylab (expression (P^{57}*"Fe" ~~ (nM))) + 
  theme_bw()+
  theme(axis.title=element_text(face = "bold", size = 14, color = "black"),
        axis.text=element_text(size = 12,face = "bold", color = "black"), 
       
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 12, face = "bold"),
        legend.position = "none") + 
  ylim(0,1.5) + 
    xlim(0,0.6)

#png("ps117_pvsd57Fe_20230617.png", width=8, height=7, units="in", res=900) 
ggdraw (x) + draw_plot (vertical_legend , x = 0.9, y= 0.15, width = 0.05, height = 0.05)  
#dev.off()
```

Figure S_protein 
```{r}
ggplot(filter (nonproteindata, grepl ("T0|T8", full_treatment_name) & grepl ("BA2|", bioassay)), 
       aes(y = protein_ug_L*3, #this is to adjust for the wrong protein calculation - need to go back and double check all of this
          x = PON_ug_L*4.78,
          color = tempfe_treatment, 
         shape = tempfe_treatment))+
      scale_color_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c('purple', 'black', 'black','darkorange2', "darkorange2")) +
   scale_shape_manual(name="",
                      breaks = c("T0", "LT_noFe", "LT_Fe", "HT_noFe", "HT_Fe"),
                      labels = c("T0", "LT -Fe", "LT +Fe", "HT -Fe", "HT +Fe"),
                      values = c(5, 1, 19,1, 19)) +
       
    geom_point(size = 4, alpha =0.65, stroke =1.5)+
   labs (x = "Total protein estimate (ug/L) PON*4.78", y = "Theoretical maximum Cd-CA (ug/L) " )+
   facet_wrap(~bioassay) +
   theme_bw()+
   theme(axis.title = element_text (size = 14, face = "bold"), 
        axis.text = element_text(size = 12, color = 'black'), 
        legend.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 10, margin = margin(r = 0.5, unit = 'cm')), 
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 14, face = "bold"))  +
   geom_abline(intercept = 0, slope = 1, color = "Blue", size = 1, linetype = 2)
  
x <- T0_T8_FenoFe(PON_ug_L*4.78, "Total protein estimate (ug/L) ")


#png("ps117_S_protein_20230617.png", width=8, height=6, units="in", res=900) 

ggdraw (x) + draw_plot (vertical_legend , x = 0.9, y= 0.15, width = 0.05, height = 0.05)  

#dev.off()


```

Figure S_refs 
```{r}
alldata <- read.csv("D:/School/PhD/Prelim/temp-nutr/LJ_prelim_data_20221206.csv")

alldata$limitation  <- paste(alldata$NO3_lim, alldata$PO4_lim)
alldata2 <- filter(alldata, !grepl("y",limitation))

library (rcartocolor)

#png("ps117_S_refs_20230617.png", width=8, height=7, units="in", res=900) 

ggplot(filter(alldata, grepl("polar",provenance)), aes(
    x = temperature,
   y = NP_MM))+
     facet_wrap(.~ type) +
  geom_point(size  = 4, alpha = 0.8, aes(color = reference))+
  
  #geom_smooth(method='lm', color  = "black")+
  
  #stat_summary(fun.data =  mean_se, geom = "errorbar", size = 1, width = 0.1) +
  

  
scale_color_carto_d()+
  theme_bw() +
    ylab (expression(N:P~~~(mol:mol))) +
    xlab (expression("Temperature (°C)")) +
    theme(axis.text.x=element_text(face = "bold", size = 12, color = "black", vjust = 0.5, hjust = 0.5),
          axis.title.x=element_text(size=16, color = "black"),
          axis.title.y=element_text(size=16, color = "black"),
          axis.text.y=element_text(face = "bold", size = 12, color = "black"),
          #legend.title = element_text(size = 10),
          legend.title = element_blank(),
          legend.text = element_text(size = 10), 
          legend.position = "bottom",
          strip.text = element_text(size = 11, face = "bold"))+
 geom_hline (aes(yintercept = 16), col = "black", lty = 2, lwd=0.6 )


#dev.off()
```

particulate metals
```{r}
x_vs_time (PMn_nM, PMn~(n*M)) + ylim(0,NA)
x_vs_time (PFe56_nM, PFe56~(n*M)) + ylim(0,NA)
x_vs_time (PFe57_nM, PFe57~(n*M)) + ylim(0,NA)
x_vs_time (PCo_nM, PCo~(n*M)) + ylim(0,NA)
x_vs_time (PNi_nM, PNi~(n*M)) + ylim(0,NA)
x_vs_time (PCu_nM, PCu~(n*M)) + ylim(0,NA)
x_vs_time (PZn_nM, PZn~(n*M)) + ylim(0,NA)
x_vs_time (PCd_nM, PCd~(n*M)) + ylim(0,NA)
```

cellcounts
```{r}
x_vs_time (alldata, daycount, (total_cells_mL), POP~nM)

#png("total_cellcount_vs_time_20220329.png", width=40*0.75, height=30*0.75, units="cm", res=600) 
x_vs_time (total_cells_mL,  Total~cells~(mL^-1)) + scale_y_log10()
#dev.off()

day6 <- filter (nonproteindata, grepl("6", daycount))

Fe_treatment_order <- c('T0', 'noFe', 'Fe')

temp_treatment_order <- c('T0', 'LT', 'HT')

ggplot(day6, aes(
                x = (factor(temperature_treatment, level = temptreatment_order)),
                y = total_cells_mL,
                color = temperature_treatment,
                shape = iron_treatment)) + 
  facet_wrap(~bioassay)+
  stat_summary(fun = mean, geom = "point", show.legend = T, size = 5, stroke = 1.8, alpha = 0.5)+
  stat_summary(fun.data =  mean_se, geom = "errorbar", show.legend = F, size = 1.3, width = 0.1)+
  
  scale_color_manual(name="Temperature",
                     breaks = c("LT","HT"),
                     labels = c("ambient", "warm"),
                     values = c('black','darkorange2', 'black')) +
    scale_shape_manual(name = "Iron",
                     breaks = c("noFe","Fe"), #this removes the T0 from the legend
                     labels = c("noFe", "Fe"),
                     values = c(1,19, 1)) + #this keeps the T0 formatting that I want
  #ylab (expression (PFe[total]~~nM)) + 
  xlab (expression ("")) +
  theme_bw()+
  theme(axis.title.x=element_text(size = 20,face = "bold", color = "black"), 
        axis.text.x=element_text(face = "bold", size = 15, color = "black", vjust = 0.5, hjust = 1),
        axis.title.y=element_text(size=20, color = "black"), 
        axis.text.y=element_text(face = "bold", size = 15, color = "black"),
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 15, face = "bold"),
        legend.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 14))
         
#dev.off()
```

cellcounts - 2023-05-17
```{r}
x_vs_time (Nano2_cells_mL/total_cells_mL, "Phaeo cell counts per mL") +scale_y_log10()

x_vs_time (total_cells_mL, "totalcells") 

x_vs_time (Nano2_volume_um3_mL, "Phaeo volume")
```

cellsize
```{r}
ggplot(filter(nonproteindata, !grepl("9", daycount)), aes(x=daycount, 
              # y=cellnumber_normalized_volume,
              #y=  (total_cell_volume_um3_mL/total_cells_mL),
              y=  total_cell_volume_um3_mL,
                color = temperature_treatment, 
                shape = iron_treatment)) + 
  facet_wrap(~bioassay)+
  #geom_point(size =3, alpha =0.5)+
  stat_summary(fun = mean, geom = "point", show.legend = T, size = 5, stroke = 1.8, alpha = 0.5)+
  #stat_summary(fun = mean, geom = "line", size = 1.1, aes(linetype = iron_treatment))+
  stat_summary(fun.data =  mean_se, geom = "errorbar", size = 1.3, width = 0.1)+
  
  scale_color_manual(name="Temperature",
                     breaks = c("LT","HT"),
                     labels = c("ambient", "warm"),
                     values = c('black', 'darkorange2', 'black')) +
  
  scale_shape_manual(name = "Iron",
                     breaks = c("noFe","Fe", "T0"), #this removes the T0 from the legend
                     labels = c("noFe", "Fe", "T0"),
                     values = c(1,19, 1)) + #this keeps the T0 formatting that I want
  #scale_y_continuous(expand = c(0,5))+
  ylab (expression ("biovolume normalized to cell count")) + 
  xlab (expression ("Day Count")) +
  theme_bw()+
  theme(axis.title.x=element_text(size = 20,face = "bold", color = "black"), 
        axis.text.x=element_text(face = "bold", size = 15, color = "black", vjust = 0.5, hjust = 1),
        axis.title.y=element_text(size=20, color = "black"), 
        axis.text.y=element_text(face = "bold", size = 15, color = "black"),
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 15, face = "bold"),
        legend.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 14))

####
cellsize <- alldata [c(14,61,63,65,67,69,71,73,75,77,79,81,83,85,87,89,91,93,95,97,99,101,103,105)]

meltedsize <- melt(cellsize, id.vars=c("full_treatment_name"))

cellsizemap <- c ("Pico2_cells_mL" = "1.18",
                  "Pico3_cells_mL" = "1.99",
                  "Pico4_cells_mL" = "2.03",
                  "Pico5_cells_mL" = "2.53",
                  "Pico6_cells_mL" = "3",
                  "Pico7_cells_mL" = "3",
                  "Nano1_cells_mL" = "3.23",
                  "Nano2_cells_mL" = "3.71",
                  "Nano3_cells_mL" = "4.31",
                  "Nano4_cells_mL" = "4.93",
                  "Nano7_cells_mL" = "5.17",
                  "Nano8_cells_mL" = "5.52",
                  "Nano9_cells_mL" = "6",
                  "Nano11_cells_mL" = "6.7",
                  "Nano12_cells_mL" = "7.06",
                  "Nano13_cells_mL" = "7.32",
                  "Nano14_cells_mL" = "7.75",
                  "Nano17_cells_mL" = "8.81",
                  "Nano19_cells_mL" = "10.36",
                  "Nano20_cells_mL" = "13.27",
                  "Nano22_cells_mL" = "15.39",
                  "Nano23_cells_mL" = "19",
                  "Nano24_cells_mL" = "19.78")

meltedsize$estimated_size <- cellsizemap [meltedsize$variable]

meltedsize <- meltedsize %>% separate(full_treatment_name, c("bioassay", "timepoint", "temperature_treatment", "iron_treatment", "replicate")) 

meltedsize_day6 <- filter(meltedsize, grepl("T6", timepoint))

meltedsize_day6 <- filter(meltedsize_day6, grepl("BA1", bioassay))

ggplot(meltedsize_day6, aes(
  #x = variable,
  x = estimated_size,
  #x=as.numeric(estimated_size), 
                y=value,
              #y=  total_cell_volume_um3_mL,
               # color = temperature_treatment))+ 
                shape = iron_treatment)) + 
  facet_wrap(~temperature_treatment, ncol=1)+
  geom_point(size =3, alpha =0.5)+
  #stat_summary(fun = mean, geom = "point", show.legend = T, size = 5, stroke = 1.8, alpha = 0.5)+
  stat_summary(fun = mean, geom = "line", size = 1.1, aes(linetype = iron_treatment))+
  #stat_summary(fun.data =  mean_se, geom = "errorbar", size = 1.3, width = 0.1)+
  
  scale_color_manual(name="Temperature",
                     breaks = c("LT","HT"),
                     labels = c("ambient", "warm"),
                     values = c('darkorange2', 'black', 'black')) +
  scale_shape_manual(name = "Iron",
                     breaks = c("noFe","Fe"), #this removes the T0 from the legend
                     labels = c("noFe", "addedFe"),
                     values = c(1,19, 1)) +
   #this keeps the T0 formatting that I want
  #scale_y_continuous(expand = c(0,5))+
 # ylab (expression ("norm relative size")) + 
  #xlab (expression ("Day Count")) +
  theme_bw()+
  theme(axis.title.x=element_text(size = 20,face = "bold", color = "black"), 
        axis.text.x=element_text(face = "bold", size = 10, color = "black" , angle = 90),
        axis.title.y=element_text(size=20, color = "black"), 
        axis.text.y=element_text(face = "bold", size = 15, color = "black"),
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 15, face = "bold"),
        legend.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 14))


#meltedsize$rounded <- round (meltedsize$value, -2)

meltedsize2 <- subset (meltedsize, rounded <200000 )

#https://aslopubs.onlinelibrary.wiley.com/doi/10.1002/lno.12023 - anderson et al. 

# https://adrienne-marshall.github.io/ggplot2_workshop/
#no big effect of temp or iron on community cell size ? 
```

macronutrients 
```{r}
x_vs_time (PON_uM, PON~(mu*M)) + ylim(0,NA)
x_vs_time (POP_nM, POP~(n*M)) + ylim(0,NA)
```

POC-chla
```{r}
alldata$PONchla <- alldata$PON_ug_L/alldata$total_chla_ug_L


ggplot(filter(alldata, grepl("9|", daycount)), aes(x=daycount, 
                y=  PONchla,
                color = temperature_treatment, 
                shape = iron_treatment)) + 
  facet_wrap(~bioassay)+
  #geom_point(size =3, alpha =0.5)+
  stat_summary(fun = mean, geom = "point", show.legend = T, size = 5, stroke = 1.8, alpha = 0.5)+
  #stat_summary(fun = mean, geom = "line", size = 1.1, aes(linetype = iron_treatment))+
  stat_summary(fun.data =  mean_se, geom = "errorbar", size = 1.3, width = 0.1)+
  
  scale_color_manual(name="Temperature",
                     breaks = c("LT","HT"),
                     labels = c("ambient", "warm"),
                     values = c('black', 'darkorange2', 'black')) +
  
  scale_shape_manual(name = "Iron",
                     breaks = c("noFe","Fe", "T0"), #this removes the T0 from the legend
                     labels = c("noFe", "Fe", "T0"),
                     values = c(1,19, 1)) + #this keeps the T0 formatting that I want
  ylab (expression ("PON:Chla (ug:ug)")) + 
  xlab (expression ("Day Count")) +
  theme_bw()+
  theme(axis.title.x=element_text(size = 20,face = "bold", color = "black"), 
        axis.text.x=element_text(face = "bold", size = 15, color = "black", vjust = 0.5, hjust = 1),
        axis.title.y=element_text(size=20, color = "black"), 
        axis.text.y=element_text(face = "bold", size = 15, color = "black"),
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 15, face = "bold"),
        legend.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 14))
```

fasciclin
```{r}
fasc <- read.csv("amCharts_fasc.csv")


 ggplot(filter (fasc, grepl ("SRF", depth)), 
                 aes(x = depth, 
                     y= Iron_5m.umol.l.,
                 size = abundanceExp)) + 
   geom_point(alpha = 0.5) 
   
  stat_summary(fun = mean, geom = "point", size = 4, stroke = 1.5, alpha = 0.5)+
  stat_summary(fun.data =  mean_se, geom = "errorbar", show.legend = F, linewidth = 1, width = 0.1)+
  scale_color_manual(name = "Temperature",
                     breaks = c("LT","HT"),
                     labels = c("LT", "HT"),
                     values = c('black', 'darkorange2','black')) +
  scale_shape_manual(name = "Iron",
                       breaks = c("Fe","noFe", "T0"), 
                       labels = c("Fe", "noFe", "T0"),
                       values = c(19,1,5)) +
  scale_x_discrete(limits = c("T0","noFe","Fe"), 
                   #labels = c("T0", "-Fe", "+Fe"))+
                  
                  labels = expression(bold("T0"), bold("-Fe")~"(T8)", bold("+Fe")~"(T8)"))+


  xlab(expression ("")) +
  ylab (substitute(yaxislabel))+ 
  facet_wrap(~bioassay)+
  theme_bw()+
  theme(axis.text=element_text(face = "bold", size = 12, color = "black"),
        axis.title.x=element_text(size = 14,face = "bold", color = "black"), 
        axis.title.y=element_text(size=16, color = "black"), 
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 12, face = "bold"),
        legend.position = "none") 

```









#######
#######
Everything after this point is experimental plots / things aren't in the paper (yet)
#######
#######

#Other things to think about 
fasciclin (supplement ? )
```{r}
fas <- read.table ("fasciclin_environmental_parameters.csv")
```

Z-score for things ? 
```{r}
silicon <- filter(norm_proteomicsdata, grepl('clust_202 |clust_19688',cluster)) clust_202 |clust_13423|clust_15873

#proteases 2022-04-28
#https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1326381/
#https://pubmed.ncbi.nlm.nih.gov/10849347/
#https://academic.oup.com/plcell/article/12/3/419/6008767
#some proteases function to stabalize the Mn cluster or keep the Mn cluster good and going. So if we have an increase in Mn cluster proteins, we'll need more proteases to maintain these proteins. 
```


#Misc 
proteins (see chucnk with same name below without a function) 2022-08-30
```{r}

protein_clusters ("centric|pennate|Bacillariop", "Plastocyanin", protein = c("Plastocyanin", "Fe-Mn SOD"), shape = c(1,4))

protein_clusters ("centric|pennate|Bacillariop", "", "Fe-Mn|Mn-clus" )

protein_clusters ("centric|pennate|Bacillariop", "", "ISIP" ) +geom_point(size =2, position= position_dodge(width=0.3) )
Mnproteins_phaeo <- protein_clusters ("Phaeo|Prymn", "", "Fe-MnSOD|Mn-clus")

Cuproteins_diatoms <- protein_clusters ("centric|pennate|Bacillariop", "", "Plastocy|Cu-Zn" )
Cuproteins_phaeo <- protein_clusters ("Phaeo|Prymn", "", "Plastocy|Cu-Zn")

Znproteins_diatoms <- protein_clusters ("centric|pennate|Bacillariop", "", "ZIP|RNA" )
Znproteins_phaeo <- protein_clusters ("Phaeo|Prymn", "", "ZIP|RNA")

Cdproteins_diatoms <- protein_clusters ("centric|pennate|Bacillariop", "", "CA" )
Cdproteins_phaeo <- protein_clusters ("Phaeo|Prymn", "", "CA")

```

legend
```{r}
x<- c(1,2,3)
y <- c("a", "b", "c")
colr <- c (-1.2,0, 1.7)
z <- data.frame(x,y, colr)

legend3b <- ggplot(z, aes (x,y, fill = colr))+
   geom_point()+
    scale_fill_gradientn(colors = c("deepskyblue2","white","darkred"),
         values = rescale(c(-1.2,0,1.7)),
         guide = "colorbar", limits=c(-1.2,1.7),
         name = "Row Z-score")+
   theme(legend.text = element_text(face = "bold", size=13),
        legend.title = element_text (face = "bold", size = 13, angle = 90), 
        legend.title.align=0.5)+
   guides(fill = guide_colourbar(title.position = "left", barwidth = 1.5, barheight = 15))
  

legend3b <- get_legend(legend3b)
legend3b <- as_ggplot(legend3b)
```

proteins heatmap
```{r}
protein_clusters_hm ("pennate", "Fe-Mn|Mn-clus|Anhyd", "BA" )
protein_clusters ("Phaeo|Prymn", "", "Fe-Mn|Mn-clus|Anhyd", "BA2")

Cuproteins_centric <- protein_clusters ("centric", "Plastocy|Cu-Zn", "BA" )
Cuproteins_pennate <- protein_clusters ("pennate", "Plastocy|Cu-Zn", "BA" )
Cuproteins_phaeo <- protein_clusters ("Phaeo|Prymn", "Plastocy|Cu-Zn", "BA")

Znproteins_diatoms <- protein_clusters ("centric|pennate|Bacillariop", "", "ZIP|RNA" )
Znproteins_phaeo <- protein_clusters ("Phaeo|Prymn", "", "ZIP|RNA")

Cdproteins_diatoms <- protein_clusters ("centric|pennate|Bacillariop", "", "CA" )
Cdproteins_phaeo <- protein_clusters ("Phaeo|Prymn", "", "CA")

```

plot
```{r}
row1 <- plot_grid(NULL, metals_Mn, NULL, Mnproteins_diatoms, NULL, Mnproteins_phaeo, 
          nrow = 1,
          rel_widths = c(0.05, 1, 0.1, 1.3, 0.05, 1.3),
          labels = c("", "A", "", "B", "", "C"), 
          align = 'h')
row2 <- plot_grid(NULL, metals_Cu, NULL, Cuproteins_diatoms, NULL, Cuproteins_phaeo, 
          nrow = 1,
          rel_widths = c(0.05, 1, 0.1, 1.3, 0.05, 1.3),
          labels = c("", "D", "", "E", "", "F"), 
          align = 'h')
row3 <- plot_grid(NULL, metals_Zn, NULL, Znproteins_diatoms, NULL, Znproteins_phaeo, 
          nrow = 1,
          rel_widths = c(0.05, 1, 0.1, 1.3, 0.05, 1.3),
          labels = c("", "G", "", "H", "", "I"), 
          align = 'h')

plot_grid(row1, row2, row3, 
          nrow = 3)
```

plot_v2
```{r}
row1 <- plot_grid(NULL, metals_Mn, NULL, Mnproteins_diatoms, NULL, Mnproteins_phaeo, 
          nrow = 1,
          rel_widths = c(0.05, 1, 0.1, 1.3, 0.05, 1.3),
          labels = c("", "A", "", "B", "", "C"), 
          align = 'h')

plot_grid(metals_Cu, Cuproteins_diatoms, Cuproteins_phaeo, 
          nrow = 3, align = 'hv', 
          rel_heights = c(1,0.5,0.5), 
          axis = 'l')


           
row3 <- plot_grid(NULL, metals_Zn, NULL, Znproteins_diatoms, NULL, Znproteins_phaeo, 
          nrow = 1,
          rel_widths = c(0.05, 1, 0.1, 1.3, 0.05, 1.3),
          labels = c("", "G", "", "H", "", "I"), 
          align = 'h')


plot_grid(row1, row2, row3, 
          nrow = 3)
```


#photophysiology
raw PAM curve fitting
```{r}
setwd("D:/School/PhD/PS117/data/photophysiology/")

rawpam <- read.csv("ps117_rawPAM.csv", header = T)
rawpam <- filter (rawpam, !grepl ("blank", iron_treatment) )

rawpam$fvfm <-  (rowMeans(rawpam[,c(47:72)]) - rowMeans(rawpam[,c(17:32)]))/ rowMeans(rawpam[,c(47:72)])

rawpam2 <- rawpam [c(143,14)]


ggplot(rawpam, aes(x=fvfm, 
                y=YII)) + 
  geom_point(size =3, alpha =0.4, aes(color = temperature_treatment, 
                shape = iron_treatment))+
  geom_abline(intercept = 0, slope = 1, color = "Blue", size = 1, linetype = 2)+
  scale_x_continuous(expand = c(0, 0), limits = c(0, 0.7)) + 
  scale_y_continuous(expand = c(0, 0), limits = c(0, 0.7))+
  #geom_smooth(method='lm', formula= y~x, color = "green")+

    scale_shape_manual(name = "Iron",
                     breaks = c("T0", "noFe","Fe"), #this removes the T0 from the legend
                     labels = c("T0","noFe", "addedFe"),
                     values = c(1, 1, 19)) + #this keeps the T0 formatting that I want
  scale_color_manual(name="Temperature",
                     breaks = c("T0","LT","HT"),
                     labels = c("T0","ambient", "warm"),
                     values = c('black', 'black', "red")) +
  ylab (expression (PAM~~F[v]/F[m])) + 
  xlab (expression (Loay~~F[v]/F[m])) +
  theme_bw()+
  theme(axis.title.x=element_text(size = 20,face = "bold", color = "black"), 
        axis.text.x=element_text(face = "bold", size = 15, color = "black", vjust = 0.5, hjust = 1),
        axis.title.y=element_text(size=20, color = "black"), 
        axis.text.y=element_text(face = "bold", size = 15, color = "black"),
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 15, face = "bold"), 
        #legend.position = "top", # c(0.04,0.88), 
        legend.title = element_text(size = 12, face = "bold"),
        legend.text = element_text(size = 12))+ 
        
   guides(color = guide_legend(override.aes = list(size= 4)))+
   guides(shape = guide_legend(override.aes = list(size= 4)))


ggplot(rawpam, aes(x=daycount, 
                y=YII,
                color = temperature_treatment, 
                shape = iron_treatment)) + 
  #geom_point(size =5, alpha =0.5)+
  stat_summary(fun = mean, geom = "point", size = 2, stroke = 2, alpha = 0.6)+
  stat_summary(fun.data =  mean_se, geom = "errorbar", size = 0.8, width = 0.2)+
  #geom_boxplot()+
  scale_shape_manual(name = "Iron",
                     breaks = c("T0", "noFe","Fe"), #this removes the T0 from the legend
                     labels = c("T0","noFe", "addedFe"),
                     values = c(1, 1, 19)) + #this keeps the T0 formatting that I want
  scale_color_manual(name="Temperature",
                     breaks = c("T0","LT","HT"),
                     labels = c("T0","ambient", "warm"),
                     values = c('black', 'black', 'red')) +
  #ylab (expression (F[v]/F[m])) + 
  xlab (expression ("Day Count")) +
  theme_bw()+
  theme(axis.title.x=element_text(size = 20,face = "bold", color = "black"), 
        axis.text.x=element_text(face = "bold", size = 15, color = "black", vjust = 0.5, hjust = 1),
        axis.title.y=element_text(size=20, color = "black"), 
        axis.text.y=element_text(face = "bold", size = 15, color = "black"),
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 15, face = "bold"), 
        #legend.position = "top", # c(0.04,0.88), 
        legend.title = element_text(size = 12, face = "bold"),
        legend.text = element_text(size = 12)) +
   guides(color = guide_legend(override.aes = list(size= 4)))+
   guides(shape = guide_legend(override.aes = list(size= 4)))
```

plots
```{r}
#photophysiology
#png("fvfm_timeseries.png", width=40*0.75, height=30*0.75, units="cm", res=600) 
ggplot(alldata, aes(x=daycount, 
                y=fvfm,
                color = temperature_treatment, 
                shape = iron_treatment)) + 
  facet_wrap(~bioassay)+
  geom_point(size = 5, alpha =0.5, stroke =1)+
  #stat_summary(fun = mean, geom = "point", size = 5, stroke = 1.8, alpha = 0.5)+
   #stat_summary(fun.data =  mean_se, geom = "errorbar", size = 1.3, width = 0.1)+
    scale_color_manual(name="Temperature",
                     breaks = c("LT","HT"),
                     labels = c("ambient", "warm"),
                     values = c('darkorange2', 'black', 'black')) +
   
    scale_shape_manual(name = "Iron",
                     breaks = c("Fe","noFe"), #this removes the T0 from the legend
                     labels = c("Fe", "noFe"),
                     values = c(19, 1, 1)) + #this keeps the T0 formatting that I want
    #scale_y_continuous(expand = c(0,5))+
  ylab (expression (F[v]/F[m])) + 
  xlab (expression ("Day Count")) +
  theme_bw()+
  theme(axis.title.x=element_text(size = 20,face = "bold", color = "black"), 
        axis.text.x=element_text(face = "bold", size = 15, color = "black", vjust = 0.5, hjust = 1),
        axis.title.y=element_text(size=20, color = "black"), 
        axis.text.y=element_text(face = "bold", size = 15, color = "black"),
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 15, face = "bold"),
        legend.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 14))

#dev.off()

#############

photophysi <- filter (alldata, grepl ("T0", full_treatment_name) )

#(DFe56_nM+DFe57_nM)

#png("fvfm_vs_iron.png", width=40*0.75, height=30*0.75, units="cm", res=600) 
ggplot(alldata, aes(x=(PFe57_nM), 
                y=fvfm,
                color = temperature_treatment, 
                shape = iron_treatment)) + 
   #stat_summary(fun = mean, geom = "point", size = 2, stroke = 2, alpha = 0.6)+
#  stat_summary(fun.data =  mean_se, geom = "errorbar", size = 0.8, width = 0.2)+
  facet_wrap(~bioassay)+
  geom_point(alpha =0.5, stroke =1, size = 6)+
  scale_color_manual(name="Temperature",
                     breaks = c("LT","HT"),
                     labels = c("ambient", "warm"),
                     values = c('darkorange2', 'black', 'black')) +
       scale_shape_manual(name = "PFe57",
                     breaks = c("Fe","noFe"), #this removes the T0 from the legend
                     labels = c("Fe", "noFe"),
                     values = c(19, 1, 1)) + #this keeps the T0 formatting that I want
  #scale_y_continuous(expand = c(0,5))+
  ylab (expression (F[v]/F[m])) + 
  xlab (expression ("Iron (nM)")) +
  theme_bw()+
  theme(axis.title.x=element_text(size = 20,face = "bold", color = "black"), 
        axis.text.x=element_text(face = "bold", size = 15, color = "black", vjust = 0.5, hjust = 1),
        axis.title.y=element_text(size=20, color = "black"), 
        axis.text.y=element_text(face = "bold", size = 15, color = "black"),
        strip.background =element_rect(fill = "white"),
        strip.text.x = element_text(size = 15, face = "bold"),
        legend.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 14)) +
       expand_limits(y = 0:0.6) 
```