Inv4m Ionome Treatment Effects
Your Name
2025-10-17
Overview
This analysis evaluates ionome responses to phosphorus starvation in maize, comparing wild-type (CTRL) and Inv4m mutant genotypes. The analysis uses spatially corrected data and focuses on:
- Main Figure (
p_main_combined): P, Zn, Ca, and S responses arranged in 4 columns - Supplementary Figure
(
p_supp_combined): Remaining minerals (Mg, Mn, K, Fe) in 4 columns - Ratio Analysis: Seed:Stover ratios for each mineral
- Interaction Analysis: Genotype × Treatment interactions
Each mineral plot uses nested facets showing: Tissue (stover/seed) > Genotype (CTRL/Inv4m), with Treatment on the x-axis.
Load Libraries
library(dplyr)
library(tidyr)
library(purrr)
library(ggplot2)
library(ggpubr)
library(ggbeeswarm)
library(ggfx)
library(rstatix)
library(ggtext)
library(tibble)
library(ggh4x)Define Plot Parameters
# Color palette for treatments
pal <- c("gold", "purple4")
# Reliable minerals analyzed
reliable <- c("Ca", "Fe", "K", "Mg", "Mn", "P", "S", "Zn")Load Spatially Corrected Data
# Read spatially corrected mineral concentrations
ionome <- read.csv("PSU2022_spatially_corrected_both_treatments.csv")
# Rename for consistency and filter for Inv4m and CTRL genotypes only
ionome <- ionome %>%
rename(Genotype = genotype) %>%
filter(Genotype %in% c("CTRL", "Inv4m")) %>%
droplevels()
ionome$Fe_stover[ionome$Fe_stover > 800] <- NA
ionome$Mn_stover[ionome$Mn_stover > 75] <- NA
# Ensure proper factor levels
ionome$Treatment <- factor(ionome$Treatment, levels = c("+P", "-P"))
ionome$Genotype <- factor(ionome$Genotype, levels = c("CTRL", "Inv4m"))
# Display sample sizes
ionome %>%
count(Genotype, Treatment)## Genotype Treatment n
## 1 CTRL +P 16
## 2 CTRL -P 16
## 3 Inv4m +P 16
## 4 Inv4m -P 16Calculate Seed:Stover Ratios
# Calculate ratios for each reliable mineral
for (ion in reliable) {
ionome[[paste0(ion, "_ratio")]] <-
ionome[[paste0(ion, "_seed")]] / ionome[[paste0(ion, "_stover")]]
}
# Create list of all mineral-tissue combinations
reliable_by_tissue <- paste(
rep(reliable, each = 2),
rep(c("seed", "stover"), times = length(reliable)),
sep = "_"
)
vars <- c(reliable_by_tissue, paste0(reliable, "_ratio"))Prepare Data for Plotting
# Pivot to long format for plotting (EXCLUDE ratio tissue)
to_plot <- ionome %>%
select(Genotype, Treatment, all_of(reliable_by_tissue)) %>%
pivot_longer(
cols = all_of(reliable_by_tissue),
names_to = "mineral_tissue",
values_to = "ppm"
) %>%
separate_wider_delim(
cols = "mineral_tissue",
names = c("mineral", "tissue"),
delim = "_"
) %>%
filter(tissue %in% c("seed", "stover")) %>%
mutate(
tissue = factor(tissue, levels = c("stover", "seed")),
Genotype = gsub("Inv4m", "<i>Inv4m</i>", Genotype),
Genotype = factor(Genotype, levels = c("CTRL", "<i>Inv4m</i>"))
)
# Display data structure
to_plot %>%
group_by(mineral, tissue, Genotype, Treatment) %>%
summarise(
n = n(),
mean_ppm = mean(ppm, na.rm = TRUE),
sd_ppm = sd(ppm, na.rm = TRUE),
.groups = "drop"
) %>%
arrange(mineral, tissue) %>%
print(n = 30)## # A tibble: 64 × 7
## mineral tissue Genotype Treatment n mean_ppm sd_ppm
## <chr> <fct> <fct> <fct> <int> <dbl> <dbl>
## 1 Ca stover CTRL +P 16 4534. 244.
## 2 Ca stover CTRL -P 16 4599. 408.
## 3 Ca stover <i>Inv4m</i> +P 16 4123. 531.
## 4 Ca stover <i>Inv4m</i> -P 16 4035. 274.
## 5 Ca seed CTRL +P 16 18.2 6.04
## 6 Ca seed CTRL -P 16 37.1 8.36
## 7 Ca seed <i>Inv4m</i> +P 16 18.0 7.46
## 8 Ca seed <i>Inv4m</i> -P 16 34.6 11.4
## 9 Fe stover CTRL +P 16 126. 75.7
## 10 Fe stover CTRL -P 16 147. 44.9
## 11 Fe stover <i>Inv4m</i> +P 16 113. 87.8
## 12 Fe stover <i>Inv4m</i> -P 16 136. 53.0
## 13 Fe seed CTRL +P 16 14.6 2.60
## 14 Fe seed CTRL -P 16 15.4 1.06
## 15 Fe seed <i>Inv4m</i> +P 16 15.0 0.967
## 16 Fe seed <i>Inv4m</i> -P 16 15.4 1.29
## 17 K stover CTRL +P 16 11127. 1298.
## 18 K stover CTRL -P 16 10601. 1873.
## 19 K stover <i>Inv4m</i> +P 16 10550. 2070.
## 20 K stover <i>Inv4m</i> -P 16 10125. 1627.
## 21 K seed CTRL +P 16 2979. 372.
## 22 K seed CTRL -P 16 3076. 177.
## 23 K seed <i>Inv4m</i> +P 16 2835. 351.
## 24 K seed <i>Inv4m</i> -P 16 2701. 140.
## 25 Mg stover CTRL +P 16 1560. 175.
## 26 Mg stover CTRL -P 16 1561. 95.0
## 27 Mg stover <i>Inv4m</i> +P 16 1589. 165.
## 28 Mg stover <i>Inv4m</i> -P 16 1543. 83.4
## 29 Mg seed CTRL +P 16 1065. 112.
## 30 Mg seed CTRL -P 16 967. 79.8
## # ℹ 34 more rowsPrepare Ratio Data for Plotting
# Pivot ratio data to long format
ratio_plot <- ionome %>%
select(Genotype, Treatment, all_of(paste0(reliable, "_ratio"))) %>%
pivot_longer(
cols = ends_with("_ratio"),
names_to = "mineral_ratio",
values_to = "ratio"
) %>%
mutate(
mineral = sub("_ratio$", "", mineral_ratio),
Genotype = gsub("Inv4m", "<i>Inv4m</i>", Genotype),
Genotype = factor(Genotype, levels = c("CTRL", "<i>Inv4m</i>"))
) %>%
select(-mineral_ratio) %>%
filter(!is.na(ratio), !is.infinite(ratio))
# Display data structure
ratio_plot %>%
group_by(mineral, Genotype, Treatment) %>%
summarise(
n = n(),
mean_ratio = mean(ratio, na.rm = TRUE),
sd_ratio = sd(ratio, na.rm = TRUE),
.groups = "drop"
) %>%
arrange(mineral) %>%
print(n = 20)## # A tibble: 32 × 6
## mineral Genotype Treatment n mean_ratio sd_ratio
## <chr> <fct> <fct> <int> <dbl> <dbl>
## 1 Ca CTRL +P 11 0.00398 0.00121
## 2 Ca CTRL -P 15 0.00812 0.00193
## 3 Ca <i>Inv4m</i> +P 13 0.00454 0.00214
## 4 Ca <i>Inv4m</i> -P 13 0.00855 0.00289
## 5 Fe CTRL +P 11 0.143 0.137
## 6 Fe CTRL -P 15 0.114 0.0404
## 7 Fe <i>Inv4m</i> +P 13 0.318 0.429
## 8 Fe <i>Inv4m</i> -P 13 0.167 0.174
## 9 K CTRL +P 11 0.277 0.0549
## 10 K CTRL -P 15 0.300 0.0538
## 11 K <i>Inv4m</i> +P 13 0.273 0.0391
## 12 K <i>Inv4m</i> -P 13 0.270 0.0473
## 13 Mg CTRL +P 11 0.679 0.0815
## 14 Mg CTRL -P 15 0.621 0.0509
## 15 Mg <i>Inv4m</i> +P 13 0.617 0.0780
## 16 Mg <i>Inv4m</i> -P 13 0.590 0.0366
## 17 Mn CTRL +P 11 0.0971 0.0134
## 18 Mn CTRL -P 15 0.0868 0.0163
## 19 Mn <i>Inv4m</i> +P 13 0.0917 0.0117
## 20 Mn <i>Inv4m</i> -P 13 0.0823 0.0133
## # ℹ 12 more rowsFit Linear Models for All Ionome Variables
# Fit linear models: mineral ~ Genotype * Treatment
effects <- lapply(vars, FUN = function(x) {
n_obs <- sum(!is.na(ionome[[x]]))
if (n_obs < 20) {
message("Skipping ", x, " - insufficient data (n=", n_obs, ")")
return(NULL)
}
form <- paste(x, "~ Genotype * Treatment")
model <- tryCatch(
lm(as.formula(form), data = ionome),
error = function(e) {
message("Model failed for ", x, ": ", e$message)
return(NULL)
}
)
if (is.null(model)) return(NULL)
out <- coef(summary(model)) %>%
as.data.frame() %>%
tibble::rownames_to_column("predictor")
colnames(out)[3] <- "Std.Error"
colnames(out)[5] <- "p.value"
out$response <- x
out %>% select(response, everything())
}) %>%
bind_rows() %>%
mutate(p.adjust = p.adjust(p.value, method = "fdr")) %>%
arrange(p.adjust)
cat("\n=== Model Summary ===\n")##
## === Model Summary ===cat("Total models fitted:", length(unique(effects$response)), "\n")## Total models fitted: 24cat("Total effects tested:", nrow(effects), "\n")## Total effects tested: 96Identify Significant Effects
# Filter for significant effects (FDR < 0.05)
significant <- effects %>%
filter(
p.adjust < 0.05,
!grepl("Intercept", predictor)
) %>%
arrange(p.adjust)
# Recode predictor names for plotting
significant <- significant %>%
mutate(
predictor = case_when(
predictor == "Treatment-P" ~ "-P",
predictor == "GenotypeInv4m" ~ "Inv4m",
predictor == "GenotypeInv4m:Treatment-P" ~ "Inv4m:-P",
TRUE ~ predictor
),
predictor = factor(predictor, levels = c("-P", "Inv4m", "Inv4m:-P"))
)
cat("\n=== Significant Effects (FDR < 0.05) ===\n")##
## === Significant Effects (FDR < 0.05) ===print(significant %>% select(response, predictor, Estimate, p.value, p.adjust))## response predictor Estimate p.value p.adjust
## 1 P_stover -P -1.592099e+03 2.236498e-26 1.130020e-25
## 2 P_ratio -P 1.994326e+00 2.740705e-20 1.252894e-19
## 3 P_seed -P -6.718244e+02 4.015387e-09 1.675988e-08
## 4 Zn_stover -P 6.854552e+00 8.781928e-08 3.242558e-07
## 5 Zn_ratio -P -2.306819e-01 1.057069e-06 3.624236e-06
## 6 Ca_seed -P 1.895800e+01 1.315143e-06 4.353576e-06
## 7 S_stover -P 1.130616e+02 1.401800e-06 4.485761e-06
## 8 Ca_ratio -P 4.138748e-03 1.347384e-05 4.172545e-05
## 9 S_stover Inv4m:-P -9.381206e+01 2.165378e-03 6.496135e-03
## 10 Mg_seed -P -9.745101e+01 2.457528e-03 7.149173e-03
## 11 Mg_seed Inv4m -9.555164e+01 3.850969e-03 1.087332e-02
## 12 Ca_stover Inv4m -4.114157e+02 5.033988e-03 1.380751e-02
## 13 S_seed -P 7.909508e+01 1.108442e-02 2.955846e-02# Count by effect type
cat("\n=== Summary by Effect Type ===\n")##
## === Summary by Effect Type ===print(table(significant$predictor))##
## -P Inv4m Inv4m:-P
## 10 2 1Focus on Genotype × Treatment Interactions
# Extract interaction effects specifically
interaction_effects <- effects %>%
filter(grepl(":", predictor)) %>%
mutate(
predictor_clean = case_when(
predictor == "GenotypeInv4m:Treatment-P" ~ "Inv4m:-P",
TRUE ~ predictor
)
) %>%
separate(
response,
into = c("mineral", "tissue"),
sep = "_",
remove = FALSE
) %>%
mutate(
signif = case_when(
p.adjust < 0.001 ~ "***",
p.adjust < 0.01 ~ "**",
p.adjust < 0.05 ~ "*",
TRUE ~ "ns"
)
) %>%
arrange(p.adjust)
cat("\n=== Genotype × Treatment Interactions ===\n")##
## === Genotype × Treatment Interactions ===cat("Total interactions tested:", nrow(interaction_effects), "\n")## Total interactions tested: 24cat("Significant interactions (FDR < 0.05):",
sum(interaction_effects$p.adjust < 0.05), "\n\n")## Significant interactions (FDR < 0.05): 1if (any(interaction_effects$p.adjust < 0.05)) {
cat("Significant interactions:\n")
interaction_effects %>%
filter(p.adjust < 0.05) %>%
select(mineral, tissue, Estimate, Std.Error, p.value, p.adjust, signif) %>%
print()
}## Significant interactions:
## mineral tissue Estimate Std.Error p.value p.adjust signif
## 1 S stover -93.81206 29.22574 0.002165378 0.006496135 **# Summary of significant interactions
cat("\nSummary of Genotype × Treatment Interactions:\n")##
## Summary of Genotype × Treatment Interactions:cat("Total mineral-tissue combinations tested:", nrow(interaction_effects), "\n")## Total mineral-tissue combinations tested: 24cat("Significant interactions (FDR < 0.05):",
sum(interaction_effects$p.adjust < 0.05), "\n\n")## Significant interactions (FDR < 0.05): 1if (any(interaction_effects$p.adjust < 0.05)) {
cat("Significant interactions:\n")
interaction_effects %>%
filter(p.adjust < 0.05) %>%
select(mineral, tissue, Estimate, p.adjust, signif) %>%
print()
}## Significant interactions:
## mineral tissue Estimate p.adjust signif
## 1 S stover -93.81206 0.006496135 **Calculate T-tests for Significance Brackets
# T-tests comparing treatments within each genotype-mineral-tissue combination
treatment.t.test <- to_plot %>%
group_by(tissue, Genotype, mineral) %>%
rstatix::t_test(ppm ~ Treatment) %>%
rstatix::adjust_pvalue(method = "fdr") %>%
rstatix::add_significance() %>%
rstatix::add_y_position(scales = "free_y", step.increase = 0.12)
# Display significant results
cat("\n=== T-test Results (FDR < 0.05) ===\n")##
## === T-test Results (FDR < 0.05) ===treatment.t.test %>%
filter(p.adj < 0.05) %>%
select(mineral, tissue, Genotype, statistic, p, p.adj, p.adj.signif) %>%
arrange(p.adj) %>%
print(n = 50)## # A tibble: 12 × 7
## mineral tissue Genotype statistic p p.adj p.adj.signif
## <chr> <fct> <fct> <dbl> <dbl> <dbl> <chr>
## 1 P stover <i>Inv4m</i> 18.6 1.15e-13 3.68e-12 ****
## 2 P stover CTRL 20.2 1.94e-12 3.10e-11 ****
## 3 P seed <i>Inv4m</i> 10.2 1.4 e- 9 1.49e- 8 ****
## 4 Ca seed CTRL -6.71 6.08e- 7 4.86e- 6 ****
## 5 S stover CTRL -6.07 1.72e- 6 1.10e- 5 ****
## 6 Zn stover CTRL -5.93 2.26e- 6 1.21e- 5 ****
## 7 P seed CTRL 5.67 2.20e- 5 1.01e- 4 ***
## 8 Zn stover <i>Inv4m</i> -4.79 4.35e- 5 1.74e- 4 ***
## 9 Ca seed <i>Inv4m</i> -4.40 2.61e- 4 9.28e- 4 ***
## 10 Mn seed <i>Inv4m</i> 3.32 3.15e- 3 1.01e- 2 *
## 11 Mg seed <i>Inv4m</i> 3.39 3.7 e- 3 1.08e- 2 *
## 12 S seed CTRL -2.61 1.73e- 2 4.61e- 2 *Calculate T-tests for Ratios
# T-tests comparing treatments within each genotype-mineral combination
ratio.t.test <- ratio_plot %>%
group_by(Genotype, mineral) %>%
filter(n() >= 6) %>%
rstatix::t_test(ratio ~ Treatment) %>%
rstatix::adjust_pvalue(method = "fdr") %>%
rstatix::add_significance()
# Add y positions separately for each Genotype-mineral combination
if (nrow(ratio.t.test) > 0) {
ratio.t.test <- ratio.t.test %>%
group_by(Genotype, mineral) %>%
group_modify(~ {
mineral_data <- ratio_plot %>%
filter(Genotype == .y$Genotype, mineral == .y$mineral)
.x %>%
rstatix::add_y_position(
data = mineral_data,
formula = ratio ~ Treatment,
step.increase = 0.12
)
}) %>%
ungroup()
# Validate y.position values
ratio.t.test <- ratio.t.test %>%
filter(
!is.na(y.position),
!is.infinite(y.position),
!is.nan(y.position)
)
}
# Display significant results
cat("\n=== Ratio T-test Results (FDR < 0.05) ===\n")##
## === Ratio T-test Results (FDR < 0.05) ===if (nrow(ratio.t.test) > 0) {
ratio.t.test %>%
filter(p.adj < 0.05) %>%
select(mineral, Genotype, statistic, p, p.adj, p.adj.signif) %>%
arrange(p.adj) %>%
print(n = 50)
} else {
cat("No valid t-test results\n")
}## # A tibble: 6 × 6
## mineral Genotype statistic p p.adj p.adj.signif
## <chr> <fct> <dbl> <dbl> <dbl> <chr>
## 1 P CTRL -16.3 4.77e-13 7.63e-12 ****
## 2 P <i>Inv4m</i> -15.5 1.58e-10 1.26e- 9 ****
## 3 Ca CTRL -6.70 6.92e- 7 3.69e- 6 ****
## 4 Zn CTRL 5.98 3.39e- 5 1.36e- 4 ***
## 5 Ca <i>Inv4m</i> -4.03 5.6 e- 4 1.79e- 3 **
## 6 Zn <i>Inv4m</i> 3.63 1.35e- 3 3.6 e- 3 **Generate Individual Mineral Plots
# Function to create plot for a single mineral
plot_mineral <- function(mineral_name, log_scale = FALSE) {
plot_data <- to_plot %>%
filter(
mineral == mineral_name,
tissue %in% c("seed", "stover")
) %>%
mutate(
tissue = factor(tissue, levels = c("stover", "seed")),
Genotype = factor(Genotype, levels = c("CTRL", "<i>Inv4m</i>")),
Treatment = factor(Treatment, levels = c("+P", "-P"))
)
# Get t-test results for this mineral
ttest_subset <- treatment.t.test %>%
filter(mineral == mineral_name)
p <- ggplot(plot_data, aes(x = Treatment, y = ppm, color = Treatment)) +
geom_boxplot(width = 0.5, linewidth = 1, outlier.shape = NA) %>%
with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
geom_quasirandom(size = 1) %>%
with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
scale_color_manual(values = pal) +
labs(
title = mineral_name,
x = NULL,
y = "ppm [mg/kg]"
)
# Add y-scale (log or linear)
if (log_scale) {
ttest_subset$y.position <- log2(ttest_subset$y.position)
p <- p +
scale_y_continuous(
trans = "log2",
breaks = c(5, 10, 25, 50, 100, 200, 400, 800, 1600, 3200, 6400),
expand = expansion(mult = c(0.05, 0.15))
) +
stat_pvalue_manual(
ttest_subset,
label = "p.adj.signif",
size = 5,
bracket.size = 0.5,
hide.ns = TRUE
) +
ggh4x::facet_nested(~ tissue + Genotype, scales = "free_y")
} else {
p <- p +
scale_y_continuous(expand = expansion(mult = c(0.05, 0.15))) +
stat_pvalue_manual(
ttest_subset,
label = "p.adj.signif",
size = 5,
bracket.size = 0.5,
hide.ns = TRUE
) +
ggh4x::facet_nested(~ tissue + Genotype, scales = "free_y")
}
# Add theme with element_markdown for italics
p <- p +
theme_classic(base_size = 14) +
theme(
legend.position = "none",
axis.text.x = element_text(color = "black", face = "bold", size = 15),
axis.title.y = element_text(face = "bold", size = 13),
plot.title = element_text(face = "bold", size = 20, hjust = 0.5),
strip.text = element_markdown(face = "bold", size = 13),
strip.background = element_blank(),
panel.spacing = unit(0.1, "lines"),
ggh4x.facet.nestline = element_line(color = "black", linewidth = 0.5)
)
return(p)
}Generate Individual Ratio Plots
# Function to create plot for a single mineral ratio
plot_ratio <- function(mineral_name) {
plot_data <- ratio_plot %>%
filter(mineral == mineral_name) %>%
mutate(
Genotype = factor(Genotype, levels = c("CTRL", "<i>Inv4m</i>")),
Treatment = factor(Treatment, levels = c("+P", "-P"))
)
# Check if we have enough data
if (nrow(plot_data) < 6) {
message("Insufficient data for ", mineral_name, " ratio plot")
return(NULL)
}
# Get t-test results for this mineral
ttest_subset <- ratio.t.test %>%
filter(mineral == mineral_name)
# Remove any rows with invalid y.position
if (nrow(ttest_subset) > 0) {
ttest_subset <- ttest_subset %>%
filter(
!is.na(y.position),
!is.infinite(y.position),
!is.nan(y.position)
)
}
# Create base plot
p <- ggplot(plot_data, aes(x = Treatment, y = ratio, color = Treatment)) +
geom_boxplot(width = 0.5, linewidth = 1, outlier.shape = NA) %>%
with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
geom_quasirandom(size = 1) %>%
with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
stat_pvalue_manual(
ttest_subset,
label = "p.adj.signif",
size = 5,
bracket.size = 0.5,
hide.ns = TRUE
) +
scale_color_manual(values = pal) +
facet_wrap(~ Genotype) +
labs(
x = NULL,
y = "Seed/Stover Ratio"
) +
scale_y_continuous(expand = expansion(mult = c(0.05, 0.15))) +
theme_classic(base_size = 14) +
theme(
legend.position = "none",
axis.text.x = element_text(color = "black", face = "bold", size = 15),
axis.title.y = element_text(color = "black", face = "bold",size = 13),
plot.title = element_blank(),
#plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
strip.text = element_markdown(face = "bold", size = 13),
strip.background = element_blank()
)
return(p)
}
# Test individual plots before combining
cat("\n=== Testing ratio plots ===\n")##
## === Testing ratio plots ===test_minerals <- c("P", "Zn", "Ca", "S")
for (m in test_minerals) {
tryCatch({
p_test <- plot_ratio(m)
if (!is.null(p_test)) {
cat(m, ": OK\n")
} else {
cat(m, ": NULL (insufficient data)\n")
}
}, error = function(e) {
cat(m, ": ERROR -", e$message, "\n")
})
}## P : OK
## Zn : OK
## Ca : OK
## S : OKArrange Main Figure (P, Zn, Ca, S)
# Create individual plots for mineral content
p_P <- plot_mineral("P")
p_Zn <- plot_mineral("Zn")
p_Ca <- plot_mineral("Ca", log_scale = TRUE)
p_S <- plot_mineral("S")
# Combine plots in 4 columns
p_main_combined <- ggarrange(
p_P,
p_Zn + ggpubr::rremove("ylab"),
p_Ca + ggpubr::rremove("ylab"),
p_S + ggpubr::rremove("ylab"),
ncol = 4,
nrow = 1,
align = "h"
)
print(p_main_combined)
Arrange Main Figure with Ratios
# Create ratio plots for main minerals with error handling
p_P_ratio <- tryCatch(plot_ratio("P"), error = function(e) {
message("Error creating P ratio plot: ", e$message)
NULL
})
p_Zn_ratio <- tryCatch(plot_ratio("Zn"), error = function(e) {
message("Error creating Zn ratio plot: ", e$message)
NULL
})
p_Ca_ratio <- tryCatch(plot_ratio("Ca"), error = function(e) {
message("Error creating Ca ratio plot: ", e$message)
NULL
})
p_S_ratio <- tryCatch(plot_ratio("S"), error = function(e) {
message("Error creating S ratio plot: ", e$message)
NULL
})
# Build plot list carefully
all_plots <- list()
all_plots[[1]] <- p_P
all_plots[[2]] <- p_Zn + rremove("ylab")
all_plots[[3]] <- p_Ca + rremove("ylab")
all_plots[[4]] <- p_S + rremove("ylab")
if (!is.null(p_P_ratio)) all_plots[[5]] <- p_P_ratio
if (!is.null(p_Zn_ratio)) all_plots[[6]] <- p_Zn_ratio + rremove("ylab")
if (!is.null(p_Ca_ratio)) all_plots[[7]] <- p_Ca_ratio + rremove("ylab")
if (!is.null(p_S_ratio)) all_plots[[8]] <- p_S_ratio + rremove("ylab")
# Only combine if we have ratio plots
if (length(all_plots) > 4) {
p_main_with_ratios <- ggarrange(
plotlist = all_plots,
ncol = 4,
nrow = 2,
align = "hv",
labels= LETTERS[1:length(all_plots)]
)
print(p_main_with_ratios)
} else {
message("No valid ratio plots for main minerals, showing content only")
print(p_main_combined)
}
Generate Supplementary Plots (Other Minerals)
# Create plots for supplementary minerals
p_Mg <- plot_mineral("Mg")
p_Fe <- plot_mineral("Fe")
p_K <- plot_mineral("K")
p_Mn <- plot_mineral("Mn", log_scale = TRUE)
# Combine supplementary plots in 4 columns
p_supp_combined <- ggarrange(
p_Mg,
p_Mn + ggpubr::rremove("ylab"),
p_K + ggpubr::rremove("ylab"),
p_Fe + ggpubr::rremove("ylab"),
ncol = 4,
nrow = 1,
align = "h"
)
print(p_supp_combined)
Arrange Supplementary Figure with Ratios
# Create ratio plots for supplementary minerals with error handling
p_Mg_ratio <- tryCatch(plot_ratio("Mg"), error = function(e) {
message("Error creating Mg ratio plot: ", e$message)
NULL
})
p_Mn_ratio <- tryCatch(plot_ratio("Mn"), error = function(e) {
message("Error creating Mn ratio plot: ", e$message)
NULL
})
p_K_ratio <- tryCatch(plot_ratio("K"), error = function(e) {
message("Error creating K ratio plot: ", e$message)
NULL
})
p_Fe_ratio <- tryCatch(plot_ratio("Fe"), error = function(e) {
message("Error creating Fe ratio plot: ", e$message)
NULL
})
# Build plot list carefully
all_supp_plots <- list()
all_supp_plots[[1]] <- p_Mg
all_supp_plots[[2]] <- p_Mn + rremove("ylab")
all_supp_plots[[3]] <- p_K + rremove("ylab")
all_supp_plots[[4]] <- p_Fe + rremove("ylab")
if (!is.null(p_Mg_ratio)) all_supp_plots[[5]] <- p_Mg_ratio + labs(title = NULL)
if (!is.null(p_Mn_ratio)) all_supp_plots[[6]] <- p_Mn_ratio + labs(title = NULL) + rremove("ylab")
if (!is.null(p_K_ratio)) all_supp_plots[[7]] <- p_K_ratio + labs(title = NULL) + rremove("ylab")
if (!is.null(p_Fe_ratio)) all_supp_plots[[8]] <- p_Fe_ratio + labs(title = NULL) + rremove("ylab")
# Only combine if we have ratio plots
if (length(all_supp_plots) > 4) {
p_supp_with_ratios <- ggarrange(
plotlist = all_supp_plots,
ncol = 4,
nrow = 2,
align = "hv",
labels = LETTERS[1:length(all_supp_plots)]
)
print(p_supp_with_ratios)
} else {
message("No valid ratio plots for supplementary minerals, showing content only")
print(p_supp_combined)
}
Save Output
# Save main combined figure (content only)
ggsave(
"ionome_main_figure.pdf",
plot = p_main_combined,
width = 12,
height = 4,
units = "in"
)
ggsave(
"ionome_main_figure.png",
plot = p_main_combined,
width = 12,
height = 4,
units = "in",
dpi = 300
)
# Save main figure with ratios (if exists)
if (exists("p_main_with_ratios")) {
ggsave(
"ionome_main_with_ratios.pdf",
plot = p_main_with_ratios,
width = 12,
height = 8,
units = "in"
)
ggsave(
"ionome_main_with_ratios.png",
plot = p_main_with_ratios,
width = 12,
height = 8,
units = "in",
dpi = 300
)
}
# Save supplementary combined figure (content only)
ggsave(
"ionome_supplementary_figure.pdf",
plot = p_supp_combined,
width = 12,
height = 4,
units = "in"
)
ggsave(
"ionome_supplementary_figure.png",
plot = p_supp_combined,
width = 12,
height = 4,
units = "in",
dpi = 300
)
# Save supplementary figure with ratios (if exists)
if (exists("p_supp_with_ratios")) {
ggsave(
"ionome_supplementary_with_ratios.pdf",
plot = p_supp_with_ratios,
width = 12,
height = 8,
units = "in"
)
ggsave(
"ionome_supplementary_with_ratios.png",
plot = p_supp_with_ratios,
width = 12,
height = 8,
units = "in",
dpi = 300
)
}
# Export CSV files
write.csv(effects, "ionome_linear_model_effects.csv", row.names = FALSE)
write.csv(significant, "ionome_significant_effects.csv", row.names = FALSE)
write.csv(interaction_effects, "ionome_interaction_effects.csv", row.names = FALSE)
write.csv(treatment.t.test, "ionome_treatment_ttests.csv", row.names = FALSE)
if (nrow(ratio.t.test) > 0) {
write.csv(ratio.t.test, "ionome_ratio_ttests.csv", row.names = FALSE)
}Pairwise Genotype Comparisons for Significant Effects
# Define genotype color palette
pal_geno <- c("CTRL" = "gold", "<i>Inv4m</i>" = "#4a0f82")
# Filter data for the two significant mineral-tissue combinations
mg_seed_data <- to_plot %>%
filter(mineral == "Mg", tissue == "seed") %>%
mutate(
Genotype = factor(Genotype, levels = c("CTRL", "<i>Inv4m</i>")),
Treatment = factor(Treatment, levels = c("+P", "-P"))
)
ca_stover_data <- to_plot %>%
filter(mineral == "Ca", tissue == "stover") %>%
mutate(
Genotype = factor(Genotype, levels = c("CTRL", "<i>Inv4m</i>")),
Treatment = factor(Treatment, levels = c("+P", "-P"))
)
# T-tests comparing genotypes within each treatment
genotype.t.test <- bind_rows(
mg_seed_data %>% mutate(mineral_tissue = "Mg_seed"),
ca_stover_data %>% mutate(mineral_tissue = "Ca_stover")
) %>%
group_by(mineral_tissue, Treatment) %>%
rstatix::t_test(ppm ~ Genotype) %>%
rstatix::adjust_pvalue(method = "fdr") %>%
rstatix::add_significance()
# Add y positions separately for each mineral-treatment combination
if (nrow(genotype.t.test) > 0) {
genotype.t.test <- genotype.t.test %>%
group_by(mineral_tissue, Treatment) %>%
group_modify(~ {
if (.y$mineral_tissue == "Mg_seed") {
test_data <- mg_seed_data %>%
filter(Treatment == .y$Treatment)
} else {
test_data <- ca_stover_data %>%
filter(Treatment == .y$Treatment)
}
.x %>%
rstatix::add_y_position(
data = test_data,
formula = ppm ~ Genotype,
step.increase = 0.12
)
}) %>%
ungroup()
# Validate y.position values
genotype.t.test <- genotype.t.test %>%
filter(
!is.na(y.position),
!is.infinite(y.position),
!is.nan(y.position)
)
}
# Split tests by mineral
mg_seed_genotype_test <- genotype.t.test %>%
filter(mineral_tissue == "Mg_seed")
ca_stover_genotype_test <- genotype.t.test %>%
filter(mineral_tissue == "Ca_stover")
# Plot for Mg seed
p_mg_seed_genotype <- ggplot(
mg_seed_data,
aes(x = Genotype, y = ppm, color = Genotype)
) +
geom_boxplot(width = 0.5, linewidth = 1, outlier.shape = NA) %>%
with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
geom_quasirandom(size = 2) %>%
with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
stat_pvalue_manual(
mg_seed_genotype_test,
label = "p.adj.signif",
size = 5,
bracket.size = 0.5,
hide.ns = TRUE
) +
scale_color_manual(values = pal_geno) +
facet_wrap(~ Treatment) +
labs(
title = "Mg (Seed)",
x = NULL,
y = "ppm [mg/kg]"
) +
scale_y_continuous(expand = expansion(mult = c(0.05, 0.15))) +
theme_classic(base_size = 14) +
theme(
legend.position = "none",
axis.text = element_text(color = "black", face = "bold"),
axis.text.x = element_markdown(angle = 0, hjust = 0.5),
axis.title.y = element_text(face = "bold", size = 14),
plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
strip.text = element_text(face = "bold", size = 11),
strip.background = element_blank()
)
# Plot for Ca stover
p_ca_stover_genotype <- ggplot(
ca_stover_data,
aes(x = Genotype, y = ppm, color = Genotype)
) +
geom_boxplot(width = 0.5, linewidth = 1, outlier.shape = NA) %>%
with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
geom_quasirandom(size = 2) %>%
with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 1) +
stat_pvalue_manual(
ca_stover_genotype_test,
label = "p.adj.signif",
size = 5,
bracket.size = 0.5,
hide.ns = TRUE
) +
scale_color_manual(values = pal_geno) +
facet_wrap(~ Treatment) +
labs(
title = "Ca (Stover)",
x = NULL,
y = "ppm [mg/kg]"
) +
theme_classic(base_size = 14) +
theme(
legend.position = "none",
axis.text = element_text(color = "black", face = "bold"),
axis.text.x = element_markdown(angle = 0, hjust = 0.5),
axis.title.y = element_text(face = "bold", size = 14),
plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
strip.text = element_text(face = "bold", size = 11),
strip.background = element_blank()
)
# Combine both plots
p_genotype_effects <- ggarrange(
p_mg_seed_genotype,
p_ca_stover_genotype + rremove("ylab"),
ncol = 2,
nrow = 1,
align = "h"
)
print(p_genotype_effects)
# Display test results
cat("\n=== Genotype Effects on Mg (Seed) ===\n")##
## === Genotype Effects on Mg (Seed) ===print(mg_seed_genotype_test)## # A tibble: 2 × 14
## mineral_tissue Treatment .y. group1 group2 n1 n2 statistic df
## <chr> <fct> <chr> <chr> <chr> <int> <int> <dbl> <dbl>
## 1 Mg_seed +P ppm CTRL <i>Inv4m</i> 11 13 2.46 16.3
## 2 Mg_seed -P ppm CTRL <i>Inv4m</i> 15 13 3.14 18.2
## # ℹ 5 more variables: p <dbl>, p.adj <dbl>, p.adj.signif <chr>,
## # y.position <dbl>, groups <named list>cat("\n=== Genotype Effects on Ca (Stover) ===\n")##
## === Genotype Effects on Ca (Stover) ===print(ca_stover_genotype_test)## # A tibble: 2 × 14
## mineral_tissue Treatment .y. group1 group2 n1 n2 statistic df
## <chr> <fct> <chr> <chr> <chr> <int> <int> <dbl> <dbl>
## 1 Ca_stover +P ppm CTRL <i>Inv4m</i> 14 16 2.78 21.7
## 2 Ca_stover -P ppm CTRL <i>Inv4m</i> 16 16 4.59 26.3
## # ℹ 5 more variables: p <dbl>, p.adj <dbl>, p.adj.signif <chr>,
## # y.position <dbl>, groups <named list># Save genotype comparison figure
ggsave(
"ionome_genotype_effects.pdf",
plot = p_genotype_effects,
width = 8,
height = 4,
units = "in"
)
ggsave(
"ionome_genotype_effects.png",
plot = p_genotype_effects,
width = 8,
height = 4,
units = "in",
dpi = 300
)
# Export genotype comparison statistics
write.csv(
genotype.t.test,
"ionome_genotype_comparison_ttests.csv",
row.names = FALSE
)Session Info
sessionInfo()## R version 4.5.1 (2025-06-13)
## Platform: x86_64-apple-darwin20
## Running under: macOS Sequoia 15.6.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggh4x_0.3.1 tibble_3.3.0 ggtext_0.1.2 rstatix_0.7.2
## [5] ggfx_1.0.2 ggbeeswarm_0.7.2 ggpubr_0.6.1 ggplot2_3.5.2
## [9] purrr_1.1.0 tidyr_1.3.1 dplyr_1.1.4
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.6 sass_0.4.10 generics_0.1.4 xml2_1.4.0
## [5] stringi_1.8.7 digest_0.6.37 magrittr_2.0.3 evaluate_1.0.4
## [9] grid_4.5.1 RColorBrewer_1.1-3 fastmap_1.2.0 jsonlite_2.0.0
## [13] backports_1.5.0 Formula_1.2-5 scales_1.4.0 textshaping_1.0.1
## [17] jquerylib_0.1.4 abind_1.4-8 cli_3.6.5 rlang_1.1.6
## [21] litedown_0.7 commonmark_2.0.0 cowplot_1.2.0 withr_3.0.2
## [25] cachem_1.1.0 yaml_2.3.10 tools_4.5.1 ggsignif_0.6.4
## [29] broom_1.0.9 vctrs_0.6.5 R6_2.6.1 lifecycle_1.0.4
## [33] magick_2.8.7 stringr_1.5.1 car_3.1-3 vipor_0.4.7
## [37] ragg_1.4.0 pkgconfig_2.0.3 beeswarm_0.4.0 pillar_1.11.0
## [41] bslib_0.9.0 gtable_0.3.6 glue_1.8.0 Rcpp_1.1.0
## [45] systemfonts_1.2.3 xfun_0.53 tidyselect_1.2.1 rstudioapi_0.17.1
## [49] knitr_1.50 farver_2.1.2 htmltools_0.5.8.1 labeling_0.4.3
## [53] rmarkdown_2.29 carData_3.0-5 compiler_4.5.1 markdown_2.0
## [57] gridtext_0.1.5