Biological Context & Analysis Goals

Research Hypothesis

Inv4m constitutively upregulates DNA replication genes (PCNA2, MCM5, others) across tissues, accelerating cell cycle in the SAM, leading to:

Larger SAM size (9% increase at 3 weeks)
Taller plants
Earlier flowering

Specific Questions

Field-to-chamber replication: Do field adult leaf DEGs replicate in chamber V3 leaf tissues (V1_Leaf, V3_Leaf, V3_S2_Leaf_Tip)?
SAM hypothesis: Are field leaf DEGs upregulated in V3_Stem_SAM (SAM-enriched tissue)?
Temperature sensitivity: Do Inv4m effects increase in magnitude under cold (highland adaptation hypothesis)? Focus on SAM and replication genes.
Root tip validation: Are PCNA2/MCM5 upregulated in V1_Root (to justify EdU/confocal experiments)?

Statistical Approach

Cell means model with explicit contrasts:

No reference levels for biological factors
Direct interpretation: each coefficient = mean expression in that condition
Symmetric treatment of all tissues
Clear temperature dependence tests

Setup

library(edgeR)

## Loading required package: limma

library(limma)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyr)
library(stringr)
library(ggplot2)
library(ggbeeswarm)
library(rstatix)

## 
## Attaching package: 'rstatix'

## The following object is masked from 'package:stats':
## 
##     filter

library(ggpubr)
library(readr)
library(pheatmap)

options(stringsAsFactors = FALSE)

Data Loading

Expression Counts

counts <- read.csv("~/Desktop/crow2020_gene_xp.csv")
gene_ids <- counts[, 1]
counts_matrix <- as.matrix(counts[, -1])
rownames(counts_matrix) <- gene_ids

cat("Expression data loaded:\n")

## Expression data loaded:

cat("  Genes:", nrow(counts_matrix), "\n")

##   Genes: 39756

cat("  Samples:", ncol(counts_matrix), "\n")

##   Samples: 346

Sample Metadata

quant_dir <- "/Volumes/DOE_CAREER/inv4m/run/results_quantify/kallisto_quant"

# Define paths
sra_file <- "/Volumes/DOE_CAREER/inv4m/crow2020/SraRunTable.csv"
sample_submission_file <- "/Volumes/DOE_CAREER/inv4m/crow2020/gc7_sample_submission.csv"

# Load data
sra_meta <- read.csv(sra_file, stringsAsFactors = FALSE, check.names = FALSE)
submission_meta <- read.csv(sample_submission_file, stringsAsFactors = FALSE)

# Standardize SRA metadata
sra_meta <- sra_meta %>%
  dplyr::rename(
    samp = `Library Name`,
    genotype = genotype,
    temp = temp,
    tissue = tissue,
    Run = Run
  )

# Join core metadata and filter to NIL lines of interest
metadata <- submission_meta %>%
  filter(old_line %in% c("PT_NIL", "Mi21_NIL")) %>%
  dplyr::select(samp, old_line, ID, source, harvest_date, tube_num_full) %>%
  inner_join(sra_meta, by = "samp") %>%
  arrange(Run)

# Create analysis factors and sampling covariates
metadata <- metadata %>%
  mutate(
    # === Biological factors ===
    # Donor
    Donor = case_when(
      grepl("^PT_", old_line) ~ "PT",
      grepl("^Mi21_", old_line) ~ "Mi21",
      TRUE ~ NA_character_
    ),
    Donor = factor(Donor, levels = c("PT", "Mi21")),
    
    # Genotype
    Genotype = factor(genotype, levels = c("B73", "Inv4m")),
    Genotype = recode(Genotype, "B73" = "CTRL", "Inv4m" = "Inv4m"),
    
    # temperature
    temp = factor(temp, levels = c("warm", "cold")),
    temp = recode(temp, "warm" = "CTRL", "cold" = "cold"),
    
    # Tissue — with leaf subtypes
    Tissue = case_when(
      tissue == "V1_Leaf_tissue"    ~ "V1_Leaf",
      tissue == "V1_Primary_root"   ~ "V1_Root",
      tissue == "V3_Leaf_tissue"    ~ "V3_Leaf",
      tissue == "V3_leaf_base"      ~ "V3_Leaf_Base",
      tissue == "V3_leaf_sheath"    ~ "V3_Leaf_Sheath",
      tissue == "V3_S2_leaf_base"   ~ "V3_S2_Leaf_Base",
      tissue == "V3_S2_leaf_tip"    ~ "V3_S2_Leaf_Tip",
      tissue == "V3_Primary_root"   ~ "V3_Root",
      tissue == "V3_Stem_SAM"       ~ "V3_SAM",
      TRUE ~ NA_character_
    ),
    Tissue = factor(Tissue, levels = c(
      "V1_Root", "V3_Root", "V1_Leaf", "V3_Leaf", 
      "V3_Leaf_Base", "V3_Leaf_Sheath", 
      "V3_S2_Leaf_Base", "V3_S2_Leaf_Tip", "V3_SAM"
    )),
    
    # === Sampling & technical covariates ===
    # Unique plant identifier (from submission metadata)
    plant_id = ID,
    
    # Collection team (source column = technician/team ID)
    collector_team = as.factor(source),
    
    # Harvest date (biological sampling day)
    harvest_date = as.Date(harvest_date, format = "%m/%d/%Y"),
    
    # Plant serial number (numeric suffix of ID; proxy for planting/sampling order)
    plant_serial_number = as.numeric(stringr::str_sub(ID, 5)),
    
    # Within-plant tissue order (0 = leaf, 1 = root, 2 = base, etc.)
    tissue_order_index = as.numeric(
      stringr::str_extract(tube_num_full, "\\.[0-9]+")
    ) %>% 
      dplyr::coalesce(0) %>% 
      abs(),
    
    # SRA upload time (proxy for library prep + sequencing batch)
    sra_upload_time = as.POSIXct(create_date, format = "%Y-%m-%dT%H:%M:%OS", tz = "UTC"),
    
    # Combined harvest batch (date + team)
    harvest_batch = interaction(harvest_date, collector_team, drop = TRUE)
  ) %>%
  # Remove samples with unmapped tissue
  filter(!is.na(Tissue))

metadata$harvest_date_f <- factor(metadata$harvest_date)

{
  cat("\nSample metadata:\n")
  cat("  Total samples:", nrow(metadata), "\n")
  cat("  Genotypes:", paste(unique(metadata$Genotype), collapse = ", "), "\n")
  cat("  temperatures:", paste(unique(metadata$temp), collapse = ", "), "\n")
  cat("  Donors:", paste(unique(metadata$Donor), collapse = ", "), "\n")
  cat("  Tissues:", paste(levels(metadata$Tissue), collapse = ", "), "\n")
  cat("  Collection teams:", paste(levels(metadata$collector_team), collapse = ", "), "\n")
  cat("  Harvest dates:", paste(sort(unique(metadata$harvest_date)), collapse = ", "), "\n")
  cat("\nTissue counts:\n")
  print(sort(table(metadata$Tissue), decreasing = TRUE))
}

## 
## Sample metadata:
##   Total samples: 346 
##   Genotypes: CTRL, Inv4m 
##   temperatures: CTRL, cold 
##   Donors: PT, Mi21 
##   Tissues: V1_Root, V3_Root, V1_Leaf, V3_Leaf, V3_Leaf_Base, V3_Leaf_Sheath, V3_S2_Leaf_Base, V3_S2_Leaf_Tip, V3_SAM 
##   Collection teams: 1, 2 
##   Harvest dates: 2017-03-07, 2017-03-08, 2017-03-09, 2017-03-13, 2017-03-14, 2017-03-15, 2017-03-16, 2017-03-17, 2017-03-24, 2017-03-28, 2017-04-27, 2017-04-28, 2017-05-02, 2017-05-03, 2017-05-04, 2017-05-12, 2017-05-15, 2017-05-16 
## 
## Tissue counts:
## 
##         V3_Root         V1_Root          V3_SAM         V1_Leaf V3_S2_Leaf_Base 
##              42              41              40              39              38 
##         V3_Leaf  V3_Leaf_Sheath    V3_Leaf_Base  V3_S2_Leaf_Tip 
##              37              37              36              36

Verify Design Balance

cat("\nGenotype × Temperature × Tissue:\n")

## 
## Genotype × Temperature × Tissue:

print(with(metadata, table(Genotype, temp, Tissue)))

## , , Tissue = V1_Root
## 
##         temp
## Genotype CTRL cold
##    CTRL    11   11
##    Inv4m   10    9
## 
## , , Tissue = V3_Root
## 
##         temp
## Genotype CTRL cold
##    CTRL    13   11
##    Inv4m   11    7
## 
## , , Tissue = V1_Leaf
## 
##         temp
## Genotype CTRL cold
##    CTRL    10   12
##    Inv4m    9    8
## 
## , , Tissue = V3_Leaf
## 
##         temp
## Genotype CTRL cold
##    CTRL    12   10
##    Inv4m    9    6
## 
## , , Tissue = V3_Leaf_Base
## 
##         temp
## Genotype CTRL cold
##    CTRL    12    6
##    Inv4m   11    7
## 
## , , Tissue = V3_Leaf_Sheath
## 
##         temp
## Genotype CTRL cold
##    CTRL    13    7
##    Inv4m   11    6
## 
## , , Tissue = V3_S2_Leaf_Base
## 
##         temp
## Genotype CTRL cold
##    CTRL    13    8
##    Inv4m    9    8
## 
## , , Tissue = V3_S2_Leaf_Tip
## 
##         temp
## Genotype CTRL cold
##    CTRL    12    8
##    Inv4m    9    7
## 
## , , Tissue = V3_SAM
## 
##         temp
## Genotype CTRL cold
##    CTRL    13    8
##    Inv4m   11    8

cat("\nDonor × Genotype (overall):\n")

## 
## Donor × Genotype (overall):

print(with(metadata, table(Donor, Genotype)))

##       Genotype
## Donor  CTRL Inv4m
##   PT    100    71
##   Mi21   90    85

Prepare Count Matrix

sample_order <- match(colnames(counts_matrix), metadata$Run)
metadata <- metadata[sample_order, ]

stopifnot(all(colnames(counts_matrix) == metadata$Run))

cat("\nSample order verified\n")

## 
## Sample order verified

DGEList and Filtering

Create DGEList with Cell Means Condition Factor

metadata$condition <- interaction(
  metadata$Tissue,
  metadata$Genotype,
  metadata$temp,
  sep = "_",
  drop = TRUE
)

y <- DGEList(counts = counts_matrix, samples = metadata)

cat("\nDGEList created\n")

## 
## DGEList created

cat("Unique conditions:", nlevels(y$samples$condition), "\n")

## Unique conditions: 36

cat("\nCondition factor levels (first 10):\n")

## 
## Condition factor levels (first 10):

print(head(levels(y$samples$condition), 10))

##  [1] "V1_Root_CTRL_CTRL"         "V3_Root_CTRL_CTRL"        
##  [3] "V1_Leaf_CTRL_CTRL"         "V3_Leaf_CTRL_CTRL"        
##  [5] "V3_Leaf_Base_CTRL_CTRL"    "V3_Leaf_Sheath_CTRL_CTRL" 
##  [7] "V3_S2_Leaf_Base_CTRL_CTRL" "V3_S2_Leaf_Tip_CTRL_CTRL" 
##  [9] "V3_SAM_CTRL_CTRL"          "V1_Root_Inv4m_CTRL"

Filter Low-Expression Genes

keep <- filterByExpr(y, group = y$samples$condition)
y <- y[keep, , keep.lib.sizes = FALSE]

cat("\nGenes retained after filtering:", nrow(y), "\n")

## 
## Genes retained after filtering: 26490

TMM Normalization

y <- calcNormFactors(y, method = "TMM")

cat("Normalization factors range:", 
    round(range(y$samples$norm.factors), 3), "\n")

## Normalization factors range: 0.512 1.608

Design Matrix: Cell Means Model

Build Design (No Reference Levels for Biological Factors)

design <- model.matrix(~ 0 + Donor + condition , data = y$samples)

colnames(design) <- gsub("^condition", "", colnames(design))
colnames(design) <- gsub("^Donor", "Donor_", colnames(design))

cat("\nDesign matrix:\n")

## 
## Design matrix:

cat("  Dimensions:", nrow(design), "samples ×", ncol(design), "coefficients\n")

##   Dimensions: 346 samples × 37 coefficients

cat("\nFirst 15 coefficient names:\n")

## 
## First 15 coefficient names:

print(head(colnames(design), 15))

##  [1] "Donor_PT"                  "Donor_Mi21"               
##  [3] "V3_Root_CTRL_CTRL"         "V1_Leaf_CTRL_CTRL"        
##  [5] "V3_Leaf_CTRL_CTRL"         "V3_Leaf_Base_CTRL_CTRL"   
##  [7] "V3_Leaf_Sheath_CTRL_CTRL"  "V3_S2_Leaf_Base_CTRL_CTRL"
##  [9] "V3_S2_Leaf_Tip_CTRL_CTRL"  "V3_SAM_CTRL_CTRL"         
## [11] "V1_Root_Inv4m_CTRL"        "V3_Root_Inv4m_CTRL"       
## [13] "V1_Leaf_Inv4m_CTRL"        "V3_Leaf_Inv4m_CTRL"       
## [15] "V3_Leaf_Base_Inv4m_CTRL"

Voom Transformation

voom_obj <- voom(y, design = design, plot = TRUE)

Contrasts: Explicit Biological Comparisons

Define Contrasts

# Verify available conditions
cat("Available condition coefficients in design matrix:\n")

## Available condition coefficients in design matrix:

cond_names <- grep("^V[0-9]", colnames(design), value = TRUE)
print(cond_names)

##  [1] "V3_Root_CTRL_CTRL"          "V1_Leaf_CTRL_CTRL"         
##  [3] "V3_Leaf_CTRL_CTRL"          "V3_Leaf_Base_CTRL_CTRL"    
##  [5] "V3_Leaf_Sheath_CTRL_CTRL"   "V3_S2_Leaf_Base_CTRL_CTRL" 
##  [7] "V3_S2_Leaf_Tip_CTRL_CTRL"   "V3_SAM_CTRL_CTRL"          
##  [9] "V1_Root_Inv4m_CTRL"         "V3_Root_Inv4m_CTRL"        
## [11] "V1_Leaf_Inv4m_CTRL"         "V3_Leaf_Inv4m_CTRL"        
## [13] "V3_Leaf_Base_Inv4m_CTRL"    "V3_Leaf_Sheath_Inv4m_CTRL" 
## [15] "V3_S2_Leaf_Base_Inv4m_CTRL" "V3_S2_Leaf_Tip_Inv4m_CTRL" 
## [17] "V3_SAM_Inv4m_CTRL"          "V1_Root_CTRL_cold"         
## [19] "V3_Root_CTRL_cold"          "V1_Leaf_CTRL_cold"         
## [21] "V3_Leaf_CTRL_cold"          "V3_Leaf_Base_CTRL_cold"    
## [23] "V3_Leaf_Sheath_CTRL_cold"   "V3_S2_Leaf_Base_CTRL_cold" 
## [25] "V3_S2_Leaf_Tip_CTRL_cold"   "V3_SAM_CTRL_cold"          
## [27] "V1_Root_Inv4m_cold"         "V3_Root_Inv4m_cold"        
## [29] "V1_Leaf_Inv4m_cold"         "V3_Leaf_Inv4m_cold"        
## [31] "V3_Leaf_Base_Inv4m_cold"    "V3_Leaf_Sheath_Inv4m_cold" 
## [33] "V3_S2_Leaf_Base_Inv4m_cold" "V3_S2_Leaf_Tip_Inv4m_cold" 
## [35] "V3_SAM_Inv4m_cold"

# Build contrasts using actual condition names
# CTRL = control temperature, cold = cold temperature
# Note: V1_Root_CTRL_CTRL was dropped as reference level

my_contrasts <- makeContrasts(
  
  # Q1 & Q2: Inv4m effect in leaf tissues at CTRL temperature
  Inv4m_V1Leaf_CTRL = V1_Leaf_Inv4m_CTRL - V1_Leaf_CTRL_CTRL,
  Inv4m_V3Leaf_CTRL = V3_Leaf_Inv4m_CTRL - V3_Leaf_CTRL_CTRL,
  Inv4m_V3Tip_CTRL = V3_S2_Leaf_Tip_Inv4m_CTRL - V3_S2_Leaf_Tip_CTRL_CTRL,
  
  # Q2: SAM hypothesis
  Inv4m_SAM_CTRL = V3_SAM_Inv4m_CTRL - V3_SAM_CTRL_CTRL,
  
  # Q3: Inv4m effect under cold temperature
  Inv4m_V1Leaf_cold = V1_Leaf_Inv4m_cold - V1_Leaf_CTRL_cold,
  Inv4m_V3Leaf_cold = V3_Leaf_Inv4m_cold - V3_Leaf_CTRL_cold,
  Inv4m_V3Tip_cold = V3_S2_Leaf_Tip_Inv4m_cold - V3_S2_Leaf_Tip_CTRL_cold,
  Inv4m_SAM_cold = V3_SAM_Inv4m_cold - V3_SAM_CTRL_cold,
  
  # Q3: Temperature dependence (cold effect - CTRL effect)
  Inv4m_V1Leaf_tempDep = 
    (V1_Leaf_Inv4m_cold - V1_Leaf_CTRL_cold) - 
    (V1_Leaf_Inv4m_CTRL - V1_Leaf_CTRL_CTRL),
  
  Inv4m_V3Leaf_tempDep = 
    (V3_Leaf_Inv4m_cold - V3_Leaf_CTRL_cold) - 
    (V3_Leaf_Inv4m_CTRL - V3_Leaf_CTRL_CTRL),
  
  Inv4m_V3Tip_tempDep = 
    (V3_S2_Leaf_Tip_Inv4m_cold - V3_S2_Leaf_Tip_CTRL_cold) - 
    (V3_S2_Leaf_Tip_Inv4m_CTRL - V3_S2_Leaf_Tip_CTRL_CTRL),
  
  Inv4m_SAM_tempDep = 
    (V3_SAM_Inv4m_cold - V3_SAM_CTRL_cold) - 
    (V3_SAM_Inv4m_CTRL - V3_SAM_CTRL_CTRL),
  
  # Q4: Root validation
  # V1_Root_CTRL_CTRL is reference level (= 0 in design)
  Inv4m_V1Root_CTRL = V1_Root_Inv4m_CTRL,
  
  Inv4m_V1Root_cold = V1_Root_Inv4m_cold - V1_Root_CTRL_cold,
  
  Inv4m_V1Root_tempDep = 
    (V1_Root_Inv4m_cold - V1_Root_CTRL_cold) - 
    V1_Root_Inv4m_CTRL,
  
  levels = design
)

cat("\nContrasts successfully defined:\n")

## 
## Contrasts successfully defined:

print(colnames(my_contrasts))

##  [1] "Inv4m_V1Leaf_CTRL"    "Inv4m_V3Leaf_CTRL"    "Inv4m_V3Tip_CTRL"    
##  [4] "Inv4m_SAM_CTRL"       "Inv4m_V1Leaf_cold"    "Inv4m_V3Leaf_cold"   
##  [7] "Inv4m_V3Tip_cold"     "Inv4m_SAM_cold"       "Inv4m_V1Leaf_tempDep"
## [10] "Inv4m_V3Leaf_tempDep" "Inv4m_V3Tip_tempDep"  "Inv4m_SAM_tempDep"   
## [13] "Inv4m_V1Root_CTRL"    "Inv4m_V1Root_cold"    "Inv4m_V1Root_tempDep"

Model Fitting

Fit Linear Models

fit <- lmFit(voom_obj, design)
fit2 <- contrasts.fit(fit, my_contrasts)
fit2 <- eBayes(fit2)

cat("\nModel fitted with", ncol(my_contrasts), "contrasts\n")

## 
## Model fitted with 15 contrasts

cat("Residual df per gene (summary):\n")

## Residual df per gene (summary):

print(summary(fit2$df.residual))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     309     309     309     309     309     309

Results Extraction

Extract TopTables for All Contrasts

extract_results <- function(fit_obj, coef_name) {
  topTable(fit_obj, coef = coef_name, number = Inf, adjust.method = "BH") %>%
    mutate(contrast = coef_name) %>%
    tibble::rownames_to_column("gene_id")
}

all_contrasts <- colnames(my_contrasts)
results_list <- lapply(all_contrasts, function(x) extract_results(fit2, x))
names(results_list) <- all_contrasts

results_combined <- bind_rows(results_list)

cat("\nResults extracted for", length(all_contrasts), "contrasts\n")

## 
## Results extracted for 15 contrasts

cat("Total rows in combined results:", nrow(results_combined), "\n")

## Total rows in combined results: 397350

Summary: Significant DEGs per Contrast

sig_counts <- results_combined %>%
  filter(adj.P.Val < 0.05) %>%
  group_by(contrast) %>%
  summarise(
    n_sig = n(),
    n_up = sum(logFC > 0),
    n_down = sum(logFC < 0)
  ) %>%
  arrange(desc(n_sig))

cat("\nSignificant DEGs (adj.P.Val < 0.05) per contrast:\n")

## 
## Significant DEGs (adj.P.Val < 0.05) per contrast:

print(sig_counts)

## # A tibble: 12 × 4
##    contrast             n_sig  n_up n_down
##    <chr>                <int> <int>  <int>
##  1 Inv4m_V1Root_cold      128    43     85
##  2 Inv4m_SAM_CTRL         100    30     70
##  3 Inv4m_V1Root_CTRL       88    32     56
##  4 Inv4m_V3Tip_CTRL        73    26     47
##  5 Inv4m_V1Leaf_CTRL       67    20     47
##  6 Inv4m_V1Leaf_cold       66    25     41
##  7 Inv4m_V3Tip_cold        64    19     45
##  8 Inv4m_SAM_cold          53    12     41
##  9 Inv4m_V3Leaf_CTRL       51    12     39
## 10 Inv4m_V3Leaf_cold       51    14     37
## 11 Inv4m_V1Leaf_tempDep     2     2      0
## 12 Inv4m_V1Root_tempDep     2     0      2

Question 1: Field-to-Chamber Replication

Load Field DEG List

# Load differential expression results
effects_df <- read.csv("~/Desktop/predictor_effects_leaf_interaction_model.csv")

# Load curated labels
curated <- read.csv("~/Desktop/selected_DEGs_curated_locus_label_2.csv") %>%
  dplyr::select(gene, locus_label)

effects_df$locus_label[effects_df$locus_label == "A0A1D6NG77"] <- NA


# Merge curated labels
effects_df <- effects_df %>%
  left_join(curated, by = "gene", suffix = c("", "_curated")) %>%
  mutate(locus_label = coalesce(locus_label_curated, locus_label)) %>%
  dplyr::select(-locus_label_curated)

## Warning in left_join(., curated, by = "gene", suffix = c("", "_curated")): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 4420 of `x` matches multiple rows in `y`.
## ℹ Row 61 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.

# Curated updates: only gene and locus_label
curated <- tibble::tribble(
  ~gene,              ~locus_label,
  # custom labels for trans network
  "Zm00001eb384770",  "zcn26",
  "Zm00001eb192330",  "roc4-1",
  "Zm00001eb193260",  "trxl2",
  "Zm00001eb194300",  "paspa",
  "Zm00001eb213980",  "dywl",
  "Zm00001eb194080",  "ugt85a2",
  "Zm00001eb112470",  "engd2",
  "Zm00001eb332450",  "roc4-2",
  "Zm00001eb193270",  "trxl5",
  "Zm00001eb066620",  "tipin",
  "Zm00001eb241410",  "etfa",
  "Zm00001eb013800",  "mrlk",
  "Zm00001eb249540",  "prpl29",
  "Zm00001eb000540",  "prd3",
  "Zm00001eb060540",  "actin-1",
  "Zm00001eb192580",  "map4k8",
  "Zm00001eb193120",  "o3l1",
  "Zm00001eb043750",  "cdrk9",
  # From original curated list
  "Zm00001eb034810",  "mgd2",
  "Zm00001eb241920",  "gpx1",
  "Zm00001eb162710",  "spx4",
  "Zm00001eb038730",  "pht7",
  "Zm00001eb370610",  "rfk1",
  "Zm00001eb386270",  "spx6",
  "Zm00001eb335670",  "sqd3",
  "Zm00001eb151650",  "pap1",
  "Zm00001eb154650",  "ppa1",
  "Zm00001eb280120",  "pfk1",
  "Zm00001eb036910",  "gpx3",
  "Zm00001eb406610",  "glk4",
  "Zm00001eb048730",  "spx2",
  "Zm00001eb361620",  "ppa2",
  "Zm00001eb297970",  "sqd2",
  "Zm00001eb347070",  "sqd1",
  "Zm00001eb144680",  "rns1",
  "Zm00001eb430590",  "nrx3",
  "Zm00001eb247580",  "ppck3",
  "Zm00001eb351780",  "ugp3",
  "Zm00001eb222510",  "pht1",
  "Zm00001eb126380",  "phos1",
  "Zm00001eb047070",  "pht2",
  "Zm00001eb041390",  "rns3",
  "Zm00001eb116580",  "spd1",
  "Zm00001eb069630",  "oct4",
  "Zm00001eb083520",  "dgd1",
  "Zm00001eb130570",  "sag21",
  "Zm00001eb239700",  "ppa3",  
  "Zm00001eb258130",  "mgd3",
  "Zm00001eb202100",  "pap14",
  "Zm00001eb144670",  "rns2",
  "Zm00001eb087720",  "pht13",
  "Zm00001eb277280",  "gst19",
  "Zm00001eb339870",  "pld16",
  "Zm00001eb289800",  "pah1",
  "Zm00001eb369590",  "nrx1",
  "Zm00001eb238670",  "pep2"
)

# Update only the locus_label using curated values
effects_df <- effects_df %>%
  left_join(curated, by = "gene", suffix = c("", "_new")) %>%
  mutate(locus_label = coalesce(locus_label_new, locus_label)) %>%
  dplyr::select(-locus_label_new)

# Replace with your actual field DEG gene IDs
field_degs <- effects_df %>% filter(is_hiconf_DEG, predictor=="Inv4m") %>% pull(gene) %>% unique()

cat("\nField DEG list loaded:\n")

## 
## Field DEG list loaded:

cat("  Total field DEGs:", length(field_degs), "\n")

##   Total field DEGs: 165

Test Replication in Chamber Leaf Tissues

leaf_contrasts_CTRL <- c("Inv4m_V1Leaf_CTRL", "Inv4m_V3Leaf_CTRL", "Inv4m_V3Tip_CTRL")

replication_results <- results_combined %>%
  filter(contrast %in% leaf_contrasts_CTRL, gene_id %in% field_degs) %>%
  select(gene_id, contrast, logFC, P.Value, adj.P.Val) %>%
  arrange(gene_id, contrast)

cat("\nField DEGs in chamber leaf tissues (warm):\n")

## 
## Field DEGs in chamber leaf tissues (warm):

# print(replication_results)

rep_summary <- replication_results %>%
  group_by(contrast) %>%
  summarise(
    n_tested = n(),
    n_sig = sum(adj.P.Val < 0.05),
    n_same_direction = sum(logFC > 0 & adj.P.Val < 0.05),
    prop_replicated = n_sig / n_tested
  )

cat("\nReplication summary:\n")

## 
## Replication summary:

print(rep_summary)

## # A tibble: 3 × 5
##   contrast          n_tested n_sig n_same_direction prop_replicated
##   <chr>                <int> <int>            <int>           <dbl>
## 1 Inv4m_V1Leaf_CTRL      144    31                8           0.215
## 2 Inv4m_V3Leaf_CTRL      144    28                8           0.194
## 3 Inv4m_V3Tip_CTRL       144    35               12           0.243

Enrichment Test: Field DEGs in Chamber Leaves

v3tip_CTRL <- results_list[["Inv4m_V3Tip_CTRL"]]
v3tip_sig <- v3tip_CTRL$gene_id[v3tip_CTRL$adj.P.Val < 0.05]

common_genes <- intersect(field_degs, rownames(voom_obj))
A <- length(intersect(field_degs, v3tip_sig))
B <- length(setdiff(field_degs, v3tip_sig))
C <- length(setdiff(v3tip_sig, field_degs))
D <- nrow(voom_obj) - (A + B + C)

fisher_mat <- matrix(c(A, B, C, D), nrow = 2,
                     dimnames = list(
                       c("In_V3Tip_sig", "Not_in_V3Tip_sig"),
                       c("In_field_list", "Not_in_field_list")
                     ))

fisher_test <- fisher.test(fisher_mat)

cat("\nFisher's exact test: Field DEGs vs V3_Tip chamber DEGs\n")

## 
## Fisher's exact test: Field DEGs vs V3_Tip chamber DEGs

print(fisher_mat)

##                  In_field_list Not_in_field_list
## In_V3Tip_sig                35                38
## Not_in_V3Tip_sig           130             26287

cat("\nOdds ratio:", round(fisher_test$estimate, 2), "\n")

## 
## Odds ratio: 186.31

cat("P-value:", signif(fisher_test$p.value, 3), "\n")

## P-value: 9.2e-59

Question 2: SAM Hypothesis

Test Field DEGs in SAM

sam_CTRL <- results_list[["Inv4m_SAM_CTRL"]]

sam_field_overlap <- sam_CTRL %>%
  filter(gene_id %in% field_degs) %>%
  select(gene_id, logFC, AveExpr, P.Value, adj.P.Val) %>%
  arrange(adj.P.Val)

cat("\nField DEGs in V3_SAM (warm):\n")

## 
## Field DEGs in V3_SAM (warm):

print(sam_field_overlap %>% pull(gene_id))

##   [1] "Zm00001eb190240" "Zm00001eb192880" "Zm00001eb193870" "Zm00001eb188220"
##   [5] "Zm00001eb189350" "Zm00001eb196780" "Zm00001eb191790" "Zm00001eb191820"
##   [9] "Zm00001eb196760" "Zm00001eb189580" "Zm00001eb194380" "Zm00001eb188750"
##  [13] "Zm00001eb192170" "Zm00001eb193300" "Zm00001eb193620" "Zm00001eb192160"
##  [17] "Zm00001eb189200" "Zm00001eb213070" "Zm00001eb187890" "Zm00001eb182850"
##  [21] "Zm00001eb188070" "Zm00001eb194700" "Zm00001eb189170" "Zm00001eb192840"
##  [25] "Zm00001eb126050" "Zm00001eb188570" "Zm00001eb190350" "Zm00001eb112470"
##  [29] "Zm00001eb195790" "Zm00001eb264460" "Zm00001eb190090" "Zm00001eb196830"
##  [33] "Zm00001eb196560" "Zm00001eb183030" "Zm00001eb188990" "Zm00001eb241410"
##  [37] "Zm00001eb190670" "Zm00001eb189710" "Zm00001eb194680" "Zm00001eb186670"
##  [41] "Zm00001eb194110" "Zm00001eb194570" "Zm00001eb192180" "Zm00001eb195180"
##  [45] "Zm00001eb041890" "Zm00001eb193130" "Zm00001eb196100" "Zm00001eb188240"
##  [49] "Zm00001eb129670" "Zm00001eb194180" "Zm00001eb007720" "Zm00001eb192750"
##  [53] "Zm00001eb192580" "Zm00001eb276970" "Zm00001eb196630" "Zm00001eb191990"
##  [57] "Zm00001eb193260" "Zm00001eb060540" "Zm00001eb188610" "Zm00001eb195370"
##  [61] "Zm00001eb186860" "Zm00001eb194020" "Zm00001eb249540" "Zm00001eb193120"
##  [65] "Zm00001eb070810" "Zm00001eb187460" "Zm00001eb195620" "Zm00001eb192630"
##  [69] "Zm00001eb196920" "Zm00001eb428740" "Zm00001eb166340" "Zm00001eb186400"
##  [73] "Zm00001eb191780" "Zm00001eb334460" "Zm00001eb157030" "Zm00001eb196200"
##  [77] "Zm00001eb191000" "Zm00001eb194300" "Zm00001eb196180" "Zm00001eb194610"
##  [81] "Zm00001eb193660" "Zm00001eb188460" "Zm00001eb371380" "Zm00001eb195800"
##  [85] "Zm00001eb248260" "Zm00001eb189920" "Zm00001eb190750" "Zm00001eb187450"
##  [89] "Zm00001eb194190" "Zm00001eb187400" "Zm00001eb013800" "Zm00001eb068500"
##  [93] "Zm00001eb195970" "Zm00001eb192330" "Zm00001eb193230" "Zm00001eb196170"
##  [97] "Zm00001eb213980" "Zm00001eb043280" "Zm00001eb186610" "Zm00001eb228100"
## [101] "Zm00001eb187030" "Zm00001eb140380" "Zm00001eb341430" "Zm00001eb192470"
## [105] "Zm00001eb187280" "Zm00001eb332450" "Zm00001eb192970" "Zm00001eb186790"
## [109] "Zm00001eb187290" "Zm00001eb122060" "Zm00001eb170590" "Zm00001eb250800"
## [113] "Zm00001eb170670" "Zm00001eb096620" "Zm00001eb192910" "Zm00001eb339770"
## [117] "Zm00001eb190450" "Zm00001eb123230" "Zm00001eb024970" "Zm00001eb187850"
## [121] "Zm00001eb189900" "Zm00001eb195600" "Zm00001eb148670" "Zm00001eb043750"
## [125] "Zm00001eb192120" "Zm00001eb196240" "Zm00001eb197370" "Zm00001eb186970"
## [129] "Zm00001eb187440" "Zm00001eb000540" "Zm00001eb195080" "Zm00001eb192350"
## [133] "Zm00001eb196600" "Zm00001eb189910" "Zm00001eb192850" "Zm00001eb192050"
## [137] "Zm00001eb251960" "Zm00001eb194270" "Zm00001eb197300" "Zm00001eb257900"
## [141] "Zm00001eb195030" "Zm00001eb188030" "Zm00001eb187430" "Zm00001eb188550"

cat("\nSignificant in SAM:", sum(sam_field_overlap$adj.P.Val < 0.05), "/", 
    nrow(sam_field_overlap), "\n")

## 
## Significant in SAM: 33 / 144

Enrichment: Field DEGs in SAM

sam_sig <- sam_CTRL$gene_id[sam_CTRL$adj.P.Val < 0.05]

A_sam <- length(intersect(field_degs, sam_sig))
B_sam <- length(setdiff(field_degs, sam_sig))
C_sam <- length(setdiff(sam_sig, field_degs))
D_sam <- nrow(voom_obj) - (A_sam + B_sam + C_sam)

fisher_mat_sam <- matrix(c(A_sam, B_sam, C_sam, D_sam), nrow = 2,
                         dimnames = list(
                           c("In_SAM_sig", "Not_in_SAM_sig"),
                           c("In_field_list", "Not_in_field_list")
                         ))

fisher_test_sam <- fisher.test(fisher_mat_sam)

cat("\nFisher's exact test: Field DEGs vs SAM DEGs\n")

## 
## Fisher's exact test: Field DEGs vs SAM DEGs

print(fisher_mat_sam)

##                In_field_list Not_in_field_list
## In_SAM_sig                33                67
## Not_in_SAM_sig           132             26258

cat("Odds ratio:", round(fisher_test_sam$estimate, 2), "\n")

## Odds ratio: 97.75

cat("P-value:", signif(fisher_test_sam$p.value, 3), "\n")

## P-value: 1.15e-48

Question 3: Temperature Sensitivity

PCNA2 & MCM5 Temperature Dependence in SAM

fork_genes <- c(pcna2="Zm00001eb192050", mcm5 ="Zm00001eb257900")

sam_contrasts <- c("Inv4m_SAM_CTRL", "Inv4m_SAM_cold", "Inv4m_SAM_tempDep")

fork_sam_temp <- results_combined %>%
  filter(contrast %in% sam_contrasts, gene_id %in% fork_genes) %>%
  select(gene_id, contrast, logFC, P.Value, adj.P.Val) %>%
  arrange(gene_id, match(contrast, sam_contrasts))

cat("\nReplication fork genes in SAM across temperatures:\n")

## 
## Replication fork genes in SAM across temperatures:

print(fork_sam_temp)

##           gene_id          contrast     logFC   P.Value adj.P.Val
## 1 Zm00001eb192050    Inv4m_SAM_CTRL 0.2189205 0.4905177 0.9999732
## 2 Zm00001eb192050    Inv4m_SAM_cold 0.4544017 0.2310523 0.9999940
## 3 Zm00001eb192050 Inv4m_SAM_tempDep 0.2354812 0.6337746 0.9998833
## 4 Zm00001eb257900    Inv4m_SAM_CTRL 0.1954433 0.5389770 0.9999732
## 5 Zm00001eb257900    Inv4m_SAM_cold 0.4416383 0.2484475 0.9999940
## 6 Zm00001eb257900 Inv4m_SAM_tempDep 0.2461950 0.6204894 0.9998833

Global Temperature Dependence Pattern

temp_dep_contrasts <- grep("tempDep$", all_contrasts, value = TRUE)

temp_dep_field <- results_combined %>%
  filter(contrast %in% temp_dep_contrasts, gene_id %in% field_degs) %>%
  select(gene_id, contrast, logFC, adj.P.Val) %>%
  arrange(contrast,adj.P.Val)

cat("\nField DEGs: Temperature dependence (cold - warm effect):\n")

## 
## Field DEGs: Temperature dependence (cold - warm effect):

#print(temp_dep_field)

temp_increase <- temp_dep_field %>%
  filter(abs(logFC) > 0.5) %>%
  group_by(contrast) %>%
  summarise(n_increased = n())

cat("\nField DEGs with increased Inv4m effect in cold:\n")

## 
## Field DEGs with increased Inv4m effect in cold:

# print(temp_increase)

Question 4: Root Validation

PCNA2 & MCM5 in V1_Root

root_contrasts <- c("Inv4m_V1Root_CTRL", "Inv4m_V1Root_cold")

fork_root <- results_combined %>%
  filter(contrast %in% root_contrasts, gene_id %in% fork_genes) %>%
  select(gene_id, contrast, logFC, P.Value, adj.P.Val) %>%
  arrange(gene_id, contrast)

cat("\nReplication fork genes in V1_Root:\n")

## 
## Replication fork genes in V1_Root:

print(fork_root)

##           gene_id          contrast       logFC    P.Value adj.P.Val
## 1 Zm00001eb192050 Inv4m_V1Root_CTRL  0.17289012 0.61592234 0.9999951
## 2 Zm00001eb192050 Inv4m_V1Root_cold -0.09012957 0.79958941 0.9627373
## 3 Zm00001eb257900 Inv4m_V1Root_CTRL -0.07495143 0.82730361 0.9999951
## 4 Zm00001eb257900 Inv4m_V1Root_cold -0.66625463 0.05818316 0.6368331

cat("\nConclusion for EdU experiment:\n")

## 
## Conclusion for EdU experiment:

if (any(fork_root$adj.P.Val < 0.05 & fork_root$logFC > 0)) {
  cat("YES - PCNA2/MCM5 show significant upregulation in roots.\n")
  cat("Proceed with EdU/confocal validation in root tips.\n")
} else {
  cat("NO - No significant upregulation detected in V1_Root.\n")
  cat("Consider alternative tissues or verify gene IDs.\n")
}

## NO - No significant upregulation detected in V1_Root.
## Consider alternative tissues or verify gene IDs.

Visualization

Volcano Plots: Key Contrasts

plot_volcano <- function(res_df, title, fdr_cut = 0.05, lfc_cut = 0.5) {
  res_df %>%
    mutate(
      sig = case_when(
        adj.P.Val < fdr_cut & abs(logFC) > lfc_cut ~ "Sig",
        TRUE ~ "NS"
      )
    ) %>%
    ggplot(aes(x = logFC, y = -log10(adj.P.Val), color = sig)) +
    geom_point(alpha = 0.6, size = 1.5) +
    scale_color_manual(values = c("Sig" = "red", "NS" = "gray60")) +
    geom_vline(xintercept = c(-lfc_cut, lfc_cut), 
               linetype = "dashed", alpha = 0.5) +
    geom_hline(yintercept = -log10(0.05), 
               linetype = "dashed", alpha = 0.5) +
    labs(title = title, x = "log2 Fold Change", y = "-log10(P-value)") +
    theme_minimal() +
    theme(legend.position = "bottom")
}

key_contrasts <- c("Inv4m_V3Tip_CTRL", "Inv4m_SAM_CTRL", 
                   "Inv4m_V1Root_CTRL", "Inv4m_SAM_tempDep")

for (ctr in key_contrasts) {
  p <- plot_volcano(results_list[[ctr]], title = ctr)
  print(p)
}

Heatmap: Field DEGs Across Conditions

common_field <- intersect(field_degs, rownames(voom_obj$E))

if (length(common_field) > 0) {
  expr_matrix <- voom_obj$E[common_field, , drop = FALSE]
  
  anno_df <- y$samples %>%
    select(Tissue, Genotype, temp) %>%
    as.data.frame()
  rownames(anno_df) <- colnames(expr_matrix)
  
  pheatmap(
    expr_matrix,
    annotation_col = anno_df,
    scale = "row",
    clustering_distance_rows = "correlation",
    clustering_distance_cols = "correlation",
    show_colnames = FALSE,
    main = "Field DEG Expression Across Crow2020 Conditions"
  )
}

Export Results

Save TopTables

output_dir <- "crow2020_reanalysis_results"
dir.create(output_dir, showWarnings = FALSE, recursive = TRUE)

for (ctr in names(results_list)) {
  write_csv(
    results_list[[ctr]],
    file.path(output_dir, paste0(ctr, "_results.csv"))
  )
}

write_csv(
  results_combined,
  file.path(output_dir, "all_contrasts_combined.csv")
)

cat("\nResults exported to:", output_dir, "\n")

## 
## Results exported to: crow2020_reanalysis_results

Save Key Gene Summary Table

key_contrasts_summary <- c(
  "Inv4m_V3Tip_CTRL", "Inv4m_V3Leaf_CTRL", "Inv4m_V1Leaf_CTRL",
  "Inv4m_SAM_CTRL", "Inv4m_SAM_cold", "Inv4m_SAM_tempDep",
  "Inv4m_V1Root_CTRL", "Inv4m_V1Root_cold"
)

key_genes_table <- results_combined %>%
  filter(contrast %in% key_contrasts_summary, gene_id %in% fork_genes) %>%
  select(gene_id, contrast, logFC, AveExpr, P.Value, adj.P.Val) %>%
  arrange(gene_id, contrast)

write_csv(
  key_genes_table,
  file.path(output_dir, "PCNA2_MCM5_summary.csv")
)

cat("\nKey genes summary saved\n")

## 
## Key genes summary saved

Gene Expression Visualization

Function to Plot Gene Expression by Temperature and Donor

#' Plot gene expression across temperatures, faceted by Donor
#'
#' Creates boxplot + beeswarm visualization comparing genotypes across
#' temperatures and donors, with statistical tests.
#'
#' @param gene_id Gene identifier (must match rownames in voom_obj$E)
#' @param tissue_name Tissue to plot (e.g., "V3_SAM")
#' @param voom_obj Voom object containing log2CPM expression data
#' @param metadata Sample metadata with columns: Run, Tissue, Genotype, Donor, temp
#'
#' @return ggplot object
plot_gene_temp_comparison_by_donor <- function(gene_id, 
                                                tissue_name,
                                                voom_obj, 
                                                metadata) {
  
  stopifnot(
    gene_id %in% rownames(voom_obj$E),
    tissue_name %in% metadata$Tissue
  )
  
  expr_data <- data.frame(
    sample = colnames(voom_obj$E),
    log2cpm = voom_obj$E[gene_id, ],
    stringsAsFactors = FALSE
  )
  
  plot_data <- metadata %>%
    filter(Tissue == tissue_name) %>%
    left_join(expr_data, by = c("Run" = "sample")) %>%
    mutate(
      temp_label = ifelse(temp == "CTRL", "warm", "cold"),
      temp_label = factor(temp_label, levels = c("warm", "cold"))
    ) %>%
    select(Run, Genotype, Donor, temp_label, log2cpm)
  
  stat_test <- plot_data %>%
    group_by(Donor, temp_label) %>%
    t_test(log2cpm ~ Genotype) %>%
    add_significance() %>%
    add_xy_position(x = "Genotype")
  
  cat("\n========================================\n")
  cat("Gene:", gene_id, "| Tissue:", tissue_name, "\n")
  cat("========================================\n")
  print(stat_test)
  
  p <- ggplot(plot_data, aes(x = Genotype, y = log2cpm, color = Genotype)) +
    geom_boxplot(width = 0.5, outlier.shape = NA, alpha = 0.7) +
    geom_beeswarm(size = 2, alpha = 0.8, cex = 2) +
    scale_color_manual(
      values = c("CTRL" = "gold", "Inv4m" = "purple4")
    ) +
    facet_grid(temp_label ~ Donor) +
    labs(
      title = paste0(gene_id, " expression in ", tissue_name),
      subtitle = "Warm vs Cold by Donor",
      x = "Genotype",
      y = "log2(CPM)"
    ) +
    stat_pvalue_manual(
      stat_test,
      label = "p = {p}",
      tip.length = 0.01,
      bracket.size = 0.5,
      size = 3.5
    ) +
    theme_classic2(base_size = 25) +
    theme(
      legend.position = "none",
      plot.title = element_text(face = "bold"),
      strip.background = element_blank(),
      panel.grid.major.x = element_blank()
    )
  
  return(p)
}

Usage Examples

V3 sampling.

JMJ

# Single gene
p1 <- plot_gene_temp_comparison_by_donor(
  gene_id = "Zm00001eb191820",
  tissue_name = "V3_S2_Leaf_Tip",
  voom_obj = voom_obj,
  metadata = y$samples
) + ggtitle("JMJ4 Zm00001eb191820  V3_S2_Leaf_Tip")

## 
## ========================================
## Gene: Zm00001eb191820 | Tissue: V3_S2_Leaf_Tip 
## ========================================
## # A tibble: 4 × 15
##   Donor temp_label .y.     group1 group2    n1    n2 statistic    df        p
##   <fct> <fct>      <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl>    <dbl>
## 1 PT    warm       log2cpm CTRL   Inv4m      6     4      6.30  7.24 0.000353
## 2 PT    cold       log2cpm CTRL   Inv4m      4     3      5.56  3.35 0.00855 
## 3 Mi21  warm       log2cpm CTRL   Inv4m      6     5      4.59  4.55 0.00741 
## 4 Mi21  cold       log2cpm CTRL   Inv4m      4     4      4.62  3.03 0.0186  
## # ℹ 5 more variables: p.signif <chr>, y.position <dbl>, groups <named list>,
## #   xmin <dbl>, xmax <dbl>

print(p1)

# Single gene
p1 <- plot_gene_temp_comparison_by_donor(
  gene_id = "Zm00001eb191820",
  tissue_name = "V3_SAM",
  voom_obj = voom_obj,
  metadata = y$samples
) + ggtitle("JMJ4 Zm00001eb191820  V3_SAM")

## 
## ========================================
## Gene: Zm00001eb191820 | Tissue: V3_SAM 
## ========================================
## # A tibble: 4 × 15
##   Donor temp_label .y.     group1 group2    n1    n2 statistic    df        p
##   <fct> <fct>      <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl>    <dbl>
## 1 PT    warm       log2cpm CTRL   Inv4m      7     5      8.34  5.21 0.000331
## 2 PT    cold       log2cpm CTRL   Inv4m      5     4      2.33  3.27 0.0948  
## 3 Mi21  warm       log2cpm CTRL   Inv4m      6     6      4.83  7.84 0.00139 
## 4 Mi21  cold       log2cpm CTRL   Inv4m      3     4      7.77  4.92 0.000606
## # ℹ 5 more variables: p.signif <chr>, y.position <dbl>, groups <named list>,
## #   xmin <dbl>, xmax <dbl>

print(p1)

# Single gene
p2 <- plot_gene_temp_comparison_by_donor(
  gene_id = "Zm00001eb191820",
  tissue_name = "V1_Root",
  voom_obj = voom_obj,
  metadata = y$samples
) + ggtitle("JMJ4 Zm00001eb191820 V1_Root")

## 
## ========================================
## Gene: Zm00001eb191820 | Tissue: V1_Root 
## ========================================
## # A tibble: 4 × 15
##   Donor temp_label .y.     group1 group2    n1    n2 statistic    df        p
##   <fct> <fct>      <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl>    <dbl>
## 1 PT    warm       log2cpm CTRL   Inv4m      6     5      5.42  4.52 0.00395 
## 2 PT    cold       log2cpm CTRL   Inv4m      6     4      4.90  3.31 0.0128  
## 3 Mi21  warm       log2cpm CTRL   Inv4m      5     5      7.35  5.75 0.000395
## 4 Mi21  cold       log2cpm CTRL   Inv4m      5     5      6.51  7.20 0.000292
## # ℹ 5 more variables: p.signif <chr>, y.position <dbl>, groups <named list>,
## #   xmin <dbl>, xmax <dbl>

print(p2)

SEC6

# Single gene
p1 <- plot_gene_temp_comparison_by_donor(
  gene_id = "Zm00001eb194380",
  tissue_name = "V3_S2_Leaf_Tip",
  voom_obj = voom_obj,
  metadata = y$samples
) + ggtitle("SEC6 Zm00001eb193790 V3_S2_Leaf_Tip")

## 
## ========================================
## Gene: Zm00001eb194380 | Tissue: V3_S2_Leaf_Tip 
## ========================================
## # A tibble: 4 × 15
##   Donor temp_label .y.     group1 group2    n1    n2 statistic    df           p
##   <fct> <fct>      <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl>       <dbl>
## 1 PT    warm       log2cpm CTRL   Inv4m      6     4    -11.9   7.66 0.00000337 
## 2 PT    cold       log2cpm CTRL   Inv4m      4     3     -6.19  2.15 0.0211     
## 3 Mi21  warm       log2cpm CTRL   Inv4m      6     5    -17.2   8.05 0.000000125
## 4 Mi21  cold       log2cpm CTRL   Inv4m      4     4    -15.1   3.79 0.000157   
## # ℹ 5 more variables: p.signif <chr>, y.position <dbl>, groups <named list>,
## #   xmin <dbl>, xmax <dbl>

print(p1)

# Single gene
p1 <- plot_gene_temp_comparison_by_donor(
  gene_id = "Zm00001eb194380",
  tissue_name = "V3_SAM",
  voom_obj = voom_obj,
  metadata = y$samples
) + ggtitle("SEC6 Zm00001eb193790 V3_SAM")

## 
## ========================================
## Gene: Zm00001eb194380 | Tissue: V3_SAM 
## ========================================
## # A tibble: 4 × 15
##   Donor temp_label .y.     group1 group2    n1    n2 statistic    df           p
##   <fct> <fct>      <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl>       <dbl>
## 1 PT    warm       log2cpm CTRL   Inv4m      7     5    -11.7   9.87 0.000000403
## 2 PT    cold       log2cpm CTRL   Inv4m      5     4     -2.47  3.64 0.0752     
## 3 Mi21  warm       log2cpm CTRL   Inv4m      6     6    -14.9   7.72 0.000000591
## 4 Mi21  cold       log2cpm CTRL   Inv4m      3     4    -21.8   4.96 0.00000409 
## # ℹ 5 more variables: p.signif <chr>, y.position <dbl>, groups <named list>,
## #   xmin <dbl>, xmax <dbl>

print(p1)

# Single gene
p2 <- plot_gene_temp_comparison_by_donor(
  gene_id = "Zm00001eb194380",
  tissue_name = "V1_Root",
  voom_obj = voom_obj,
  metadata = y$samples
) + ggtitle("SEC6 Zm00001eb194380 V1_Root")

## 
## ========================================
## Gene: Zm00001eb194380 | Tissue: V1_Root 
## ========================================
## # A tibble: 4 × 15
##   Donor temp_label .y.     group1 group2    n1    n2 statistic    df          p
##   <fct> <fct>      <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl>      <dbl>
## 1 PT    warm       log2cpm CTRL   Inv4m      6     5     -8.20  7.54 0.0000511 
## 2 PT    cold       log2cpm CTRL   Inv4m      6     4     -7.23  7.94 0.000093  
## 3 Mi21  warm       log2cpm CTRL   Inv4m      5     5     -7.73  6.75 0.000137  
## 4 Mi21  cold       log2cpm CTRL   Inv4m      5     5    -10.2   8.00 0.00000749
## # ℹ 5 more variables: p.signif <chr>, y.position <dbl>, groups <named list>,
## #   xmin <dbl>, xmax <dbl>

print(p2)

PCNA2

# Single gene
p2 <- plot_gene_temp_comparison_by_donor(
  gene_id = "Zm00001eb192050",
  tissue_name = "V3_S2_Leaf_Tip",
  voom_obj = voom_obj,
  metadata = y$samples
) + ggtitle("PCNA2 Zm00001eb192050 V3_S2_Leaf_Tip")

## 
## ========================================
## Gene: Zm00001eb192050 | Tissue: V3_S2_Leaf_Tip 
## ========================================
## # A tibble: 4 × 15
##   Donor temp_label .y.     group1 group2    n1    n2 statistic    df     p
##   <fct> <fct>      <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl>
## 1 PT    warm       log2cpm CTRL   Inv4m      6     4     1.41   4.25 0.228
## 2 PT    cold       log2cpm CTRL   Inv4m      4     3    -1.74   3.16 0.176
## 3 Mi21  warm       log2cpm CTRL   Inv4m      6     5    -0.757  8.63 0.469
## 4 Mi21  cold       log2cpm CTRL   Inv4m      4     4    -1.77   4.33 0.146
## # ℹ 5 more variables: p.signif <chr>, y.position <dbl>, groups <named list>,
## #   xmin <dbl>, xmax <dbl>

print(p2)

# Single gene
p2 <- plot_gene_temp_comparison_by_donor(
  gene_id = "Zm00001eb192050",
  tissue_name = "V3_SAM",
  voom_obj = voom_obj,
  metadata = y$samples
) + ggtitle("PCNA2 Zm00001eb192050 V3_SAM")

## 
## ========================================
## Gene: Zm00001eb192050 | Tissue: V3_SAM 
## ========================================
## # A tibble: 4 × 15
##   Donor temp_label .y.     group1 group2    n1    n2 statistic    df     p
##   <fct> <fct>      <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl>
## 1 PT    warm       log2cpm CTRL   Inv4m      7     5     0.322  5.39 0.759
## 2 PT    cold       log2cpm CTRL   Inv4m      5     4    -0.448  4.49 0.675
## 3 Mi21  warm       log2cpm CTRL   Inv4m      6     6    -1.26   5.54 0.257
## 4 Mi21  cold       log2cpm CTRL   Inv4m      3     4    -0.921  2.19 0.447
## # ℹ 5 more variables: p.signif <chr>, y.position <dbl>, groups <named list>,
## #   xmin <dbl>, xmax <dbl>

print(p2)

# Single gene
p3 <- plot_gene_temp_comparison_by_donor(
  gene_id = "Zm00001eb192050",
  tissue_name = "V1_Root",
  voom_obj = voom_obj,
  metadata = y$samples
) + ggtitle("PCNA2 Zm00001eb192050 V1_Root")

## 
## ========================================
## Gene: Zm00001eb192050 | Tissue: V1_Root 
## ========================================
## # A tibble: 4 × 15
##   Donor temp_label .y.     group1 group2    n1    n2 statistic    df     p
##   <fct> <fct>      <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl>
## 1 PT    warm       log2cpm CTRL   Inv4m      6     5    -0.237  8.46 0.819
## 2 PT    cold       log2cpm CTRL   Inv4m      6     4     1.27   3.06 0.293
## 3 Mi21  warm       log2cpm CTRL   Inv4m      5     5    -0.986  6.23 0.361
## 4 Mi21  cold       log2cpm CTRL   Inv4m      5     5    -1.86   5.13 0.121
## # ℹ 5 more variables: p.signif <chr>, y.position <dbl>, groups <named list>,
## #   xmin <dbl>, xmax <dbl>

print(p3)

MCM5

# Single gene
p3 <- plot_gene_temp_comparison_by_donor(
  gene_id = "Zm00001eb257900",
  tissue_name = "V3_S2_Leaf_Tip",
  voom_obj = voom_obj,
  metadata = y$samples
) + ggtitle("MCM5 Zm00001eb257900 V3_S2_Leaf_Tip")

## 
## ========================================
## Gene: Zm00001eb257900 | Tissue: V3_S2_Leaf_Tip 
## ========================================
## # A tibble: 4 × 15
##   Donor temp_label .y.     group1 group2    n1    n2 statistic    df     p
##   <fct> <fct>      <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl>
## 1 PT    warm       log2cpm CTRL   Inv4m      6     4     1.58   6.53 0.161
## 2 PT    cold       log2cpm CTRL   Inv4m      4     3    -0.997  4.79 0.366
## 3 Mi21  warm       log2cpm CTRL   Inv4m      6     5    -1.08   8.93 0.31 
## 4 Mi21  cold       log2cpm CTRL   Inv4m      4     4    -1.01   3.15 0.383
## # ℹ 5 more variables: p.signif <chr>, y.position <dbl>, groups <named list>,
## #   xmin <dbl>, xmax <dbl>

print(p3)

# Single gene
p3 <- plot_gene_temp_comparison_by_donor(
  gene_id = "Zm00001eb257900",
  tissue_name = "V3_SAM",
  voom_obj = voom_obj,
  metadata = y$samples
) + ggtitle("MCM5 Zm00001eb257900 V3_SAM")

## 
## ========================================
## Gene: Zm00001eb257900 | Tissue: V3_SAM 
## ========================================
## # A tibble: 4 × 15
##   Donor temp_label .y.     group1 group2    n1    n2 statistic    df     p
##   <fct> <fct>      <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl>
## 1 PT    warm       log2cpm CTRL   Inv4m      7     5    -0.173  5.07 0.869
## 2 PT    cold       log2cpm CTRL   Inv4m      5     4    -0.721  4.47 0.507
## 3 Mi21  warm       log2cpm CTRL   Inv4m      6     6    -1.19   8.64 0.266
## 4 Mi21  cold       log2cpm CTRL   Inv4m      3     4    -0.546  2.12 0.637
## # ℹ 5 more variables: p.signif <chr>, y.position <dbl>, groups <named list>,
## #   xmin <dbl>, xmax <dbl>

print(p3)

# Single gene
p3 <- plot_gene_temp_comparison_by_donor(
  gene_id = "Zm00001eb257900",
  tissue_name = "V1_Root",
  voom_obj = voom_obj,
  metadata = y$samples
) + ggtitle("MCM5 Zm00001eb257900 V1_Root")

## 
## ========================================
## Gene: Zm00001eb257900 | Tissue: V1_Root 
## ========================================
## # A tibble: 4 × 15
##   Donor temp_label .y.     group1 group2    n1    n2 statistic    df     p
##   <fct> <fct>      <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl>
## 1 PT    warm       log2cpm CTRL   Inv4m      6     5     0.132  8.90 0.898
## 2 PT    cold       log2cpm CTRL   Inv4m      6     4     1.37   3.09 0.262
## 3 Mi21  warm       log2cpm CTRL   Inv4m      5     5    -0.351  5.95 0.738
## 4 Mi21  cold       log2cpm CTRL   Inv4m      5     5     0.701  7.20 0.506
## # ℹ 5 more variables: p.signif <chr>, y.position <dbl>, groups <named list>,
## #   xmin <dbl>, xmax <dbl>

print(p3)

ATG8/Mediator/Fibrillarin-2

# Single gene
p3 <- plot_gene_temp_comparison_by_donor(
  gene_id = "Zm00001eb191090",
  tissue_name = "V3_S2_Leaf_Tip",
  voom_obj = voom_obj,
  metadata = y$samples
) + ggtitle("ATG8 Zm00001eb191090 V3_S2_Leaf_Tip")

## 
## ========================================
## Gene: Zm00001eb191090 | Tissue: V3_S2_Leaf_Tip 
## ========================================
## # A tibble: 4 × 15
##   Donor temp_label .y.     group1 group2    n1    n2 statistic    df      p
##   <fct> <fct>      <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl>  <dbl>
## 1 PT    warm       log2cpm CTRL   Inv4m      6     4   -2.12    7.99 0.0667
## 2 PT    cold       log2cpm CTRL   Inv4m      4     3   -1.53    3.16 0.219 
## 3 Mi21  warm       log2cpm CTRL   Inv4m      6     5   -0.0834  8.64 0.935 
## 4 Mi21  cold       log2cpm CTRL   Inv4m      4     4   -0.354   5.38 0.737 
## # ℹ 5 more variables: p.signif <chr>, y.position <dbl>, groups <named list>,
## #   xmin <dbl>, xmax <dbl>

print(p3)

# Single gene
p3 <- plot_gene_temp_comparison_by_donor(
  gene_id = "Zm00001eb191090",
  tissue_name = "V3_SAM",
  voom_obj = voom_obj,
  metadata = y$samples
) + ggtitle("ATG8 Zm00001eb191090 V3_SAM")

## 
## ========================================
## Gene: Zm00001eb191090 | Tissue: V3_SAM 
## ========================================
## # A tibble: 4 × 15
##   Donor temp_label .y.     group1 group2    n1    n2 statistic    df      p
##   <fct> <fct>      <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl>  <dbl>
## 1 PT    warm       log2cpm CTRL   Inv4m      7     5    -1.75   8.58 0.115 
## 2 PT    cold       log2cpm CTRL   Inv4m      5     4    -0.908  5.93 0.399 
## 3 Mi21  warm       log2cpm CTRL   Inv4m      6     6    -1.65   5.79 0.151 
## 4 Mi21  cold       log2cpm CTRL   Inv4m      3     4    -2.78   2.88 0.0724
## # ℹ 5 more variables: p.signif <chr>, y.position <dbl>, groups <named list>,
## #   xmin <dbl>, xmax <dbl>

print(p3)

# Single gene
p3 <- plot_gene_temp_comparison_by_donor(
  gene_id = "Zm00001eb191090",
  tissue_name = "V1_Root",
  voom_obj = voom_obj,
  metadata = y$samples
) + ggtitle("ATG8 Zm00001eb191090 V1_Root")

## 
## ========================================
## Gene: Zm00001eb191090 | Tissue: V1_Root 
## ========================================
## # A tibble: 4 × 15
##   Donor temp_label .y.     group1 group2    n1    n2 statistic    df       p
##   <fct> <fct>      <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl>   <dbl>
## 1 PT    warm       log2cpm CTRL   Inv4m      6     5    -1.33   9.00 0.217  
## 2 PT    cold       log2cpm CTRL   Inv4m      6     4     0.469  3.89 0.664  
## 3 Mi21  warm       log2cpm CTRL   Inv4m      5     5    -0.846  5.17 0.435  
## 4 Mi21  cold       log2cpm CTRL   Inv4m      5     5    -3.97   5.89 0.00766
## # ℹ 5 more variables: p.signif <chr>, y.position <dbl>, groups <named list>,
## #   xmin <dbl>, xmax <dbl>

print(p3)

# # Multiple genes
# genes_to_plot <- c("Zm00001eb191820", "PCNA2", "MCM5")
# 
# plots <- lapply(genes_to_plot, function(gene) {
#   if (gene %in% rownames(voom_obj$E)) {
#     plot_gene_temp_comparison_by_donor(
#       gene_id = gene,
#       tissue_name = "V3_SAM",
#       voom_obj = voom_obj,
#       metadata = y$samples
#     )
#   }
# })
# 
# plots <- plots[!sapply(plots, is.null)]
# 
# wrap_plots(plots, ncol = 1)

# # Same gene across tissues
# tissues <- c("V3_SAM", "V3_S2_Leaf_Tip", "V1_Root")
# 
# tissue_plots <- lapply(tissues, function(tissue) {
#   plot_gene_temp_comparison_by_donor(
#     gene_id = "Zm00001eb191820",
#     tissue_name = tissue,
#     voom_obj = voom_obj,
#     metadata = y$samples
#   )
# })
# 
# wrap_plots(tissue_plots, ncol = 1)

Session Info

sessionInfo()

## R version 4.5.1 (2025-06-13)
## Platform: x86_64-apple-darwin20
## Running under: macOS Sequoia 15.6.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] pheatmap_1.0.13  readr_2.1.6      ggpubr_0.6.2     rstatix_0.7.3   
##  [5] ggbeeswarm_0.7.2 ggplot2_4.0.1    stringr_1.6.0    tidyr_1.3.1     
##  [9] dplyr_1.1.4      edgeR_4.8.0      limma_3.66.0    
## 
## loaded via a namespace (and not attached):
##  [1] utf8_1.2.6         sass_0.4.10        generics_0.1.4     stringi_1.8.7     
##  [5] lattice_0.22-7     hms_1.1.4          digest_0.6.39      magrittr_2.0.4    
##  [9] evaluate_1.0.5     grid_4.5.1         RColorBrewer_1.1-3 fastmap_1.2.0     
## [13] jsonlite_2.0.0     backports_1.5.0    Formula_1.2-5      purrr_1.2.0       
## [17] scales_1.4.0       jquerylib_0.1.4    abind_1.4-8        cli_3.6.5         
## [21] crayon_1.5.3       rlang_1.1.6        bit64_4.6.0-1      withr_3.0.2       
## [25] cachem_1.1.0       yaml_2.3.10        parallel_4.5.1     tools_4.5.1       
## [29] tzdb_0.5.0         ggsignif_0.6.4     locfit_1.5-9.12    broom_1.0.10      
## [33] png_0.1-8          vctrs_0.6.5        R6_2.6.1           lifecycle_1.0.4   
## [37] bit_4.6.0          car_3.1-3          vroom_1.6.6        vipor_0.4.7       
## [41] pkgconfig_2.0.3    beeswarm_0.4.0     pillar_1.11.1      bslib_0.9.0       
## [45] gtable_0.3.6       glue_1.8.0         statmod_1.5.1      xfun_0.54         
## [49] tibble_3.3.0       tidyselect_1.2.1   rstudioapi_0.17.1  knitr_1.50        
## [53] farver_2.1.2       htmltools_0.5.8.1  labeling_0.4.3     rmarkdown_2.30    
## [57] carData_3.0-5      compiler_4.5.1     S7_0.2.1

Inv4m field vs growth chamber transcriptome: Crow2020 Reanalysis II

Fausto Rodríguez-Zapata

2025-12-04

Biological Context & Analysis Goals

Research Hypothesis

Specific Questions

Statistical Approach

Setup

Data Loading

Expression Counts

Sample Metadata

Verify Design Balance

Prepare Count Matrix

DGEList and Filtering

Create DGEList with Cell Means Condition Factor

Filter Low-Expression Genes

TMM Normalization

Design Matrix: Cell Means Model

Build Design (No Reference Levels for Biological Factors)

Voom Transformation

Contrasts: Explicit Biological Comparisons

Define Contrasts

Model Fitting

Fit Linear Models

Results Extraction

Extract TopTables for All Contrasts

Summary: Significant DEGs per Contrast

Question 1: Field-to-Chamber Replication

Load Field DEG List

Test Replication in Chamber Leaf Tissues

Enrichment Test: Field DEGs in Chamber Leaves

Question 2: SAM Hypothesis

Test Field DEGs in SAM

Enrichment: Field DEGs in SAM

Question 3: Temperature Sensitivity

PCNA2 & MCM5 Temperature Dependence in SAM

Global Temperature Dependence Pattern

Question 4: Root Validation

PCNA2 & MCM5 in V1_Root

Visualization

Volcano Plots: Key Contrasts

Heatmap: Field DEGs Across Conditions

Export Results

Save TopTables

Save Key Gene Summary Table

Gene Expression Visualization

Function to Plot Gene Expression by Temperature and Donor

Usage Examples

JMJ

SEC6

PCNA2

MCM5

ATG8/Mediator/Fibrillarin-2

Session Info