This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Problem with lfcShrink() in differential gene expression analysis.

I am getting different number of dysregulated genes when performing lfcSrink() on the DESeq run. Sample size: I have 2 groups High(n=16) and Low(n=8). Following is my code for creating the DESeq object:

# Creating DESeq dataset
dds_risk <- DESeqDataSetFromMatrix(
  countData = count.mat %>% 
    column_to_rownames("Geneid"),
  colData = sample.pdata %>% column_to_rownames("Sample.Name"),
  design = ~ Risk.group
)

# Setting reference
dds_risk$Risk.group <- relevel(dds_risk$Risk.group, ref = "Low")

# Run DESeq (quiet)
dds_risk_run <- DESeq(dds_risk)

Now, when I am performing the lfcShrink() using following code I get 159 upregulated and 265 downregulated genes

High_low.result <- results(dds_risk_run, alpha = 0.05, lfcThreshold = 0.585,
                           contrast = c("Risk.group", "High", "Low"))
High_low.resultShrink <- lfcShrink(dds_risk_run,
                                   coef = "Risk.group_High_vs_Low", 
                                   type = "apeglm")

But when I use the following code I am getting only 8 upregulated genes.

High_low.result <- results(dds_risk_run, alpha = 0.05, lfcThreshold = 0.585,
                           contrast = c("Risk.group", "High", "Low"))
High_low.resultShrink <- lfcShrink(dds_risk_run,
                                   coef = "Risk.group_High_vs_Low", 
                                   type = "apeglm",
                                   res = High_low.result)

The only difference in those steps are that I am supplying the results along with DESeq run object and coef in the later method. Which as per documentation should give same result when the results are not supplied, but in my case I am getting different results. Can anybody explain me why this is happening?

lfcsrink dge deseq2 apeglm

1 answer

It is expected. In the first chunk the results() is run internally by lfcShrink() (so with default settings), and the default for lfcThreshold in results() is 0. In the second, you explicitely pass the results() but with a much more stringent lfcThreshold of .585. Hence, in the second you get fewer DEGs.

Personal note on best practices, I would only use the thresholding option to prioritize DEGs in case you have a lot (hundreds, thousands) of them. Here, you probably want to keep defaults to at least have some genes to work with. Making sense out of only 8 genes is hard I guess.

Thanks for the help. It makes sense now.

Log in to answer this question.