Overview

This notebook analyzes shoot apical meristem (SAM) morphology in Inv4m and control (B73) plants. SAMs were cleared and imaged using differential interference contrast (DIC) microscopy from two-week-old seedlings grown in controlled environment chambers.

Measurements:

  • Height (h): Distance from the base to the apex of the SAM (μm)
  • Width: Maximum diameter of the SAM (μm)
  • Radius (r): Half of the width (μm)

Derived metrics (paraboloid model):

  • h/r ratio: Height-to-radius ratio (dimensionless)
  • Shape factor (h/r²): Indicates elongation; higher values = more elongated (nm⁻¹)
  • Arc length, surface area, volume: Geometric properties assuming paraboloid shape

Key hypothesis: Inv4m plants have elongated SAMs compared to B73 controls, consistent with their accelerated flowering and increased plant height phenotypes.

Load Libraries

library(dplyr)
library(tidyr)
library(purrr)
library(ggplot2)
library(ggfx)
library(ggtext)
library(ggpubr)
library(ggbeeswarm)
library(rstatix)
library(broom)
library(knitr)

Load and Process Data

# Data file: SAM measurements from DIC microscopy
# Height and Width measured in micrometers (μm)
data_file <- file.path(paths$data, "MeristemMesurements_B73-INV_December2024.csv")

# Check if file exists
if (!file.exists(data_file)) {
  stop("Data file not found: ", data_file,
       "\nPlease copy MeristemMesurements_B73-INV_December2024.csv to the data directory.")
}

microscopy_raw <- read.csv(data_file) %>%
  dplyr::select(-starts_with("X"))  # Remove any empty columns

cat("Raw data dimensions:", dim(microscopy_raw), "\n")
## Raw data dimensions: 38 6
cat("Columns:", paste(colnames(microscopy_raw), collapse = ", "), "\n")
## Columns: Plant., ID, Genotype, Width, Height, SurfaceArea
# Calculate derived metrics using paraboloid approximation of SAM shape
microscopy <- microscopy_raw %>%
  mutate(
    # Basic derived measurements
    radius = Width / 2,
    hr_ratio = Height / radius,
    shape_factor = 1000 * Height / (radius^2),  # Units: nm⁻¹ (1000 factor converts μm to nm)

    # Paraboloid model calculations
    arc_length = sqrt(radius^2 + 4 * Height^2 +
                      (radius^2 / (2 * Height)) + asinh(2 * Height / radius)),
    area_section = (4/3) * (Height * radius),
    area_surface = pi * radius * ((radius^2 + 4 * Height^2)^(3/2) - radius^3) / (6 * Height^2),
    volume = pi * (Height * radius^2) / 2
  ) %>%
  # Standardize genotype labels
  mutate(
    Genotype = case_when(
      grepl("B73", Genotype, ignore.case = TRUE) ~ "CTRL",
      grepl("INV", Genotype, ignore.case = TRUE) ~ "Inv4m",
      TRUE ~ Genotype
    ),
    Genotype = factor(Genotype, levels = c("CTRL", "Inv4m"))
  )

cat("\nProcessed data dimensions:", dim(microscopy), "\n")
## 
## Processed data dimensions: 38 13
cat("Genotype counts:\n")
## Genotype counts:
print(table(microscopy$Genotype))
## 
##  CTRL Inv4m 
##    18    20

Summary Statistics

# Define columns for analysis
primary_metrics <- c("Height", "radius", "hr_ratio", "shape_factor")

# Calculate summary statistics by genotype
summary_stats <- microscopy %>%
  group_by(Genotype) %>%
  summarise(
    n = n(),
    across(
      all_of(primary_metrics),
      list(
        mean = ~ mean(.x, na.rm = TRUE),
        sd = ~ sd(.x, na.rm = TRUE),
        sem = ~ sd(.x, na.rm = TRUE) / sqrt(sum(!is.na(.x))),
        median = ~ median(.x, na.rm = TRUE)
      ),
      .names = "{.col}_{.fn}"
    ),
    .groups = "drop"
  )

# Display summary
summary_stats %>%
  select(Genotype, n, starts_with("Height"), starts_with("radius")) %>%
  kable(digits = 2, caption = "SAM Height and Radius Summary Statistics")
SAM Height and Radius Summary Statistics
Genotype n Height_mean Height_sd Height_sem Height_median radius_mean radius_sd radius_sem radius_median
CTRL 18 111.23 17.71 4.18 112.64 62.74 4.52 1.06 63.61
Inv4m 20 120.88 16.17 3.62 121.66 62.77 2.07 0.46 62.64 
summary_stats %>%
  select(Genotype, n, starts_with("hr_ratio"), starts_with("shape_factor")) %>%
  kable(digits = 3, caption = "SAM Shape Metrics Summary Statistics")
SAM Shape Metrics Summary Statistics
Genotype n hr_ratio_mean hr_ratio_sd hr_ratio_sem hr_ratio_median shape_factor_mean shape_factor_sd shape_factor_sem shape_factor_median
CTRL 18 1.767 0.209 0.049 1.757 28.212 3.043 0.717 28.078
Inv4m 20 1.926 0.241 0.054 1.956 30.706 3.881 0.868 31.045

Statistical Analysis: Linear Model Effects

We fit a linear model for each SAM metric to estimate the effect of Inv4m genotype compared to CTRL. This provides effect sizes with standard errors and p-values comparable to a two-sample t-test.

# Fit linear models and extract Inv4m effects
lm_results <- map_df(primary_metrics, function(trait) {
  formula <- as.formula(paste(trait, "~ Genotype"))
  model <- lm(formula, data = microscopy)

  tidy(model) %>%
    filter(term == "GenotypeInv4m") %>%
    mutate(
      trait = trait,
      # Calculate 95% CI
      ci_lower = estimate - qt(0.975, model$df.residual) * std.error,
      ci_upper = estimate + qt(0.975, model$df.residual) * std.error,
      # Cohen's d effect size
      cohens_d = estimate / sd(microscopy[[trait]], na.rm = TRUE)
    ) %>%
    select(trait, estimate, std.error, statistic, p.value, ci_lower, ci_upper, cohens_d)
})

# Add units and format for display
lm_results_formatted <- lm_results %>%
  mutate(
    units = case_when(
      trait == "Height" ~ "μm",
      trait == "radius" ~ "μm",
      trait == "hr_ratio" ~ "",
      trait == "shape_factor" ~ "nm⁻¹"
    ),
    effect_str = case_when(
      trait %in% c("Height", "radius") ~ sprintf("%.2f ± %.2f %s", estimate, std.error, units),
      trait == "hr_ratio" ~ sprintf("%.3f ± %.3f", estimate, std.error),
      trait == "shape_factor" ~ sprintf("%.2f ± %.2f %s", estimate, std.error, units)
    ),
    p_str = ifelse(p.value < 0.001, sprintf("%.2e", p.value), sprintf("%.4f", p.value)),
    significance = case_when(
      p.value < 0.0001 ~ "****",
      p.value < 0.001 ~ "***",
      p.value < 0.01 ~ "**",
      p.value < 0.05 ~ "*",
      TRUE ~ "ns"
    )
  )

cat("\n=== Inv4m Effect on SAM Morphology (lm coefficients) ===\n\n")
## 
## === Inv4m Effect on SAM Morphology (lm coefficients) ===
lm_results_formatted %>%
  select(trait, effect_str, p_str, significance, cohens_d) %>%
  kable(
    col.names = c("Trait", "Effect (Inv4m - CTRL)", "p-value", "Sig.", "Cohen's d"),
    digits = 3,
    caption = "Effect of Inv4m genotype on SAM morphology metrics"
  )
Effect of Inv4m genotype on SAM morphology metrics
Trait Effect (Inv4m - CTRL) p-value Sig. Cohen’s d
Height 9.65 ± 5.50 μm 0.0876 ns 0.555
radius 0.03 ± 1.12 μm 0.9779 ns 0.009
hr_ratio 0.158 ± 0.074 0.0382 * 0.667
shape_factor 2.49 ± 1.14 nm⁻¹ 0.0354 * 0.677

Statistical Analysis: One-Tailed t-Tests

Given our directional hypothesis that Inv4m plants have larger/more elongated SAMs, we also perform one-tailed t-tests.

# Prepare data for t-tests
to_plot <- microscopy %>%
  pivot_longer(
    cols = all_of(primary_metrics),
    names_to = "trait",
    values_to = "value"
  ) %>%
  mutate(
    trait = factor(trait, levels = primary_metrics),
    Genotype_label = case_when(
      Genotype == "Inv4m" ~ "<i>Inv4m</i>",
      Genotype == "CTRL" ~ "CTRL"
    ),
    Genotype_label = factor(Genotype_label, levels = c("CTRL", "<i>Inv4m</i>"))
  )

# One-tailed t-tests (Inv4m > CTRL)
t_test_results <- to_plot %>%
  group_by(trait) %>%
  t_test(
    value ~ Genotype_label,
    ref.group = "CTRL",
    alternative = "less",  # Tests if CTRL < Inv4m
    detailed = TRUE
  ) %>%
  add_significance("p") %>%
  add_y_position(scales = "free_y")

# Adjust y positions for plotting
t_test_results$y.position <- 1.05 * t_test_results$y.position

cat("\n=== One-Tailed t-Tests (H₀: μ_Inv4m ≤ μ_CTRL) ===\n\n")
## 
## === One-Tailed t-Tests (H₀: μ_Inv4m ≤ μ_CTRL) ===
t_test_results %>%
  select(trait, estimate1, estimate2, statistic, p, p.signif) %>%
  mutate(
    diff = estimate2 - estimate1,
    pct_change = 100 * diff / estimate1
  ) %>%
  kable(
    col.names = c("Trait", "Mean CTRL", "Mean Inv4m", "t", "p", "Sig.", "Diff", "% Change"),
    digits = 3,
    caption = "One-tailed t-tests comparing SAM metrics between genotypes"
  )
One-tailed t-tests comparing SAM metrics between genotypes
Trait Mean CTRL Mean Inv4m t p Sig. Diff % Change
Height 111.229 120.881 -1.748 0.045 * 9.652 8.678
radius 62.739 62.770 -0.027 0.489 ns 0.031 0.050
hr_ratio 1.767 1.926 -2.169 0.018 * 0.158 8.966
shape_factor 28.212 30.706 -2.215 0.017 * 2.494 8.838

Visualization

SAM Morphology Boxplots (Raw Values)

# Color palette matching other figures
pal <- c(CTRL = "gold", `<i>Inv4m</i>` = "purple4")

# Theme matching ionome/phenotype plots
pheno_theme <- theme_classic2(base_size = 20) +
  theme(
    plot.title = element_markdown(hjust = 0.5, face = "bold"),
    axis.title.y = element_markdown(),
    axis.title.x = element_blank(),
    axis.text.x = element_markdown(color = "black", size = 20),
    strip.background = element_blank(),
    strip.text = element_markdown(size = 20),
    legend.position = "none"
  )

# Y-axis labels with units
yunit_labels <- c(
  Height = "Height<br><i>h</i> [μm]",
  radius = "Radius<br><i>r</i> [μm]",
  hr_ratio = "Ratio<br><i>h/r</i>",
  shape_factor = "Shape factor<br><i>h/r<sup>2</sup></i> [nm<sup>-1</sup>]"
)
# Function to create individual trait plots
plot_sam_trait <- function(trait_name, data, t_results) {
  data %>%
    filter(trait == trait_name) %>%
    ggplot(aes(x = Genotype_label, y = value, color = Genotype_label)) +
    ggtitle(gsub("_", " ", trait_name)) +
    ylab(yunit_labels[trait_name]) +
    geom_boxplot(width = 0.25, linewidth = 1, alpha = 0) %>%
    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(
      t_results %>% filter(trait == trait_name),
      label = "p.signif",
      size = 10,
      bracket.size = 0.8,
      hide.ns = TRUE
    ) +
    scale_color_manual(values = pal) +
    scale_y_continuous(expand = expansion(mult = c(0.05, 0.15))) +
    pheno_theme +
    theme(plot.margin = unit(c(1, 0.5, 1, 0.5), "lines"))
}
# Create plots for primary metrics
sam_plots <- lapply(primary_metrics, plot_sam_trait, data = to_plot, t_results = t_test_results)

# Arrange in single row
sam_figure <- ggarrange(
  plotlist = sam_plots,
  nrow = 1,
  ncol = 4,
  align = "h"
)

print(sam_figure)

Export Results

# Export effect sizes for results text
effects_export <- lm_results_formatted %>%
  select(trait, estimate, std.error, p.value, ci_lower, ci_upper, cohens_d, units)

write.csv(
  effects_export,
  file.path(paths$intermediate, "SAM_morphology_effects.csv"),
  row.names = FALSE
)

# Export summary statistics
write.csv(
  summary_stats,
  file.path(paths$intermediate, "SAM_morphology_summary_stats.csv"),
  row.names = FALSE
)

# Save figure
ggsave(
  file.path(paths$figures, "SAM_morphology_boxplots.pdf"),
  sam_figure,
  width = 12,
  height = 4
)

ggsave(
  file.path(paths$figures, "SAM_morphology_boxplots.png"),
  sam_figure,
  width = 12,
  height = 4,
  dpi = 300
)

cat("\nExported files:\n")
## 
## Exported files:
cat("- ", file.path(paths$intermediate, "SAM_morphology_effects.csv"), "\n")
## -  /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/inversion_paper/intermediate/SAM_morphology_effects.csv
cat("- ", file.path(paths$intermediate, "SAM_morphology_summary_stats.csv"), "\n")
## -  /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/inversion_paper/intermediate/SAM_morphology_summary_stats.csv
cat("- ", file.path(paths$figures, "SAM_morphology_boxplots.pdf"), "\n")
## -  /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/inversion_paper/figures/SAM_morphology_boxplots.pdf
cat("- ", file.path(paths$figures, "SAM_morphology_boxplots.png"), "\n")
## -  /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/inversion_paper/figures/SAM_morphology_boxplots.png

Results Summary for Manuscript

# Generate formatted text for results section
cat("\n=== Formatted Results for Manuscript ===\n\n")
## 
## === Formatted Results for Manuscript ===
for (i in 1:nrow(lm_results_formatted)) {
  row <- lm_results_formatted[i, ]
  trait_name <- switch(row$trait,
    "Height" = "SAM height",
    "radius" = "SAM radius",
    "hr_ratio" = "height-to-radius ratio (h/r)",
    "shape_factor" = "shape factor (h/r²)"
  )

  cat(sprintf(
    "%s: %.2f ± %.2f %s (p = %s)\n",
    trait_name,
    row$estimate,
    row$std.error,
    row$units,
    row$p_str
  ))
}
## SAM height: 9.65 ± 5.50 μm (p = 0.0876)
## SAM radius: 0.03 ± 1.12 μm (p = 0.9779)
## height-to-radius ratio (h/r): 0.16 ± 0.07  (p = 0.0382)
## shape factor (h/r²): 2.49 ± 1.14 nm⁻¹ (p = 0.0354)
# Percent changes
cat("\n=== Percent Changes ===\n\n")
## 
## === Percent Changes ===
pct_changes <- t_test_results %>%
  mutate(
    pct_change = 100 * (estimate2 - estimate1) / estimate1
  ) %>%
  select(trait, pct_change)

for (i in 1:nrow(pct_changes)) {
  cat(sprintf("%s: %.1f%% increase\n", pct_changes$trait[i], pct_changes$pct_change[i]))
}
## Height: 8.7% increase
## radius: 0.0% increase
## hr_ratio: 9.0% increase
## shape_factor: 8.8% increase

Session Information

sessionInfo()
## R version 4.5.2 (2025-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sequoia 15.7.3
## 
## Matrix products: default
## BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/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] knitr_1.51       broom_1.0.11     rstatix_0.7.3    ggbeeswarm_0.7.3
##  [5] ggpubr_0.6.2     ggtext_0.1.2     ggfx_1.0.3       ggplot2_4.0.1   
##  [9] purrr_1.2.0      tidyr_1.3.2      dplyr_1.1.4      here_1.0.2      
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10        generics_0.1.4     xml2_1.5.1         stringi_1.8.7     
##  [5] digest_0.6.39      magrittr_2.0.4     evaluate_1.0.5     grid_4.5.2        
##  [9] RColorBrewer_1.1-3 fastmap_1.2.0      rprojroot_2.1.1    jsonlite_2.0.0    
## [13] backports_1.5.0    Formula_1.2-5      scales_1.4.0       textshaping_1.0.4 
## [17] jquerylib_0.1.4    abind_1.4-8        cli_3.6.5          rlang_1.1.6       
## [21] litedown_0.9       commonmark_2.0.0   cowplot_1.2.0      withr_3.0.2       
## [25] cachem_1.1.0       yaml_2.3.12        otel_0.2.0         tools_4.5.2       
## [29] ggsignif_0.6.4     vctrs_0.6.5        R6_2.6.1           lifecycle_1.0.4   
## [33] magick_2.9.0       stringr_1.6.0      car_3.1-3          vipor_0.4.7       
## [37] ragg_1.5.0         beeswarm_0.4.0     pkgconfig_2.0.3    pillar_1.11.1     
## [41] bslib_0.9.0        gtable_0.3.6       glue_1.8.0         Rcpp_1.1.0        
## [45] systemfonts_1.3.1  xfun_0.55          tibble_3.3.0       tidyselect_1.2.1  
## [49] farver_2.1.2       htmltools_0.5.9    labeling_0.4.3     carData_3.0-5     
## [53] rmarkdown_2.30     compiler_4.5.2     S7_0.2.1           markdown_2.0      
## [57] gridtext_0.1.5