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

Hello michael love,

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, n.sv=3)
## 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
resultsNames(deseq.wald)
**1] "Intercept"     "SV1"           "SV2"           "SV3"
[5] "groupA"   "groupB"   "groupC" "groupD"**
res.wald <- results(deseq.wald,contrast = c("group", "B", "A"))
## order
res.wald <- as.data.frame(res.wald)[order(as.data.frame(res.wald)[6]),]
## how many significant tags?
sum(res.wald$padj < 0.05, na.rm=TRUE) **0** **## Using LRT test** deseq.lrt <- DESeq(object = ddssva, test = "LRT", reduced = ~1) ## get contrast resultsNames(deseq.lrt) **[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 <- as.data.frame(res.lrt)[order(as.data.frame(res.lrt)[6]),] ## how many significant tags? sum(res.lrt$padj < 0.05, na.rm=TRUE)
**345**

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 ??

rnaseq sva deseq2 • 609 views
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 ?
0
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().