Hello,

I am looking for advice about how to correctly use the "coef" argument to lfcshrink when I am using DESeq2 with a grouping variable. I am analyzing an RNAseq dataset with multiple factors and variables using DESeq2:

• SurgeryGroup (injury versus sham)
• TimePoint (1 versus 2)
• Sex (M versus F)

My exploratory analysis of the data indicated that there is (surprisingly) very little effect of Sex so I control for it in my model and focus on the other variables. I am interested in asking the following questions:

1. What genes are differentially expressed in injury versus sham at timepoint 1?

2. What genes are differentially expressed in injury versus sham at timepoint 2?

To this end, I have created a grouping variable with levels corresponding to the 4 possible combinations of SurgeryGroup and TimePoint. I can then call results and shrink the LFCs by specifying my contrasts of interest.

> dds$group <- factor(paste0(dds$SurgeryGroup, dds$TimePoint))
> design(dds) <- ~ Sex + group
> dds <- DESeq(dds)
> resultsNames(dds)
[1] "Intercept"                "Sex_M_vs_F"               "group_sham1_vs_sham2"
[4] "group_injury2_vs_sham2" "group_injury1_vs_sham2"

> resultsTimePoint1 <- results(dds, contrast=c("group", "injury1", "sham1"), alpha=0.05)
> resultsTimePoint1 <- lfcShrink(dds, contrast=c("group", "injury1", "sham1"), res=resultsTimePoint1)

However, I'm interested in using "apeglm" to shrink the LFCs, which requires the "coef" argument to be used instead of "contrast" for lfcshrink. Is it possible to do this with a grouping variable as I have done? Perhaps I am just confused about how to specify the correct contrasts I am interested in by its coefficient, so if that's the case, I'd appreciate advice about how to do this.

Thanks for your consideration!

Here is some example code.

You'd just add sex to the beginning of the design to adapt to your code, ~sex + time + time:condition.

> dds <- makeExampleDESeqDataSet()
> dds$condition
 [1] A A A A A A B B B B B B
Levels: A B
> dds$time <- factor(rep(c(1,2,1,2),each=3))

> design(dds) <- ~time + time:condition

> dds <- DESeq(dds, quiet=TRUE)
> resultsNames(dds)
[1] "Intercept"        "time_2_vs_1"      "time1.conditionB" "time2.conditionB"

> res1 <- lfcShrink(dds, coef=3, type="apeglm")
> res2 <- lfcShrink(dds, coef=4, type="apeglm")

These last two could equivalently have been:

> res1 <- lfcShrink(dds, coef="time1.conditionB", type="apeglm")
> res2 <- lfcShrink(dds, coef="time2.conditionB", type="apeglm")

These two are the condition B vs A effect for time 1 and for time 2, respectively.