01B_LH-SIP_Data_Reduction_Correction.Rmd

---
title: "01_LH-SIP_Data_Reduction"
author: "Tristan Caro"
date: "`r Sys.Date()`"
output:
    html_document:
    code_folding: show
    df_print: paged
    number_sections: yes
    toc: yes
    toc_depth: 3
    toc_float: yes
editor_options:
  chunk_output_type: console
---

# Setup

## Clear the environment

```{r}
rm(list=ls())
```


```{r setup, include = FALSE}
# global knitting options for code rendering
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>")

# global knitting options for automatic saving of all plots as .png and .pdf
knitr::opts_chunk$set(
  dev = c("png", "pdf"),
  dev.args = list(pdf = list(encoding = "WinAnsi", useDingbats = FALSE)),
  fig.keep = "all",
  fig.path = file.path("fig_output", paste0(gsub("\\.[Rr]md", "", knitr::current_input()), "_"))
)
```

## Load packages

```{r, message=FALSE, warning=FALSE}
# Libraries
library(tidyverse)    # CRAN v1.3.1
library(forcats)      # CRAN v0.5.1 
library(isoreader)    # [::NA/NA] v1.3.0 # CRAN v1.3.0
library(isoprocessor) # [github::isoverse/isoprocessor] v0.6.7
library(isotopia)     # [github::isoverse/isotopia] v0.5.8 
library(readxl)       # CRAN v1.3.1 
library(ggsci)        # CRAN v2.9 
library(ggborderline) # [github::wurli/ggborderline] v0.1.0
library(ggrepel)      # CRAN v0.9.1 
library(latex2exp)    # CRAN v0.5.0
library(ggdist)       # CRAN v3.0.0 
library(ggsignif)     # CRAN v0.6.2
library(ggridges)     # CRAN v0.5.3
library(cowplot)      # CRAN v1.1.1
library(isotopia)
```

## Load sourced functions
These functions are found in the `libs/` directory.
```{r}
# Sourced Functions
source(file.path("libs", "visualization.R"))              # Visualization scripts
source(file.path("libs", "chromeleon.R"))                 # Chromeleon Reader
source(file.path("libs", "error_prop.R"))                 # Error Propagation
source(file.path("libs", "calculate_turnover.R"))         # Turnover Calculator
source(file.path("libs", "d2H_to_2F.R"))                  # Convert delta to frac abund
source(file.path("libs", "rename_FAs.R"))                 # Rename FAs to correctly ID'd FAs
var_to_str = function(v) {return(deparse(substitute(v)))} # Convert variable name to string
```




# Load data

## Load lipid d2H assimilation efficiency data


```{r}
d2H_assim <- read_xlsx("data/d2H_data_isotope_measurements.xlsx", sheet = 'lipids')
```


## Load F_L 2F Data

```{r}
F_L_2F <- read_csv("data/Label_2H_Data/dD_cleaned.csv") %>% 
  group_by(soil) %>% 
  summarize(F_label_ppm_mn = mean(F_label_ppm),
            F_label_ppm_se = sd(F_label_ppm)) %>% 
  mutate(F_label_mn = F_label_ppm_mn / 10000,
         F_label_se = F_label_ppm_se / 10000) %>% 
  mutate(terrain = case_when(
    soil == "Conifer" ~ "Gordon Gulch Conifer Forest",
    soil == "Tundra" ~ "Niwot Ridge Tundra",
    soil == "Grassland" ~ "Marshall Mesa Grassland"
  ))

write_csv(F_L_2F, file = "data/Label_2H_Data/dD_summarized.csv")
```

## Load Compound nC and nH Data sheet

```{r}
compound_nC_nH <- read_xlsx("data/compound_nC_nH.xlsx")
```

## Load Phthalic Acid IRMS Data

> NOTE: For this analysis, we use the corrected PAME dataset: `data/PAME_TAC_jan_2022.xlsx`. The original script `01_LH-SIP_Data_Reduction.rmd` is retained for reference.

```{r}
PAME_summarized <- readxl::read_excel("data/PAME_TAC_jan_2022.xlsx") |> 
  mutate(
    # known isotopic composition of the 4 non-exchangeable Hs of phthalic acid
    pa_d2H = -81.9,
    pa_d2H_err = 1.2,
    # mass balance calculation of methanol methyl group H d2H
    meoh_d2H = 10/6 * pame_mean - 4/6 * pa_d2H,
    meoh_d2H_err = sqrt((10/6 * pame_sd)^2 + (4/6 * pa_d2H_err)^2),
    # convert delta d2H to atom percent at2H
    at2H_mn = d2H_to_2F(meoh_d2H)
  )

# uncomment this to write to drive:
# writexl::write_xlsx(PAME_summarized, path = "data/PAME_TAC_jan_2022.xlsx")
write_rds(PAME_summarized, file = "cache/PAME_summarized.rds")
```



## Load IRMS Data

### Load Metadata and Problematic Run Data

```{r}
mtda <- read_csv(file.path("data", "sample_metadata_2021.csv"))
sample_names <- mtda$id1

problematic_runs <- read_xlsx(file.path("data", "problematic_runs.xlsx"))
problematic_analyses <- problematic_runs$Analysis
problematic_analyses_num <- parse_number(problematic_analyses)
```

### Load Peak Map

```{r}
# this information is often maintained in a csv or Excel file
peak_map <- 
  # initial peak map:
  # readxl::read_excel(file.path("data","IRMS", "peak_map_manual_general.xlsx"))
  # corrected peak map
  readxl::read_excel(file.path("data", "peak_maps_resolved.xlsx"))
```

### Read raw IRMS data files

```{r, warning=FALSE}
# Set file path(s) to data files, folders or rds collections 
# can be multiple folders or mix of folders and files
# Isoverse will iteratively search subfolders. Huzzah! ^_^
data_path <- file.path("data", "IRMS", "raw_data")
ref_ratio <- get_standard("2H") %>% as.numeric()

# read files
iso_files_raw <- 
  # path to data files
  data_path %>% 
  # read data files in parallel for fast read
  iso_read_continuous_flow() %>%
  # filter out files with read errors (e.g. from aborted analysis)
  iso_filter_files_with_problems()
```

### Process file info & peak table

```{r}
# process IRMS file information
iso_files <- iso_files_raw %>% 
  # rename key file info columns
  iso_rename_file_info(analysis = Analysis, id1 = `Identifier 1`, id2 = `Identifier 2`) %>% 
  # parse text info into numbers
  iso_parse_file_info(number = analysis) %>% 
  # process other file information that is specific to the naming conventions
  # of this particular sequence
  iso_mutate_file_info(
    # what is the type of each analysis?
    type = case_when(
      str_detect(id1, "[Zz]ero")      ~ "on_off",
      str_detect(id1, "H3")           ~ "H3_factor",
      str_detect(id1, "F8")           ~ "F8_std",
      str_detect(id1, "F9")           ~ "F9_std",
      TRUE                            ~ "sample"
    ),
    # what was the concentration? (assuming Preparation = concentration or volume)
    concentration = 
      ifelse(type == "std",  
             str_extract(Preparation, "[0-9.]+ ?ng( per |/)uL") %>% 
               parse_number() %>% iso_double_with_units("ng/uL"),
             NA),
    # what folder are the data files in? (assuming folder = sequence)
    folder = basename(dirname(file_path))
  ) %>% 
  # focus only on the relevant file info, discarding the rest
  iso_select_file_info(
    folder, analysis, file_datetime, id1, type, concentration
  ) %>% 
  # add in additional sample metadata (could be any info)
  # note: this would typically be stored in / read from a csv or excel file
  iso_add_file_info(
    read_csv(file.path("data", "sample_metadata_2021.csv")),
    join_by = "id1"
  )

# set peak table from vendor data table with default isodat template
iso_files <- iso_set_peak_table_from_isodat_vendor_data_table(iso_files) %>% 
  # convert units from mV to V for amplitudes and area
  iso_convert_peak_table_units(V = mV, Vs = mVs)

# focus on sample files
sample_files <- iso_filter_files(iso_files, type == "sample")

# Include only sample files
sample_files <- iso_filter_files(sample_files, id1 %in% sample_names)

# EXCLUDE CRC SAMPLES FROM THIS ANALYSIS
sample_files <- iso_filter_files(sample_files, str_detect(id1, "cri", negate = TRUE))

sample_files <- iso_filter_files(sample_files, str_detect(id1, "cr0", negate = TRUE))

# EXCLUDE RUN WE'VE DEEMED PROBLEMATIC
sample_files <- iso_filter_files(sample_files, !analysis %in% problematic_analyses_num)



```

Show file information

```{r}
# display file information
iso_files %>% 
  iso_get_file_info() %>% select(-file_id, -folder) %>% 
  iso_make_units_explicit() %>% knitr::kable()
```


### Example chromatogram

```{r "example_chromatograms", fig.width=8, fig.height=8}
# plot an example chromatogram

# sample_files[10] %>% # choosing arbitrary file to plot
# iso_plot_continuous_flow_data(
#   # select data and aesthetics
#   data = c(2), color = id1, panel = id1,
#   # zoom in on time interval
#   time_interval = c(750, 4000),
#   # peak labels for all peaks > 2V
#   peak_bounds = TRUE,
#   peak_marker = FALSE,
#   peak_label = iso_format(rt),
#   #peak_label_size = 3,
#   peak_label_filter = analysis == 5685 & amp2 > 1
# ) + scale_color_npg() + 
#   theme_classic() +
#   theme(strip.background = element_blank())
```

## Load memory-corrected IRMS data!

```{r}
LH_SIP_memcorr <- readxl::read_excel("data/samples_data_memory_corrected_20221228.xlsx") %>% 
# Trim it down
  select(-c(type, # or `standard` if using original mc dataset 
            #is_std_peak,
            folder, 
            file_datetime, 
            id2, 
            #individual,
            #injection_volume,
            #gc_method,
            #map_id, 
            gc_ramp, 
            group, 
            #peak_type,
            #ref_nr,
            #calib_peak, 
            peak_nr,
            #is_ref,
            rt_start,
            rt, 
            rt_end, 
            n_overlapping,
            n_matches,
            #use_in_calib,
            #d2H_in_calib, 
            #d2H_resid, 
            d2H_calib_points
            )
         ) %>% 
  filter(!id1 %in% c("cr0f0x", "cr0f0y", "cr0f0z",
                     "cr0f3x", "cr0f3y", "cr0f3z",
                     "cr0f7x", "cr0f7y", "cr0f7z",
                     "cri0", "cri3", "cri7")) %>% 
  rename(sample_id = id1,
         Analysis = analysis) %>% 
  select(Analysis, sample_id, compound, area2, area3, note,
          below_area_range, 
          calibrated_d2H_without_area,
          calibrated_d2H_without_area_se,
          calibrated_d2H_without_memory, 
          calibrated_d2H_without_memory_se,
          calibrated_d2H_with_memory_correction, 
          calibrated_d2H_with_memory_correction_se) %>% 
  mutate(Analysis = as.numeric(Analysis)) %>% 
  # Convert Delta Values in Permil (vs. VSMOW) to fractional abundance (at%)
  mutate(
    # Without area correction:
    calibrated_at2H_without_area = d2H_to_2F(calibrated_d2H_without_area),
    calibrated_at2H_without_area_se = d2H_to_2F(calibrated_d2H_without_area_se),
    # With Area, without memory correction:
    calibrated_at2H_without_memory = d2H_to_2F(calibrated_d2H_without_memory),
    calibrated_at2H_without_memory_se = d2H_to_2F(calibrated_d2H_without_memory_se),
    # With memory correction:
    calibrated_at2H_with_memory_correction = d2H_to_2F(calibrated_d2H_with_memory_correction), 
    calibrated_at2H_with_memory_correction_se = d2H_to_2F(calibrated_d2H_with_memory_correction_se)
  )
```

# Data Reduction

### Reduce data and get linear model

```{r}
# Find Fatty acid analytes manually

fas_in_d2H_assim <- d2H_assim %>% 
  filter(str_detect(Analyte, "acid") |
           str_detect(Analyte, "ate")) %>% 
  pull(Analyte) %>% 
  unique()

d2H_assim_lm <- d2H_assim %>% 
  # Select columns of interest
  select(
    org_id,
    exp_id,
    Analyte,
    `Water dD`,
    `Lipid dD`
  ) %>% 
  # Filter to only look at FAs
  filter(Analyte %in% fas_in_d2H_assim) %>% 
  drop_na() %>% 
  # Group by organism, experiment, analyte
  group_by(
    org_id,
    exp_id,
    Analyte
  ) %>% 
  # Filter out samples with three or fewer observations
  # (so that we can do a linear model)
  filter(n() > 3) %>% 
  tidyr::nest(
    data = c(`Water dD`, `Lipid dD`)
  ) %>% 
  mutate(
    fit = purrr::map(data, ~lm(`Lipid dD` ~ `Water dD`, data = .x)),
    estimates = purrr::map(fit, broom::tidy),
    summary = purrr::map(fit, broom::glance)
  ) %>% tidyr::unnest(estimates)

# Simplify by only looking at the slopes
# This is what gets added to the peaks_mapped dataframes!
d2H_assim_lm_summary <- d2H_assim_lm %>% 
  ungroup() %>% 
  filter(term == "`Water dD`") %>% 
  filter(estimate >= 0,
         estimate <= 1) %>% 
  summarize(sa = sd(estimate, na.rm = TRUE),
            a_mn = mean(estimate, na.rm = TRUE))
```


## IRMS Peak Mapping

```{r}
sample_peak_table <- sample_files %>%
  iso_set_peak_table_from_isodat_vendor_data_table() %>%
  iso_get_peak_table() %>% 
  mutate(Analysis = substr(file_id, 1, 7),
         sample_id = substr(file_id, 10,14))
```

### Generate Isoverse peak map
```{r}
# Generate an isoverse peak map
peaks_mapped <- sample_peak_table %>%
  iso_map_peaks(peak_map, map_id = Analysis) %>%
  filter(!is.na(compound))
```

### Remove readychecks, blanks, etc. from analysis
```{r}
# Filter based on what is in the sample_names metadata
peaks_mapped <- peaks_mapped %>% 
  filter(sample_id %in% sample_names)

# Filter out runs that are on our problematic_runs list
peaks_mapped <- peaks_mapped %>% 
  filter(!Analysis %in% problematic_runs)
```

### Check problematic peak assignments
```{r}
problem_peaks <- iso_get_problematic_peak_mappings(peaks_mapped)
```

### Export summary of unresolved peaks
```{r}
# Spit out peaks that need to be resolved
peaks_mapped %>%
  ungroup() %>%
  select(Analysis, compound, rt) %>%
  mutate(rt = round(rt, digits=0)) %>%
  pivot_wider(
    names_from = Analysis, 
    names_prefix = "rt:",
    values_from = rt, 
    id_cols = compound, 
    values_fn = function(x) paste(x, collapse = "; ")
  ) %>%
  openxlsx::write.xlsx(file.path("data", "IRMS", "peak_maps_resolve.xlsx"))
```

## Metadata and cleanup

### Assign incubation params
```{r}
# Add Some Parameters
peaks_mapped <- peaks_mapped %>% 
  group_by(sample_id) %>% 
  mutate(inc_time_d = as.numeric(substr(sample_id, 4,4)),
         inc_time_d_str = paste(inc_time_d, "Days"),
         f_label = 0.005,
         t_series_id = paste0(substr(sample_id, 1,1), 
                              substr(sample_id, 5,5)),
         soil_id = substr(t_series_id, 1,1),
         replicate_id = substr(t_series_id, 2,2),
         terrain = case_when(
           substr(sample_id, 1, 1) == "t" ~ "Niwot Ridge Tundra",
           substr(sample_id, 1 ,1) == "c" ~ "Gordon Gulch Conifer Forest",
           substr(sample_id, 1 ,1) == "m" ~ "Gordon Gulch Meadow",
           substr(sample_id, 1 ,1) == "g" ~ "Marshall Mesa Grassland"))
# 
# inc_time_zero <- peaks_mapped %>% 
#   filter(inc_time_d == 0) %>% 
#   select(sample_id, compound, at2H, t_series_id) %>% 
#   mutate(f_start = at2H) %>% 
#   select(-at2H)
# 
#   

  
```

### Join memory-corrected data

```{r}
peaks_mapped_mc <- peaks_mapped %>% 
# Need to remove "BF" from Analysis number
  mutate(Analysis = parse_number(Analysis)) %>%
  select(-c(area2, area3)) %>% 
  # Join
  left_join(
    LH_SIP_memcorr, 
    by = c("Analysis", "sample_id", "compound"))
```


### Extract zero timepoint values
These SIP measurements compare an enriched isotopic value to a starting isotopic composition. To have a time `t = 0` value, we need to define the starting point of our incubation.

```{r}
zero_peaks <- peaks_mapped_mc %>%
  filter(inc_time_d_str == "0 Days") %>%
  group_by(t_series_id, compound) %>%
  mutate(n_analyical_reps = n()) %>%
  ungroup() %>%
  # exclude ambiguous peaks
  filter(!is_ambiguous)

zero_peaks_averaged <- zero_peaks %>% 
  group_by(soil_id, compound) %>% 
  summarise(
    # at2H at zero (no mc)
    at2H_mn_zero = mean(calibrated_at2H_without_memory,
                        na.rm = TRUE),
    # Memory corrected value
    at2H_mn_zero_mc = mean(calibrated_at2H_with_memory_correction, 
                           na.rm = TRUE),
    # d2H at zero (no mc)
    d2H_mn_zero = mean(calibrated_d2H_without_memory, 
                       na.rm = TRUE),
    # Memory corrected value
    d2H_mn_zero_mc = mean(calibrated_d2H_with_memory_correction, 
                          na.rm = TRUE))
  
# Cache the zero-timepoint dataframe
zero_peaks_averaged %>% 
  openxlsx::write.xlsx(
    file.path("cache", 
              "zero_peaks_averaged.xlsx"),
    overwrite = TRUE)


# This is the dataset we will use for generation time calculations
peaks_mapped_with_zeros <- 
  peaks_mapped_mc %>%
  left_join(
    zero_peaks_averaged,
    by = c("soil_id", "compound")
  )
stopifnot(nrow(peaks_mapped_with_zeros) == nrow(peaks_mapped_mc))


peaks_mapped_with_zeros <- peaks_mapped_with_zeros %>% 
  mutate(sample_id = substr(file_id, 10,14))
```


## Apply corrections

### Correct for PAME

Correcting for the isotopic composition of the derivatization agent by measuring phthalic acid methyl ester (PAME).

```{r}

peaks_mapped_with_zeros_FID_mc_pame <- peaks_mapped_with_zeros %>% 
  left_join(compound_nC_nH, by = "compound") %>% 
  mutate(`2F_Me` = PAME_summarized %>% pull(at2H_mn),
         # define calibrated 2F
         `2F_FAME` = calibrated_at2H_without_memory,
         `2F_alk` = (`2F_FAME` - `2F_Me`*x_me) / x_alk,
         # define calibrated 2F mc
         `2F_FAME_mc` = calibrated_at2H_with_memory_correction,
         `2F_alk_mc` = (`2F_FAME_mc` - `2F_Me`*x_me) / x_alk)
```

### Separate analytes and standards, calculate F_0 standard error

```{r}
# For 0-3,0-7 timepoints
# Analytes
peaks_mapped_with_zeros_analytes <- peaks_mapped_with_zeros_FID_mc_pame %>% 
  # Don't want our standard compounds to be used for generation time calculations! x_x
  filter(str_detect(compound, "STD", negate = TRUE)) %>% 
  # Also don't want our yeast-extract stimulated condition in our analyte pool
  filter(str_detect(sample_id, "t0t", negate = TRUE)) %>% 
  # Add in label strengths
  left_join(F_L_2F, by = "terrain")

# Standards
peaks_mapped_with_zeros_standards<- peaks_mapped_with_zeros_FID_mc_pame %>% 
  filter(str_detect(compound, "STD"))

# 23:0 and 21:0 PC standard errors at inc time zero
# Determines our standard error in FAME measurements at F_0
sF_0 <- peaks_mapped_with_zeros_standards %>% 
  # Select only 21 and 23 PC standards at initial time point
  filter(compound %in% c("21:0 (STD)"),
         inc_time_d == 0) %>%
  # Calculate mean and se of all 21:0 and 23:0 standard at2H values
  ungroup() %>% 
  summarize(F_0_at2H_mn = mean(at2H),
            F_0_at2H_se = sd(at2H))


# Add sF_0 to analyte data

peaks_mapped_with_zeros_analytes <- peaks_mapped_with_zeros_analytes %>% 
  mutate(F_0_at2H_se = sF_0 %>% pull(F_0_at2H_se))

```

### Eliminate extraneous columns

```{r}
# Reduce the 0d --> 3d --> 7d incubation data
peaks_mapped_with_zeros_analytes_reduced <- peaks_mapped_with_zeros_analytes %>%
  # Add assimilation efficiency data
  mutate(
    a = d2H_assim_lm_summary$a_mn,
    sa = d2H_assim_lm_summary$sa
  ) %>% 
  # REPAIR out-of-area t=0 datapoints:
  mutate(
    calibrated_at2H_with_memory_correction =
      case_when(
        below_area_range == TRUE & inc_time_d == 0 ~ at2H_mn_zero,
        TRUE ~ calibrated_at2H_with_memory_correction
        ),
    below_area_range = 
      case_when(
        below_area_range == TRUE & inc_time_d == 0 ~ FALSE,
        below_area_range == TRUE ~ TRUE,
        below_area_range == FALSE ~ FALSE
      )
  ) %>% 
  # REMOVE peaks that are outside area range of accurate correction!
  filter(below_area_range == FALSE) %>% 
  # Column reduction
  select(Analysis,
         sample_id,
         terrain,
         soil,
         inc_time_d,
         inc_time_d_str,
         compound,
         #area_FID,
         #rel_area_FID,
         rt,
         d2H_mn_zero,
         d2H,
         at2H_mn_zero,
         at2H_mn_zero_mc, # MC
         at2H, # not-corrected value
         calibrated_at2H_without_area, # F_T, no area corr
         calibrated_at2H_without_memory, # F_T, no memcorr
         calibrated_at2H_with_memory_correction, # F_T, MC!
         area2,
         area3,
         #area_FID,
         #compound_ug_per_g_soil,
         a, # assimilation efficiency
         `2F_Me`,
         `2F_FAME`,
         `2F_FAME_mc`, # MC
         `2F_alk`,
         `2F_alk_mc`,  # MC
         F_label_mn,
         # Uncertainty terms:
         F_label_se, # error in the label strength (pipetting error)
         F_0_at2H_se, # error in the initial F_0 (instrument error)
         calibrated_at2H_without_area_se, # error in F_T, no area corr
         calibrated_at2H_without_memory_se, # error in F_T, no memcorr
         calibrated_at2H_with_memory_correction_se, # error in F_T, MC!
         sa, # error in assimilation efficiency (se of literature values)
  ) %>% 
  # Turn all unit dbls into unitless
  mutate(
    across(
      where(is.numeric),
      as.numeric
    )
  ) %>% 
  # Rename FAs to match correctly ID'd analytes
  rename_FAs()

```

```{r}
LH_SIP_memcorr <- LH_SIP_memcorr %>% rename_FAs()
```



# Export
```{r}
peaks_mapped_with_zeros_analytes_reduced |> writexl::write_xlsx(path = "data/LH_SIP_PAME_CORR.xlsx")
```