PSU2025 — Spatial analysis (statgenHTP + SpATS) with GDD integration
Fausto Rodríguez Zapata
October 01, 2025
Purpose
This R Markdown adapts the working PSU2025 statgenHTP + SpATS
workflow to the PSU2025 data structure and integrates GDD (Growing
Degree Days) calculations. 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_pheno2.csv"
meta_file <- "~/Desktop/PSU2025_metadata.csv"
gdd_lookup_file <- "gdd_lookup_table.csv"
planting_date <- mdy("5/20/2025")
environment_name <- "PSU2025"0b. Load GDD lookup table
gdd_at_phenology <- read_csv(gdd_lookup_file) %>%
filter(environment == environment_name)
cat("Loaded GDD lookup for", environment_name, "with", nrow(gdd_at_phenology), "dates\n")## Loaded GDD lookup for PSU2025 with 129 datescat("GDD range:", min(gdd_at_phenology$cumulative_gdd_from_planting), "to",
max(gdd_at_phenology$cumulative_gdd_from_planting), "°C\n")## GDD range: 2.285 to 1386.835 °C1. Read + tidy (produce plot-level table matching PSU2025 column names)
pheno <- read_csv(pheno_file) %>%
clean_names() %>%
rename_with(~ str_remove(., "_cm$"), ends_with("_cm"))
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,
BL = leaf_length,
BW = leaf_width
) %>%
# Calculate Estimated Blade Area (EBA) = 0.75 * BL * BW FIRST
mutate(EBA = 0.75 * BL * BW) %>%
filter(inv4m %in% c("INV4M", "CTRL")) %>%
select(rowid, nested_plantid, donor, inv4m,
anthesis_date, silking_date, DTA, DTS, PH, BL, BW, EBA) %>%
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 %>%
mutate(
genotype = as.factor(paste(donor, inv4m, sep = ".")),
environment = environment_name
) %>%
group_by(plot_id = rowid, rep, plot_row = x, plot_col = y, donor, inv4m_gt = inv4m, genotype) %>%
left_join(
gdd_at_phenology %>%
rename(anthesis_date = date, GDDA = cumulative_gdd_from_planting),
by = c("environment", "anthesis_date")
) %>%
left_join(
gdd_at_phenology %>%
rename(silking_date = date, GDDS = cumulative_gdd_from_planting),
by = c("environment", "silking_date")
) %>%
select(plot_id, rep, plot_row, plot_col, donor, inv4m_gt, genotype,
DTA, DTS, GDDA, GDDS, PH, BL, BW, EBA) %>%
summarise(
DTA = mean(DTA, na.rm = TRUE),
DTS = mean(DTS, na.rm = TRUE),
GDDA = mean(GDDA, na.rm = TRUE),
GDDS = mean(GDDS, na.rm = TRUE),
PH = mean(PH, na.rm = TRUE),
BL = mean(BL, na.rm = TRUE),
BW = mean(BW, na.rm = TRUE),
EBA = mean(EBA, na.rm = TRUE),
.groups = "drop"
)
cat("Produced plot-level data with dimensions:", dim(plot_data), "\n")## Produced plot-level data with dimensions: 64 15cat("GDD columns added: GDDA (anthesis), GDDS (silking)\n")## GDD columns added: GDDA (anthesis), GDDS (silking)cat("Plants with GDDA:", sum(!is.na(plot_data$GDDA)), "\n")## Plants with GDDA: 64cat("Plants with GDDS:", sum(!is.na(plot_data$GDDS)), "\n")## Plants with GDDS: 64head(plot_data)## # A tibble: 6 × 15
## plot_id rep plot_row plot_col donor inv4m_gt genotype DTA DTS GDDA
## <dbl> <dbl> <dbl> <dbl> <chr> <chr> <fct> <dbl> <dbl> <dbl>
## 1 1352 3 9 12 TMEX INV4M TMEX.INV4M 78.8 77.1 884.
## 2 1353 3 8 12 TMEX CTRL TMEX.CTRL 79.7 80.1 894.
## 3 1354 3 7 12 MI21 CTRL MI21.CTRL 76 75.9 849.
## 4 1355 3 6 12 MI21 INV4M MI21.INV4M 77.8 76.4 871.
## 5 1356 1 5 12 MI21 INV4M MI21.INV4M 78 77.2 873.
## 6 1357 1 4 12 TMEX CTRL TMEX.CTRL 79 80.2 886.
## # ℹ 5 more variables: GDDS <dbl>, PH <dbl>, BL <dbl>, BW <dbl>, EBA <dbl>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", "GDDA", "GDDS", "PH","EBA", "BL", "BW")
TP <- createTimePoints(
dat = single_time,
experimentName = environment_name,
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
## GDDA GDDS PH BL BW EBA timePoint
## 1 884.4133 862.1119 256.8000 77.00000 8.500000 488.8750 1970-01-01 00:00:01
## 2 894.4925 901.3164 235.8333 72.14286 8.416667 467.6875 1970-01-01 00:00:01
## 3 848.5421 846.4093 240.8571 70.85714 8.285714 441.6964 1970-01-01 00:00:01
## 4 870.9217 852.2781 252.5714 78.42857 8.625000 517.8750 1970-01-01 00:00:01
## 5 872.8213 863.4963 250.5000 70.00000 7.687500 405.8036 1970-01-01 00:00:01
## 6 885.7167 902.0294 226.7500 74.14286 8.333333 486.4375 1970-01-01 00:00:01
## 7 852.7164 862.3986 231.5714 73.71429 8.857143 490.2857 1970-01-01 00:00:01
## 8 875.3167 883.7671 259.1667 74.85714 8.785714 494.7321 1970-01-01 00:00:01
## 9 851.7400 867.5186 232.5714 73.50000 8.642857 470.7500 1970-01-01 00:00:01
## 10 885.9693 887.3856 227.2857 80.14286 8.857143 533.1429 1970-01-01 00:00:01
## 11 877.8950 883.0830 228.8000 75.60000 9.000000 510.0750 1970-01-01 00:00:01
## 12 919.4040 927.2260 245.8000 77.60000 8.600000 499.9500 1970-01-01 00:00:01
## 13 891.8900 920.6606 251.4286 71.28571 8.687500 462.5893 1970-01-01 00:00:01
## 14 915.4179 935.5929 226.8571 74.28571 8.166667 468.7500 1970-01-01 00:00:01
## 15 865.2438 862.1013 228.2500 75.75000 8.937500 508.9688 1970-01-01 00:00:01
## 16 860.2250 858.4350 235.3333 69.50000 7.833333 410.8750 1970-01-01 00:00:01
## 17 903.1894 914.3175 229.7143 78.14286 9.083333 549.3750 1970-01-01 00:00:01
## 18 879.3150 891.6250 247.8000 73.80000 8.900000 493.2750 1970-01-01 00:00:01
## 19 875.1164 881.6144 232.4286 72.00000 8.500000 463.5536 1970-01-01 00:00:01
## 20 875.3444 884.9563 224.2500 78.00000 9.187500 538.4531 1970-01-01 00:00:01
## 21 865.6770 878.4571 223.5000 76.66667 8.916667 504.4500 1970-01-01 00:00:01
## 22 885.6920 917.5020 245.6000 75.20000 9.400000 530.0250 1970-01-01 00:00:01
## 23 845.3750 854.7057 234.4286 71.14286 8.785714 471.0000 1970-01-01 00:00:01
## 24 852.3828 864.3556 226.0000 78.00000 9.000000 527.5000 1970-01-01 00:00:01
## 25 887.0081 925.6969 226.5714 81.14286 8.857143 539.5179 1970-01-01 00:00:01
## 26 868.9283 873.2121 254.8571 72.85714 9.428571 518.0893 1970-01-01 00:00:01
## 27 902.0917 906.7071 238.3333 76.80000 9.083333 535.8750 1970-01-01 00:00:01
## 28 862.3279 860.3307 231.5714 71.42857 8.571429 464.1964 1970-01-01 00:00:01
## 29 910.4700 924.4779 227.1429 82.14286 9.083333 561.0000 1970-01-01 00:00:01
## 30 869.7775 887.0919 266.5000 76.00000 8.500000 490.5000 1970-01-01 00:00:01
## 31 873.0030 875.3090 241.0000 72.60000 8.300000 452.8500 1970-01-01 00:00:01
## 32 853.2089 861.6161 235.0000 76.44444 8.888889 510.2917 1970-01-01 00:00:01
## 33 866.3725 868.9865 265.3000 79.50000 9.050000 540.3000 1970-01-01 00:00:01
## 34 870.0350 883.8888 238.1250 76.12500 9.187500 524.7188 1970-01-01 00:00:01
## 35 870.1728 871.7261 233.6667 75.00000 8.277778 468.6667 1970-01-01 00:00:01
## 36 907.8494 900.5322 243.6667 71.66667 8.722222 470.5833 1970-01-01 00:00:01
## 37 876.8405 884.4195 240.8000 77.50000 8.950000 520.5000 1970-01-01 00:00:01
## 38 890.0764 891.9244 231.7143 76.00000 9.083333 517.4375 1970-01-01 00:00:01
## 39 860.4150 879.3825 254.4000 70.33333 9.500000 527.1000 1970-01-01 00:00:01
## 40 865.0450 877.9156 233.7500 69.50000 8.187500 427.1719 1970-01-01 00:00:01
## 41 885.5531 903.7211 255.3750 77.50000 9.000000 531.1875 1970-01-01 00:00:01
## 42 900.7550 919.7708 225.0000 73.66667 8.916667 495.4375 1970-01-01 00:00:01
## 43 875.8806 871.5728 234.8889 76.00000 8.611111 492.1250 1970-01-01 00:00:01
## 44 893.3350 888.9617 233.6667 65.88889 8.000000 409.6875 1970-01-01 00:00:01
## 45 904.6738 922.2350 261.6000 78.80000 9.100000 537.9750 1970-01-01 00:00:01
## 46 899.7829 912.4436 233.1429 77.00000 9.166667 528.2500 1970-01-01 00:00:01
## 47 874.5240 886.8860 236.6000 77.50000 9.300000 541.3875 1970-01-01 00:00:01
## 48 874.1605 867.0760 244.6000 71.00000 8.944444 478.9167 1970-01-01 00:00:01
## 49 898.9229 919.0414 245.8571 74.66667 9.000000 506.7500 1970-01-01 00:00:01
## 50 885.6336 889.0507 236.2857 78.50000 9.500000 557.8750 1970-01-01 00:00:01
## 51 881.5606 888.5922 234.3333 73.88889 9.062500 502.3594 1970-01-01 00:00:01
## 52 904.0763 913.7300 239.0000 67.25000 8.333333 418.5000 1970-01-01 00:00:01
## 53 879.4644 885.5913 238.1250 77.00000 9.250000 536.0156 1970-01-01 00:00:01
## 54 895.2283 922.4522 236.7778 77.44444 8.722222 511.7083 1970-01-01 00:00:01
## 55 885.6336 914.1100 253.5714 76.85714 9.428571 543.8571 1970-01-01 00:00:01
## 56 890.9136 901.0156 240.3333 75.75000 8.812500 504.8906 1970-01-01 00:00:01
## 57 911.7080 916.8410 228.2000 79.60000 9.700000 580.2750 1970-01-01 00:00:01
## 58 891.4121 893.2850 227.5714 74.85714 9.357143 524.4107 1970-01-01 00:00:01
## 59 868.3790 860.2230 244.2000 71.40000 8.800000 473.3250 1970-01-01 00:00:01
## 60 898.0450 922.6380 266.5000 76.80000 9.500000 546.3000 1970-01-01 00:00:01
## 61 906.6039 929.5520 243.5556 81.40000 9.055556 567.7083 1970-01-01 00:00:01
## 62 873.1408 875.7550 257.6667 77.00000 9.142857 535.4375 1970-01-01 00:00:01
## 63 882.4617 888.1117 244.5556 78.00000 9.277778 544.2500 1970-01-01 00:00:01
## 64 889.5006 899.3644 226.0000 73.37500 9.125000 517.9286 1970-01-01 00:00:01
## plotId rowNum colNum
## 1 1352 9 12
## 2 1353 8 12
## 3 1354 7 12
## 4 1355 6 12
## 5 1356 5 12
## 6 1357 4 12
## 7 1358 3 12
## 8 1359 2 12
## 9 1362 2 13
## 10 1363 3 13
## 11 1364 4 13
## 12 1365 5 13
## 13 1366 6 13
## 14 1367 7 13
## 15 1368 8 13
## 16 1369 9 13
## 17 1412 9 14
## 18 1413 8 14
## 19 1414 7 14
## 20 1415 6 14
## 21 1416 5 14
## 22 1417 4 14
## 23 1418 3 14
## 24 1419 2 14
## 25 1422 2 15
## 26 1423 3 15
## 27 1424 4 15
## 28 1425 5 15
## 29 1426 6 15
## 30 1427 7 15
## 31 1428 8 15
## 32 1429 9 15
## 33 1472 9 16
## 34 1473 8 16
## 35 1474 7 16
## 36 1475 6 16
## 37 1476 5 16
## 38 1477 4 16
## 39 1478 3 16
## 40 1479 2 16
## 41 1482 2 17
## 42 1483 3 17
## 43 1484 4 17
## 44 1485 5 17
## 45 1486 6 17
## 46 1487 7 17
## 47 1488 8 17
## 48 1489 9 17
## 49 1532 9 18
## 50 1533 8 18
## 51 1534 7 18
## 52 1535 6 18
## 53 1536 5 18
## 54 1537 4 18
## 55 1538 3 18
## 56 1539 2 18
## 57 1542 2 19
## 58 1543 3 19
## 59 1544 4 19
## 60 1545 5 19
## 61 1546 6 19
## 62 1547 7 19
## 63 1548 8 19
## 64 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")
}
}
cat("Successfully fitted models for:", paste(names(all_models), collapse = ", "), "\n")## Successfully fitted models for: DTA, DTS, GDDA, GDDS, PH, EBA, BL, BW4. Combine corrected data
if (length(corrected_list) > 0) {
sp_corrected <- lapply(
corrected_list,
function(corr_df) {
corr_df %>%
inner_join(TP$`1970-01-01 00:00:01` %>%
select(plotId, repId), by = "plotId") %>%
select(timeNumber, timePoint, genotype, rowId, colId, plotId, repId)
}
) %>%
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")
cat("Corrected traits:", paste(intersect(traits, colnames(sp_corrected)),
collapse = ", "), "\n")
print(head(sp_corrected))
} else {
message("No models successfully fitted.")
}## Final corrected dimensions: 64 16
## Corrected traits: DTA, DTS, GDDA, GDDS, PH, EBA, BL, BW
## timeNumber timePoint donor genotype rowId colId plotId repId
## 1 1 1970-01-01 00:00:01 TMEX INV4M 9 12 1352 3
## 2 1 1970-01-01 00:00:01 TMEX CTRL 8 12 1353 3
## 3 1 1970-01-01 00:00:01 MI21 CTRL 7 12 1354 3
## 4 1 1970-01-01 00:00:01 MI21 INV4M 6 12 1355 3
## 5 1 1970-01-01 00:00:01 MI21 INV4M 5 12 1356 1
## 6 1 1970-01-01 00:00:01 TMEX CTRL 4 12 1357 1
## DTA DTS GDDA GDDS PH EBA BL BW
## 1 79.07806 78.91915 886.9483 884.9339 253.8939 529.1286 78.31431 9.068398
## 2 80.11322 81.56636 899.7659 919.2897 230.8841 501.5207 73.25060 8.787715
## 3 76.48460 76.89870 854.4195 859.3479 234.3971 480.5153 72.04539 8.861557
## 4 77.94059 76.83674 872.3086 857.2655 245.8155 542.9155 78.93808 8.959852
## 5 78.48991 77.96787 879.0295 872.3249 244.8992 435.3163 70.78156 8.151510
## 6 79.64875 81.05198 893.9541 912.1468 223.1626 506.6444 74.59847 8.647286plot(data= sp_corrected %>% filter(donor=="MI21"), GDDA ~ DTA)
5. Spatial diagnostic plots
if (length(all_models) > 0) {
for (trait in names(all_models)) {
cat("\n### Spatial plots for", trait, "\n")
plot(all_models[[trait]], timePoints = 1)
plot(all_models[[trait]],
timePoints = 1,
plotType = "spatial",
spaTrend = "raw")
}
}##
## ### Spatial plots for DTA
##
## ### Spatial plots for DTS
##
## ### Spatial plots for GDDA
##
## ### Spatial plots for GDDS
##
## ### Spatial plots for PH
##
## ### Spatial plots for EBA
##
## ### Spatial plots for BL
##
## ### Spatial plots for BW
6. Statistical comparisons and visualization
if (exists("sp_corrected") && nrow(sp_corrected) > 0) {
usable_phenotypes <- traits
corrected_phenotypes <- intersect(usable_phenotypes, colnames(sp_corrected))
if (length(corrected_phenotypes) > 0) {
comparison_genotypes <- unique(sp_corrected$genotype)
cat("Available genotypes:", paste(comparison_genotypes, collapse = ", "), "\n")
donors <- unique(sp_corrected$donor)
cat("Available donors:", paste(donors, collapse = ", "), "\n")
htest <- sp_corrected %>%
pivot_longer(
cols = all_of(corrected_phenotypes),
names_to = "trait",
values_to = "value"
) %>%
filter(!is.na(value)) %>%
group_by(donor, trait) %>%
filter(n_distinct(genotype) == 2) %>%
t_test(value ~ genotype) %>%
adjust_pvalue(method = "fdr") %>%
add_y_position(scales = "free") %>%
add_significance()
print(htest)
pheno_theme2 <- theme_classic2(base_size = 22) +
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 = 18, face = "bold", color = "black"),
strip.background = element_blank(),
strip.text = element_markdown(size = 18, face = "bold"),
legend.position = "none"
)
trait_plots <- list()
for (current_trait in corrected_phenotypes) {
cat("\n=== Processing trait:", current_trait, "===\n")
plot_data_trait <- sp_corrected %>%
select(donor, genotype, all_of(current_trait)) %>%
rename(value = all_of(current_trait)) %>%
filter(!is.na(value))
test_to_plot_trait <- htest %>%
filter(trait == current_trait)
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 = 8,
bracket.size = 0.8
) +
scale_y_continuous(expand = expansion(mult = c(0.1, 0.15))) +
facet_wrap(~ donor, ncol = 2) +
pheno_theme2
trait_plots[[current_trait]] <- trait_plot
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)
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)
}
if (length(trait_plots) > 0) {
cat("\n=== Creating combined plot ===\n")
combined_plot <- ggarrange(
plotlist = trait_plots,
ncol = 2,
nrow = 4,
common.legend = TRUE,
legend = "bottom"
) %>%
annotate_figure(
top = text_grob(paste(environment_name,
"— Spatially Corrected Phenotypes"),
color = "black",
face = "bold",
size = 18)
)
print(combined_plot)
}
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: 16 × 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 BL value CTRL INV4M 16 16 3.13 24.0 4.52e- 3 1.45e- 2
## 2 MI21 BW value CTRL INV4M 16 16 4.32 28.1 1.75e- 4 9.33e- 4
## 3 MI21 DTA value CTRL INV4M 16 16 -1.92 24.7 6.63e- 2 1.06e- 1
## 4 MI21 DTS value CTRL INV4M 16 16 -0.584 25.8 5.64e- 1 6.02e- 1
## 5 MI21 EBA value CTRL INV4M 16 16 3.56 23.4 1.63e- 3 6.52e- 3
## 6 MI21 GDDA value CTRL INV4M 16 16 -1.96 24.8 6.19e- 2 1.06e- 1
## 7 MI21 GDDS value CTRL INV4M 16 16 -0.517 26.1 6.09e- 1 6.09e- 1
## 8 MI21 PH value CTRL INV4M 16 16 -6.30 27.9 8.47e- 7 6.78e- 6
## 9 TMEX BL value CTRL INV4M 16 16 1.96 29.9 5.93e- 2 1.06e- 1
## 10 TMEX BW value CTRL INV4M 16 16 -1.40 29.1 1.73e- 1 2.13e- 1
## 11 TMEX DTA value CTRL INV4M 16 16 2.08 30.0 4.6 e- 2 1.06e- 1
## 12 TMEX DTS value CTRL INV4M 16 16 1.52 30.0 1.39e- 1 1.85e- 1
## 13 TMEX EBA value CTRL INV4M 16 16 1.24 28.0 2.24e- 1 2.56e- 1
## 14 TMEX GDDA value CTRL INV4M 16 16 2.05 30.0 4.88e- 2 1.06e- 1
## 15 TMEX GDDS value CTRL INV4M 16 16 1.57 30.0 1.27e- 1 1.85e- 1
## 16 TMEX PH value CTRL INV4M 16 16 -17.7 28.7 5.28e-17 8.45e-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 MI21 -1.92 0.0663 0.106 ns
## 2 TMEX 2.08 0.046 0.106 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.185 ns
## 2 MI21 -0.584 0.564 0.602 ns
##
## === Processing trait: GDDA ===
##
## --- Summary for GDDA ---
## # 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 870. 877. 7.75 12.7
## 2 TMEX 16 16 894. 885. 12.8 12.7
## Statistical test results for GDDA :
## # A tibble: 2 × 5
## donor statistic p p.adj p.adj.signif
## <chr> <dbl> <dbl> <dbl> <chr>
## 1 MI21 -1.96 0.0619 0.106 ns
## 2 TMEX 2.05 0.0488 0.106 ns
##
## === Processing trait: GDDS ===
##
## --- Summary for GDDS ---
## # 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 876. 878. 10.0 15.1
## 2 TMEX 16 16 908. 899. 17.1 17.6
## Statistical test results for GDDS :
## # A tibble: 2 × 5
## donor statistic p p.adj p.adj.signif
## <chr> <dbl> <dbl> <dbl> <chr>
## 1 TMEX 1.57 0.127 0.185 ns
## 2 MI21 -0.517 0.609 0.609 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 8.45e-16 ****
## 2 MI21 -6.30 8.47e- 7 6.78e- 6 ****
##
## === Processing trait: EBA ===
##
## --- Summary for EBA ---
## # 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 507. 465. 22.7 40.9
## 2 TMEX 16 16 525. 516. 23.3 17.6
## Statistical test results for EBA :
## # A tibble: 2 × 5
## donor statistic p p.adj p.adj.signif
## <chr> <dbl> <dbl> <dbl> <chr>
## 1 MI21 3.56 0.00163 0.00652 **
## 2 TMEX 1.24 0.224 0.256 ns
##
## === Processing trait: BL ===
##
## --- Summary for BL ---
## # 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 75.3 72.1 2.09 3.61
## 2 TMEX 16 16 77.4 75.6 2.63 2.51
## Statistical test results for BL :
## # A tibble: 2 × 5
## donor statistic p p.adj p.adj.signif
## <chr> <dbl> <dbl> <dbl> <chr>
## 1 MI21 3.13 0.00452 0.0145 *
## 2 TMEX 1.96 0.0593 0.106 ns
##
## === Processing trait: BW ===
##
## --- Summary for BW ---
## # 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 8.94 8.53 0.23 0.302
## 2 TMEX 16 16 8.93 9.03 0.224 0.188
## Statistical test results for BW :
## # A tibble: 2 × 5
## donor statistic p p.adj p.adj.signif
## <chr> <dbl> <dbl> <dbl> <chr>
## 1 MI21 4.32 0.000175 0.000933 ***
## 2 TMEX -1.40 0.173 0.213 ns
##
## === Creating combined plot ===
##
## === Overall Summary Across All Traits and Donors ===
## # A tibble: 8 × 5
## trait n_donors_tested n_donors_significant min_pvalue max_pvalue
## <chr> <int> <int> <dbl> <dbl>
## 1 PH 2 2 8.45e-16 0.00000678
## 2 BW 2 1 9.33e- 4 0.213
## 3 EBA 2 1 6.52e- 3 0.256
## 4 BL 2 1 1.45e- 2 0.106
## 5 DTA 2 0 1.06e- 1 0.106
## 6 GDDA 2 0 1.06e- 1 0.106
## 7 GDDS 2 0 1.85e- 1 0.609
## 8 DTS 2 0 1.85e- 1 0.6027. Export standardized dataset for GxE analysis
if (exists("sp_corrected") && nrow(sp_corrected) > 0) {
psu2025_standardized <- sp_corrected %>%
mutate(
genotype = case_when(
genotype == "INV4M" ~ "Inv4m",
genotype == "CTRL" ~ "CTRL",
TRUE ~ genotype
),
environment = environment_name
) %>%
select(
plotId,
genotype,
environment,
donor,
repId,
all_of(intersect(traits, colnames(sp_corrected)))
) %>%
filter(genotype %in% c("CTRL", "Inv4m")) %>%
filter(complete.cases(.))
output_filename <- paste0(environment_name, "_spatially_corrected_phenotypes.csv")
write_csv(psu2025_standardized, output_filename)
cat("Exported standardized dataset for GxE analysis:\n")
cat("File:", output_filename, "\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")
cat("Traits:", paste(names(psu2025_standardized)[6:ncol(psu2025_standardized)],
collapse = ", "), "\n")
cat("\nSample of exported data:\n")
print(head(psu2025_standardized))
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 13
## Genotypes: Inv4m, CTRL
## Donors: TMEX, MI21
## Traits: DTA, DTS, GDDA, GDDS, PH, EBA, BL, BW
##
## Sample of exported data:
## plotId genotype environment donor repId DTA DTS GDDA GDDS
## 1 1352 Inv4m PSU2025 TMEX 3 79.07806 78.91915 886.9483 884.9339
## 2 1353 CTRL PSU2025 TMEX 3 80.11322 81.56636 899.7659 919.2897
## 3 1354 CTRL PSU2025 MI21 3 76.48460 76.89870 854.4195 859.3479
## 4 1355 Inv4m PSU2025 MI21 3 77.94059 76.83674 872.3086 857.2655
## 5 1356 Inv4m PSU2025 MI21 1 78.48991 77.96787 879.0295 872.3249
## 6 1357 CTRL PSU2025 TMEX 1 79.64875 81.05198 893.9541 912.1468
## PH EBA BL BW
## 1 253.8939 529.1286 78.31431 9.068398
## 2 230.8841 501.5207 73.25060 8.787715
## 3 234.3971 480.5153 72.04539 8.861557
## 4 245.8155 542.9155 78.93808 8.959852
## 5 244.8992 435.3163 70.78156 8.151510
## 6 223.1626 506.6444 74.59847 8.647286
##
## 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/New_York
## 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