I am currently conducting a differential expression analysis of small RNA-seq data using the DESeq2 package in R, and I have encountered unexpected behavior when using the lfcShrink() function to adjust fold-changes (FC). Specifically, I have noticed that in some cases, the lfcShrink() function appears to inflate the FC values.

I have an experimental design with 5 groups: a control group and 4 from different conditions. I filter the matrix of counts for each unique sequence so that I keep those sequences that have at least 5 counts in more than half of the samples in at least 1 group.

dds <- DESeqDataSetFromMatrix(countData = counts_table_final,
                                  colData = sampletable,
                                  design = ~ Group)
dds$Group <- relevel(dds$Group, ref ="C")
dds_2 <- DESeq(dds)

res_test <- results(dds_2, name ="Group_A2_vs_C", lfcThreshold = 0.585, alpha=0.05)
shrunk_test <- lfcShrink(dds_2, coef="Group_A2_vs_C", res=res_test)

res_test_tb <- res_test %>%
  data.frame() %>%
  rownames_to_column(var="seq") %>% 
  as_tibble()

res_test_tb$lfcShrunk <- shrunk_test$log2FoldChange
res_test_tb$lfcShrunkSE <- shrunk_test$lfcSE

In general, I see that the correction works as expected, especially by decreasing very extreme LFCs. But if I inspect the sequences I see that strange things happen, as some sequences (aprox 13%) present an LFC after shrinkage in absolute value higher than the original LFC. It happens mostly in sequences with low baseMean, that are NS and/or that we cannot annotate as endogenous sRNA (it is important to note that we perform the annotation after DE, and therefore many of those sequences may have no biological relevance).

Although in many cases the difference is not very large, I also observe some extreme cases, e.g., a sequence going from 5.9 to 17.37 LFC after adjustment (A2 vs C), here the counts:

In most sequences I have plotted where I see such radical changes they present very low baseMean and 0 counts in the controls. When the baseMean is higher and the controls present counts, the increase in LFC Shrunken with respect to the original LFC is more subtle.

Although these are mostly sequences of little relevance to my study, I would like to understand why this occurs.

I have carefully reviewed my data and the documentation for the DESeq2 package, but I have not been able to identify the underlying cause of this phenomenon. I suspect that it may be related to the distribution of expression values and the characteristics of the data.

The actual LFC is not defined here, it is some positive value vs 0, so +Inf. So the value provided by any method is a bit of a stand-in for: "cannot be estimated without deeper sequencing".