Question: One-way hypothesis test in limma
0
gravatar for s1437643
18 months ago by
s14376430
s14376430 wrote:

Is it possible to test for fold-change in only one direction (LFC > 0) using the treat function from limma? Similar to using the altHypothesis argument in the results function from DESeq2?

I tried myself to introduce some changes to the code in treat to perform this:

lfc <- 0
se <- stdev.unscaled*sqrt(fit$s2.post)
tstat <- (coefficients-lfc)/se
fit$t <- array(0,dim(coefficients),dimnames=dimnames(coefficients))
fit$p.value <- pt(tstat, df=df.total, lower.tail=FALSE)
tstat <- pmax(tstat,0)
fit$t <- tstat
fit$treat.lfc <- lfc

This seems correct, but wanted to double-check incase a method was already implemented.

limma • 442 views
ADD COMMENTlink modified 18 months ago by Gordon Smyth38k • written 18 months ago by s14376430
Answer: One-way hypothesis test in limma
2
gravatar for Aaron Lun
18 months ago by
Aaron Lun25k
Cambridge, United Kingdom
Aaron Lun25k wrote:

That seems fine to me.

Keep in mind that the usual two-sided test in treat() considers the null hypothesis where the true log-fold changes are distributed at -lfc and lfc with 50% probability each. Your one-sided test considers a different null hypothesis where the true log-fold change is lfc alone. This is probably the sensible thing to do, but it means that the nulls are not exactly comparable, which is something to be aware of.

Also, I would avoid hacking into treat() like you've done here. If you know enough to compute your own test statistics, then you know enough to form your own output data.frame from fit. Perhaps there is an argument for adding a one-sided test option to limma functions, but this tends to be a very bespoke setting - I have only ever done this on occasion for ChIP-seq analysis, and never for RNA-seq.

ADD COMMENTlink modified 18 months ago • written 18 months ago by Aaron Lun25k

Many thanks, Aaron.

ADD REPLYlink written 18 months ago by s14376430
Answer: One-way hypothesis test in limma
0
gravatar for Gordon Smyth
18 months ago by
Gordon Smyth38k
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Gordon Smyth38k wrote:

It would be easy to implement a one-sided test in treat(), pretty much as you've done. The two-sided test was much more difficult to derive! When the threshold lfc > 0, the two-sided p-value can't be derived from the one-sided p-values or vice-versa.

I haven't implemented one-sided tests in treat() because I have yet to see a real experiment for which that would be scientifically the right thing to do. I don't want to give users rope to hang themselves with.

ADD COMMENTlink modified 18 months ago • written 18 months ago by Gordon Smyth38k

So if one wanted to implement one-sided tests (LFC > LT and LFC < -LT seperately), with an lfc threshold of LT, would the code look like this?

# testing LFC > 0 (condition1 > condition2 for the given contrast, same code as the original post)
lfc <- log2(2)#arbitrary lfc threshold
se <- stdev.unscaled*sqrt(fit$s2.post)
tstat <- (coefficients-lfc)/se          
fit$t <- array(0,dim(coefficients),dimnames=dimnames(coefficients))
fit$p.value <- pt(tstat, df=df.total, lower.tail=FALSE)
tstat <- pmax(tstat,0)
fit$t <- tstat
fit$treat.lfc <- lfc
# testing LFC < 0 (condition1 < condition2 for the given contrast)
lfc <- log2(2) #arbitrary lfc threshold
se <- stdev.unscaled*sqrt(fit$s2.post)
tstat <- ((-coefficients)-lfc)/se          
fit$t <- array(0,dim(coefficients),dimnames=dimnames(coefficients))
fit$p.value <- pt(tstat, df=df.total, lower.tail=FALSE)
tstat <- pmax(tstat,0)
fit$t <- tstat
fit$treat.lfc <- lfc

P.S. I am interested in the one-sided p-values because I will use them as input for GO enrichment, so I want to know when GO terms are up vs. down regulated for each comparison

ADD REPLYlink written 8 weeks ago by sidreed340
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 319 users visited in the last hour