I'm having a bit of trouble understanding how `lfcThreshold`

parameter of the `results`

function of `DESeq2`

affects the p-value. Imposing a stricter `lfcThreshold`

reduces the number of significant results, which makes sense, however, it does so above and beyond what I would have expected.

From my data:

lfc <- 0 thrHuh <- results(dds_huh7, contrast = c("group", "infected_3", "uninfected_3"), alpha = alpha, lfcThreshold = lfc) summary(thrHuh) out of 72860 with nonzero total read count adjusted p-value < 0.05 LFC > 0 (up) : 1545, 2.1% LFC < 0 (down) : 1278, 1.8% outliers [1] : 39, 0.054% low counts [2] : 42376, 58% (mean count < 6) sig <- subset(thrHuh, thrHuh$padj < 0.05) nrow(subset(sig, sig$log2FoldChange >= 2 | sig$log2FoldChange <= -2)) [1] 495

My interpretation of this is that I have 495 (significant) genes that have a greater than 4-fold change. My (obviously incorrect) understanding is that if I set `lfcThreshold <- 2`

I would get the same results, but I don't:

lfc <- 2 thrHuh <- results(dds_huh7, contrast = c("group", "infected_3", "uninfected_3"), alpha = alpha, lfcThreshold = lfc) summary(thrHuh) out of 72860 with nonzero total read count adjusted p-value < 0.05 LFC > 0 (up) : 21, 0.029% LFC < 0 (down) : 3, 0.0041% outliers [1] : 39, 0.054% low counts [2] : 8476, 12% (mean count < 0)

I've looked at the DESEq2 paper, the vignette, and the workflow, and cannot figure why one method gives 495 and the other 24. Is it simply because there are fewer tests being run in the latter case, and this is reflected in the adjusted p value? If I only want to consider genes that are a) significant and b) have at least a 4-fold change or more, which is the better one to use?

Thanks,

Russ

edits: formatting mistakes...

Much clearer. Thanks for the detailed explanation! I will lower the lfcThreshold - was just using 2 as it provided an extreme illustration of the effect.

Thanks James.

And to connect to the DESeq2 paper, this is the relevant paragraph, under the section '

Hypothesis tests with thresholds on effect size':A common procedure is to disregard genes whose estimated LFC

βiris below some threshold, |βir|≤θ. However, this approach loses the benefit of an easily interpretable FDR, as the reportedPvalue and adjustedPvalue still correspond to the test ofzeroLFC. It is therefore desirable to include the threshold in the statistical testing procedure directly, i.e., not to filter post hoc on a reported fold-changeestimate, but rather to evaluate statistically directly whether there is sufficient evidence that the LFC is above the chosen threshold.Very nice explanation, I have just one question. If I have understood you, to apply somekind of LFC filter it must be done by the addition of lfcThreshold option. In this way we conserve the real meaning of the padj result of our FDR test, but if we use:

sig_Wild_LA <- Wild_LA[which(Wild_LA$log2FoldChange > 1.99999), ] (example)

My subset of DGE now has no sense padj values.

Is this correct?

I have this doubt because I found a lot of tutorial including this filtering step and not even one using the lfcThreshold option.

Yes, you've understood correctly. There is a lot of bad advice out there, and one of the most common bits of bad advice is filtering based on fold change separately from significance testing. If you filter fold changes in this way, the FDR estimation has no knowledge of your fold change filtering criterion, so there's no way that the reported FDR could possibly have taken your filtering criterion into account. Depending on the data and the specific way in which you implement this flawed double filtering strategy, you could end up with FDR values that either greatly overestimate or underestimate the true significance.

Thank you Ryan. I'm not surprised, I was working with RTqPCR and data analysis is a blackhole even for a technique with more than 15 years-old...

RNAseq papers usually has a M&M quite schematic and I have seen a lot of people using software in a wrong way... without any problem to cross "referee-barrier".

Thank you again!