1. Disentangling Stress Tolerance: Interaction Modeling

In plant biotechnology, simply identifying genes induced by drought across all plants is insufficient. Both tolerant and sensitive plants wilt under severe water deficit and activate core survival responses. The key biological objective is identifying cultivar-specific tolerance mechanisms: genes that respond differently to drought in the tolerant potato cultivar compared to the sensitive cultivar.

This requires a two-factor interaction model: \[ \log_2(\mu_{ij}) = \beta_0 + \beta_G \cdot \text{Genotype} + \beta_D \cdot \text{Drought} + \beta_{G:D} \cdot (\text{Genotype} \times \text{Drought}) \]

2. Implementing Likelihood Ratio Tests (LRT) in R

We fit a full model including the interaction term, and compare it against a reduced model lacking the interaction using an ANOVA-like Likelihood Ratio Test:

library(DESeq2)

# Set up full factorial design
design(dds) <- ~ Genotype + Drought_Stage + Genotype:Drought_Stage

# Execute Likelihood Ratio Test (LRT)
# The reduced model removes the interaction term: Genotype:Drought_Stage
dds_lrt <- DESeq(dds, test = "LRT", reduced = ~ Genotype + Drought_Stage)

# Extract genes with significant interaction p-values (FDR < 0.01)
res_interaction <- results(dds_lrt, alpha = 0.01)
sig_interaction_genes <- subset(res_interaction, padj < 0.01)

cat(sprintf("Identified %d genes displaying significant Genotype × Drought interaction!\n", nrow(sig_interaction_genes)))

# Flag genes into a comprehensive publication-ready results table
results_annotated <- as.data.frame(res_interaction) %>%
  rownames_to_column(var = "Gene_ID") %>%
  mutate(
    Significance = case_when(
      padj < 0.01 & log2FoldChange >= 1.5  ~ "Upregulated in Tolerant Interaction",
      padj < 0.01 & log2FoldChange <= -1.5 ~ "Downregulated in Tolerant Interaction",
      padj < 0.05                          ~ "Modest Interaction Effect",
      TRUE                                 ~ "Not Significant"
    )
  )

table(results_annotated$Significance)