How to account for solvent controls in DESeq2 when comparing follicular vs luteal phase?

Question

[DareDevil]

I have an RNA-seq dataset with the following samples (genes as rows, samples as columns):

AA - FP SC T1   (Follicular phase, solvent control, replicate 1)
AB - FP T1      (Follicular phase, treated, replicate 1)
AC - FP SC T2   (Follicular phase, solvent control, replicate 2)
AD - FP T2      (Follicular phase, treated, replicate 2)
AG - LP SC T1   (Luteal phase, solvent control, replicate 1)
AE - LP T1      (Luteal phase, treated, replicate 1)
AF - LP SC T2   (Luteal phase, solvent control, replicate 2)
AH - LP T2      (Luteal phase, treated, replicate 2)

So I have two biological replicates (T1, T2) for each phase (FP, LP), and each replicate has both a solvent control (SC) and a treated sample. I want to perform differential expression analysis between FP and LP, but I want to remove the effect of solvent controls so that I only see biological differences between the phases.

Possible approaches I considered

1. Subtract SC from treated counts within each replicate before running DESeq2 (e.g., AB–AA = FP_T1_corrected).: This reduces to 4 corrected samples: FP_T1, FP_T2, LP_T1, LP_T2. But here I am concerned about whether the subtraction can create negative counts, and only 2 replicates per group remain.

2. Model SC directly in DESeq2 using an interaction term: Metadata includes phase (FP/LP), condition (SC/Treated), and replicate (T1/T2).

Design: ~ replicate + phase * condition.

Then use the interaction term (phaseLP.conditionTreated) to test whether LP vs FP differences hold after adjusting for SC.

What is the best practice in this situation?

Should I subtract SC counts before DESeq2 (Approach 1)?

Or should I keep all samples and model SC in the design formula with an interaction term (Approach 2)?

Are there other recommended ways to handle this setup?

Thanks for your help

Answer 1

[swbarnes2]

You generally can't correct for replicates in such a tiny experiment like this, so I'd drop that element.

But yes, if you want to do an a/b - c/d comparison, interactions is how you do that. You don't muck with the counts.

[DareDevil]

I have checked DEGs between Solvent Control (SC) of FP and LP and there are noticebale differences in terms of DEGs.

# Checking SC in FP and LP
# Read featureCounts table
counts = read.table("counts.txt", header=TRUE, sep="\t", row.names=1)

# Metadata
targets = data.frame(ID = c("YHDAA","YHDAB","YHDAC","YHDAD","YHDAE","YHDAF","YHDAG","YHDAH"),
                     Name = c("FP_SC_T1","FP_T1","FP_SC_T2","FP_T2", "LP_SC_T1","LP_T1","LP_SC_T2","LP_T2"),
                     phase = c("FP","FP","FP","FP","LP","LP","LP","LP"),
                     Group = c("FP_SC","FP","FP_SC","FP","LP_SC","LP","LP_SC","LP"),
                     condition = c("SC","Treated","SC","Treated","SC","Treated","SC","Treated"))

#converting 'NA' values to 0
counts[is.na(counts)] <- 0
#make sample_info and counts in same order
counts = counts[, targets$ID]

# Remove genes with 0 counts
keep_genes <- rowSums(counts(dds)) > 0
dds <- dds[ keep_genes, ]
# genes to have more than a sum total of 10 reads of support in all samples
dds = dds[rowSums(counts(dds)) > 10, ]

dds = DESeq(dds, betaPrior=FALSE, parallel=TRUE, BPPARAM=SnowParam(5))
res = results(dds, contrast = c("Group", "FP_SC", "LP_SC"), parallel = TRUE)

This has identified quite a good number of DEGs.

Next approach,

#D ESeq2 object with interaction design
dds = DESeqDataSetFromMatrix(countData = counts, colData = sample_info, design = ~ phase * condition)
# Equivalent to: ~ phase + condition + phase:condition

How do I calculate the res here by negating the SC effect?

res <- results(dds, name = "phaseLP.conditionTreated")