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

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

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.
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

rnaseq replicates deseq2 edger

• 4,373 views • link

written 11 months ago by DareDevil

1 answer

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.

• 0 views • link

written 11 months ago by swbarnes2

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")

• 0 views • link

written 11 months ago by DareDevil