Detecting differentially expressed genes with foldchange >= 2 and FDR < 0.05

Question

I'm using edgeR for differential analysis. Using glmTreat function I'm detecting differentially expressed genes between Tumor and Normal. I have set the arguments like below:

tr <- glmTreat(fit, contrast=contrast.matrix, lfc=2)
tab2 <- topTags(tr,n=Inf, adjust.method = "BH")
keep <- tab2$table$FDR <= 0.05

This gave me differentially expressed genes with logFC > 2 at FDR < 0.05 and genes with logFC < 2 at FDR < 0.05. [There are few genes with logFC 1.9, 1.89 also]

I have also tried with lfc=log2(2).

tr <- glmTreat(fit, contrast=contrast.matrix, lfc=log2(2))
tab2 <- topTags(tr,n=Inf, adjust.method = "BH")
keep <- tab2$table$FDR <= 0.05

To my surprise I see DE genes with logFC >=1, FDR < 0.05 and also logFC <=1 at FDR < 0.05.

Here is an example of DE genes I also found when I used lfc = log2(2)

GeneSymbols logFC   unshrunk.logFC  logCPM  PValue  FDR
DPYD-IT1    1.14493675  1.560389373 0.67468395  0.002972841 0.018704644
LINC00945   1.144288473 1.525191847 0.693845996 0.00651175  0.038064312
ITPKB-AS1   1.14046752  1.484991704 0.703887505 0.006620218 0.038585862
AC105206.1  1.122495266 1.4871905   0.699167758 0.006545299 0.038194633
LINC02066   1.120660287 1.558458881 0.66016102  0.004344743 0.026424298
AC090985.1  1.10235045  1.58968763  0.64343419  0.002745089 0.017426607
KCNJ6-AS1   1.100399688 1.717028637 0.613324799 0.002224039 0.01438451
LINC01494   1.090270138 1.490899364 0.676937984 0.006770608 0.039348027
AC099541.1  1.080957092 1.53817727  0.651356989 0.008449298 0.048356352
AL133396.2  1.012984207 1.512530231 0.63130643  0.008090259 0.046411928
AC011363.1  1.01207753  1.499708945 0.635740622 0.007226791 0.041797256

I'm very confused with these foldchange cutoffs.

Which one should I use to select differentially expressed genes based on fold change >=2 and FDR < 0.05?

Answer 1

Fold-change >= 2 is the same as logFC (log2(fold-change)) >= 1, so your example is doing exactly what you want. logFC is generally easier to think in and work with than fold-change, since (aside from computational reasons) then increases and decreases in expression differ only in sign and are, thus, easier to compare in magnitude (e.g., it's easier to tell that logFCs of 1.5 and -1.5 are the same magnitude change than it is that 3 and 0.333 are the same magnitude change in the fold-change).

For prioritizing the results, typically you would then rank these filtered results by fold-change and relevance (either given background knowledge or for genes with known function).

Answer 2

The edgeR authors recommend that you use a relatively low logFC threshold for glmTreat such as lfc=log2(1.2). A lfc value as high as lfc=log2(2) is seldom required and only would only be appropriate for datasets with really large amounts of DE. The first threshold you used of lfc=2 is equivalent to lfc=log2(4) and is too high for any dataset.

There's a few subtleties here:

1. First, there's the difference between logFC and FC as explained by Devon.

2. The lfc in glmTreat is not a fold-change cutoff! It is instead a testing threshold. All the significant genes will have actual logFCs well above the threshold.

3. The lfc threshold applies to the raw fold-changes (unshrunk.logFC in the output), not to the shrunk versions (logFC in the table).

4. The intent of glmTreat is to give a better gene ranking than would be provided by FDR or fold-change, so please don't re-sort the results by fold-change -- that would tend to defeat the purpose of the method.

I understand the documentation for glmTreat is a bit light, but you can get an idea of how to use it from our workflows.

• https://bioconductor.org/packages/release/workflows/vignettes/RnaSeqGeneEdgeRQL/inst/doc/edgeRQL.html

• https://bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow.html

I understand the simple fold-change cutoff of FC > 2 with FDR < 0.05 is often used in the literature, but the edgeR authors do not recommend it and hence we don't provide an automatic way to get it from the edgeR pipeline. For RNA-seq, fold-change cutoffs tend to prioritize genes with low expression and low counts. The glmTreat approach is more subtle but we think gives better results, both in terms of biological relevance and statistical correctness.