PSU2025 — Spatial analysis (statgenHTP + SpATS) — adapted from CLY2025 workflow
Auto-adapted script
2025-09-15
Purpose
This R Markdown adapts the working CLY2025 statgenHTP + SpATS
workflow to the PSU2025 data structure. It preserves the same pipeline:
create plot-level data, build a TimePoints object, fit
SpATS models using fitModels(), extract spatially corrected
values with getCorrected(), then run donor-specific
comparisons (INV4M vs CTRL).
0. Configuration / file paths
pheno_file <- "~/Desktop/PSU2025_pheno.csv"
meta_file <- "~/Desktop/PSU2025_metadata.csv"
planting_date <- mdy("5/20/2025")1. Read + tidy (produce plot-level table matching CLY2025 column names)
pheno <- read_csv(pheno_file) %>% clean_names()
metadata <- read_csv(meta_file) %>% clean_names()
field_data <- pheno %>%
separate(row_plant, into = c("rowid", "nested_plantid"), sep = "-", remove = FALSE) %>%
mutate(
rowid = as.integer(rowid),
donor = case_when(
str_detect(note, "TMEX") ~ "TMEX",
str_detect(note, "MI21") ~ "MI21",
str_detect(note, "B73") ~ "B73",
TRUE ~ "Unknown"
),
inv4m = case_when(
str_detect(note, "INV4M") ~ "INV4M",
str_detect(note, "CTRL") ~ "CTRL",
str_detect(note, "B73") ~ "B73",
TRUE ~ "Unknown"
),
anthesis_date = mdy(anthesis),
silking_date = mdy(silking),
DTA = as.numeric(anthesis_date - planting_date),
DTS = as.numeric(silking_date - planting_date),
PH = actual_height
) %>%
filter(inv4m %in% c("INV4M", "CTRL")) %>%
select(rowid, nested_plantid, donor, inv4m, DTA, DTS, PH) %>%
inner_join(metadata %>% select(rowid, x=x_row, y=y_range), by = "rowid")
field_data <- field_data %>%
mutate(
rep_row = x - min(x) + 1,
rep_col = y - min(y) + 1,
rep = case_when(
rep_row <= 4 & rep_col <= 4 ~ 1,
rep_row <= 4 & rep_col > 4 ~ 2,
rep_row > 4 & rep_col <= 4 ~ 3,
rep_row > 4 & rep_col > 4 ~ 4,
TRUE ~ NA_integer_
)
)
plot_data <- field_data %>%
group_by(plot_id = rowid, rep, plot_row = x, plot_col = y, donor, inv4m_gt = inv4m) %>%
summarise(
DTA = mean(DTA, na.rm = TRUE),
DTS = mean(DTS, na.rm = TRUE),
PH = mean(PH, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(genotype = as.factor(paste(donor, inv4m_gt, sep = "."))) %>%
select(plot_id, rep, plot_row, plot_col, donor, inv4m_gt, genotype, DTA, DTS, PH)
cat("Produced plot-level data with dimensions:", dim(plot_data), "\n")## Produced plot-level data with dimensions: 64 10head(plot_data)## # A tibble: 6 × 10
## plot_id rep plot_row plot_col donor inv4m_gt genotype DTA DTS PH
## <dbl> <dbl> <dbl> <dbl> <chr> <chr> <fct> <dbl> <dbl> <dbl>
## 1 1352 3 9 12 TMEX INV4M TMEX.INV4M 78.8 77.1 257.
## 2 1353 3 8 12 TMEX CTRL TMEX.CTRL 79.7 80.1 236.
## 3 1354 3 7 12 MI21 CTRL MI21.CTRL 76 75.9 241.
## 4 1355 3 6 12 MI21 INV4M MI21.INV4M 77.8 76.4 253.
## 5 1356 1 5 12 MI21 INV4M MI21.INV4M 78 77.2 250.
## 6 1357 1 4 12 TMEX CTRL TMEX.CTRL 79 80.2 227.2. Build TimePoints object
single_time <- plot_data %>%
mutate(
time = 1,
plotId = as.character(plot_id),
rep = as.integer(rep),
plot_row = as.integer(plot_row),
plot_col = as.integer(plot_col)
)
traits <- c("DTA", "DTS", "PH")
TP <- createTimePoints(
dat = single_time,
experimentName = "PSU2025",
genotype = "genotype",
timePoint = "time",
plotId = "plotId",
repId = "rep",
rowNum = "plot_row",
colNum = "plot_col",
addCheck = FALSE
)
TP## $`1970-01-01 00:00:01`
## plot_id repId rowId colId donor inv4m_gt genotype DTA DTS
## 1 1352 3 9 12 TMEX INV4M TMEX.INV4M 78.83333 77.12500
## 2 1353 3 8 12 TMEX CTRL TMEX.CTRL 79.66667 80.14286
## 3 1354 3 7 12 MI21 CTRL MI21.CTRL 76.00000 75.85714
## 4 1355 3 6 12 MI21 INV4M MI21.INV4M 77.83333 76.37500
## 5 1356 1 5 12 MI21 INV4M MI21.INV4M 78.00000 77.25000
## 6 1357 1 4 12 TMEX CTRL TMEX.CTRL 79.00000 80.25000
## 7 1358 1 3 12 MI21 CTRL MI21.CTRL 76.42857 77.14286
## 8 1359 1 2 12 TMEX INV4M TMEX.INV4M 78.16667 78.85714
## 9 1362 1 2 13 MI21 INV4M MI21.INV4M 76.28571 77.57143
## 10 1363 1 3 13 TMEX CTRL TMEX.CTRL 79.00000 79.12500
## 11 1364 1 4 13 MI21 CTRL MI21.CTRL 78.40000 78.80000
## 12 1365 1 5 13 TMEX INV4M TMEX.INV4M 81.60000 82.20000
## 13 1366 3 6 13 TMEX INV4M TMEX.INV4M 79.42857 81.62500
## 14 1367 3 7 13 TMEX CTRL TMEX.CTRL 81.28571 82.85714
## 15 1368 3 8 13 MI21 CTRL MI21.CTRL 77.37500 77.12500
## 16 1369 3 9 13 MI21 INV4M MI21.INV4M 77.00000 76.83333
## 17 1412 3 9 14 TMEX CTRL TMEX.CTRL 80.37500 81.25000
## 18 1413 3 8 14 TMEX INV4M TMEX.INV4M 78.50000 79.50000
## 19 1414 3 7 14 MI21 INV4M MI21.INV4M 78.14286 78.66667
## 20 1415 3 6 14 MI21 CTRL MI21.CTRL 78.12500 78.87500
## 21 1416 1 5 14 TMEX CTRL TMEX.CTRL 77.40000 78.42857
## 22 1417 1 4 14 TMEX INV4M TMEX.INV4M 79.00000 81.60000
## 23 1418 1 3 14 MI21 INV4M MI21.INV4M 75.71429 76.57143
## 24 1419 1 2 14 MI21 CTRL MI21.CTRL 76.33333 77.33333
## 25 1422 1 2 15 TMEX CTRL TMEX.CTRL 79.12500 82.12500
## 26 1423 1 3 15 TMEX INV4M TMEX.INV4M 77.66667 78.00000
## 27 1424 1 4 15 MI21 INV4M MI21.INV4M 80.33333 80.71429
## 28 1425 1 5 15 MI21 CTRL MI21.CTRL 77.14286 77.00000
## 29 1426 3 6 15 TMEX CTRL TMEX.CTRL 81.00000 82.00000
## 30 1427 3 7 15 TMEX INV4M TMEX.INV4M 77.75000 79.12500
## 31 1428 3 8 15 MI21 INV4M MI21.INV4M 78.00000 78.20000
## 32 1429 3 9 15 MI21 CTRL MI21.CTRL 76.44444 77.11111
## 33 1472 4 9 16 TMEX INV4M TMEX.INV4M 77.40000 77.80000
## 34 1473 4 8 16 TMEX CTRL TMEX.CTRL 77.75000 78.87500
## 35 1474 4 7 16 MI21 CTRL MI21.CTRL 77.77778 77.88889
## 36 1475 4 6 16 MI21 INV4M MI21.INV4M 80.75000 80.22222
## 37 1476 2 5 16 MI21 CTRL MI21.CTRL 78.30000 78.90000
## 38 1477 2 4 16 TMEX CTRL TMEX.CTRL 79.28571 79.50000
## 39 1478 2 3 16 TMEX INV4M TMEX.INV4M 77.00000 78.50000
## 40 1479 2 2 16 MI21 INV4M MI21.INV4M 77.37500 78.37500
## 41 1482 2 2 17 TMEX INV4M TMEX.INV4M 79.00000 80.44444
## 42 1483 2 3 17 TMEX CTRL TMEX.CTRL 80.16667 81.66667
## 43 1484 2 4 17 MI21 CTRL MI21.CTRL 78.22222 77.88889
## 44 1485 2 5 17 MI21 INV4M MI21.INV4M 79.55556 79.22222
## 45 1486 4 6 17 TMEX INV4M TMEX.INV4M 80.50000 81.80000
## 46 1487 4 7 17 TMEX CTRL TMEX.CTRL 80.14286 81.14286
## 47 1488 4 8 17 MI21 CTRL MI21.CTRL 78.10000 79.10000
## 48 1489 4 9 17 MI21 INV4M MI21.INV4M 78.10000 77.50000
## 49 1532 4 9 18 TMEX INV4M TMEX.INV4M 80.00000 81.57143
## 50 1533 4 8 18 TMEX CTRL TMEX.CTRL 79.00000 79.28571
## 51 1534 4 7 18 MI21 CTRL MI21.CTRL 78.66667 79.22222
## 52 1535 4 6 18 MI21 INV4M MI21.INV4M 80.50000 81.25000
## 53 1536 2 5 18 MI21 CTRL MI21.CTRL 78.50000 79.00000
## 54 1537 2 4 18 TMEX CTRL TMEX.CTRL 79.77778 81.88889
## 55 1538 2 3 18 TMEX INV4M TMEX.INV4M 79.00000 81.28571
## 56 1539 2 2 18 MI21 INV4M MI21.INV4M 79.42857 80.25000
## 57 1542 2 2 19 TMEX CTRL TMEX.CTRL 81.00000 81.40000
## 58 1543 2 3 19 MI21 CTRL MI21.CTRL 79.42857 79.57143
## 59 1544 2 4 19 MI21 INV4M MI21.INV4M 77.60000 77.00000
## 60 1545 2 5 19 TMEX INV4M TMEX.INV4M 80.00000 82.00000
## 61 1546 4 6 19 TMEX CTRL TMEX.CTRL 80.66667 82.40000
## 62 1547 4 7 19 TMEX INV4M TMEX.INV4M 78.00000 78.28571
## 63 1548 4 8 19 MI21 INV4M MI21.INV4M 78.77778 79.22222
## 64 1549 4 9 19 MI21 CTRL MI21.CTRL 79.33333 80.11111
## PH timePoint plotId rowNum colNum
## 1 256.8000 1970-01-01 00:00:01 1352 9 12
## 2 235.8333 1970-01-01 00:00:01 1353 8 12
## 3 240.8571 1970-01-01 00:00:01 1354 7 12
## 4 252.5714 1970-01-01 00:00:01 1355 6 12
## 5 250.5000 1970-01-01 00:00:01 1356 5 12
## 6 226.7500 1970-01-01 00:00:01 1357 4 12
## 7 231.5714 1970-01-01 00:00:01 1358 3 12
## 8 259.1667 1970-01-01 00:00:01 1359 2 12
## 9 232.5714 1970-01-01 00:00:01 1362 2 13
## 10 227.2857 1970-01-01 00:00:01 1363 3 13
## 11 228.8000 1970-01-01 00:00:01 1364 4 13
## 12 245.8000 1970-01-01 00:00:01 1365 5 13
## 13 251.4286 1970-01-01 00:00:01 1366 6 13
## 14 226.8571 1970-01-01 00:00:01 1367 7 13
## 15 228.2500 1970-01-01 00:00:01 1368 8 13
## 16 235.3333 1970-01-01 00:00:01 1369 9 13
## 17 229.7143 1970-01-01 00:00:01 1412 9 14
## 18 247.8000 1970-01-01 00:00:01 1413 8 14
## 19 232.4286 1970-01-01 00:00:01 1414 7 14
## 20 224.2500 1970-01-01 00:00:01 1415 6 14
## 21 223.5000 1970-01-01 00:00:01 1416 5 14
## 22 245.6000 1970-01-01 00:00:01 1417 4 14
## 23 234.4286 1970-01-01 00:00:01 1418 3 14
## 24 226.0000 1970-01-01 00:00:01 1419 2 14
## 25 226.5714 1970-01-01 00:00:01 1422 2 15
## 26 254.8571 1970-01-01 00:00:01 1423 3 15
## 27 238.3333 1970-01-01 00:00:01 1424 4 15
## 28 231.5714 1970-01-01 00:00:01 1425 5 15
## 29 227.1429 1970-01-01 00:00:01 1426 6 15
## 30 266.5000 1970-01-01 00:00:01 1427 7 15
## 31 241.0000 1970-01-01 00:00:01 1428 8 15
## 32 235.0000 1970-01-01 00:00:01 1429 9 15
## 33 265.3000 1970-01-01 00:00:01 1472 9 16
## 34 238.1250 1970-01-01 00:00:01 1473 8 16
## 35 233.6667 1970-01-01 00:00:01 1474 7 16
## 36 243.6667 1970-01-01 00:00:01 1475 6 16
## 37 240.8000 1970-01-01 00:00:01 1476 5 16
## 38 231.7143 1970-01-01 00:00:01 1477 4 16
## 39 254.4000 1970-01-01 00:00:01 1478 3 16
## 40 233.7500 1970-01-01 00:00:01 1479 2 16
## 41 255.3750 1970-01-01 00:00:01 1482 2 17
## 42 225.0000 1970-01-01 00:00:01 1483 3 17
## 43 234.8889 1970-01-01 00:00:01 1484 4 17
## 44 233.6667 1970-01-01 00:00:01 1485 5 17
## 45 261.6000 1970-01-01 00:00:01 1486 6 17
## 46 233.1429 1970-01-01 00:00:01 1487 7 17
## 47 236.6000 1970-01-01 00:00:01 1488 8 17
## 48 244.6000 1970-01-01 00:00:01 1489 9 17
## 49 245.8571 1970-01-01 00:00:01 1532 9 18
## 50 236.2857 1970-01-01 00:00:01 1533 8 18
## 51 234.3333 1970-01-01 00:00:01 1534 7 18
## 52 239.0000 1970-01-01 00:00:01 1535 6 18
## 53 238.1250 1970-01-01 00:00:01 1536 5 18
## 54 236.7778 1970-01-01 00:00:01 1537 4 18
## 55 253.5714 1970-01-01 00:00:01 1538 3 18
## 56 240.3333 1970-01-01 00:00:01 1539 2 18
## 57 228.2000 1970-01-01 00:00:01 1542 2 19
## 58 227.5714 1970-01-01 00:00:01 1543 3 19
## 59 244.2000 1970-01-01 00:00:01 1544 4 19
## 60 266.5000 1970-01-01 00:00:01 1545 5 19
## 61 243.5556 1970-01-01 00:00:01 1546 6 19
## 62 257.6667 1970-01-01 00:00:01 1547 7 19
## 63 244.5556 1970-01-01 00:00:01 1548 8 19
## 64 226.0000 1970-01-01 00:00:01 1549 9 19
##
## attr(,"experimentName")
## [1] "PSU2025"
## attr(,"timePoints")
## timeNumber timePoint
## 1 1 1970-01-01 00:00:01
## attr(,"plotLimObs")
## character(0)
## attr(,"class")
## [1] "TP" "list"3. Fit SpATS models
corrected_list <- list()
fit_info <- list()
all_models <- list()
for (phen in traits) {
no_nas <- single_time %>% filter(!is.na(.data[[phen]]))
if (nrow(no_nas) < 10) {
message("Insufficient data for ", phen, " (n=", nrow(no_nas), ")")
next
}
fit_try <- try(
fitModels(
TP = TP,
trait = phen,
timePoints = 1,
what = "fixed"
),
silent = TRUE
)
if (inherits(fit_try, "try-error")) {
message("Fit failed for ", phen)
next
}
all_models[[phen]] <- fit_try
corr <- try(getCorrected(fit_try), silent = TRUE)
if (!inherits(corr, "try-error") && !is.null(corr)) {
corr <- corr %>% select(-phen)
corrected_list[[phen]] <- corr
fit_info[[phen]] <- list(n_obs = nrow(no_nas), status = "success")
}
}4. Combine corrected data
if (length(corrected_list) > 0) {
sp_corrected <- lapply(
corrected_list,
function(corr_df) {
corr_df %>%
select(timeNumber, timePoint, genotype, rowId, colId, plotId)
}
) %>%
bind_rows() %>%
arrange(plotId) %>%
distinct()
for (phen in names(corrected_list)) {
phen_corr <- paste0(phen, "_corr")
corr_df <- corrected_list[[phen]] %>%
select(plotId, all_of(phen_corr))
sp_corrected <- left_join(sp_corrected, corr_df, by = "plotId")
}
colnames(sp_corrected) <- gsub("_corr", "", colnames(sp_corrected))
sp_corrected <- sp_corrected %>%
separate(genotype, c("donor","genotype"), remove = FALSE)
cat("Final corrected dimensions:", dim(sp_corrected), "\n")
print(head(sp_corrected))
} else {
message("No models successfully fitted.")
}## Final corrected dimensions: 64 10
## timeNumber timePoint donor genotype rowId colId plotId DTA
## 1 1 1970-01-01 00:00:01 TMEX INV4M 9 12 1352 79.07806
## 2 1 1970-01-01 00:00:01 TMEX CTRL 8 12 1353 80.11322
## 3 1 1970-01-01 00:00:01 MI21 CTRL 7 12 1354 76.48460
## 4 1 1970-01-01 00:00:01 MI21 INV4M 6 12 1355 77.94059
## 5 1 1970-01-01 00:00:01 MI21 INV4M 5 12 1356 78.48991
## 6 1 1970-01-01 00:00:01 TMEX CTRL 4 12 1357 79.64875
## DTS PH
## 1 78.91915 253.8939
## 2 81.56636 230.8841
## 3 76.89870 234.3971
## 4 76.83674 245.8155
## 5 77.96787 244.8992
## 6 81.05198 223.16265. Spatial diagnostic plots
if (length(all_models) > 0) {
for (trait in names(all_models)) {
cat("\n### Spatial plots for", trait, "\n")
# Model diagnostic plots
plot(all_models[[trait]], timePoints = 1)
# Spatial trend plots
plot(all_models[[trait]],
timePoints = 1,
plotType = "spatial",
spaTrend = "raw")
}
}##
## ### Spatial plots for DTA
##
## ### Spatial plots for DTS
##
## ### Spatial plots for PH
6. Statistical comparisons and visualization
if (exists("sp_corrected") && nrow(sp_corrected) > 0) {
# Get corrected phenotype names (exclude metadata columns)
usable_phenotypes <- traits
corrected_phenotypes <- intersect(usable_phenotypes, colnames(sp_corrected))
if (length(corrected_phenotypes) > 0) {
# Define comparison groups
comparison_genotypes <- unique(sp_corrected$genotype)
cat("Available genotypes:", paste(comparison_genotypes, collapse = ", "), "\n")
# Get unique donors
donors <- unique(sp_corrected$donor)
cat("Available donors:", paste(donors, collapse = ", "), "\n")
# T-test comparing genotypes with FDR correction (per donor)
htest <- sp_corrected %>%
pivot_longer(
cols = all_of(corrected_phenotypes),
names_to = "trait",
values_to = "value"
) %>%
filter(!is.na(value)) %>%
group_by(donor, trait) %>%
# Only perform t-test if we have exactly 2 genotype levels
filter(n_distinct(genotype) == 2) %>%
t_test(value ~ genotype) %>%
adjust_pvalue(method = "fdr") %>%
add_y_position(scales = "free") %>%
add_significance()
print(htest)
# Plot theme
pheno_theme2 <- theme_classic2(base_size = 16) +
theme(
plot.title = element_markdown(hjust = 0.5, face = "bold"),
axis.title.y = element_markdown(),
axis.title.x = element_blank(),
axis.text.x = element_text(size = 14, face = "bold", color = "black"),
strip.background = element_blank(),
strip.text = element_markdown(size = 14, face = "bold"),
legend.position = "none"
)
# Create individual plots for each phenotype
trait_plots <- list()
for (current_trait in corrected_phenotypes) {
cat("\n=== Processing trait:", current_trait, "===\n")
# Prepare data for this trait
plot_data_trait <- sp_corrected %>%
select(donor, genotype, all_of(current_trait)) %>%
rename(value = all_of(current_trait)) %>%
filter(!is.na(value))
# Get statistical tests for this trait
test_to_plot_trait <- htest %>%
filter(trait == current_trait)
# Create trait-specific plot with donors as facets
trait_plot <- ggplot(plot_data_trait, aes(x = genotype, y = value, color = genotype)) +
geom_boxplot(width = 0.25, linewidth = 1.5, alpha = 0) +
geom_quasirandom(size = 2.5) +
scale_color_manual(values = c("CTRL" = "gold", "INV4M" = "purple4")) +
labs(
title = current_trait,
y = paste("Corrected", current_trait),
x = "Genotype"
) +
stat_pvalue_manual(
test_to_plot_trait,
tip.length = 0.01,
step.height = 0.08,
size = 4,
bracket.size = 0.8
) +
scale_y_continuous(expand = expansion(mult = c(0.1, 0.15))) +
facet_wrap(~ donor, scales = "free_y", ncol = 2) +
pheno_theme2
# Store the plot
trait_plots[[current_trait]] <- trait_plot
# Print summary statistics for this trait
cat("\n--- Summary for", current_trait, "---\n")
summary_stats_trait <- plot_data_trait %>%
group_by(donor, genotype) %>%
summarise(
n = n(),
mean = round(mean(value, na.rm = TRUE), 3),
sd = round(sd(value, na.rm = TRUE), 3),
.groups = "drop"
) %>%
pivot_wider(
names_from = genotype,
values_from = c(n, mean, sd),
names_sep = "_"
)
print(summary_stats_trait)
# Print significance results for this trait
sig_results_trait <- test_to_plot_trait %>%
select(donor, statistic, p, p.adj, p.adj.signif) %>%
arrange(p.adj)
cat("Statistical test results for", current_trait, ":\n")
print(sig_results_trait)
}
# Combine all plots using ggarrange
if (length(trait_plots) > 0) {
cat("\n=== Creating combined plot ===\n")
combined_plot <- ggarrange(
plotlist = trait_plots,
ncol = 3,
nrow = 1,
common.legend = TRUE,
legend = "bottom"
) %>%
annotate_figure(
top = text_grob("PSU2025 — Spatially Corrected Phenotypes",
color = "black",
face = "bold",
size = 18)
)
print(combined_plot)
}
# Overall summary across all traits and donors
cat("\n=== Overall Summary Across All Traits and Donors ===\n")
alpha <- 0.05
overall_sig <- htest %>%
mutate(significant = p.adj < alpha) %>%
group_by(trait) %>%
summarise(
n_donors_tested = n(),
n_donors_significant = sum(significant),
min_pvalue = min(p.adj),
max_pvalue = max(p.adj),
.groups = "drop"
) %>%
arrange(desc(n_donors_significant), min_pvalue)
print(overall_sig)
} else {
message("No corrected phenotypes available for analysis")
}
} else {
message("No spatially corrected data available")
}## Available genotypes: INV4M, CTRL
## Available donors: TMEX, MI21
## # A tibble: 6 × 14
## donor trait .y. group1 group2 n1 n2 statistic df p p.adj
## <chr> <chr> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 MI21 DTA value CTRL INV4M 16 16 -1.92 24.7 6.63e- 2 9.94e- 2
## 2 MI21 DTS value CTRL INV4M 16 16 -0.584 25.8 5.64e- 1 5.64e- 1
## 3 MI21 PH value CTRL INV4M 16 16 -6.30 27.9 8.47e- 7 2.54e- 6
## 4 TMEX DTA value CTRL INV4M 16 16 2.08 30.0 4.6 e- 2 9.2 e- 2
## 5 TMEX DTS value CTRL INV4M 16 16 1.52 30.0 1.39e- 1 1.67e- 1
## 6 TMEX PH value CTRL INV4M 16 16 -17.7 28.7 5.28e-17 3.17e-16
## # ℹ 3 more variables: y.position <dbl>, groups <named list>, p.adj.signif <chr>
##
## === Processing trait: DTA ===
##
## --- Summary for DTA ---
## # A tibble: 2 × 7
## donor n_CTRL n_INV4M mean_CTRL mean_INV4M sd_CTRL sd_INV4M
## <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 MI21 16 16 77.8 78.3 0.623 1.03
## 2 TMEX 16 16 79.7 78.9 1.01 0.996
## Statistical test results for DTA :
## # A tibble: 2 × 5
## donor statistic p p.adj p.adj.signif
## <chr> <dbl> <dbl> <dbl> <chr>
## 1 TMEX 2.08 0.046 0.092 ns
## 2 MI21 -1.92 0.0663 0.0994 ns
##
## === Processing trait: DTS ===
##
## --- Summary for DTS ---
## # A tibble: 2 × 7
## donor n_CTRL n_INV4M mean_CTRL mean_INV4M sd_CTRL sd_INV4M
## <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 MI21 16 16 78.2 78.4 0.792 1.21
## 2 TMEX 16 16 80.8 80.0 1.34 1.39
## Statistical test results for DTS :
## # A tibble: 2 × 5
## donor statistic p p.adj p.adj.signif
## <chr> <dbl> <dbl> <dbl> <chr>
## 1 TMEX 1.52 0.139 0.167 ns
## 2 MI21 -0.584 0.564 0.564 ns
##
## === Processing trait: PH ===
##
## --- Summary for PH ---
## # A tibble: 2 × 7
## donor n_CTRL n_INV4M mean_CTRL mean_INV4M sd_CTRL sd_INV4M
## <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 MI21 16 16 232. 240. 2.94 3.90
## 2 TMEX 16 16 231. 256. 3.51 4.36
## Statistical test results for PH :
## # A tibble: 2 × 5
## donor statistic p p.adj p.adj.signif
## <chr> <dbl> <dbl> <dbl> <chr>
## 1 TMEX -17.7 5.28e-17 3.17e-16 ****
## 2 MI21 -6.30 8.47e- 7 2.54e- 6 ****
##
## === Creating combined plot ===
##
## === Overall Summary Across All Traits and Donors ===
## # A tibble: 3 × 5
## trait n_donors_tested n_donors_significant min_pvalue max_pvalue
## <chr> <int> <int> <dbl> <dbl>
## 1 PH 2 2 3.17e-16 0.00000254
## 2 DTA 2 0 9.2 e- 2 0.0994
## 3 DTS 2 0 1.67e- 1 0.5647. Export standardized dataset for GxE analysis
if (exists("sp_corrected") && nrow(sp_corrected) > 0) {
# Create standardized dataset for GxE analysis
psu2025_standardized <- sp_corrected %>%
mutate(
# Standardize genotype names to match CLY2025 format
genotype = case_when(
genotype == "INV4M" ~ "Inv4m",
genotype == "CTRL" ~ "CTRL",
TRUE ~ genotype
),
# Add environment identifier
environment = "PSU2025",
# Ensure proper column names for compatibility
repId = 1 # Add if rep info available, otherwise set to 1
) %>%
# Select only essential columns for GxE analysis
select(
plotId,
genotype,
environment,
donor,
repId,
all_of(traits) # DTA, DTS, PH
) %>%
# Filter to ensure we have the genotypes of interest
filter(genotype %in% c("CTRL", "Inv4m")) %>%
# Remove any rows with missing phenotype data
filter(complete.cases(.))
# Export the standardized dataset
write.csv(psu2025_standardized,
"PSU2025_spatially_corrected_phenotypes.csv",
row.names = FALSE)
cat("Exported standardized dataset for GxE analysis:\n")
cat("File: PSU2025_spatially_corrected_phenotypes.csv\n")
cat("Dimensions:", dim(psu2025_standardized), "\n")
cat("Genotypes:", paste(unique(psu2025_standardized$genotype), collapse = ", "), "\n")
cat("Donors:", paste(unique(psu2025_standardized$donor), collapse = ", "), "\n")
# Display sample of the exported data
cat("\nSample of exported data:\n")
print(head(psu2025_standardized))
# Summary statistics for verification
cat("\nSample sizes by donor and genotype:\n")
print(table(psu2025_standardized$donor, psu2025_standardized$genotype))
} else {
message("No spatially corrected data available for export")
}## Exported standardized dataset for GxE analysis:
## File: PSU2025_spatially_corrected_phenotypes.csv
## Dimensions: 64 8
## Genotypes: Inv4m, CTRL
## Donors: TMEX, MI21
##
## Sample of exported data:
## plotId genotype environment donor repId DTA DTS PH
## 1 1352 Inv4m PSU2025 TMEX 1 79.07806 78.91915 253.8939
## 2 1353 CTRL PSU2025 TMEX 1 80.11322 81.56636 230.8841
## 3 1354 CTRL PSU2025 MI21 1 76.48460 76.89870 234.3971
## 4 1355 Inv4m PSU2025 MI21 1 77.94059 76.83674 245.8155
## 5 1356 Inv4m PSU2025 MI21 1 78.48991 77.96787 244.8992
## 6 1357 CTRL PSU2025 TMEX 1 79.64875 81.05198 223.1626
##
## Sample sizes by donor and genotype:
##
## CTRL Inv4m
## MI21 16 16
## TMEX 16 168. Session information
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/Bogota
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggbeeswarm_0.7.2 ggtext_0.1.2 ggpubr_0.6.1 rstatix_0.7.2
## [5] broom_1.0.9 statgenHTP_1.0.9.1 janitor_2.2.1 lubridate_1.9.4
## [9] forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4 purrr_1.1.0
## [13] readr_2.1.5 tidyr_1.3.1 tibble_3.3.0 ggplot2_3.5.2
## [17] tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] dotCall64_1.2 gtable_0.3.6 beeswarm_0.4.0 spam_2.11-1
## [5] xfun_0.53 bslib_0.9.0 tzdb_0.5.0 vctrs_0.6.5
## [9] tools_4.5.1 generics_0.1.4 parallel_4.5.1 pkgconfig_2.0.3
## [13] SpATS_1.0-19 data.table_1.17.8 RColorBrewer_1.1-3 lifecycle_1.0.4
## [17] compiler_4.5.1 farver_2.1.2 ggforce_0.5.0 carData_3.0-5
## [21] snakecase_0.11.1 litedown_0.7 vipor_0.4.7 htmltools_0.5.8.1
## [25] sass_0.4.10 yaml_2.3.10 Formula_1.2-5 crayon_1.5.3
## [29] pillar_1.11.0 car_3.1-3 jquerylib_0.1.4 MASS_7.3-65
## [33] cachem_1.1.0 abind_1.4-8 commonmark_2.0.0 tidyselect_1.2.1
## [37] digest_0.6.37 stringi_1.8.7 labeling_0.4.3 cowplot_1.2.0
## [41] polyclip_1.10-7 fastmap_1.2.0 grid_4.5.1 cli_3.6.5
## [45] magrittr_2.0.3 utf8_1.2.6 withr_3.0.2 scales_1.4.0
## [49] backports_1.5.0 bit64_4.6.0-1 timechange_0.3.0 rmarkdown_2.29
## [53] bit_4.6.0 gridExtra_2.3 ggsignif_0.6.4 hms_1.1.3
## [57] evaluate_1.0.4 knitr_1.50 markdown_2.0 rlang_1.1.6
## [61] gridtext_0.1.5 Rcpp_1.1.0 glue_1.8.0 tweenr_2.0.3
## [65] LMMsolver_1.0.11 xml2_1.4.0 rstudioapi_0.17.1 vroom_1.6.5
## [69] jsonlite_2.0.0 R6_2.6.1