160519_run2_calibration.Rmd

---
title: 'Run 2: Oman 2015 NO$_3$ $\delta^{15}N$ and $\delta^{18}O$'
subtitle: "Source file: 160519_run2_calibration.Rmd"
date: "`r Sys.Date()`"
output:
  html_document: 
    css: stylesheet.css
    fig_caption: yes
    number_sections: yes
    df_print: paged
    toc: yes
    toc_float: true
    toc_depth: 3
    code_folding: show
editor_options:
  chunk_output_type: console
---


```{r setup, include=FALSE}
library(isoreader)
library(isoprocessor)
library(isorunN2O)
library(tidyverse)
library(latex2exp)
library(readxl)
library(ggplot2)
library(plotly)
library(stringr)
library(knitr)

opts_knit$set(root.dir = ".", output_dir = file.path("docs"))
source(file.path("scripts", "functions.R"))
```


The dissolved $NO_3$ isotopic composition was measured using the denitrifier method coupled to GC-IRMS as described in the methods section of the manuscript. The following code steps through the calibration of measured isotopic values. The fluid samples processed here in this document were collected in Oman in the 2015 field season. 

Run 2 was based on acidified samples (for IC analysis) that were acidified using hydrochloric acid to a pH of ~3. At this low pH, all nitrite would have reacted away during storage so the measured sample was pure nitrate rather than $NO_x$.


This processing script uses isoreader version `r packageVersion("isoreader")`, isorunN2O version `r packageVersion("isorunN2O")`, and isoprocessor version `r packageVersion("isoprocessor")`.


## Load Raw Data
### To create .rds file from folder of .dxf files for fast loading
```{r eval = FALSE}
#Load raw data from .dxf files
  folder <- "160519_OM_run2"
  raw_data <-
    file.path("data", "raw", folder) %>%
    iso_read_continuous_flow()
  #look for problematic files
 if (iso_has_problems(raw_data)) {
   message("are these problems okay?")
    print(iso_get_problems_summary(raw_data))
  }

#omit any files with problems and export to an .rds file
  raw_data %>%
    iso_omit_files_with_problems() %>%
    iso_save("data/raw/run2_raw") 
```

### Load raw data from .rds file
```{r "load run"}
#load entire run
raw_data <- iso_read_continuous_flow(file.path("data", "raw", "run2_raw.cf.rds"))
```

### Parse file information
```{r "retrieve data"}
df.raw <- raw_data %>% 
  # aggregate peak data
  iso_get_vendor_data_table(include_file_info = c(file_path, file_datetime, Row, `Identifier 1`, `Identifier 2`, Analysis)) %>% 
  # select N2O peak (select_N2O_peak requires file field)
  rename(file = file_id) %>% 
  select_N2O_peak(c(310, 370)) %>% 
  # extract relevant file information
	 # focus on the columns we care about:
  rename(d45 = `d 45N2O/44N2O`, d46 = `d 46N2O/44N2O`, area = `Intensity All`) %>%
  mutate(folder = basename(dirname(file_path)), run_number = parse_integer(Row), analysis = as.numeric(sub(pattern = "MAT253", replacement = "", x = Analysis)),
         name = `Identifier 1`, volume = parse_number(`Identifier 2`)) %>% 
  # select relevant columns
  select(analysis, name,
         folder, date = file_datetime, run_number, volume,
         area, d45, d46) %>% 
  # extract categories from names
  mutate(category = extract_word(name, capture_n = 1, include_dash = TRUE),
    #add big picture categories
    is_standard = str_detect(category, "(USGS|IAEA)"),
    is_sample = category %in% c("OM15"),
    year = case_when(
    	category == "OM15" ~ "2015"),
    #samples with NR are nitrate only
    compound = "NO3", 
    volume = ifelse(category == "N2O"| category == "Conditioner", NA, volume/1000), #convert volume to mL
    id = sub("NR", "", name)) %>%
	#label conditioners
  change_category(run_number < 6, "Conditioner") %>%
  change_category(grepl("HPW", name), "Background") %>%
  change_category(grepl("ICDSTD", name), "control") %>% #select only relevant categories from metadata for downstream processing
  filter(category %in% c("Background", "control", "N2O", "Conditioner", "OM15", "OM16", "USGS-34","IAEA-NO3"))

df.raw
```

## Summary of the raw data
$\delta^{45}N_{2}O$ and $\delta^{46}N_{2}O$ in $‰$ are represented as d45 and d46 respectively in the data table. Raw d45 and d46 peak areas and standard deviations are reported. Two potassium nitrate isotopic standards (USGS-34 and IAEA) dissolved in HPW (high purity water) in two concentrations (30 $\mu$M and 250 $\mu$M) were spaced every 8-10 samples for calibration of  $\delta^{15}N$, $\delta^{18}O$, and nitrate concentration. Vials of 20 nmoles of injected N2O were also evenly spaced every 8 samples + standards to monitor drift.
```{r "raw data summary"}
#table of average and standard deviations of peak areas, d45, and d46 for standards and controls
df.raw %>% group_by(category, name) %>% 
 iso_summarize_data_table(area, d45, d46, cutoff = 3) 
```


## First look at standards
With the exception of run_number 91 and 102, standards look uniform in area and d45. Standards are plotted by subrun. Subrun 1 is samples that contain nitrate and potentially nitrite; subrun 2 contains samples and standards that underwent nitrite removal procedures.
```{r "plot standards", warning= FALSE, plotly=TRUE}
plot_raw_standards <- df.raw %>%
	#select just isotopic and N2O standards
	filter(category %in% c("IAEA-NO3", "USGS-34", "N2O")) %>%
  #plot with panels with ggplot  
  iso_plot_data(
  x = run_number, y = d45, points = TRUE,
  # shape = 21 with fill makes data points with black borders
  shape = 21, fill = name, size = area,
  panel = category ~ .
)

#change to colorblind friendly palette and add axis labels
plot_raw_standards + 
  	scale_fill_manual(values = cbPalette) + 
	  scale_color_manual(values = cbPalette) +
    ylab(latex2exp::TeX("$\\delta^{45}N_{2}O$")) +
    xlab("run number") 

#plot interactive version
# ggplotly(plot_raw_standards, dynamicTicks = TRUE)
```


## Corrections

### Drift correction
Evaluate drift throughout run duration. Linear drift is observed; correction applied.
```{r}
df.drift <- df.raw %>% evaluate_drift(d45, d46, correct = TRUE, method = "lm", correct_with = category %in% c("USGS-34", "IAEA-NO3", "N2O"), plot = TRUE)
#linear calibration applied; new columns d45.drift, d46.drift, p.drift created
```

### $^{17}$O correction
Correction for $^{17}$O isobaric interference in the mass 45 peak according to Kaiser et al. (2008) to derive raw $\delta^{15}N$ and $\delta^{18}O$ values.
```{r}
df.O17 <- df.drift %>%
  correct_N2O_for_17O(d45.drift, d46.drift) %>% #apply the 17O corrections to the drift-corrected d45.drift and d46.drift columns
	#corrected columns are d15.raw and 18.raw; p.17Ocor column added with correction parameters
  select_columns(-d45, -d45.drift, -d46, -d46.drift) 
```


## Concentration Calibration with Isoprocessor

### Load Standards
```{r}
#load in concentration and isotopic standards from excel file
stds_isotopes <- read_excel(file.path("data", "standards", "standards_run2.xlsx"), sheet = "isotopes")
stds_concs <- read_excel(file.path("data", "standards", "standards_run2.xlsx"), sheet = "concentrations")
kable(stds_isotopes)
kable(stds_concs)
```

### Generate concentration calibration
```{r "add standards", warning = FALSE}
df_calibp <- 
  df.O17 %>% 
  filter(category != "excluded") %>% #filter out excluded categories and conditioners
  filter(category != "Conditioner") %>%
  select(-p.17Ocor, -p.drift, -date) %>% #removed nested parameter columns
  # add standards
  iso_add_standards(stds_isotopes, match_by = category, is_std_peak = is_iso_std) %>% #add isotope standards
  iso_add_standards(stds_concs, match_by = name, is_std_peak = is_conc_std)  #add concentration standards
```

```{r "concentration calibration", warning=FALSE}
df_calib <- df_calibp %>% 
	 mutate(amount.nmol = true_conc.uM * volume) %>% #calculate nmoles from concentration
  iso_prepare_for_calibration() %>% 
  iso_generate_calibration(
    model = lm(area ~ amount.nmol - 1), #linear model where area is the response variable; the - 1 removes the implied intercept term
    calibration = "area",
    use_in_calib = is_conc_std) 

# check for problems
df_calib %>% iso_get_problematic_calibrations(calibration = "area")
```

### Inspect calibration
Residuals ~2% for most concentration standards
```{r warning = FALSE}
# parameter overview table
df_calib %>% 
  iso_get_calibration_parameters(calibration = "area")
```

```{r warning=FALSE}
# visualization of residuals
conc_resid <- df_calib %>% 
  iso_get_calibration_data() %>% 
  filter(is_conc_std) %>%
	filter(category != "excluded") %>% #filter for only standards used in calibration
  mutate(
    `Var: rel. residual area [%]` = 100*area_resid/area #calculate % relative residual area
  ) %>%
	ggplot(
		aes(x = run_number, y =`Var: rel. residual area [%]`, color = name)) +
			geom_point(size = 4) +
			theme_figure() +
	scale_color_manual(values = cbPalette) +
			xlab("Run Number") +
			ylab("Relative Residual Area (%)") + 
	guides(color = guide_legend(title = NULL)) 
	
ggplotly(conc_resid) 
```

### Apply calibration
```{r warning = FALSE}
# apply concentration calibration 
df_w_conc <- df_calib %>% 
  iso_apply_calibration(
    predict = amount.nmol,
    calibration = "area",
    calculate_error = TRUE, #standard error using Wald method
    predict_range = c(0, 20) #prediction range from standards
  ) %>% 
  iso_get_calibration_data()

df_w_conc
```

### Check data
A concentration control was not used in calibration, but is calculated to be 162.9 uM +/- 1.66 uM. The known value is 161 uM
```{r}
# visualize (note that most are technically out of range of the calibration which is only defined for the narrow signal observed in standards; the rest are extrapolations beyond the calibration model)
plot_nmoles <- df_w_conc %>%
	filter(category != "excluded") %>%
	ggplot(
	#define global plot aesthetics 
		aes(x = run_number, y = amount.nmol_pred, color = category, label = name)) +
	scale_color_manual(values = cbPalette) + 
	xlab("Run number") +
	ylab("Predicted nmoles of N") +
	theme_figure(grid = TRUE, legend = TRUE, text_size = 20) +
	guides(color = guide_legend(title = NULL))

#ggplot version with check or if calibration was ok
plot_nmoles_ggplot <- plot_nmoles + 
	geom_point(aes(shape = df_w_conc$area_calib_ok, fill = df_w_conc$category, color = df_w_conc$category)) + 
	scale_shape_manual(values = c(21, 22)) + 
	scale_fill_manual(values = cbPalette) + 
	theme_figure(grid = TRUE, legend = TRUE, text_size = 20) + 
	guides(shape = guide_legend(title = "Calibration OK"), fill = guide_legend(title = NULL))
plot_nmoles_ggplot

#interactive plot to see individual sample names from plot above	
ggplotly((plot_nmoles + geom_point()), tooltip=c("name", "amount.nmol_pred", "run_number"))  
```

```{r "exclude samples with <5 nmoles N"}
#exclude samples that have <5 nmoles predicted N
df_w_conc <- df_w_conc %>%
	filter(
		if_else(is_sample == TRUE, amount.nmol_pred > 5, amount.nmol_pred > 0)
	)
df_w_conc
```

```{r "calculate concentration"}
# calculate concentration
df_w_conc <- 
  df_w_conc %>% 
  mutate(conc_uM =  amount.nmol_pred / volume,
         conc_uM_err = amount.nmol_pred_se / volume)

# check amount and concentration data
df_w_conc %>% filter(category != "N2O") %>%
  group_by(category, true_conc.uM) %>% 
  iso_summarize_data_table(
    conc_uM,
    amount.nmol_pred, 
    amount.nmol_pred_se, 
    cutoff = 3) 
```


## Delta calibration with Isoprocessor

### Generate calibration for N
```{r warning = FALSE}
df_calibN <- df_w_conc %>% 
  iso_prepare_for_calibration() %>% 
  iso_generate_calibration(
    model = c(
      # generate calibration using different regression models 
      Nlinear = lm(d15.raw ~ true_d15N), #simple linear
      Nwith_conc = lm(d15.raw ~ true_d15N + conc_uM), #multiple regression with concentration
      Nwith_volume = lm(d15.raw ~ true_d15N + volume) #multiple regression with volume
    ),
    calibration = "d15N",
    use_in_calib = is_iso_std) 
```

### Inspect calibration for N
```{r}
# parameter overview table
df_calibN %>% 
  iso_get_calibration_parameters(calibration = "d15N")

#visualization of residuals for each N isotope calibration
N_residuals <- df_calibN %>%
	iso_get_calibration_data() %>%
	filter(is_iso_std) %>%
	filter(category != "excluded") %>% #filter for only standards used in calibration
	ggplot(
		aes(x = run_number, y = d15N_resid, color = d15N_calib)) +
			geom_point(size = 4) +
			theme_figure() +
	scale_color_manual(values = cbPalette) +
			xlab("Run Number") +
			ylab("Residuals") + 
	guides(color = guide_legend(title = NULL)) 
	
ggplotly(N_residuals) 

#multiple regressions using concentration or volume are not significantly better than the simple linear calibration model
```

### Generate calibration for O
```{r warning = FALSE}
df_calibO <- df_w_conc %>% 
  iso_prepare_for_calibration() %>% 
  iso_generate_calibration(
    model = c(
      # generate calibration using different regression models
      Olinear = lm(d18.raw ~ true_d18O), #linear
      Owith_conc = lm(d18.raw ~ true_d18O + conc_uM), #multiple regression with concentration
      Owith_volume = lm(d18.raw ~ true_d18O + volume) #multiple regression with volume
    ),
    calibration = "d18O",
    use_in_calib = is_iso_std) 
```

### Inspect calibration for O
```{r}
# parameter overview table
df_calibO %>% 
  iso_get_calibration_parameters(calibration = "d18O")
  
#visualization of residuals for each N isotope calibration
O_residuals <- df_calibO %>%
	iso_get_calibration_data() %>%
	filter(is_iso_std) %>%
	filter(category != "excluded") %>% #filter for only standards used in calibration
	ggplot(
		aes(x = run_number, y = d18O_resid, color = d18O_calib)) +
			geom_point(size = 4) +
			theme_figure() +
	scale_color_manual(values = cbPalette) +
			xlab("Run Number") +
			ylab("Residuals") + 
	guides(color = guide_legend(title = NULL)) 
	
ggplotly(O_residuals) 

#multiple regression using concentration is better than the simple linear calibration model# visualization of residuals
```

### Apply calibration
Applying a multiple regression calibration with volume for N and O
```{r warning = FALSE}
# apply linear calibration for N 
df_finalN <- df_calibN %>%
filter(d15N_calib == "Nlinear")%>% 
  iso_apply_calibration(
    predict = true_d15N,
    calibration = "d15N",
    calculate_error = TRUE
  ) %>%
  iso_get_calibration_data() %>% iso_remove_list_columns()
 
  
# apply multiple regression calibration with volume for O
df_finalO <- df_calibO %>% 
  filter(d18O_calib == "Owith_volume") %>% 
  iso_apply_calibration(
    predict = true_d18O,
    calibration = "d18O",
    calculate_error = TRUE
  ) %>%
 iso_get_calibration_data() %>%
  iso_remove_list_columns()

#join O and N calibrated data into final dataframe, df_final
 df_final <- left_join(df_finalN, df_finalO)
 
 df_final
```

### Check data
```{r}
# summary of calibrated isotopic data for major categories
df_final %>% 
	filter(category != "N2O") %>%
  group_by(category, true_d15N, true_d18O) %>% 
  iso_summarize_data_table(
    d15N = true_d15N_pred, 
    d15N_se = true_d15N_pred_se,
    d18O = true_d18O_pred, 
    d18O_se = true_d18O_pred_se,
    cutoff = 3) 
```

### Visualize Calibrated Nitrogen Isotopic Data
```{r "Visualize Calibrated Nitrogen Isotopic Data"}
#visualization of calibrated N isotopic data
N_plot <- df_final %>%
	filter(category != "N2O") %>%
	filter(category != "excluded") %>% #remove excluded samples and standards from plot
	ggplot(
		aes(x = run_number, y = true_d15N_pred, color = category, label = name, shape = d15N_calib_ok)) +
			geom_point(size = 4) +
		  geom_errorbar(aes(ymax = true_d15N_pred + true_d15N_pred_se, ymin = true_d15N_pred - true_d15N_pred_se)) +
			theme_figure() +
	scale_color_manual(values = cbPalette) +
			xlab("Run Number") +
	guides(color = guide_legend(title = NULL), shape = guide_legend(title = NULL)) 

#ggplot version with check for if calibration was ok
plot_Nisotopes_ggplot <- N_plot + 
	scale_fill_manual(values = cbPalette) + 
	theme_figure(grid = TRUE, legend = TRUE, text_size = 20) + 
	ylab(latex2exp::TeX("$\\delta^{15}N$")) +
	guides(shape = guide_legend(title = "Calibration OK"), fill = guide_legend(title = NULL), color = guide_legend(title = NULL))

plot_Nisotopes_ggplot
```

### Visualize Calibrated Oxygen Isotopic Data
```{r "Visualize Calibrated Oxygen Isotopic Data"}
#visualization of calibrated O isotopic data
O_plot <- df_final %>%
		filter(category != "N2O") %>% 
	filter(category != "excluded") %>% #remove excluded samples and standards from plot
	ggplot(
		aes(x = run_number, y = true_d18O_pred, color = category, label = name, shape = d18O_calib_ok)) +
			geom_point(size = 4) +
		  geom_errorbar(aes(ymax = true_d18O_pred + true_d18O_pred_se, ymin = true_d18O_pred - true_d18O_pred_se)) +
			theme_figure() +
	scale_color_manual(values = cbPalette) +
			xlab("Run Number")  +
guides(color = guide_legend(title = NULL), shape = guide_legend(title = NULL)) 
 

#ggplot version with check for if the calibration was ok
plot_Oisotopes_ggplot <- O_plot + 
	scale_fill_manual(values = cbPalette) + 
	theme_figure(grid = TRUE, legend = TRUE, text_size = 20) + 
	ylab(latex2exp::TeX("$\\delta^{18}O$")) +
	guides(shape = guide_legend(title = "Calibration OK"), fill = guide_legend(title = NULL), color = guide_legend(title = NULL))

plot_Oisotopes_ggplot
```

The first part of this run was direct injections of sampling volume into denitrifier vials since the resuspension medium was well buffered, no issues were expected with the low pH for low injections. The second part of the run was combined with base addition to neutralize acidified samples. However, the volume dependent discrepancy (higher *measured* concentration with base addition) between the first and second part of the run suggests base contamination with residual nitrate and only the first part of the run is considered in downstream data processing.
```{r}
df_final <- df_final %>% filter(analysis < 126215)
```


## Summary of calibrated data
```{r}
df_final <- df_final %>%
  mutate(
    conc = conc_uM,
    conc_err = conc_uM_err,
    d15N = true_d15N_pred,
    d15N_err = true_d15N_pred_se,
    d18O = true_d18O_pred,
    d18O_err = true_d18O_pred_se
  )  

df_final %>% group_by(category, name) %>% 
  iso_summarize_data_table(cutoff = 1, d15N, d15N_err, d18O, d18O_err, conc, conc_err)
```


## Export Samples Data
```{r, include=FALSE}
df_final %>%
  filter(category %in% c("OM15")) %>%
  select(analysis, run_number, volume, area, name, category, conc, conc_err, d15N, d15N_err, d18O, d18O_err, compound) %>%
write.csv(file = file.path("data", "calibrated", "run2_calib_data_export.csv"))

df_final %>% 
 filter(category %in% c("OM15")) %>% 
  saveRDS(file = file.path("data", "calibrated", "df_final_run2"))
# saveRDS
``` 

## Export standards data
```{r, include=FALSE}
df_final %>%
  filter(category %in% c("USGS-34", "IAEA-NO3")) %>%
  select(analysis, run_number, volume, area, name, category, conc, conc_err, d15N, d15N_err, d18O, d18O_err, compound) %>%
write.csv(file = file.path("data", "calibrated", "run2_stds_calib_data_export.csv"))

df_final %>% 
 filter(category %in% c("USGS-34", "IAEA-NO3")) %>%
  saveRDS(file = file.path("data", "calibrated", "df_stds_run2"))
# saveRDS
```