Question: why my resultsNames(deseq) is not the same as model.matrix ?
gravatar for unique379
3.4 years ago by
unique3790 wrote:

Hello michael love,

could you please answer from my cross post

anyway i am repeating here:

I tried with full model to remove my batch from sva package then ran deseq2 wald and LRT test to see the differences but got very surprising results for the same contrast. How it possible that, Wald Test gave me tags 0 at padj <0.05 while LRT test identified 345 tags at the same level of cut-off. This is really i dont understand why it happened, could you please explain me why or where i am making mistake ?? Here is my code:

dds <- DESeqDataSetFromMatrix(countData= filtered, colData=DataFrame(group), design = ~group)
dds <- estimateSizeFactors(dds)
dat <- counts(dds, normalized=TRUE)
idx <- rowMeans(dat) > 1
dat <- dat[idx,]
## log2 transform in order to use sva instead svaseq
dat <- log2(dat)
## full model
mod <- model.matrix(~group, colData(dds))
## taking first group as intersept
mod0 <- model.matrix(~1, colData(dds))
sva.dat <- sva(dat, mod, mod0,
## batch corrected dds object
ddssva <- dds
ddssva$SV1 <- sva.dat$sv[,1]
ddssva$SV2 <- sva.dat$sv[,2]
ddssva$SV3 <- sva.dat$sv[,3]
design(ddssva) <- as.formula(~ SV1 + SV2 + SV3 + group)

## Using Wald test (default)
deseq.wald <- DESeq(object = ddssva)
## get contrast
**1] "Intercept"     "SV1"           "SV2"           "SV3"          
[5] "groupA"   "groupB"   "groupC" "groupD"**  
res.wald <- results(deseq.wald,contrast = c("group", "B", "A"))
## order
res.wald <-[order([6]),]
## how many significant tags?
sum(res.wald$padj < 0.05, na.rm=TRUE)

**## Using LRT test**
deseq.lrt <- DESeq(object = ddssva, test = "LRT", reduced = ~1)
## get contrast
**[1] "Intercept"                "SV1"                     
[3] "SV2"                      "SV3"                     
[5] "group_B_vs_A"   "group_C_vs_A"
[7] "group_D_vs_A"**  
res.lrt <- results(deseq.lrt,contrast = c("group", "B", "A"))
## order
res.lrt <-[order([6]),]
## how many significant tags?
sum(res.lrt$padj < 0.05, na.rm=TRUE)


I understand that there is problem in contrast of resultsNames(deseq.wald) and resultsNames(deseq.lrt), but why ?? using two different tests made different different contrasts ??

Thank once again for your kind reply.

rnaseq sva deseq2 • 609 views
ADD COMMENTlink modified 3.4 years ago by Michael Love26k • written 3.4 years ago by unique3790
Answer: why my resultsNames(deseq) is not the same as model.matrix ?
gravatar for Michael Love
3.4 years ago by
Michael Love26k
United States
Michael Love26k wrote:

These are not the same statistical tests. 

The Wald test you have specified is testing for genes with differences in B vs A controlling for differences due to all the other variables in the design.

The LRT you have specified is, for which genes are any of the terms in the model (SV 1 - 3 and all levels of group) significant for explaining differences in gene expression. 

The use of contrast in the last results() call is just changing the log2FoldChange column, not the p-values.

See ?results about how the LRT p-values are calculated.

If you want to specify a Wald contrast after fitting an LRT model, you need to set test="Wald" in your call results(). 

ADD COMMENTlink modified 3.4 years ago • written 3.4 years ago by Michael Love26k
Please log in to add an answer.


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