Hello everyone, I've been working on an RNAseq dataset, where I infected cells with virus and extracted RNA at different timepoints (3hpi, 6hpi, 9hpi, 24hpi) after infection to study DEG changes over time. I compare the 3hpi, 6hpi and 9hpi samples all to a 0hpi control sample, but although I set the p.value cutoff as 0.05 when I visualise the DEGs with Volcanoplot, it seems the cutoff is different for each of the contrasts. Is that normal?

Here is my R code (from after the lmFit() function) for the first thrre timepoints. For 24hpi the code is essentially the same, except that I compare 24hpi to a mock infected sample 24hpi, since this sample and the first three are from two separate experiments, but with the same virus and cells:

contr <- makeContrasts("3hpi" = group3-group0, "6hpi" = group6-group0, "9hpi" = group9-group0, levels = mm)
fit2 <- contrasts.fit(fit, contr)
fit3 <- limma::eBayes(fit2)
Results <- decideTests(fit3, method= "separate", adjust.method="BH", p.value=0.05, lfc = log2(1.5))

plotSA(fit3, main = "Final Model: Mean-variance Trend") #Final Model: Mean-variance Trend

summary(Results)

head(fit3)


#The Results object because it is not a data frame with named columns. 
#Instead, you can use square brackets ([ ]) to subset the column.
table(Results[, '9hpi'])

unique(volcano_data$Significant)

#replace with whatever timepoint
comparison_name <- "9hpi"

# Convert levels in Results to characters
Results_char <- as.character(Results[, comparison_name])

# Create a data frame with log2 fold changes and -log10(p-values)
volcano_data <- data.frame(
  logFC = fit3$coefficients[, comparison_name],
  logPvalue = -log10(fit3$p.value[, comparison_name]),
  Significant = as.character(Results_char) # Convert to character
)

# Set the desired p-value cutoff
p_value_cutoff <- 0.05

# Create a volcano plot using ggplot2
ggplot(volcano_data, aes(x = logFC, y = logPvalue, color = Significant)) +
  geom_point(alpha = 0.5) +
  scale_color_manual(values = c("-1" = "#377DC5", "0" = "grey", "1" = "#FF5733"), name = "Significant") +
  theme_minimal() +
  labs(
    title = "0hpi vs 9hpi",
    x = "log2 Fold Change",
    y = "-log10(p-value)"
  )+
  # Add horizontal line for p-value cutoff
  geom_hline(yintercept = -log10(p_value_cutoff), linetype = "dashed", color = "gray")

Also, here the resulting volcanoplots

I would be incredibly grateful if anyone knew whether those seemingly differing p-value cutoffs are normal when doing separate comparisons like this.

Cheers and thanks in advance for reading!