Dear Community,

i would like to ask you an important question about the interpretation of the topTable function output results. Specifically, i know that this is not a general statistics blog, but i checked the argument confint and i used after a specific implementation with limma, i used confint=0.95 in order to return confidence intervals for logFCs.

In detail, here is a small output of some selected genes (after i have subsetted my topTable):

> head(significant, 20)
             GENE_SYMBOL     logFC   adj.P.Val      MAD      CI.L       CI.R
A_23_P114903       HSPA6  3.595423 0.048103409 3.691256  1.729956  5.4608897
A_24_P245379    SERPINB2  2.910139 0.027037437 2.955397  1.676020  4.1442576
A_23_P161698        MMP3  2.581726 0.022127328 2.692339  1.569174  3.5942793
A_23_P66241         MT1M  2.857147 0.030223818 2.517837  1.601280  4.1130140
A_23_P206724        MT1E  2.574222 0.023364607 2.464120  1.535268  3.6131756
A_32_P87013        CXCL8  3.531625 0.015656484 2.316576  2.304604  4.7586457
A_24_P125096        MT1X  2.528568 0.017441204 2.314788  1.622389  3.4347482
A_23_P37983         MT1B  2.421996 0.018654318 2.282815  1.531431  3.3125618
A_23_P206707        MT1G  2.395756 0.025898285 2.259628  1.395022  3.3964891
A_23_P71037          IL6  2.135955 0.037906167 2.179081  1.119542  3.1523685
A_23_P427703        MT1L  2.364336 0.017614832 2.136130  1.514988  3.2136838
A_23_P163782      MT1HL1  2.284161 0.020399185 2.112673  1.414799  3.1535226
A_23_P315364       CXCL2  2.426432 0.006127387 1.995657  1.900835  2.9520286
A_23_P414343        MT1H  2.380918 0.014319566 1.995214  1.598750  3.1630865
A_23_P365738         ARC  2.102056 0.030560036 1.993690  1.170771  3.0333421
A_23_P1691          MMP1  2.154745 0.020770068 1.982411  1.329721  2.9797693
A_23_P108842       DUSP2  1.972652 0.014964790 1.973428  1.306712  2.6385916
A_23_P54840         MT1A  1.974725 0.024914831 1.839428  1.158456  2.7909943
A_23_P15727       FKBP10 -1.908220 0.040797553 1.724086 -2.840906 -0.9755334
A_24_P251764       CXCL3  1.938549 0.007849170 1.667274  1.424668  2.4524294

Thus, how i can interpret and "evaluate" the returned confidence intervals about a specific gene with a specific logFC ? That for instance, for the first gene, HSPA6 which has a "relatively" big fold change, is more "significant" due to the fact that both CI.L & CI.R >1 ? Or even the case that one of these is >1 ? as here this specific gene is upregulated ? Or my approach to this matter is completely wrong ? For instance, if a gene above with a significant p-value-i.e FDR < 0.05--and a logFC of -0.5, had CI.R=-0.4 & CI.L=-0.8, which is the "evaluation" of this example ?

Thank you,

Konstantinos

The confidence interval gives you a measure of how precisely the log-fold change is estimated. The tighter the interval, the more precise the estimate. This can be useful if the actual value of the log-fold change estimate is going to be used for further work. For example, I often use the confidence intervals to check if near-zero log-fold changes are being precisely estimated; this allows me to identify genes that are likely to be non-DE, which is hard to do with conventional significance testing. Confidence intervals may also be helpful for visualizing this estimation uncertainty in plots. For your results, the CIs are basically saying that the true log-fold change is probably quite large, given that the interval is quite a fair distance from zero.

However, if you want to interpret significance of DE, you're better off looking at the p-value. The p-value calculation already accounts for the estimation precision by using the standard error of the log-fold change in the moderated t-test. The p-values can also be adjusted for multiple testing, whereas that's not done (as easily) for the CIs. Of course, all of these things are interrelated. If you have tight CIs, you're more likely to get a significant result. If you have a large log-fold change, then having a wider CI doesn't matter so much as the evidence against the null is still strong.

For your example, the interpretation would be something like: "The gene is significant because it has an adjusted p-value below the 5% threshold. Also, the log-fold change estimate of -0.5 is reasonably precise as the confidence interval is fairly small."

It's a little unfortunate that you've deleted the P.Value column from the topTable, because the CI relates directly to the unadjusted p-value.

If the p-value is > 0.05, then the CI will overlap 0, i.e., CI.L<0 and CI.R>0. (Note that logFC=0 corresponds to no DE.)

If the p-value is < 0.05, then the CI will be entirely above zero or entirely below zero.

Your gene HSPA6 has a large logFC but also a high variability. If you were focusing on this gene as being of a priori interest, then the 95% confidence interval for the log2-fold-change is (1.73,  5.46). Another way to say the same thing is that the log2-fold change is 3.60 +- 1.86.

If I unlog the logFC and CI values for HSPA6, then the estimated fold change is 2^3.60 = 12.1, and the range of feasible values for the fold change is from as low as 2^1.73 = 3.32 to as high as 2^5.46 = 44.0.

PS. You seem to be ranking genes by fold change in the table in your post. The limma documents recommends against that, see help("topTable"). If you want to rank genes by importance. it is usually better to stick to p-value, possibly after using treat.