This is not an easy thing to do, as we cannot work with the log-fold changes directly as we would in limma
. Rather, we need to test against the upper and lower log-fold change thresholds of the TOST through the offsets.
Assume we have a DGEList
named y
for which we have estimated the dispersions, a design
matrix and a coefficient of interest coef
for our DE comparison. We could then test against the null hypothesis that the coefficient is set at the upper TOST threshold. In this case, we use a 3-fold change as the threshold.
orig.offs <- getOffset(y)
y$offset <- orig.offs + log(3)*design[,coef]
fit.up <- glmFit(y, design)
lrt.up <- glmLRT(fit.up, coef=coef)
We do the same for the lower TOST threshold.
y$offset <- orig.offs - log(3)*design[,coef]
fit.down <- glmFit(y, design)
lrt.down <- glmLRT(fit.down, coef=coef)
We then compute the one-sided p-values for each direction, and intersect them to get those with log-fold changes that lie between the upper and lower TOST boundaries. Those that are outside get a p-value of 1.
p.up <- ifelse(lrt.up$table$logFC < 0, lrt.up$table$PValue/2, 1)
p.down <- ifelse(lrt.down$table$logFC > 0, lrt.down$table$PValue/2, 1)
p.tost <- pmax(p.up, p.down)
So, as you can see, it's not a pleasant process, and it's pretty conservative. You'll probably need to use fairly wide TOST thresholds to get any significance. Obviously, if you're planning to apply this genome-wide, you'll have to do a BH correction on p.tost
.
Thanks Aaron for the awesome solution! And just to make sure: if I need to define contrasts, I would just need to edit the offset of one of the two mates, right?
I'm not entirely sure what you mean by the mates, but the problem becomes more complicated if you have a contrast vector instead of simply dropping a coefficient. I'd suggest reparameterizing your design matrix such that your DE comparison can be expressed by just dropping a coefficient. This should be easy enough to do for most comparisons. Otherwise, you'll have to do something like that in
glmLRT
, which involves some complicated stuff with QR decompositions to get the null design matrix.Cool! I'll change my design matrix as you suggest then. Thanks a lot for your suggestions!