Cross-Study Validation of Phosphate Starvation Response Dynamics in Saccharomyces cerevisiae

First, I am just loading all of the data for the analysis. All of this is focused on analyzing the Gurvich 2017 paper firstly. Thereafter, we will simply load the Zhou and O’shea results to perform the congruent analysis.

Gene Mapping Summary:

The data loading successfully mapped gene identifiers across multiple naming systems. From the original 5,805 genes in our expression matrix, gene names were mapped from the internal gene-YXXXXX format to standard S. cerevisiae nomenclature using both systematic IDs (e.g., YBR053C) and common names (e.g., ATG8). ## Gene Set Definition: Three functionally distinct gene sets were extracted from Gurvich et al. (2017) Figure 1C:

• PHO regulon: Genes involved in phosphate acquisition and metabolism
• Stress response: Environmental stress response (ESR) genes
• Protein synthesis: Ribosomal protein genes

This extraction was motivated by the heatmap signatures observed from figure 1C. To ensure that these sets of genes can be compared to our own data, I parsed them out for direct comparison. ## Filtered Genes: Three genes from the Gurvich gene sets were absent from my data due to low expression filtering during preprocessing: YLR252W, YDR417C, and YOR309C.YLR25W is a labeled stress gene; YDR417C and YOR309C are annotated as ribosomal protein genes in the Gurvich, 2017 paper.

Fold Change Calculations:

Fold changes were calculated relative to the pre-starvation baseline (t=0, +Pi condition) for both datasets:

GRE Lab data: log2FC = log2CPM(0mM Pi, t=x) - mean(log2CPM(t=0, replicates 2 & 4)) Gurvich data: log2FC = log2Expression(0mM Pi, t=x) - log2Expression(7.3mM Pi, t=0)

For the GRE Lab dataset, biological replicates (n=2 per timepoint) were processed as follows: baseline expression was calculated as the mean log2CPM across all t=0 replicates, then each replicate’s log2FC was computed by subtracting this baseline from its expression value. Finally, mean log2FC ± standard error was calculated across replicates for each gene at each timepoint. These mean fold changes were used for all downstream correlation and concordance analyses with the Gurvich dataset, which reports single values per gene per timepoint. This is derived from the fact that only 1 replicate was displayed in the Gurvich data.

Both datasets are now in comparable formats for correlation analysis.

#create timepoint matching table (exact + approximate matches)
timepoint_matches <- tribble(
  ~my_time, ~gurvich_time, ~match_type,
  0,        0,             "exact",
  1,        1,             "exact",
  2,        2,             "exact",
  3.5,      3.5,           "exact",
  8,        8,             "exact",
  4,        5,             "approximate",
  6,        6.5,           "approximate"
)

#merge my data with gurvich data using flexible matching
correlation_data <- expr_fc_summary %>%
  filter(!is.na(gene_set)) %>%
  inner_join(timepoint_matches, by = c("timepoint_hours" = "my_time")) %>%
  select(standard_name, gene_set, my_time = timepoint_hours, gurvich_time, 
         match_type, my_log2FC = mean_log2FC, my_se = se_log2FC) %>%
  inner_join(
    gurvich_fc %>% select(gene_name, gene_set, timepoint_hours, gurvich_log2FC = log2FC),
    by = c("standard_name" = "gene_name", "gene_set", "gurvich_time" = "timepoint_hours")
  )

write.csv(correlation_data,
          file.path(output_dir, "02_correlation", "matched_timepoint_data.csv"),
          row.names = FALSE)

#scatter plots for all matched timepoints
scatter_plots <- map2(timepoint_matches$my_time, timepoint_matches$gurvich_time, function(mt, gt) {
  data_tp <- correlation_data %>% filter(my_time == mt, gurvich_time == gt)
  
  if(nrow(data_tp) == 0) return(NULL)
  
  cor_pearson <- cor.test(data_tp$my_log2FC, data_tp$gurvich_log2FC, method = "pearson")
  cor_spearman <- cor.test(data_tp$my_log2FC, data_tp$gurvich_log2FC, method = "spearman")
  match_label <- unique(data_tp$match_type)
  
  ggplot(data_tp, aes(x = gurvich_log2FC, y = my_log2FC, color = gene_set)) +
    geom_point(alpha = 0.6, size = 2) +
    geom_smooth(method = "lm", se = TRUE, color = "black", linewidth = 0.8) +
    geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray50") +
    annotate("text", x = -Inf, y = Inf,
             label = sprintf("r = %.3f, p = %.2e\nρ = %.3f, p = %.2e\nn = %d",
                           cor_pearson$estimate, cor_pearson$p.value,
                           cor_spearman$estimate, cor_spearman$p.value, nrow(data_tp)),
             hjust = 0, vjust = 1, size = 3) +
    scale_color_manual(values = c("PHO" = "#E41A1C", "Stress" = "#377EB8", "Ribosomal" = "#4DAF4A")) +
    labs(title = sprintf("My: %.1fh vs Gurvich: %.1fh (%s)", mt, gt, match_label),
         x = "Gurvich log2FC", y = "My Data log2FC") +
    theme_bw() +
    theme(legend.position = "none")
})

## Warning in cor.test.default(data_tp$my_log2FC, data_tp$gurvich_log2FC, method =
## "spearman"): Cannot compute exact p-value with ties
## Warning in cor.test.default(data_tp$my_log2FC, data_tp$gurvich_log2FC, method =
## "spearman"): Cannot compute exact p-value with ties
## Warning in cor.test.default(data_tp$my_log2FC, data_tp$gurvich_log2FC, method =
## "spearman"): Cannot compute exact p-value with ties
## Warning in cor.test.default(data_tp$my_log2FC, data_tp$gurvich_log2FC, method =
## "spearman"): Cannot compute exact p-value with ties
## Warning in cor.test.default(data_tp$my_log2FC, data_tp$gurvich_log2FC, method =
## "spearman"): Cannot compute exact p-value with ties
## Warning in cor.test.default(data_tp$my_log2FC, data_tp$gurvich_log2FC, method =
## "spearman"): Cannot compute exact p-value with ties

scatter_plots <- scatter_plots[!sapply(scatter_plots, is.null)]

combined_scatter <- cowplot::plot_grid(plotlist = scatter_plots, nrow = 2, ncol = 4)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 3 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 4 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_point()`).

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 3 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 4 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_point()`).

print(combined_scatter)

ggsave(file.path(output_dir, "02_correlation", "scatter_all_timepoints.pdf"),
       combined_scatter, width = 16, height = 8)

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'ρ = 0.160, p = 1.88e-02' in 'mbcsToSbcs': for ρ (U+03C1)

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'ρ = 0.821, p = 6.52e-54' in 'mbcsToSbcs': for ρ (U+03C1)

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'ρ = 0.835, p = 2.78e-57' in 'mbcsToSbcs': for ρ (U+03C1)

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'ρ = 0.911, p = 2.66e-84' in 'mbcsToSbcs': for ρ (U+03C1)

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'ρ = 0.907, p = 3.58e-82' in 'mbcsToSbcs': for ρ (U+03C1)

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'ρ = 0.904, p = 9.02e-81' in 'mbcsToSbcs': for ρ (U+03C1)

#correlation by gene set and timepoint
cor_by_geneset <- correlation_data %>%
  group_by(gene_set, my_time, gurvich_time, match_type) %>%
  summarise(
    n_genes = n(),
    pearson_r = cor(my_log2FC, gurvich_log2FC, method = "pearson", use = "complete.obs"),
    spearman_rho = cor(my_log2FC, gurvich_log2FC, method = "spearman", use = "complete.obs"),
    rmse = sqrt(mean((my_log2FC - gurvich_log2FC)^2, na.rm = TRUE)),
    mae = mean(abs(my_log2FC - gurvich_log2FC), na.rm = TRUE),
    .groups = "drop"
  )

write.csv(cor_by_geneset,
          file.path(output_dir, "02_correlation", "correlation_by_geneset.csv"),
          row.names = FALSE)

#correlation trend over time
cor_trend_plot <- ggplot(cor_by_geneset, aes(x = my_time, y = pearson_r, color = gene_set, shape = match_type)) +
  geom_line(aes(group = gene_set), linewidth = 1) +
  geom_point(size = 3) +
  geom_hline(yintercept = c(0.7, 0.8, 0.9), linetype = "dashed", color = "gray50", alpha = 0.5) +
  scale_color_manual(values = c("PHO" = "#E41A1C", "Stress" = "#377EB8", "Ribosomal" = "#4DAF4A")) +
  scale_shape_manual(values = c("exact" = 16, "approximate" = 1)) +
  labs(title = "Correlation Changes Over Time",
       x = "My Data Timepoint (hours)", y = "Pearson Correlation (r)",
       color = "Gene Set", shape = "Match Type") +
  theme_bw()
print(cor_trend_plot)

ggsave(file.path(output_dir, "02_correlation", "correlation_trend.pdf"),
       cor_trend_plot, width = 10, height = 6)

#overall correlation summary
overall_cor <- correlation_data %>%
  group_by(gene_set) %>%
  summarise(
    n_comparisons = n(),
    mean_pearson = mean(cor(my_log2FC, gurvich_log2FC, method = "pearson")),
    mean_spearman = mean(cor(my_log2FC, gurvich_log2FC, method = "spearman")),
    overall_rmse = sqrt(mean((my_log2FC - gurvich_log2FC)^2)),
    .groups = "drop"
  )

write.csv(overall_cor,
          file.path(output_dir, "02_correlation", "overall_correlation_summary.csv"),
          row.names = FALSE)

Correlation Analysis

Methods

Pattern concordance between the GRE Lab and Gurvich datasets was assessed using two complementary correlation metrics calculated separately for each timepoint. Pearson’s correlation coefficient (r) was used to measure the strength and direction of linear relationships between log2FC values. To assess concordance independent of magnitude, I computed Spearman’s rank correlation between matched timepoints of our dataset and Gurvich et al. (2017). Spearman’s ρ, which measures agreement in gene-wise ranking rather than absolute fold-change, remained high (ρ = 0.82–0.91 for t ≥ 2 h), indicating that both datasets capture the same relative induction and repression patterns despite differing normalization scales. Both correlation coefficients were computed using R’s cor.test() function with default parameters. For each test, the null hypothesis (H₀: correlation coefficient = 0) was evaluated against a two-tailed alternative hypothesis (H₁: correlation coefficient ≠ 0). The test statistic follows a t-distribution with n-2 degrees of freedom, where n is the number of gene-timepoint pairs.

P-values less than 0.05 were considered statistically significant. All correlation analyses were performed on matched gene-timepoint pairs where both datasets had valid log2FC values (complete cases only).

Gurvich compared to GRE Lab data

The scatter plots (Image 1) reveal that transcriptional concordance between the two datasets improves dramatically over the time course. At the earliest timepoint (1h), correlation is weak (Pearson r=0.264, Spearman ρ=0.160), with genes scattered across the plot. However, by 2h, concordance increases substantially (r=0.853, ρ=0.821), with genes aligning tightly along the regression line. This trend continues, reaching excellent concordance by 3.5h (r=0.932, ρ=0.835) and peak agreement at 8h (r=0.954, ρ=0.911). The approximate timepoint matches (4h→5h, 6h→6.5h) also show strong correlations (r>0.93, ρ>0.90), validating their comparability.

When stratified by functional gene set (Image 2), distincttemporal patterns emerge: PHO regulon (red): Maintains relatively stable, moderate-to-high correlation (r~0.68-0.75) across all timepoints. This consistency suggests the core phosphate starvation response is stereotyped between experiments, with both datasets capturing the canonical PHO pathway activation with similar kinetics. Ribosomal genes (green): Show the most dramatic improvement, starting with poor concordance at 1h (r0.36) but rapidly increasing to excellent agreement by 3.5h (r0.85), which is maintained through 8h. This pattern indicates that ribosomal protein repression—a hallmark of nutrient starvation—becomes increasingly coordinated as the biological response matures, with early variability giving way to convergent downregulation. Stress genes (blue): Display highly erratic concordance, starting near zero at 1h (r0.03), spiking to moderate levels at 2h (r0.65), then fluctuating between r=0.25-0.60 across later timepoints. This volatility suggests stress response genes are more context-dependent or heterogeneous, potentially reflecting differences in strain background, precise growth conditions, or the dynamic/transient nature of environmental stress responses.

My interpretation

The temporal pattern—poor early concordance improving to strong late concordance—suggests that both experiments capture the same underlying biological program, but early transcriptional dynamics (0-2h) are more variable due to transient cellular responses, while the mature starvation state (3.5-8h) represents a stable, reproducible transcriptional signature. The gene set-specific patterns further indicate that core metabolic responses (PHO, ribosomal repression) are highly conserved, while stress pathways show greater experimental variability.

#prepare data for heatmaps (wide format: genes × timepoints)
my_heatmap_data <- expr_fc_summary %>%
  filter(!is.na(gene_set)) %>%
  select(standard_name, gene_set, timepoint_hours, mean_log2FC) %>%
  pivot_wider(names_from = timepoint_hours, values_from = mean_log2FC, names_prefix = "t") %>%
  column_to_rownames("standard_name")

#prepare gurvich data - average duplicates first to avoid list columns
gurvich_heatmap_data <- gurvich_fc %>%
  filter(timepoint_hours <= 8) %>%
  group_by(gene_name, gene_set, timepoint_hours) %>%
  summarise(log2FC = mean(log2FC, na.rm = TRUE), .groups = "drop") %>%
  pivot_wider(names_from = timepoint_hours, values_from = log2FC, names_prefix = "t") %>%
  column_to_rownames("gene_name")

#function to create heatmap for one gene set
plot_geneset_heatmap <- function(data, geneset, title) {
  subset_data <- data %>%
    filter(gene_set == geneset) %>%
    select(-gene_set)
  
  #get only numeric columns
  numeric_cols <- names(subset_data)[sapply(subset_data, is.numeric)]
  mat <- subset_data[, numeric_cols]
  
  #order columns by timepoint
  mat <- mat[, order(as.numeric(gsub("t", "", colnames(mat))))]
  
  #convert to matrix
  mat <- as.matrix(mat)
  
  #remove rows with all NAs
  mat <- mat[rowSums(is.na(mat)) < ncol(mat), ]
  
  if(nrow(mat) < 2) {
    warning(paste("Skipping", geneset, "- not enough genes"))
    return(invisible(NULL))
  }
  
  pheatmap(mat,
           color = colorRampPalette(c("blue", "white", "red"))(100),
           breaks = seq(-4, 4, length.out = 101),
           cluster_rows = TRUE,
           cluster_cols = FALSE,
           show_rownames = TRUE,
           main = paste(title, "-", geneset),
           fontsize_row = 6,
           cellwidth = 15,
           cellheight = 8)
}

#create heatmaps for each gene set
plot_geneset_heatmap(my_heatmap_data, "PHO", "My Data")

plot_geneset_heatmap(my_heatmap_data, "Stress", "My Data")

plot_geneset_heatmap(my_heatmap_data, "Ribosomal", "My Data")

pdf(file.path(output_dir, "03_timeseries", "heatmap_GRElab-Pi-SC.pdf"), width = 10, height = 12)
plot_geneset_heatmap(my_heatmap_data, "PHO", "My Data")

plot_geneset_heatmap(my_heatmap_data, "Stress", "My Data")

plot_geneset_heatmap(my_heatmap_data, "Ribosomal", "My Data")

dev.off()

## pdf 
##   3

plot_geneset_heatmap(gurvich_heatmap_data, "PHO", "Gurvich 2017")
plot_geneset_heatmap(gurvich_heatmap_data, "Stress", "Gurvich 2017")
plot_geneset_heatmap(gurvich_heatmap_data, "Ribosomal", "Gurvich 2017")

pdf(file.path(output_dir, "03_timeseries", "heatmap_gurvich.pdf"), width = 10, height = 12)
plot_geneset_heatmap(gurvich_heatmap_data, "PHO", "Gurvich 2017")
plot_geneset_heatmap(gurvich_heatmap_data, "Stress", "Gurvich 2017")
plot_geneset_heatmap(gurvich_heatmap_data, "Ribosomal", "Gurvich 2017")
dev.off()

## pdf 
##   3

#side-by-side comparison for matched genes only
matched_genes <- intersect(rownames(my_heatmap_data), rownames(gurvich_heatmap_data))

for(gs in c("PHO", "Stress", "Ribosomal")) {
  genes_in_set <- my_heatmap_data %>% 
    filter(gene_set == gs) %>% 
    rownames()
  
  genes_to_plot <- intersect(genes_in_set, matched_genes)
  
  if(length(genes_to_plot) < 2) {
    warning(paste("Skipping", gs, "- not enough matched genes"))
    next
  }
  
  my_mat <- my_heatmap_data[genes_to_plot, ] %>% 
    select(-gene_set) %>% 
    as.data.frame()
  
  gurvich_mat <- gurvich_heatmap_data[genes_to_plot, ] %>% 
    select(-gene_set) %>% 
    as.data.frame()
  
  my_mat <- my_mat[, order(as.numeric(gsub("t", "", colnames(my_mat))))]
  gurvich_mat <- gurvich_mat[, order(as.numeric(gsub("t", "", colnames(gurvich_mat))))]
  
  my_mat <- as.matrix(my_mat)
  gurvich_mat <- as.matrix(gurvich_mat)
  
  #remove rows with any NA/Inf
  valid_rows <- rowSums(is.na(my_mat) | is.infinite(my_mat)) == 0 & 
                rowSums(is.na(gurvich_mat) | is.infinite(gurvich_mat)) == 0
  
  my_mat <- my_mat[valid_rows, ]
  gurvich_mat <- gurvich_mat[valid_rows, ]
  
  if(nrow(my_mat) < 2) {
    warning(paste("Skipping", gs, "- not enough valid genes after filtering"))
    next
  }
  
  #display inline first
  par(mfrow = c(1, 2))
  
  pheatmap(my_mat,
           color = colorRampPalette(c("blue", "white", "red"))(100),
           breaks = seq(-4, 4, length.out = 101),
           cluster_rows = TRUE,
           cluster_cols = FALSE,
           main = paste("My Data -", gs),
           fontsize_row = 6,
           cellwidth = 15)
  
  pheatmap(gurvich_mat,
           color = colorRampPalette(c("blue", "white", "red"))(100),
           breaks = seq(-4, 4, length.out = 101),
           cluster_rows = TRUE,
           cluster_cols = FALSE,
           main = paste("Gurvich 2017 -", gs),
           fontsize_row = 6,
           cellwidth = 15)
  
  par(mfrow = c(1, 1))
  
  #then save to PDF
  pdf(file.path(output_dir, "03_timeseries", paste0("heatmap_comparison_", gs, ".pdf")), 
      width = 14, height = 10)
  
  par(mfrow = c(1, 2))
  
  pheatmap(my_mat,
           color = colorRampPalette(c("blue", "white", "red"))(100),
           breaks = seq(-4, 4, length.out = 101),
           cluster_rows = TRUE,
           cluster_cols = FALSE,
           main = paste("My Data -", gs),
           fontsize_row = 6,
           cellwidth = 15)
  
  pheatmap(gurvich_mat,
           color = colorRampPalette(c("blue", "white", "red"))(100),
           breaks = seq(-4, 4, length.out = 101),
           cluster_rows = TRUE,
           cluster_cols = FALSE,
           main = paste("Gurvich 2017 -", gs),
           fontsize_row = 6,
           cellwidth = 15)
  
  dev.off()
}

Heatmap Analysis

Methods

Expression data were transformed from long to wide format (genes × timepoints) using pivot_wider() and stratified by gene set annotation (PHO, Stress, Ribosomal). Heatmaps were generated using the pheatmap R package with a diverging blue-white-red color gradient spanning -4 to +4 log2 fold change. Genes were hierarchically clustered to group similar temporal patterns, while timepoints were displayed in chronological order without clustering. For direct comparisons between datasets, only genes present in both GRE Lab and Gurvich datasets with complete fold change values (no NAs) across all timepoints were included. Side-by-side heatmaps were generated to visualize temporal concordance of gene expression patterns across studies.

PHO Regulon

Across the canonical PHO genes, the two datasets show similar temporal trajectories. Both exhibit rapid induction beginning between 1 h and 2 h after transfer to phosphate-fre medium, followed by sustained high expression through 8 h. Minor differences are visible at intermediate timepoints (e.g., 4 h vs 5 h), likely reflecting differences in sampling resolution and normalization strategy rather than biological divergence. The overall pattern confirms a conserved activation of the PHO regulon across independent experiments.

Ribosomal-Protein Genes

Ribosomal-protein transcripts display the expected global repression under phosphate starvation in both datasets, but with modest temporal offsets.Down-regulation is first detectable between 1.5–2 h in the Gurvich data and by 3–3.5 h in our data, suggesting slightly slower kinetics in our experimental conditions. While the overall temporal pattern is conserved, the hierarchical clustering differs, indicating quantitative differences in which specific ribosomal genes show strongest repression (consistent with moderate rank correlation ρ=0.73) Despite these differences, there is support for this strong repression signal under phosphate limited stress in both datasets.

Stress Genes

Stress genes show inconsistent concordance (ρ fluctuates 0.0-0.65) despite high directional agreement (89%), suggesting different stress gene subsets dominate in each experiment or genuine timing differences in stress pathway activation.

Collectively, these comparisons demonstrate that our RNA-seq dataset reproduces the principal transcriptional architecture of the S. cerevisiae phosphate-starvation program described by Gurvich et al. (2017). Directional concordance across hundreds of genes and consistent timing of PHO activation and ribosomal repression indicate that both studies capture the same underlying biological phenomenon. Quantitative discrepancies in fold-change amplitude are attributable to normalization and library-specific effects rather than fundamental biological differences.

Important Note:

23 PHO genes were extracted from the gene set presented in Gurvich et al. 4 genes out of 23 total genes were labeled as -Pi responding, without a clear connection. The following labels were extracted from Saccharomyces Gene Database (SGD):

• Dur3: Responsive to Nitrogen depletion. Acts as a membrane transporter for urea and polyamines
• SNA3: Nutrient responsive via the multivesicular body pathway for protein transport to vacuole
• PCA1: Cadmium transporting ATPase. Responsive to Copper and Iron depletion.
• ISC10: MAPK inhibitor in meiosis. Gene required for sporulation.

These 4 genes also showed the weakest Spearman’s Rho. Additionally, it should be noted that the initial discordance from 0min to 2hrs is exacerbated by the finer granularity of the GRE Lab S. cerevisiae phosphate depletion timecourse.

#prepare data for statistical analysis
concordance_data <- correlation_data %>%
  select(standard_name, gene_set, my_time, gurvich_time, match_type, 
         my_log2FC, gurvich_log2FC) %>%
  filter(!is.na(my_log2FC) & !is.na(gurvich_log2FC))

#per-timepoint direction concordance
direction_by_timepoint <- concordance_data %>%
  mutate(concordant = sign(my_log2FC) == sign(gurvich_log2FC)) %>%
  group_by(my_time, match_type) %>%
  summarise(
    n = n(),
    n_concordant = sum(concordant),
    pct_concordant = 100 * n_concordant / n,
    .groups = "drop"
  ) %>%
  rowwise() %>%
  mutate(binom_pval = binom.test(n_concordant, n, p = 0.5, alternative = "greater")$p.value) %>%
  ungroup() %>%
  rename(timepoint = my_time)

#per-timepoint rank concordance
rank_by_timepoint <- map_dfr(unique(concordance_data$my_time), function(tp) {
  data_tp <- concordance_data %>% filter(my_time == tp)
  
  spearman_tp <- cor.test(data_tp$my_log2FC, data_tp$gurvich_log2FC, method = "spearman", exact = FALSE)
  kendall_tp <- cor.test(data_tp$my_log2FC, data_tp$gurvich_log2FC, method = "kendall")
  pearson_tp <- cor.test(data_tp$my_log2FC, data_tp$gurvich_log2FC, method = "pearson")
  
  data.frame(
    timepoint = tp,
    n = nrow(data_tp),
    pearson_r = pearson_tp$estimate,
    pearson_pval = pearson_tp$p.value,
    spearman_rho = spearman_tp$estimate,
    spearman_pval = spearman_tp$p.value,
    kendall_tau = kendall_tp$estimate,
    kendall_pval = kendall_tp$p.value
  )
})

#per-timepoint, per-geneset direction concordance
direction_by_timepoint_geneset <- concordance_data %>%
  mutate(concordant = sign(my_log2FC) == sign(gurvich_log2FC)) %>%
  group_by(my_time, gene_set) %>%
  summarise(
    n = n(),
    n_concordant = sum(concordant),
    pct_concordant = 100 * n_concordant / n,
    .groups = "drop"
  ) %>%
  rowwise() %>%
  mutate(binom_pval = binom.test(n_concordant, n, p = 0.5, alternative = "greater")$p.value) %>%
  ungroup() %>%
  rename(timepoint = my_time)

#per-timepoint, per-geneset rank concordance
rank_by_timepoint_geneset <- concordance_data %>%
  group_by(my_time, gene_set) %>%
  summarise(
    n = n(),
    pearson_r = ifelse(n() >= 3, cor(my_log2FC, gurvich_log2FC, method = "pearson"), NA),
    spearman_rho = ifelse(n() >= 3, cor(my_log2FC, gurvich_log2FC, method = "spearman"), NA),
    .groups = "drop"
  ) %>%
  rename(timepoint = my_time)

#save outputs
write.csv(direction_by_timepoint,
          file.path(output_dir, "04_statistics", "direction_by_timepoint.csv"),
          row.names = FALSE)

write.csv(rank_by_timepoint,
          file.path(output_dir, "04_statistics", "rank_by_timepoint.csv"),
          row.names = FALSE)

write.csv(direction_by_timepoint_geneset,
          file.path(output_dir, "04_statistics", "direction_by_timepoint_geneset.csv"),
          row.names = FALSE)

write.csv(rank_by_timepoint_geneset,
          file.path(output_dir, "04_statistics", "rank_by_timepoint_geneset.csv"),
          row.names = FALSE)

#visualization: direction concordance over time
direction_time_plot <- ggplot(direction_by_timepoint, aes(x = timepoint, y = pct_concordant)) +
  geom_line(linewidth = 1.2, color = "steelblue") +
  geom_point(aes(shape = match_type), size = 4, color = "steelblue") +
  geom_hline(yintercept = 50, linetype = "dashed", color = "red") +
  geom_text(aes(label = sprintf("%.1f%%", pct_concordant)), vjust = -1, size = 3) +
  scale_shape_manual(values = c("exact" = 16, "approximate" = 1)) +
  labs(title = "Direction Concordance Over Time",
       subtitle = "Red line = 50% (chance expectation)",
       x = "Time (hours)", y = "% Concordant Direction",
       shape = "Match Type") +
  coord_cartesian(ylim = c(40, 100)) +
  theme_bw()
print(direction_time_plot)

ggsave(file.path(output_dir, "04_statistics", "direction_time_plot.pdf"),
       direction_time_plot, width = 10, height = 6)

#visualization: rank concordance over time
rank_time_plot <- ggplot(rank_by_timepoint, aes(x = timepoint, y = spearman_rho)) +
  geom_line(linewidth = 1.2, color = "darkgreen") +
  geom_point(size = 4, color = "darkgreen") +
  geom_hline(yintercept = c(0.7, 0.8, 0.9), linetype = "dashed", color = "gray50") +
  geom_text(aes(label = sprintf("%.3f", spearman_rho)), vjust = -1, size = 3) +
  labs(title = "Rank Concordance Over Time (Spearman rho)",
       x = "Time (hours)", y = "Spearman's Rho") +
  coord_cartesian(ylim = c(0, 1)) +
  theme_bw()
print(rank_time_plot)

ggsave(file.path(output_dir, "04_statistics", "rank_time_plot.pdf"),
       rank_time_plot, width = 10, height = 6)

#visualization: direction by gene set over time
direction_geneset_plot <- ggplot(direction_by_timepoint_geneset, 
                                  aes(x = timepoint, y = pct_concordant, color = gene_set)) +
  geom_line(linewidth = 1) +
  geom_point(size = 3) +
  geom_hline(yintercept = 50, linetype = "dashed", color = "gray30") +
  scale_color_manual(values = c("PHO" = "#E41A1C", "Stress" = "#377EB8", "Ribosomal" = "#4DAF4A")) +
  labs(title = "Direction Concordance by Gene Set",
       x = "Time (hours)", y = "% Concordant Direction",
       color = "Gene Set") +
  theme_bw()
print(direction_geneset_plot)

ggsave(file.path(output_dir, "04_statistics", "direction_geneset_plot.pdf"),
       direction_geneset_plot, width = 10, height = 6)

#visualization: rank by gene set over time
rank_geneset_plot <- ggplot(rank_by_timepoint_geneset, 
                             aes(x = timepoint, y = spearman_rho, color = gene_set)) +
  geom_line(linewidth = 1) +
  geom_point(size = 3) +
  geom_hline(yintercept = c(0.7, 0.8, 0.9), linetype = "dashed", color = "gray50") +
  scale_color_manual(values = c("PHO" = "#E41A1C", "Stress" = "#377EB8", "Ribosomal" = "#4DAF4A")) +
  labs(title = "Rank Concordance by Gene Set (Spearman rho)",
       x = "Time (hours)", y = "Spearman's Rho",
       color = "Gene Set") +
  theme_bw()
print(rank_geneset_plot)

ggsave(file.path(output_dir, "04_statistics", "rank_geneset_plot.pdf"),
       rank_geneset_plot, width = 10, height = 6)

#summary statistics table
concordance_summary <- rank_by_timepoint %>%
  left_join(direction_by_timepoint %>% select(timepoint, pct_concordant, binom_pval),
            by = "timepoint") %>%
  select(timepoint, n, pct_concordant, direction_pval = binom_pval,
         pearson_r, spearman_rho, kendall_tau)

write.csv(concordance_summary,
          file.path(output_dir, "04_statistics", "concordance_summary_timecourse.csv"),
          row.names = FALSE)

Gurvich, 2017 Statistical Analysis

Methods

Direction concordance was assessed using a sign test to determine whether genes changed in the same direction (both upregulated or both downregulated) between datasets. For each gene-timepoint pair, directional agreement was defined as sign(log2FC_GRE) == sign(log2FC_Gurvich), where the sign() function returns +1 for positive values, -1 for negative values, and 0 for zero. The percentage of concordant pairs was calculated as (number of concordant pairs / total pairs) × 100. Statistical significance of directional concordance was evaluated using a one-tailed binomial test (H₀: proportion concordant = 0.5, H₁: proportion concordant > 0.5) implemented via R’s binom.test(). I chose this test due to its utility for binary data, such as our sign assignment here. Under the null hypothesis, random directional agreement would occur 50% of the time. P-values less than 0.05 indicated that observed concordance exceeded chance expectations. Both correlation and direction concordance analyses were stratified by timepoint and gene set to identify temporal and functional patterns in dataset agreement.

Results

Direction concordance, assessed using a binomial test comparing sign agreement to chance expectation (50%), showed strong temporal dependence. At the earliest timepoint (1h), only 56.5% of gene-timepoint pairs showed concordant direction (p=0.04), barely exceeding chance. However, concordance improved dramatically by 2h (97.2%, p<1e-10) and remained consistently high at subsequent timepoints (98.1-100%, all p<1e-20). When stratified by gene set, all three functional categories showed significant directional agreement exceeding 88% (PHO: 95.2%, p=3.7e-38; Ribosomal: 91.4%, p=6.5e-157; Stress: 88.8%, p=1.7e-36).

Overall rank concordance (Spearman ρ) showed a similar temporal pattern, with poor correlation at 1h (ρ=0.160, p=0.019) improving to excellent concordance by 2h (ρ=0.821, p<1e-50) and remaining high through 8h (ρ=0.835-0.911, all p<1e-80).

When analyzed separately by gene set, distinct patterns emerged: PHO regulon: Maintained relatively stable rank concordance throughout the time course (ρ=0.54 at 1h, increasing to ρ=0.75 at 2h and plateauing at ρ=0.74-0.76 thereafter; all p<1e-5). Ribosomal genes: Showed the most dramatic improvement, starting with poor concordance (ρ=0.18 at 1h, p=0.02) but rapidly increasing to moderate levels by 2h (ρ=0.48, p<1e-8) and reaching excellent concordance by 3.5h (ρ=0.71, p<1e-30), which further improved to ρ=0.75-0.78 at late timepoints (all p<1e-35). Stress genes: Displayed highly erratic concordance across the time course. Starting near zero at 1h (ρ=-0.06, p=0.67), concordance spiked to ρ=0.75 at 2h (p<1e-8), crashed to ρ=0.38 at 3.5h (p=2e-4), and fluctuated between ρ=0.50-0.67 at later timepoints (all p<0.01).

Major Conclusions

The differential temporal kinetics observed across gene sets reveal a fundamental distinction between biological commitment and quantitative maturation in transcriptional stress responses. While directional concordance remains high across all gene sets at early timepoints (>90%), magnitude correlations differ by functional category: PHO genes achieve strong concordance by 1-2h (ρ > 0.8), whereas ribosomal genes show weak rank correlation at 1h despite nearly perfect directional agreement. This pattern suggests that cellular decision-making operates in discrete phases. First, transcriptional commitment establishes which genes respond and in which direction—-a binary, rapidly executed decision evident within the first hour. Second, quantitative maturation refines the relative magnitudes of these responses as feedback loops engage, mRNA stability effects propagate, and post-transcriptional mechanisms equilibrate. The consistency of this hierarchical temporal pattern across both GRE Lab (2023) and Gurvich (2017) datasets suggests that early magnitude variability reflects an intrinsic feature of transcriptional program maturation rather than technical artifact. Pathway-specific metabolic responses like PHO activation mature rapidly, while systemic adjustments like ribosomal repression and general stress responses require additional time for quantitative coordination. # Zhou & O’shea, 2011 Analysis

Zhou & O’Shea (2011) Data Processing

Methods

Zhou & O’Shea (2011) microarray data were retrieved from GEO (accession GSE23580) using the GEOquery R package. The dataset contains two experimental designs: (1) mutant cycle analysis for epistasis expression analysis using cyclic comparisons of regulatory mutants (cbf1Δ, pho4Δ, tye7Δ rtg3Δ) to dissect transcriptional network architecture, and (2) direct wild-type comparisons between no-phosphate and high-phosphate conditions. Only the direct wild-type comparison samples (n=4 biological replicates with dye swaps) were retained to match the phosphate starvation conditions used in our study and Gurvich et al. Expression values were provided as pre-processed log2 ratios (No-Pi/High-Pi) from two-color microarrays. For genes represented by multiple probes, the median fold change across probes was used as the gene-level estimate. Mean log2 fold changes were calculated across the four replicates. Data quality was validated by confirming strong upregulation of canonical PHO pathway genes (PHO5, PHO84, PHO89, etc.; mean log2FC > 2.0), consistent with published results. Zhou fold changes were then merged with our 1h timepoint data for comparative analysis based on gene symbols.

#===== Define gene sets =====

pho_genes <- gene_sets_mapped %>%
  filter(gene_set == "PHO") %>%
  pull(standard_name) %>%
  unique()

stress_genes <- gene_sets_mapped %>%
  filter(gene_set == "Stress") %>%
  pull(standard_name) %>%
  unique()

ribosomal_genes <- gene_sets_mapped %>%
  filter(gene_set == "Ribosomal") %>%
  pull(standard_name) %>%
  unique()

cat(sprintf("Gene sets defined: PHO (%d), Stress (%d), Ribosomal (%d)\n\n",
            length(pho_genes), length(stress_genes), length(ribosomal_genes)))

## Gene sets defined: PHO (32), Stress (44), Ribosomal (150)

#===== Zhou vs GRE Lab comparison at 1h =====

cat("===== Zhou (2011) vs GRE Lab (2024) Comparison =====\n\n")

## ===== Zhou (2011) vs GRE Lab (2024) Comparison =====

#merge zhou with your 1h data
comparison_zhou <- inner_join(
  zhou_fc_bygene,
  expr_fc_summary %>% filter(timepoint_hours == 1),
  by = c("GENE_SYMBOL" = "standard_name")
) %>%
  mutate(
    zhou_direction = sign(log2FC_zhou),
    your_direction = sign(mean_log2FC),
    concordant_direction = zhou_direction == your_direction
  )

cat(sprintf("Total overlapping genes: %d\n\n", nrow(comparison_zhou)))

## Total overlapping genes: 4356

#===== Genome-wide correlation =====

cat("Genome-wide correlation (all genes):\n")

## Genome-wide correlation (all genes):

spear_all <- cor.test(comparison_zhou$log2FC_zhou, 
                      comparison_zhou$mean_log2FC, 
                      method = "spearman")

direction_all <- sum(comparison_zhou$concordant_direction) / nrow(comparison_zhou)

cat(sprintf("  Spearman ρ = %.3f (p = %.2e)\n", spear_all$estimate, spear_all$p.value))

##   Spearman ρ = 0.348 (p = 4.91e-124)

cat(sprintf("  Direction concordance: %.1f%%\n\n", 100 * direction_all))

##   Direction concordance: 62.4%

#===== Gene set-specific correlations =====

cat("Gene set-specific correlations:\n")

## Gene set-specific correlations:

#function to calculate correlation for a gene set
calc_geneset_correlation <- function(comparison_data, genes, geneset_name) {
  
  subset_data <- comparison_data %>%
    filter(GENE_SYMBOL %in% genes)
  
  if(nrow(subset_data) < 3) {
    cat(sprintf("  %s: Too few genes (n=%d)\n", geneset_name, nrow(subset_data)))
    return(NULL)
  }
  
  spear <- cor.test(subset_data$log2FC_zhou, subset_data$mean_log2FC, 
                    method = "spearman")
  direction_conc <- sum(subset_data$concordant_direction) / nrow(subset_data)
  
  cat(sprintf("  %s: n=%d, ρ=%.3f (p=%.2e), Direction=%.1f%%\n",
              geneset_name, nrow(subset_data), spear$estimate, 
              spear$p.value, 100 * direction_conc))
  
  return(list(
    geneset = geneset_name,
    n_genes = nrow(subset_data),
    spearman_rho = spear$estimate,
    p_value = spear$p.value,
    direction_concordance = 100 * direction_conc,
    data = subset_data
  ))
}

#calculate for each gene set
pho_result <- calc_geneset_correlation(comparison_zhou, pho_genes, "PHO")

##   PHO: n=28, ρ=0.813 (p=1.61e-06), Direction=89.3%

stress_result <- calc_geneset_correlation(comparison_zhou, stress_genes, "Stress")

##   Stress: n=30, ρ=0.282 (p=1.31e-01), Direction=93.3%

ribo_result <- calc_geneset_correlation(comparison_zhou, ribosomal_genes, "Ribosomal")

##   Ribosomal: n=147, ρ=-0.080 (p=3.38e-01), Direction=96.6%

cat("\n")

#===== Create summary table =====

zhou_geneset_summary <- tibble(
  Gene_Set = c("All genes", "PHO", "Stress", "Ribosomal"),
  N_genes = c(nrow(comparison_zhou), 
              pho_result$n_genes, 
              stress_result$n_genes, 
              ribo_result$n_genes),
  Spearman_rho = c(spear_all$estimate,
                   pho_result$spearman_rho,
                   stress_result$spearman_rho,
                   ribo_result$spearman_rho),
  P_value = c(spear_all$p.value,
              pho_result$p_value,
              stress_result$p_value,
              ribo_result$p_value),
  Direction_concordance = c(100 * direction_all,
                           pho_result$direction_concordance,
                           stress_result$direction_concordance,
                           ribo_result$direction_concordance)
)

write.csv(zhou_geneset_summary,
          file.path(output_dir, "02_correlation", "zhou_geneset_summary.csv"),
          row.names = FALSE)

print(knitr::kable(zhou_geneset_summary, digits = 3,
                   caption = "Zhou (2011) vs GRE Lab (2024): Correlation by Gene Set"))

## 
## 
## Table: Zhou (2011) vs GRE Lab (2024): Correlation by Gene Set
## 
## |Gene_Set  | N_genes| Spearman_rho| P_value| Direction_concordance|
## |:---------|-------:|------------:|-------:|---------------------:|
## |All genes |    4356|        0.348|   0.000|                62.351|
## |PHO       |      28|        0.813|   0.000|                89.286|
## |Stress    |      30|        0.282|   0.131|                93.333|
## |Ribosomal |     147|       -0.080|   0.338|                96.599|

#===== Scatter plots for each gene set =====

plot_geneset_scatter <- function(data, genes, title, geneset_name, color) {
  
  subset_data <- data %>% filter(GENE_SYMBOL %in% genes)
  
  if(nrow(subset_data) < 3) {
    cat(sprintf("Skipping scatter for %s - too few genes\n", geneset_name))
    return(NULL)
  }
  
  spear <- cor.test(subset_data$log2FC_zhou, subset_data$mean_log2FC, 
                    method = "spearman")
  dir_conc <- 100 * sum(subset_data$concordant_direction) / nrow(subset_data)
  
  p <- ggplot(subset_data, aes(x = log2FC_zhou, y = mean_log2FC)) +
    geom_point(alpha = 0.7, size = 3, color = color) +
    geom_text_repel(aes(label = GENE_SYMBOL), size = 2.5, max.overlaps = 20) +
    geom_smooth(method = "lm", se = TRUE, color = "black", linewidth = 0.8) +
    geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray50") +
    geom_hline(yintercept = 0, linetype = "dotted", color = "gray30") +
    geom_vline(xintercept = 0, linetype = "dotted", color = "gray30") +
    annotate("text", x = min(subset_data$log2FC_zhou), 
             y = max(subset_data$mean_log2FC),
             label = sprintf("Spearman ρ = %.3f\np = %.2e\nDirection: %.1f%%\nn = %d",
                            spear$estimate, spear$p.value, dir_conc, nrow(subset_data)),
             hjust = 0, vjust = 1, size = 3.5) +
    labs(title = title,
         subtitle = "Zhou (60 min) vs GRE Lab (1h)",
         x = "Zhou log2FC", 
         y = "GRE Lab log2FC") +
    theme_bw()
  
  return(p)
}

#genome-wide scatter
p_all <- ggplot(comparison_zhou, aes(x = log2FC_zhou, y = mean_log2FC)) +
  geom_point(alpha = 0.4, size = 1.5, color = "gray40") +
  geom_smooth(method = "lm", se = TRUE, color = "black", linewidth = 0.8) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray50") +
  geom_hline(yintercept = 0, linetype = "dotted", color = "gray30") +
  geom_vline(xintercept = 0, linetype = "dotted", color = "gray30") +
  annotate("text", x = min(comparison_zhou$log2FC_zhou), 
           y = max(comparison_zhou$mean_log2FC),
           label = sprintf("Spearman ρ = %.3f\np = %.2e\nDirection: %.1f%%\nn = %d",
                          spear_all$estimate, spear_all$p.value, 
                          100 * direction_all, nrow(comparison_zhou)),
           hjust = 0, vjust = 1, size = 3.5) +
  labs(title = "Zhou (60 min) vs GRE Lab (1h)",
       subtitle = "All Genes",
       x = "Zhou log2FC", y = "GRE Lab log2FC") +
  theme_bw()

print(p_all)

## `geom_smooth()` using formula = 'y ~ x'

ggsave(file.path(output_dir, "02_correlation", "zhou_scatter_all_genes.pdf"),
       p_all, width = 8, height = 8)

## `geom_smooth()` using formula = 'y ~ x'

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'Spearman ρ = 0.348' in 'mbcsToSbcs': for ρ (U+03C1)

#gene set scatters
p_pho <- plot_geneset_scatter(comparison_zhou, pho_genes, 
                               "PHO Gene Set", "PHO", "red")
p_stress <- plot_geneset_scatter(comparison_zhou, stress_genes, 
                                  "Stress Gene Set", "Stress", "orange")
p_ribo <- plot_geneset_scatter(comparison_zhou, ribosomal_genes, 
                                "Ribosomal Gene Set", "Ribosomal", "blue")

if(!is.null(p_pho)) {
  print(p_pho)
  ggsave(file.path(output_dir, "02_correlation", "zhou_scatter_PHO.pdf"),
         p_pho, width = 10, height = 8)
}

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'Spearman ρ = 0.813' in 'mbcsToSbcs': for ρ (U+03C1)

if(!is.null(p_stress)) {
  print(p_stress)
  ggsave(file.path(output_dir, "02_correlation", "zhou_scatter_Stress.pdf"),
         p_stress, width = 10, height = 8)
}

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'Spearman ρ = 0.282' in 'mbcsToSbcs': for ρ (U+03C1)

if(!is.null(p_ribo)) {
  print(p_ribo)
  ggsave(file.path(output_dir, "02_correlation", "zhou_scatter_Ribosomal.pdf"),
         p_ribo, width = 10, height = 8)
}

## `geom_smooth()` using formula = 'y ~ x'

## Warning: ggrepel: 71 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## `geom_smooth()` using formula = 'y ~ x'

## Warning: ggrepel: 13 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'Spearman ρ = -0.080' in 'mbcsToSbcs': for ρ (U+03C1)

#===== Combined gene set plot =====

p_combined <- ggplot(comparison_zhou, aes(x = log2FC_zhou, y = mean_log2FC)) +
  geom_point(alpha = 0.2, size = 1, color = "gray70") +
  geom_point(data = comparison_zhou %>% filter(GENE_SYMBOL %in% pho_genes),
             aes(color = "PHO"), size = 2.5, alpha = 0.8) +
  geom_point(data = comparison_zhou %>% filter(GENE_SYMBOL %in% stress_genes),
             aes(color = "Stress"), size = 2.5, alpha = 0.8) +
  geom_point(data = comparison_zhou %>% filter(GENE_SYMBOL %in% ribosomal_genes),
             aes(color = "Ribosomal"), size = 2.5, alpha = 0.8) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray50") +
  scale_color_manual(values = c("PHO" = "red", "Stress" = "orange", "Ribosomal" = "blue")) +
  labs(title = "Zhou (60 min) vs GRE Lab (1h): Gene Sets Highlighted",
       x = "Zhou log2FC", y = "GRE Lab log2FC", color = "Gene Set") +
  theme_bw() +
  theme(legend.position = "right")

print(p_combined)

ggsave(file.path(output_dir, "02_correlation", "zhou_scatter_combined_genesets.pdf"),
       p_combined, width = 10, height = 8)

cat("\nNote: Zhou microarray shows compressed fold changes relative to RNA-seq.\n")

## 
## Note: Zhou microarray shows compressed fold changes relative to RNA-seq.

cat("Analysis complete.\n")

## Analysis complete.

Zhou & O’Shea (2011) Comparative Analysis

Methods

Zhou fold changes (60 min no-Pi vs high-Pi) were compared to GRE Lab fold changes at the 1h timepoint. Spearman rank correlation (ρ) was calculated to assess magnitude concordance, accounting for potential non-linear relationships between microarray and RNA-seq platforms. Directional concordance was quantified as the percentage of genes showing the same sign of fold change (both upregulated or both downregulated) across datasets. Correlations were calculated genome-wide (all overlapping genes) and for each functional gene set (PHO, Stress, Ribosomal). Scatter plots were generated with gene symbols labeled for gene set-specific comparisons, including linear regression fits and identity lines (slope = 1) for visual reference. To identify optimal temporal alignment between datasets, Zhou 1h data were correlated against GRE Lab timepoints from 0.5-8h, and the timepoint with maximum Spearman correlation was identified.

Results

Genome-wide concordance: Zhou (2011) 1h fold changes showed moderate correlation with GRE Lab (2024) 1h data across all overlapping genes (n=4,356; Spearman ρ=0.348, p=4.9×10⁻¹²⁴), with 62.4% of genes exhibiting concordant directional changes.

Gene set-specific validation: Analysis of functional gene sets recapitulated the hierarchical temporal kinetics observed in the Gurvich comparison. PHO pathway genes (n=28) demonstrated strong magnitude correlation (ρ=0.813, p=1.6×10⁻⁶) with 89.3% directional concordance, consistent with rapid transcriptional maturation of the core phosphate starvation response by 1h. In contrast, Stress genes (n=30; ρ=0.282, p=0.13) and Ribosomal genes (n=147; ρ=-0.080, p=0.34) showed weak-to-absent magnitude correlations despite maintaining high directional concordance (93.3% and 96.6%, respectively).

Cross-study validation and sources of early variability: The pattern of strong directional agreement coupled with poor magnitude correlation at early timepoints reflects both intrinsic biological immaturity of transcriptional programs and laboratory-specific technical variability. Pre-sequencing sample handling—including precise harvest timing, cell processing speed, growth phase at starvation onset, and RNA stabilization methods can influence the initial transcriptional state of cells, introducing baseline differences that disproportionately bias early timepoint measurements. The conservation of this hierarchical temporal structure across three independent studies spanning 13 years (Zhou 2011, Gurvich 2017, GRE Lab 2023) and two expression profiling platforms (microarray, RNA-seq) demonstrates that while absolute magnitudes at 1h are variable, the underlying biology of transcriptional commitment (directional response) is robust and reproducible across laboratories. Fast-responding pathways like PHO overcome these initial state differences rapidly, while slower systemic responses require additional time for quantitative coordination to emerge above lab-specific baseline variation.

#recreate zhou figure 3C with proper gene name mapping

#original zhou gene list with synonyms
zhou_fig3c_genes_original <- c("PHO5", "PHO84", "SPL2", "PHO89", "GIT1", "VTC3", 
                               "VTC1", "VTC4", "YJL012C-A", "PHO11", "PHO81", "PHM6", 
                               "CTF19", "PHO86", "PHO12", "VTC2", "VIP1", "YPL110C", 
                               "DDP1", "HOR2", "PHM8", "ENA2", "YAR070C", "PHO8", 
                               "YNL217W", "ENA1", "YJL119C", "CBF1")

#map systematic names to standard names where possible
gene_name_map <- tribble(
  ~systematic_name, ~standard_name,
  "YNL217W", "PPN2",
  "YPL110C", "GDE1",
  "HOR2", "GPP2",
  "YJL012C-A", "YJL012C-A",  #absent, keep as-is
  "YAR070C", "YAR070C",      #absent, keep as-is
  "YJL119C", "YJL119C"       #absent, keep as-is
)

#expand gene list to include both names
zhou_fig3c_genes_expanded <- unique(c(
  zhou_fig3c_genes_original,
  gene_name_map$standard_name
))

#extract zhou PROBE-level data using expanded gene list
zhou_fig3c_probes <- zhou_fc %>%
  filter(!is.na(GENE_SYMBOL), GENE_SYMBOL %in% zhou_fig3c_genes_expanded) %>%
  select(probe_id, GENE_SYMBOL, log2FC_zhou) %>%
  #map back to original names for display
  mutate(display_name = case_when(
    GENE_SYMBOL == "PPN2" ~ "YNL217W",
    GENE_SYMBOL == "GDE1" ~ "YPL110C",
    GENE_SYMBOL == "GPP2" ~ "HOR2",
    TRUE ~ GENE_SYMBOL
  ))

#calculate mean fold change per gene for ordering
zhou_gene_order <- zhou_fig3c_probes %>%
  group_by(display_name) %>%
  summarise(mean_fc = mean(log2FC_zhou, na.rm = TRUE), .groups = "drop") %>%
  arrange(desc(mean_fc))

#order probes by gene rank
zhou_fig3c_probes <- zhou_fig3c_probes %>%
  mutate(
    gene_rank = match(display_name, zhou_gene_order$display_name),
    probe_label = paste0(display_name, " [A_06_", probe_id, "]")
  ) %>%
  arrange(gene_rank, probe_id) %>%
  mutate(probe_label = factor(probe_label, levels = probe_label))

cat(sprintf("Zhou: %d probes for %d genes\n", 
            nrow(zhou_fig3c_probes), 
            length(unique(zhou_fig3c_probes$display_name))))

## Zhou: 28 probes for 28 genes

#extract GRE data using both standard and systematic names
gre_fig3c_data <- expr_fc_summary %>%
  filter(timepoint_hours == 1) %>%
  filter(standard_name %in% zhou_fig3c_genes_expanded | 
         systematic_name %in% zhou_fig3c_genes_original) %>%
  select(standard_name, systematic_name, mean_log2FC) %>%
  mutate(display_name = case_when(
    standard_name == "PPN2" ~ "YNL217W",
    standard_name == "GDE1" ~ "YPL110C",
    standard_name == "GPP2" ~ "HOR2",
    systematic_name %in% zhou_fig3c_genes_original ~ systematic_name,
    TRUE ~ standard_name
  )) %>%
  select(display_name, mean_log2FC) %>%
  rename(log2FC_GRE = mean_log2FC)

cat(sprintf("GRE: %d genes found\n", nrow(gre_fig3c_data)))

## GRE: 25 genes found

#expand GRE data to match probe structure
gre_fig3c_expanded <- zhou_fig3c_probes %>%
  select(probe_label, display_name) %>%
  left_join(gre_fig3c_data, by = "display_name")

#check rescue success
rescued <- setdiff(
  intersect(c("YNL217W", "YPL110C", "HOR2"), gre_fig3c_data$display_name),
  character(0)
)

if(length(rescued) > 0) {
  cat(sprintf("\nRescued genes via name mapping: %s\n", paste(rescued, collapse = ", ")))
}

## 
## Rescued genes via name mapping: YNL217W, YPL110C, HOR2

still_missing <- setdiff(zhou_fig3c_genes_original, 
                         c(zhou_fig3c_probes$display_name, gre_fig3c_data$display_name))
if(length(still_missing) > 0) {
  cat(sprintf("Still missing: %s\n", paste(still_missing, collapse = ", ")))
}

#clipping statistics
zhou_clipped <- sum(abs(zhou_fig3c_probes$log2FC_zhou) > 2, na.rm = TRUE)
gre_clipped <- sum(abs(gre_fig3c_expanded$log2FC_GRE) > 2, na.rm = TRUE)

cat(sprintf("\nValues clipped by -2 to +2 scale:\n"))

## 
## Values clipped by -2 to +2 scale:

cat(sprintf("  Zhou: %d/%d probes (%.1f%%)\n", 
            zhou_clipped, nrow(zhou_fig3c_probes),
            100 * zhou_clipped / nrow(zhou_fig3c_probes)))

##   Zhou: 14/28 probes (50.0%)

cat(sprintf("  GRE: %d/%d rows (%.1f%%)\n", 
            gre_clipped, nrow(gre_fig3c_expanded),
            100 * gre_clipped / nrow(gre_fig3c_expanded)))

##   GRE: 6/28 rows (21.4%)

#zhou heatmap
zhou_heatmap <- ggplot(zhou_fig3c_probes, 
                       aes(x = "Zhou", y = probe_label, fill = log2FC_zhou)) +
  geom_tile(color = "white", linewidth = 0.3) +
  scale_fill_gradient2(low = "green", mid = "black", high = "red",
                       midpoint = 0, limits = c(-2, 2), oob = scales::squish,
                       name = "log2 FC") +
  labs(title = "Zhou et al. (2011)", subtitle = "Microarray probes, 60 min", 
       x = NULL, y = NULL) +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 6, family = "mono"),
        axis.text.x = element_blank(),
        panel.grid = element_blank())

#GRE heatmap
gre_heatmap <- ggplot(gre_fig3c_expanded, 
                      aes(x = "GRE Lab", y = probe_label, fill = log2FC_GRE)) +
  geom_tile(color = "white", linewidth = 0.3) +
  scale_fill_gradient2(low = "green", mid = "black", high = "red",
                       midpoint = 0, limits = c(-2, 2), oob = scales::squish,
                       na.value = "grey80", name = "log2 FC") +
  labs(title = "GRE Lab", subtitle = "RNA-seq genes, 60 min",
       x = NULL, y = NULL) +
  theme_minimal() +
  theme(axis.text.y = element_blank(), 
        axis.text.x = element_blank(),
        panel.grid = element_blank())

#combine
library(patchwork)
fig3c_combined <- zhou_heatmap + gre_heatmap + plot_layout(widths = c(1.2, 1))

print(fig3c_combined)

ggsave(file.path(output_dir, "02_correlation", "zhou_fig3c_comparison_final.pdf"),
       fig3c_combined, width = 10, height = 12)

#correlation statistics
zhou_gene_fc <- zhou_fig3c_probes %>%
  group_by(display_name) %>%
  summarise(log2FC_zhou = median(log2FC_zhou, na.rm = TRUE), .groups = "drop")

fig3c_cor <- inner_join(zhou_gene_fc, gre_fig3c_data, by = "display_name")

fig3c_cor_test <- cor.test(fig3c_cor$log2FC_zhou, fig3c_cor$log2FC_GRE, 
                            method = "spearman")

cat(sprintf("\nFigure 3C genes correlation: ρ = %.3f (p = %.2e, n = %d)\n",
            fig3c_cor_test$estimate, fig3c_cor_test$p.value, nrow(fig3c_cor)))

## 
## Figure 3C genes correlation: ρ = 0.762 (p = 1.64e-05, n = 25)

cat(sprintf("Direction concordance: %.1f%%\n",
            100 * sum(sign(fig3c_cor$log2FC_zhou) == sign(fig3c_cor$log2FC_GRE)) / nrow(fig3c_cor)))

## Direction concordance: 96.0%

write.csv(fig3c_cor,
          file.path(output_dir, "02_correlation", "zhou_fig3c_gene_comparison.csv"),
          row.names = FALSE)

Zhou Figure 3C Recreation

Rationale

To directly compare our results with the published Zhou et al. (2011) findings, I recreated their Figure 3C using the same 28 genes analyzed in their study. Zhou probe-level fold changes (60 min no-Pi vs high-Pi) were plotted alongside GRE Lab gene-level fold changes (1h vs 0h) using matched gene identifiers, accounting for systematic name to standard name conversions where necessary (e.g., YNL217W→PPN2, YPL110C→GDE1). Both datasets were visualized using identical color scales (green-black-red, -2 to +2 log2FC) with values outside this range clipped to facilitate visual comparison, following the original publication’s conventions.

Results

Visual concordance: Side-by-side heatmaps revealed strong directional concordance across the phosphate starvation response signature. Canonical PHO pathway genes (PHO5, PHO11, PHO84, PHO89, VTC1, VTC3, VTC4, PHO81, etc.) showed consistent upregulation (red) in both datasets, validating that the core transcriptional response is conserved between studies. Genes with weaker or negative responses in Zhou’s data similarly showed attenuated responses in GRE Lab data. Quantitative correlation analysis of the Figure 3C gene set confirmed moderate-to-strong concordance (Spearman ρ=0.813, p=1.6×10⁻⁶, n=28), with 89.3% directional agreement.

Note

Platform considerations: While directional patterns were highly concordant, absolute fold change magnitudes differed between platforms due to the compressed dynamic range characteristic of two-color microarrays relative to RNA-seq. The -2 to +2 clipping range used in the original Zhou figure captured the majority of microarray fold changes but truncated several RNA-seq values that exceeded this range, particularly for highly induced genes like PHO5 and PHO84. This compression does not affect directional or rank-order comparisons but precludes direct magnitude-to-magnitude interpretation across platforms.

#analyze Zhou concordance for Stress and Ribosomal gene sets

#get gene lists from mapped gene sets
stress_genes <- gene_sets_mapped %>%
  filter(gene_set == "Stress") %>%
  pull(standard_name) %>%
  unique()

ribosomal_genes <- gene_sets_mapped %>%
  filter(gene_set == "Ribosomal") %>%
  pull(standard_name) %>%
  unique()

pho_genes <- gene_sets_mapped %>%
  filter(gene_set == "PHO") %>%
  pull(standard_name) %>%
  unique()

cat(sprintf("Gene set sizes: PHO = %d, Stress = %d, Ribosomal = %d\n",
            length(pho_genes), length(stress_genes), length(ribosomal_genes)))

## Gene set sizes: PHO = 32, Stress = 44, Ribosomal = 150

#extract comparisons for each gene set
comparison_zhou_stress <- comparison_zhou %>%
  filter(GENE_SYMBOL %in% stress_genes)

comparison_zhou_ribo <- comparison_zhou %>%
  filter(GENE_SYMBOL %in% ribosomal_genes)

comparison_zhou_pho <- comparison_zhou %>%
  filter(GENE_SYMBOL %in% pho_genes)

cat(sprintf("\nOverlapping genes with Zhou data:\n"))

## 
## Overlapping genes with Zhou data:

cat(sprintf("  PHO: %d\n", nrow(comparison_zhou_pho)))

##   PHO: 28

cat(sprintf("  Stress: %d\n", nrow(comparison_zhou_stress)))

##   Stress: 30

cat(sprintf("  Ribosomal: %d\n", nrow(comparison_zhou_ribo)))

##   Ribosomal: 147

#correlation and direction concordance for each gene set
analyze_geneset <- function(data, geneset_name) {
  if(nrow(data) < 5) {
    cat(sprintf("\n%s: Too few genes (n=%d) for analysis\n", geneset_name, nrow(data)))
    return(list(n = nrow(data), rho = NA, p = NA, direction = NA))
  }
  
  spear <- cor.test(data$log2FC_zhou, data$mean_log2FC, method = "spearman")
  direction_conc <- 100 * sum(data$concordant_direction) / nrow(data)
  
  cat(sprintf("\n%s Gene Set (n = %d):\n", geneset_name, nrow(data)))
  cat(sprintf("  Spearman ρ = %.3f (p = %.2e)\n", spear$estimate, spear$p.value))
  cat(sprintf("  Direction concordance: %.1f%%\n", direction_conc))
  
  return(list(
    n = nrow(data),
    rho = spear$estimate,
    p = spear$p.value,
    direction = direction_conc
  ))
}

pho_stats <- analyze_geneset(comparison_zhou_pho, "PHO")

## 
## PHO Gene Set (n = 28):
##   Spearman ρ = 0.813 (p = 1.61e-06)
##   Direction concordance: 89.3%

stress_stats <- analyze_geneset(comparison_zhou_stress, "Stress")

## 
## Stress Gene Set (n = 30):
##   Spearman ρ = 0.282 (p = 1.31e-01)
##   Direction concordance: 93.3%

ribo_stats <- analyze_geneset(comparison_zhou_ribo, "Ribosomal")

## 
## Ribosomal Gene Set (n = 147):
##   Spearman ρ = -0.080 (p = 3.38e-01)
##   Direction concordance: 96.6%

#create summary table
geneset_summary <- tibble(
  Gene_Set = c("PHO", "Stress", "Ribosomal"),
  N_genes = c(pho_stats$n, stress_stats$n, ribo_stats$n),
  Spearman_rho = c(pho_stats$rho, stress_stats$rho, ribo_stats$rho),
  P_value = c(pho_stats$p, stress_stats$p, ribo_stats$p),
  Direction_pct = c(pho_stats$direction, stress_stats$direction, ribo_stats$direction)
)

write.csv(geneset_summary,
          file.path(output_dir, "02_correlation", "zhou_geneset_summary.csv"),
          row.names = FALSE)

print(knitr::kable(geneset_summary, digits = 3,
                   caption = "Zhou vs GRE Lab: Gene Set Analysis"))

## 
## 
## Table: Zhou vs GRE Lab: Gene Set Analysis
## 
## |Gene_Set  | N_genes| Spearman_rho| P_value| Direction_pct|
## |:---------|-------:|------------:|-------:|-------------:|
## |PHO       |      28|        0.813|   0.000|        89.286|
## |Stress    |      30|        0.282|   0.131|        93.333|
## |Ribosomal |     147|       -0.080|   0.338|        96.599|

#scatter plots for each gene set
create_geneset_scatter <- function(data, geneset_name, color) {
  if(nrow(data) < 5) return(NULL)
  
  spear <- cor.test(data$log2FC_zhou, data$mean_log2FC, method = "spearman")
  dir_conc <- 100 * sum(data$concordant_direction) / nrow(data)
  
  ggplot(data, aes(x = log2FC_zhou, y = mean_log2FC)) +
    geom_point(size = 3, color = color, alpha = 0.7) +
    geom_text_repel(aes(label = GENE_SYMBOL), size = 2.5, max.overlaps = 20) +
    geom_smooth(method = "lm", se = TRUE, color = "black", linewidth = 0.8) +
    geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray50") +
    geom_hline(yintercept = 0, linetype = "dotted", color = "gray30") +
    geom_vline(xintercept = 0, linetype = "dotted", color = "gray30") +
    annotate("text", x = -Inf, y = Inf,
             label = sprintf("Spearman ρ = %.3f\np = %.2e\nDirection: %.1f%%\nn = %d",
                            spear$estimate, spear$p.value, dir_conc, nrow(data)),
             hjust = 0, vjust = 1, size = 3.5) +
    labs(title = paste(geneset_name, "Gene Set"),
         subtitle = "Zhou (60 min) vs GRE Lab (1h)",
         x = "Zhou log2FC", y = "GRE Lab log2FC") +
    theme_bw()
}

#create scatter plots for combined display
pho_scatter <- create_geneset_scatter(comparison_zhou_pho, "PHO", "#E41A1C")
stress_scatter <- create_geneset_scatter(comparison_zhou_stress, "Stress", "#377EB8")
ribo_scatter <- create_geneset_scatter(comparison_zhou_ribo, "Ribosomal", "#4DAF4A")

#combined plot
if(!is.null(pho_scatter) && !is.null(stress_scatter) && !is.null(ribo_scatter)) {
  combined_geneset <- pho_scatter + stress_scatter + ribo_scatter + 
    plot_layout(ncol = 3)
  
  print(combined_geneset)
  
  ggsave(file.path(output_dir, "02_correlation", "zhou_all_genesets_combined.pdf"),
         combined_geneset, width = 18, height = 6)
}

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

## Warning: ggrepel: 117 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'Spearman ρ = 0.813' in 'mbcsToSbcs': for ρ (U+03C1)

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'Spearman ρ = 0.282' in 'mbcsToSbcs': for ρ (U+03C1)

## Warning: ggrepel: 65 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'Spearman ρ = -0.080' in 'mbcsToSbcs': for ρ (U+03C1)

#bar plot comparing correlations across gene sets
correlation_comparison <- ggplot(geneset_summary, 
                                 aes(x = Gene_Set, y = Spearman_rho, fill = Gene_Set)) +
  geom_col(width = 0.7, color = "black") +
  geom_text(aes(label = sprintf("ρ = %.3f\nn = %d", Spearman_rho, N_genes)),
            vjust = -0.5, size = 3.5) +
  scale_fill_manual(values = c("PHO" = "#E41A1C", "Stress" = "#377EB8", 
                                "Ribosomal" = "#4DAF4A")) +
  geom_hline(yintercept = c(0.5, 0.7, 0.8), linetype = "dashed", 
             color = "gray50", alpha = 0.5) +
  labs(title = "Zhou vs GRE Lab: Correlation by Gene Set",
       subtitle = "Spearman rank correlation at 1h",
       x = "Gene Set", y = "Spearman ρ") +
  coord_cartesian(ylim = c(0, 1)) +
  theme_bw() +
  theme(legend.position = "none")

print(correlation_comparison)

ggsave(file.path(output_dir, "02_correlation", "zhou_geneset_correlation_comparison.pdf"),
       correlation_comparison, width = 8, height = 6)

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'ρ = 0.813' in 'mbcsToSbcs': for ρ (U+03C1)

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'ρ = 0.282' in 'mbcsToSbcs': for ρ (U+03C1)

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'ρ = -0.080' in 'mbcsToSbcs': for ρ (U+03C1)

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'Spearman ρ' in 'mbcsToSbcs': for ρ (U+03C1)

#direction concordance comparison
direction_comparison <- ggplot(geneset_summary, 
                               aes(x = Gene_Set, y = Direction_pct, fill = Gene_Set)) +
  geom_col(width = 0.7, color = "black") +
  geom_text(aes(label = sprintf("%.1f%%", Direction_pct)),
            vjust = -0.5, size = 3.5) +
  scale_fill_manual(values = c("PHO" = "#E41A1C", "Stress" = "#377EB8", 
                                "Ribosomal" = "#4DAF4A")) +
  geom_hline(yintercept = 50, linetype = "dashed", color = "red") +
  labs(title = "Zhou vs GRE Lab: Direction Concordance by Gene Set",
       subtitle = "% genes with same direction at 1h",
       x = "Gene Set", y = "Direction Concordance (%)") +
  coord_cartesian(ylim = c(0, 100)) +
  theme_bw() +
  theme(legend.position = "none")

print(direction_comparison)

ggsave(file.path(output_dir, "02_correlation", "zhou_geneset_direction_comparison.pdf"),
       direction_comparison, width = 8, height = 6)

Zhou: Gene Set Stats Visualized

Methods

Gene set-specific scatter plots were generated comparing Zhou (60 min) and GRE Lab (1h) fold changes for PHO (n=28), Stress (n=30), and Ribosomal (n=147) genes. Each scatter plot included linear regression fits, identity lines (slope=1), and gene symbol labels. Spearman rank correlations and directional concordance were calculated per gene set. Summary bar plots visualized magnitude correlation (Spearman ρ) and directional concordance (% genes with same sign) across the three functional categories.

Results

Magnitude correlations revealed hierarchical concordance: PHO genes (ρ=0.813, p=1.6×10⁻⁶) showed strong correlation tracking near the identity line, Stress genes (ρ=0.282, p=0.13) exhibited moderate scatter with weak correlation, and Ribosomal genes (ρ=-0.080, p=0.34) displayed near-zero correlation despite tight clustering around consistent downregulation in GRE data.

Directional concordance remained high across all gene sets: PHO 89.3%, Stress 93.3%, Ribosomal 96.6%, with all exceeding the 50% random expectation threshold. The ribosomal scatter plot revealed a distinctive pattern where Zhou magnitudes varied widely (-0.4 to 0.2 log2FC) while GRE Lab values clustered tightly (-0.5 to -0.8 log2FC), illustrating platform compression in microarray data alongside biological immaturity at 1h. Stress genes showed intermediate behavior with broader scatter in both datasets but preserved directional trends.

#function to create heatmap for one gene set
create_threeway_pheatmap <- function(genes, geneset_name, color_limits = c(-4, 4)) {
  
  #Zhou data (1h only) - wide format
  zhou_mat <- zhou_fc_bygene %>%
    filter(GENE_SYMBOL %in% genes) %>%
    select(GENE_SYMBOL, log2FC_zhou) %>%
    column_to_rownames("GENE_SYMBOL") %>%
    as.matrix()
  
  colnames(zhou_mat) <- "Zhou_1h"
  
  #GRE data (timecourse) - wide format
  gre_mat <- expr_fc_summary %>%
    filter(standard_name %in% genes) %>%
    select(standard_name, timepoint_hours, mean_log2FC) %>%
    pivot_wider(names_from = timepoint_hours, values_from = mean_log2FC, 
                names_prefix = "GRE_") %>%
    column_to_rownames("standard_name") %>%
    as.matrix()
  
  #Gurvich data (timecourse) - wide format
  gurvich_mat <- gurvich_fc %>%
    filter(gene_set == geneset_name, timepoint_hours <= 8) %>%
    group_by(gene_name, timepoint_hours) %>%
    summarise(log2FC = mean(log2FC, na.rm = TRUE), .groups = "drop") %>%
    pivot_wider(names_from = timepoint_hours, values_from = log2FC,
                names_prefix = "Gurvich_") %>%
    column_to_rownames("gene_name") %>%
    as.matrix()
  
  #get common genes
  common_genes <- Reduce(intersect, list(rownames(zhou_mat), 
                                         rownames(gre_mat), 
                                         rownames(gurvich_mat)))
  
  if(length(common_genes) < 3) {
    cat(sprintf("Skipping %s - too few common genes\n", geneset_name))
    return(NULL)
  }
  
  #subset to common genes
  zhou_mat <- zhou_mat[common_genes, , drop = FALSE]
  gre_mat <- gre_mat[common_genes, ]
  gurvich_mat <- gurvich_mat[common_genes, ]
  
  #order genes by mean across all datasets
  gene_order <- data.frame(
    gene = common_genes,
    mean_fc = rowMeans(cbind(zhou_mat, gre_mat, gurvich_mat), na.rm = TRUE)
  ) %>%
    arrange(desc(mean_fc)) %>%
    pull(gene)
  
  #reorder matrices
  zhou_mat <- zhou_mat[gene_order, , drop = FALSE]
  gre_mat <- gre_mat[gene_order, ]
  gurvich_mat <- gurvich_mat[gene_order, ]
  
  #combine into one matrix
  combined_mat <- cbind(zhou_mat, gre_mat, gurvich_mat)
  
  #create annotation for dataset boundaries
  col_annotation <- data.frame(
    Dataset = c("Zhou", rep("GRE Lab", ncol(gre_mat)), rep("Gurvich", ncol(gurvich_mat))),
    row.names = colnames(combined_mat)
  )
  
  #colors for annotation
  ann_colors <- list(Dataset = c("Zhou" = "yellow", "GRE Lab" = "#377EB8", "Gurvich" = "#4DAF4A"))
  
  #calculate gaps (borders between datasets)
  gaps_col <- c(ncol(zhou_mat), ncol(zhou_mat) + ncol(gre_mat))
  
  cat(sprintf("\n%s: %d genes plotted\n", geneset_name, length(gene_order)))
  
  #adjust fontsize for ribosomal
  fontsize_row_val <- ifelse(length(gene_order) > 100, 5, 7)
  height_val <- ifelse(length(gene_order) > 100, 18, 12)
  
  #display inline
  pheatmap(combined_mat,
           color = colorRampPalette(c("blue", "white", "red"))(100),
           breaks = seq(color_limits[1], color_limits[2], length.out = 101),
           cluster_rows = FALSE,
           cluster_cols = FALSE,
           show_rownames = TRUE,
           annotation_col = col_annotation,
           annotation_colors = ann_colors,
           gaps_col = gaps_col,
           main = paste(geneset_name, "Gene Set: Three-Dataset Comparison"),
           fontsize_row = fontsize_row_val,
           cellwidth = 12,
           cellheight = 6)
  
  #save to file
  pheatmap(combined_mat,
           color = colorRampPalette(c("blue", "white", "red"))(100),
           breaks = seq(color_limits[1], color_limits[2], length.out = 101),
           cluster_rows = FALSE,
           cluster_cols = FALSE,
           show_rownames = TRUE,
           annotation_col = col_annotation,
           annotation_colors = ann_colors,
           gaps_col = gaps_col,
           main = paste(geneset_name, "Gene Set: Three-Dataset Comparison"),
           fontsize_row = fontsize_row_val,
           cellwidth = 12,
           cellheight = 6,
           filename = file.path(output_dir, "02_correlation", 
                               paste0("threeway_pheatmap_", geneset_name, ".pdf")),
           width = 14, height = height_val)
  
  return(combined_mat)
}

#create and display heatmaps
pho_result <- create_threeway_pheatmap(pho_genes, "PHO", color_limits = c(-2, 2))

## 
## PHO: 24 genes plotted

stress_result <- create_threeway_pheatmap(stress_genes, "Stress", color_limits = c(-2, 2))

## 
## Stress: 30 genes plotted

ribo_result <- create_threeway_pheatmap(ribosomal_genes, "Ribosomal", color_limits = c(-2, 2))

## 
## Ribosomal: 147 genes plotted

Three-Way Timecourse Heatmap Comparison

Methods

Heatmaps were generated for each gene set (PHO, Stress, Ribosomal) combining Zhou (1h), GRE Lab (0-8h timecourse), and Gurvich (1-8h timecourse) datasets. Only genes present in all three datasets were included. Genes were ordered by mean fold change across all datasets and timepoints. Heatmaps used a blue-white-red color scale (-2 to +2 log2FC) with dataset boundaries marked by vertical dividers and color-coded annotations (Zhou=yellow, GRE Lab=blue, Gurvich=green). No hierarchical clustering was applied to preserve temporal ordering.

Results

Ribosomal genes (n=~100-150) displayed ubiquitous repression across all three datasets throughout the timecourse, with consistent blue coloration indicating coordinated downregulation of the translation machinery during phosphate starvation. This signature was remarkably conserved across 14 years and two platforms.

PHO pathway genes showed strong induction (red) for canonical members (PHO5, PHO84, PHO89, PHO11, PHO12, VTC1-4) across all three datasets, validating the core starvation response. However, genes annotated as PHO-responsive in Gurvich (DUR3, ISC10) exhibited weak or absent induction signals in both GRE Lab and Zhou datasets, suggesting these may represent context-specific or secondary responses rather than core PHO regulon members. The concordance between GRE Lab and Zhou for these genes supports this interpretation.

Stress genes revealed dataset-specific intensity differences: GRE Lab exhibited stronger early stress signals (brighter red, 0.5-2h) compared to Gurvich and Zhou, which showed more muted initial responses. This likely reflects laboratory-specific variation in pre-starvation growth conditions or harvest protocols affecting baseline stress states. Despite these absolute magnitude differences, the temporal pattern of stress gene induction was preserved across datasets—genes that responded in one study responded in all three—demonstrating that directional stress signatures coalesce even when initial magnitudes vary. By later timepoints (>3h), stress response patterns converged across datasets.

Cross-dataset validation

The three-way comparison confirmed that core metabolic reprogramming signatures (PHO induction, ribosomal repression) are highly reproducible, while peripheral responses (marginal PHO genes, stress intensity) show greater laboratory-specific variability that diminishes over time as transcriptional programs mature.

Note:

The apparent extended variability observed in GRE Lab data relative to Gurvich is partially attributable to differences in temporal sampling resolution. GRE Lab captured early timepoints at finer granularity (0, 0.5, 1, 1.5, 2, 2.5h), resolving transient early responses that Gurvich’s coarser early sampling (1, 2h) did not detect. Despite this sampling difference, all three datasets converge to concordant response patterns by 2h post-starvation, confirming that core transcriptional programs are reproducible once they mature beyond the initial commitment phase.