# load packages
# libraries
library(isoorbi) # load isoorbi R package
library(dplyr) # for mutating data frames
library(stringr) # for extarcting text information
library(ggplot2) # for data visualization
library(cowplot) # arrange multipanel plotsFlow injection
Setup
Using R version 4.5.0 (2025-04-11) , tidyverse version 2.0.0, and isoorbi version 1.5.1.
Data
Optimal tuning
flow_injections_optimal <-
"data/nitrate/optimal_tuning" |>
orbi_find_isox() |>
orbi_read_isox() |>
orbi_simplify_isox() |>
# get injection number and sample name from filename
mutate(
injection = filename |> str_extract("\\d+$") |> as.integer(),
sample_name = ifelse(injection %% 2 == 0, "USGS35", "USGS32")
) |>
# focus on injections with signal
filter(injection > 1, injection < 15) |>
# define data block
orbi_define_block_for_flow_injection(start_time.min = 0.5, end_time.min = 6.5) |>
# flag extreme TICxIT values
orbi_flag_outliers(agc_fold_cutoff = 2)Bad tuning
flow_injections_bad <-
"data/nitrate/bad_tuning" |>
orbi_find_isox() |>
orbi_read_isox() |>
orbi_simplify_isox() |>
# get injection number and sample name from filename
mutate(
injection = filename |> str_extract("\\d+$") |> as.integer(),
sample_name = ifelse(injection %% 2 == 0, "USGS35", "USGS32")
) |>
# focus on injections with signal
filter(injection > 1, injection < 15) |>
# define data block
orbi_define_block_for_flow_injection(start_time.min = 0.5, end_time.min = 6.5) |>
# flag extreme TICxIT values
orbi_flag_outliers(agc_fold_cutoff = 2)Calculations
Ratios
# basepeak & ratios
flow_injections_optimal_w_basepeak <-
flow_injections_optimal |>
orbi_define_basepeak("M0")
flow_injections_bad_w_basepeak <-
flow_injections_bad |>
orbi_define_basepeak("M0")
# ratio summaries
flow_injections_optimal_summary <-
flow_injections_optimal_w_basepeak |>
orbi_summarize_results(ratio_method = "sum") |>
mutate(is_ref = sample_name == "USGS35")
flow_injections_bad_summary <-
flow_injections_bad_w_basepeak |>
orbi_summarize_results(ratio_method = "sum") |>
mutate(is_ref = sample_name == "USGS35")Deltas
# known ratios vs Air/VSMOWM
known_ratios <- tibble(
isotopocule = c("15N", "17O", "18O"),
reference = c("Air-N2", "VSMOW", "VSMOW"),
ratio_known = c(0.003676, 0.001140, 0.006016)
)
delta_USGS35 <- tibble(
isotopocule = c("15N", "17O", "18O"),
delta_known = c(+2.7, +50.35, +57.5)
)
delta_USGS32 <- tibble(
isotopocule = c("15N", "17O", "18O"),
delta_known = c(+180.0, +13.2, +25.7)
)
# delta ratios
ratios <-
bind_rows(
flow_injections_optimal_summary |> mutate(tuning = "optimal"),
flow_injections_bad_summary |> mutate(tuning = "bad")
) |>
select(tuning, is_ref, injection, isotopocule, basepeak, ratio, ratio_sem)
# calculate deltas by bracketing
deltas <-
ratios |>
# bracketing
filter(!is_ref) |>
left_join(
ratios |>
filter(is_ref) |>
rename(
ref_before = injection, ref_before_ratio = ratio,
ref_before_ratio_sem = ratio_sem
) |> select(-is_ref),
by = c("isotopocule", "tuning"),
relationship = "many-to-many") |>
filter(injection == ref_before + 1L) |>
left_join(
ratios |>
filter(is_ref) |>
rename(
ref_after = injection, ref_after_ratio = ratio,
ref_after_ratio_sem = ratio_sem
) |> select(-is_ref),
by = c("isotopocule", "tuning"),
relationship = "many-to-many") |>
filter(injection == ref_after - 1L) |>
mutate(sample_name = "USGS32") |>
# delta calculations
left_join(rename(delta_USGS35, delta_ref = delta_known), by = "isotopocule") |>
mutate(
# average the bracketing standard
ref_ratio = 0.5 * (ref_before_ratio + ref_after_ratio),
ref_ratio_sem = 0.5 * sqrt(ref_before_ratio_sem^2 + ref_after_ratio_sem^2),
# calculate deltas
delta = ratio/ref_ratio * (delta_ref + 1000) - 1000,
delta_sem = (delta + 1000) * sqrt( (ratio_sem/ratio)^2 + (ref_ratio_sem/ref_ratio)^2)
) |>
select(-is_ref) |>
arrange(isotopocule, injection)
deltas |> orbi_export_data_to_excel("output/nitrate_data_all.xlsx")
# deltas summary
deltas_summary <-
deltas |>
left_join(known_ratios, by = "isotopocule") |>
mutate(ratio_corr = (delta/1000 + 1) * ratio_known) |>
group_by(tuning, sample_name, isotopocule) |>
summarize(
isotopocule_ratio = sprintf("%s/%s", isotopocule[1], basepeak[1]),
n = n(),
ratio_raw_mean = mean(ratio),
ratio_raw_sdev_rel = sd(ratio) / ratio_raw_mean,
ratio_corrected_mean = mean(ratio_corr),
ratio_corrected_sdev_rel = sd(ratio_corr) / ratio_corrected_mean,
delta_mean = mean(delta),
delta_sdev = sd(delta),
.groups = "drop"
) |>
# add expected values
left_join(rename(delta_USGS32, delta_expected = delta_known), by = "isotopocule") |>
left_join(rename(known_ratios, reference_ratio = ratio_known), by = "isotopocule") |>
mutate(ratio_expected = (delta_expected/1000 + 1) * reference_ratio) |>
select(-isotopocule)Tables
Table 2
deltas_summary |>
filter(tuning == "optimal") |>
select(
sample_name, isotopocule_ratio, ratio_expected, n,
starts_with("ratio_raw"), starts_with("ratio_corrected")
) |>
orbi_export_data_to_excel("output/table2_nitrate_ratios.xlsx")Table 3
deltas_summary |>
filter(tuning == "optimal") |>
select(
isotopocule_ratio, reference, reference_ratio,
ratio_corrected_mean, delta_expected, n, delta_mean, delta_sdev
) |>
orbi_export_data_to_excel("output/table3_nitrate_deltas.xlsx")Table 4
deltas_summary |>
filter(tuning == "bad") |>
select(
isotopocule_ratio, reference, reference_ratio,
ratio_corrected_mean, delta_expected, n, delta_mean, delta_sdev
) |>
orbi_export_data_to_excel("output/table4_nitrate_deltas_bad_tuning.xlsx")Figures
Bonus Figure: raw ratios
flow_injections_optimal_w_basepeak |> orbi_plot_raw_data(y = ratio) 
Figure 12 panel a (EIC trace)
fig_a <- flow_injections_bad |>
orbi_plot_raw_data(
isotopocule = "M0",
# FIME: should be intensity!
y = ions.incremental,
y_scale = "log",
color = sample_name,
show_outliers = FALSE
) +
# customize intensity plot
facet_grid(. ~ injection, scales = "free_x", space = "free_x") +
#coord_cartesian(ylim = c(3.1e7, 1.2e9)) +
labs(y = "EIC nitrate / arb. unit") +
theme(
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank()
)
fig_a
Figure 12 panel b (ratio trace)
fig_b <-
flow_injections_bad_w_basepeak |>
orbi_plot_raw_data(
isotopocule = "18O",
y = ratio,
color = sample_name,
show_outliers = FALSE
) +
# customize ratio plot
facet_grid(. ~ injection, scales = "free_x", space = "free_x") +
theme(
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank()
)
fig_b
Figure 12 panel c (ratios)
fig_c <-
flow_injections_bad_summary |>
filter(isotopocule == "18O") |>
# plot
ggplot() +
aes(
x = injection,
y = ratio, ymin = ratio - ratio_sem, ymax = ratio + ratio_sem,
color = sample_name, shape = sample_name
) +
geom_pointrange() +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(21:25)) +
scale_y_continuous(breaks = scales::pretty_breaks(5)) +
scale_x_continuous(
breaks = 0:20,
labels = function(x) { x[x %%2 > 0] <- ""; x }
) +
# theme definitions
theme_bw() +
theme(
text = element_text(size = 16),
panel.grid = element_blank(),
legend.position = c(0.8, 0.5),
legend.background = element_rect(fill = "gray90", color = "black"),
legend.key = element_blank()
) +
# labels
labs(
y = expression("ratio ("^18*O/M0*")"),
color = "sample", shape = "sample"
)
fig_c
Figure 12 panel d (delta values)
fig_d <- fig_c %+%
filter(deltas, tuning == "bad", isotopocule == "18O") %+%
aes(y = delta, ymin = delta - delta_sem, ymax = delta + delta_sem) +
scale_y_continuous(breaks = 25:35, labels = function(x) paste0(x, "\U2030")) +
labs(y = expression(delta^18*O~"(USGS32 vs. VSMOW)"))
fig_d
Figure 12 combined
fig_cd <-
plot_grid(
fig_c,
fig_d + theme(legend.position = "none"),
align = "h", axis = "tb", nrow = 1
)
fig_ab <-
plot_grid(
fig_a + theme(legend.position = "none"),
fig_b + theme(legend.position = "none"),
align = "v", axis = "lr", ncol = 1
)
plot_grid(fig_ab, fig_cd, ncol = 1, rel_heights = c(2, 1))