Hello michael love,

could you please answer from my cross post https://www.biostars.org/p/200157/

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

Thank once again for your kind reply.

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().