# Author: Shirley Pu
# Code used for the paper "Examining Environmental Matrix Effects on Quantitative Non-Targeted Analysis Estimates of Per- and Polyfluoroalkyl Substances"
# Contains code used for supporting information

# Set working directory
#setwd("~/Shirley/scripts/PFAS Water Matrices")

# Load libraries
library(tidyverse)
library(investr)
library(ggridges)
library(lme4)
library(lmerTest)

# Source custom functions
source("Pu_et_al_Env_Matrices_functions.R")

# Read in data
load("Matrices_Data_cleaned_final.RData")

# Read in results
load("Pu_et_al_qNTA_Env_Matrices_final_results.RData")

# Unique media (includes media that aren't used for final quantitation approach comparison)
uniq_media <- Matrices_Data_cleaned %>% 
  pull('Media') %>% 
  unique()

## Comparing RF distributions of original reconstitution method (Reconstitution Method 1) and more nonpolar reconstitution method (Reconstitution Method 2) for DW, Sludge, STD
reconstitution_RF_plot_data <- Matrices_Data_cleaned %>%
  filter(Media %in% c("DW","DW95","Sludge","Sludge95","Solvent","Solvent95")) %>%
  mutate(Media_Name = gsub("95","",Media) %>%
           factor(levels=c("Solvent","DW","Sludge")),
         `Reconstitution Method` = case_when(grepl("95",Media) ~ "1", TRUE ~ "2"))

# Figure S1
ggplot(reconstitution_RF_plot_data,aes(x=RF_uM,y=reorder(Chem,-RT),color=`Reconstitution Method`,fill=`Reconstitution Method`))+
  geom_density_ridges(scale=2,alpha=.5,bandwidth=0.1)+
  scale_x_log10(breaks = scales::trans_breaks("log10", function(x) 10^x),
                labels = scales::trans_format("log10", scales::math_format(10^.x)))+
  ylab("Chemical")+
  xlab(expression(paste("Response Factor (", mu, "M"^"-1",")")))+
  theme_bw(base_size=22)+
  facet_wrap(facets = vars(Media_Name))
# Replace fill with grey and white instead of red and blue?

# Calibration curve parameters
cal_curve_params <- data.frame()
# Calibration curves
# Nested loop for making plots (loop over Media, loop over Chem in that Media)
for(i in uniq_media){
  uniq_media_rows<-Matrices_Data_cleaned %>% filter(Media==i)
  uniq_media_chem<-unique(uniq_media_rows %>% pull('Chem'))
  for(j in uniq_media_chem){
    uniq_media_chem_rows<-uniq_media_rows %>% filter(Chem==j)
    plot_lm<-lm(data=uniq_media_chem_rows,LogArea~LogConc_uM)
    # get prediction interval
    pred_interval_log<-predict(plot_lm, interval = "prediction", level = 0.95) %>% #log space pred intervals
      as_tibble() %>%
      rename(fit_log = fit,
             upr_log = upr,
             lwr_log = lwr)
    # make the plot
    plot_data<-cbind(uniq_media_chem_rows,pred_interval_log)
    cf<-round(coef(plot_lm),2)
    ggtitle_eq<-paste0("Log\u2080\u2081(Area)",cf[1],ifelse(sign(cf[2])==1,"+","-"),abs(cf[2]),"Log\u2080\u2081(Conc) \u0370M")
    ggtitle_r2<-paste0("\n R\u00b2 = ",formatC(summary(plot_lm)$r.squared,digits=2))
    ggplot(data=plot_data,aes(x=LogConc_uM))+
      theme_minimal()+
      geom_point(aes(y = LogArea))+ #measurements
      geom_smooth(aes(y = LogArea), method = lm, se = TRUE)+ #linear fit with 95% confidence interval
      geom_line(aes(y = upr_log), linetype = 2, color = "red") + #upr and lwr prediction intervals from predict()
      geom_line(aes(y = lwr_log), linetype = 2, color = "red") +
      xlab(expression(paste("Log"[10]*"(Conc) ",mu,"M")))+
      ylab(expression("Log"[10]*"Area"))+
      ggtitle(paste("Calibration curve for",i,j),paste(ggtitle_eq,ggtitle_r2))
    # save plot - just for matrices used for final quantitation approach comparison
    if(i %in% c("Solvent","DW","Sludge")){
      ggsave(paste(gsub(":","",j),i,"Cal Curve uM.png"),bg='white')
    }
    # store calibration curve slope and Rsquared
    cal_curve_params <- rbind(cal_curve_params,data.frame(Media=i,Chem=j,slope=coef(plot_lm)[2],R2=summary(plot_lm)$r.squared))
  }
}
rownames(cal_curve_params) <- c()

# Histograms of slope, R-squared, looking only at DW, Sludge, STD data (Reconstitution Methods 1 and 2)
cal_curve_params_RM <- cal_curve_params %>% 
  filter(Media %in% c("DW","DW95","Sludge","Sludge95","Solvent","Solvent95")) %>%
  mutate(RM = case_when(grepl("95",Media) ~ "Reconstitution Method 1",
                        TRUE ~ "Reconstitution Method 2"),
         Media = gsub("95","",Media) %>%
           factor(levels=c("Solvent","DW","Sludge")))

# Figure S2
ggplot(cal_curve_params_RM,aes(slope))+
  geom_histogram()+
  xlab("Slope")+
  ylab("Number of curves")+
  theme_bw(base_size=24)+
  facet_grid(rows = vars(Media),cols = vars(RM),switch="y")+
  theme(panel.spacing.x=unit(28,"pt"))+ 
  scale_y_continuous(breaks=c(2,4,6,8),position="right")

ggplot(cal_curve_params_RM,aes(R_squared))+
  geom_histogram()+
  xlab(expression(R^2))+
  ylab("Number of curves")+
  theme_bw(base_size=24)+
  facet_grid(rows = vars(Media),cols = vars(RM),switch="y")+
  theme(panel.spacing.x=unit(28,"pt"))+ 
  scale_y_continuous(position = "right")

# Figure S3 generated in Pu_et_al_Env_Matrices_Supplementary_File_2.R

# Just save cal curve parameters for Solvent, DW, Sludge
cal_curve_params <- cal_curve_params %>% filter(Media %in% c("Solvent","DW","Sludge"))
save(cal_curve_params,file="Pu_et_al_Matrices_Env_Matrices_cal_curve_params.RData")

# Figure S4
# Plot LogRF vs LogIE for all matrices on one plot
ggplot(Matrices_Data_IE %>%
         rename(Matrix=Media),aes(x=logIE_predicted_LOO,y=LogRF_uM,color=Matrix,shape=Matrix))+
  geom_point(aes(x=logIE_predicted_LOO,y=LogRF_uM))+
  xlab(expression('Log'[10]*'(IE)'))+
  ylab(expression('Log'[10]*'(RF)'))+
  ggtitle(paste0("Log(RF) vs. Log(IE)"))+
  theme_bw(base_size=22)

# PFBA has much lower logIE than all the other compounds
# Check whether it has high leverage (the hat matrix has leverage for each data point)
# Note: hatvalues() are conditional on estimated variance-covariance of the random effects
RF_vs_IE_per_matrix_hat <- vector(mode="list",length=3L)
for(i in 1:length(uniq_media_RM2)){
  curr_RF_vs_IE_per_matrix <- lme4::lmer(LogRF_uM ~ logIE_predicted_LOO + (1|Chem),Matrices_Data_IE %>% filter(Media==uniq_media_RM2[i]))
  RF_vs_IE_per_matrix_hat[[i]] <- data.frame(model.frame(curr_RF_vs_IE_per_matrix),Leverage=hatvalues(curr_RF_vs_IE_per_matrix ))
}

# can we facet this? (would need to store the hat results as a column and matrix as a column)
for(i in 1:length(RF_vs_IE_per_matrix_hat)){
  ggplot(RF_vs_IE_per_matrix_hat[[i]],aes(x=logIE_predicted_LOO,y=Leverage))+
    geom_point()+
    scale_y_continuous(breaks=c(0.04,0.05,0.06,0.07,0.08),
                       limits=c(0.04,0.08))+
    xlab(expression('Log'[10]*'(IE)'))+
    ylab("Leverage")+
    ggtitle(paste0("Leverage vs. LogIE ",uniq_media_RM2[i]))+
    theme_bw(base_size=22)+
    theme(axis.text=element_text(size=22))
  ggsave(paste0("Pu_et_al_Env_Matrices_Leverage_vs_LogIE_",uniq_media_RM2[i],".png"),width=6,height=4,units="in")
}

# Use just DW, Sludge, and Solvent RM2 for the full quantitative approach comparison
Matrices_Data_cleaned_RM2 <- Matrices_Data_cleaned %>% 
  filter(Media %in% c("DW","Sludge","Solvent")) %>%
  mutate(Media = factor(Media,levels=c("Solvent","DW","Sludge")),
         Chem = gsub("FTS42","4:2 FTS",Chem) %>%
           gsub("FTS62","6:2 FTS",.) %>%
           gsub("FTS82","8:2 FTS",.))

uniq_media_RM2 <- unique(Matrices_Data_cleaned_RM2 %>% pull('Media') %>% as.character())

## PERFORMANCE METRICS ACROSS CHEMICALS USING EACH qNTA APPROACH IN EACH MATRIX
all_approaches_results_perchemmedia <- all_approaches_results_long %>% 
  group_by(Chem,Media,approach) %>%
  summarize_at(vars("AQ","AAQ","CLFR"), ~ median(abs(.x),na.rm=TRUE))

all_approaches_results_perchemmedia_RP <- all_approaches_results_long %>% 
  group_by(Chem,Media,approach) %>%
  summarize_at(vars("either_exceeded"), ~ 1 - (sum(.x[!is.na(.x)])/length(.x[!is.na(.x)]))) %>%
  rename_with(~gsub("either_exceeded","ORP",.x))

# Add a column with retention time (RT) for ordering the chemicals
perchemmedia_RT <- Matrices_Data_cleaned_RM2 %>% 
  group_by(Chem,Media) %>% 
  summarize_at("RT",median)

# Cumulative distribution plots
# Figure S5
# AAQ
AAQ_ECDF_data <- all_approaches_results_long %>%
  drop_na() %>% 
  group_by(approach) %>%
  mutate(rank = rank(AAQ)/length(AAQ)) 

ggplot(AAQ_ECDF_data %>% 
         mutate(approach = gsub("cal_curve","Calibration Curve",approach) %>%
                  gsub("expert_boot","Expert Bootstrap",.) %>%
                  gsub("global_boot", "Global Bootstrap",.) %>%
                  gsub("IE_regression","IE regression",.)), aes(x = AAQ, y = rank*100, group = approach, color = approach)) +
  #geom_line(size = 1)+
  geom_point()+
  scale_x_continuous(n.breaks = 6)+
  scale_y_continuous(n.breaks = 10)+
  scale_color_manual(limits=c("Calibration Curve","Expert Bootstrap","Global Bootstrap","IE regression"),
                     values=c("orange","blue","black","purple"),
                     labels=c("CC",expression(BRF[expert]),expression(BRF[global]),expression(IE[reg])),
                     name="Legend")+
  facet_grid(cols = vars(Media))+
  theme_bw(base_size = 20)+
  xlab("AAQ")+
  ylab("Cumulative Percentile")+
  theme(legend.position = c(0.9, 0.4),
        legend.title = element_blank(),
        legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid'))

# CLFR
CLFR_ECDF_data <- all_approaches_results_long %>%
  drop_na() %>% 
  group_by(Media,approach) %>%
  mutate(rank = rank(CLFR)/length(CLFR))

ggplot(CLFR_ECDF_data %>% 
         mutate(approach = gsub("cal_curve","Calibration Curve",approach) %>%
                  gsub("expert_boot","Expert Bootstrap",.) %>%
                  gsub("global_boot", "Global Bootstrap",.) %>%
                  gsub("IE_regression","IE regression",.)), aes(x = CLFR, y = rank*100, group = approach, color = approach)) +
  #geom_line(size = 1)+
  geom_point()+
  scale_x_continuous(n.breaks = 7)+
  scale_y_continuous(n.breaks = 10)+
  scale_color_manual(limits=c("Calibration Curve","Expert Bootstrap","Global Bootstrap","IE regression"),
                     values=c("orange","blue","black","purple"),
                     labels=c("CC",expression(BRF[expert]),expression(BRF[global]),expression(IE[reg])),
                     name="Legend")+
  facet_grid(cols = vars(Media))+
  theme_bw(base_size = 20)+
  xlab("CLFR")+
  ylab("Cumulative Percentile")+
  theme(legend.position = c(0.9, 0.4),
        legend.title = element_blank(),
        legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid'))

# Figure S6
# Random per-chemical plot
random_per_chem <- data.frame()
DW_data <- all_approaches_results_long %>% filter(Media == "DW")
uniq_chem_DW <- DW_data %>% pull('Chem') %>% unique()

set.seed(1)
for(i in 1:length(uniq_chem_DW)){
  per_chem_rows <- DW_data %>% filter(Chem == uniq_chem_DW[i]) %>% drop_na()
  random_per_chem <- rbind(random_per_chem,per_chem_rows[sample(nrow(per_chem_rows),size=1),])
}

random_per_chem_plot<-dplyr::left_join(random_per_chem,all_approaches_results_wide)

# Facet by Media
random_per_chem_facet_plot <- dplyr::left_join(random_per_chem %>% select(Conc,Chem,Rep),all_approaches_results_wide) %>% 
  mutate(Media = factor(Media,levels=c("Solvent","DW","Sludge")),
         Chem = gsub("FTS42","4:2 FTS",Chem) %>%
           gsub("FTS62","6:2 FTS",.) %>%
           gsub("FTS82","8:2 FTS",.))
ggplot(data=random_per_chem_facet_plot %>% filter(Chem != "PFBA"))+
  geom_point(aes(x=reorder(Chem,Conc),y=estimate_ng_g_cal_curve,color="Calibration Curve"),position=position_nudge(x=-0.2))+
  geom_linerange(aes(x=reorder(Chem,Conc),ymin=lower_ng_g_cal_curve,ymax=upper_ng_g_cal_curve),color="orange",position=position_nudge(x=-0.2))+
  geom_point(aes(x=reorder(Chem,Conc),y=estimate_ng_g_expert_boot,color="Expert surrogates bootstrap"))+
  geom_linerange(aes(x=reorder(Chem,Conc),ymin=lower_ng_g_expert_boot,ymax=upper_ng_g_expert_boot),color="blue")+
  geom_point(aes(x=reorder(Chem,Conc),y=estimate_ng_g_global_boot,color="Global surrogates bootstrap"),position=position_nudge(x=0.2))+
  geom_linerange(aes(x=reorder(Chem,Conc),ymin=lower_ng_g_global_boot,ymax=upper_ng_g_global_boot),position=position_nudge(x=0.2),color="black")+
  geom_point(aes(x=reorder(Chem,Conc),y=estimate_ng_g_IE_regression,color="IE Regression"),position=position_nudge(x=0.4))+
  geom_linerange(aes(x=reorder(Chem,Conc),ymin=lower_ng_g_IE_regression,ymax=upper_ng_g_IE_regression),position=position_nudge(x=0.4),color="purple")+
  geom_segment(aes(x=as.numeric(reorder(Chem,Conc))-0.5,xend=as.numeric(reorder(Chem,Conc))+0.5,y=Conc,yend=Conc),color="red")+
  scale_y_log10(breaks = scales::trans_breaks("log10", function(y) 10^y),
                labels = scales::trans_format("log10", scales::math_format(10^.x)))+
  theme_bw(base_size = 14)+
  theme(axis.text.x=element_text(angle=90,vjust=0.5,hjust=1),text=element_text(size=14))+
  ggtitle("Concentration predictions with 95% confidence intervals")+
  ylab("Concentration (ng/mL)")+
  xlab("Chemical")+
  scale_color_manual(limits=c("Calibration Curve","Expert surrogates bootstrap","Global surrogates bootstrap","IE Regression","True concentration"),
                     values=c("orange","blue","black","purple","red"),
                     labels=c("CC",expression(BRF[expert]),expression(BRF[global]),expression(IE[reg]),"True concentration"),
                     name="Legend")+
  guides(color=guide_legend(override.aes=list(linetype=c(0,0,0,0,0))))+
  theme(legend.position = c(0.9, 0.15),
        legend.title = element_blank(),
        legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid'))+
  facet_grid(cols=vars(Media))

# Figure S7
# AQ
ggplot(all_approaches_results_long %>% 
         mutate(approach = gsub("cal_curve","Calibration Curve",approach) %>%
                  gsub("expert_boot","Expert Bootstrap",.) %>%
                  gsub("global_boot", "Global Bootstrap",.) %>%
                  gsub("IE_regression","IE regression",.) %>%
                  factor(levels=c("Calibration Curve","Expert Bootstrap","Global Bootstrap","IE regression"),
                         labels=c("CC",expression(BRF[expert]),expression(BRF[global]),expression(IE[reg])))), aes(approach, AQ)) + 
  geom_boxplot(outlier.shape=NA) + #avoid plotting outliers twice
  geom_jitter(position=position_jitter(width=.2,height=0)) +
  scale_x_discrete(labels = ~parse(text = .x))+
  scale_y_log10(breaks = scales::trans_breaks("log10", function(y) 10^y),
                labels = scales::trans_format("log10", scales::math_format(10^.x)))+
  
  theme_bw() +
  facet_grid(cols = vars(Media)) +
  xlab("Approach")+
  ylab("Accuracy Quotient (AQ)") +
  theme(axis.text.x=element_text(angle=90,vjust=0.5,hjust=1),text=element_text(size=18)) +
  theme_bw()

# Graphical abstract heatmaps
all_approaches_results_perchemmedia_heatmap <- all_approaches_results_perchemmedia %>% 
  dplyr::left_join(all_approaches_results_perchemmedia_RP) %>%
  mutate(Chem = gsub("FTS42","4:2 FTS",Chem) %>%
           gsub("FTS62","6:2 FTS",.) %>%
           gsub("FTS82","8:2 FTS",.)) %>%
  dplyr::left_join(perchemmedia_RT) %>%
  mutate(across(AQ:ORP,~round(.,digits=2))) %>%
  arrange(RT)

# Median by media, approach for AQ, AAQ, CLFR
all_approaches_perapproachmedia_heatmap <- all_approaches_results_long %>%
  group_by(Media,approach) %>%
  summarize_at(vars("AQ","AAQ","CLFR"), ~ median(abs(.x),na.rm=TRUE)) %>%
  dplyr::left_join(all_approaches_results_long %>% 
                     group_by(Media,approach) %>%
                     summarize_at(vars("either_exceeded"), ~ 1 - (sum(.x[!is.na(.x)])/length(.x[!is.na(.x)]))) %>%
                     rename_with(~gsub("either_exceeded","ORP",.x))) %>%
  mutate(Media = factor(Media,levels=c("Solvent","DW","Sludge")))

all_approaches_perapproachmedia_heatmap <- all_approaches_perapproachmedia_heatmap %>% 
  mutate(approach = gsub("cal_curve","Calibration Curve",approach) %>%
           gsub("expert_boot","Expert Bootstrap",.) %>%
           gsub("global_boot", "Global Bootstrap",.) %>%
           gsub("IE_regression","IE regression",.) %>%
           factor(levels=c("IE regression","Global Bootstrap","Expert Bootstrap","Calibration Curve"),
                  labels=c(expression(IE[reg]),expression(BRF[global]),expression(BRF[expert]),"CC")))

# Normalizing AQ, AAQ, CLFR, ORP values
# Don't include calibration curve approach, just compare qNTA approaches
all_approaches_perqNTAmedia_heatmap_norm <- all_approaches_perapproachmedia_heatmap %>%
  filter(approach != "CC") %>%
  ungroup() %>%
  mutate(AQ_norm = (AQ - min(AQ))/(max(AQ)-min(AQ)),
         AAQ_norm = (AAQ - min(AAQ))/(max(AAQ)-min(AAQ)),
         CLFR_norm = (CLFR - min(CLFR))/(max(CLFR)-min(CLFR)),
         ORP_norm = (ORP - min(ORP))/(max(ORP)-min(ORP))) %>%
  mutate(AAQ_norm_rev = 1-AAQ_norm,
         CLFR_norm_rev = 1-CLFR_norm)

# Per media x per approach
# AAQ
ggplot(all_approaches_perqNTAmedia_heatmap_norm, aes(Media, approach, fill=AAQ_norm_rev)) + 
  geom_tile(colour = "white", na.rm = TRUE) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  scale_fill_gradient("AAQ_norm_rev",low="lemonchiffon",high="seagreen") +
  #geom_text(aes(label=paste(AAQ)), size=4)+
  scale_y_discrete(labels = ~parse(text = .x))+
  theme(strip.background = element_rect(fill="white"),legend.position="None") +
  xlab("Matrix")+
  ylab("Approach")

# CLFR
ggplot(all_approaches_perqNTAmedia_heatmap_norm, aes(Media, approach,  fill=CLFR_norm_rev)) + 
  geom_tile(colour = "white", na.rm = TRUE) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  scale_fill_gradient("CLFR_norm_rev",low="lemonchiffon",high="seagreen") +
  #geom_text(aes(label=paste(CLFR)), size=4)+
  scale_y_discrete(labels = ~parse(text = .x))+
  theme(strip.background = element_rect(fill="white"),legend.position="None") +
  xlab("Matrix")+
  ylab("Approach")

# ORP
ggplot(all_approaches_perqNTAmedia_heatmap_norm, aes(Media, approach, fill=ORP_norm)) + 
  geom_tile(colour = "white", na.rm = TRUE) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  scale_fill_gradient("ORP_norm",low="lemonchiffon",high="seagreen") +
  #geom_text(aes(label=paste(ORP)), size=4)+
  scale_y_discrete(labels = ~parse(text = .x))+
  theme(strip.background = element_rect(fill="white"),legend.position="None") +
  xlab("Matrix")+
  ylab("Approach")

perchem_metrics <- all_approaches_results_long %>%
  group_by(Chem,Media,approach) %>%
  summarize_at(vars("AQ","AAQ","CLFR"), ~ median(abs(.x),na.rm=TRUE)) %>%
  dplyr::left_join(all_approaches_results_long %>% 
                     group_by(Media,approach) %>%
                     summarize_at(vars("either_exceeded"), ~ 1 - (sum(.x[!is.na(.x)])/length(.x[!is.na(.x)]))) %>%
                     rename_with(~gsub("either_exceeded","ORP",.x))) %>%
  mutate(Media = factor(Media,levels=c("Solvent","DW","Sludge")))
