Question

Limma vs deseq2 results

1

Entering edit mode

k___ ▴ 30

@ed926d96

Last seen 15 hours ago

United Kingdom

Hi there,

I am trying to use limma for some rna seq data and model subjects as a random effect as I have paired samples. I have already done analysis in DESeq2 and got 400 (up and down) significant DGEs. When repeating in limma I only get 19 DGEs. I'm not sure where I am going wrong in limma - my code is below

counts <- DGE(counts)
condition <- factor(designTable$condition, levels=c("normal","tumour"))
design <- model.matrix(~ condition)
voom <- voom(counts, design)
dup <- duplicateCorrelation(voom, design, block=designTable$subjects)
voom <- voom(counts,design, block=designTable$subjects, correlation = dup$consensus)
dup <- duplicateCorrelation(voom, design, block=designTable$subjects)
fit <- lmFit(voom,design, block=designTable$subjects, correlation = dup$consensus)
fit <- eBayes(fit,robust =T) 
summary(decideTests(fit))

       (Intercept) conditiontumour
Down          3439                 4
NotSig        3161             18940
Up           12359                15

I don't understand how they are so different, any help would be great

limma RNASeq R • 4.8k views

ADD COMMENT • link 23 months ago k___ ▴ 30

score 6 · Answer 1 · 2022-05-30

6

Entering edit mode

Gordon Smyth 50k

@gordon-smyth

Last seen 50 minutes ago

WEHI, Melbourne, Australia

There is no guarantee that different packages will give the same results, but here is it clear that you have done different analyses in the two packages since DESeq2 has no equivalent to duplicateCorrelation.

If you have properly paired samples, then you should be using design <- model.matrix(~subjects + condition) and not using duplicateCorrelation. See Section 9.4.1 of the limma User's Guide.

You also need to filter unexpressed genes using filterByExpr before running voom.

There is no Bioconductor function called DGE. I am not sure whether the first line of code in your question is a typo or you have used a private function.

ADD COMMENT • link 23 months ago Gordon Smyth 50k

0

Entering edit mode

Hi, yes sorry I used the function DGEList(). I changed it to how you suggested but I am still confused as limma's results show no DEGs for tumour vs normal, am I going wrong somewhere? My code is below

> designTable
          subjects tissue
file1      A       tumour
file2      A       normal
file3      B       tumour
file4      B       normal
file5      C       tumour
file6      C       normal

#deseq2

summary(res)
out of 18991 with nonzero total read count
adjusted p-value < 0.05
LFC > 0 (up)       : 157, 0.83%
LFC < 0 (down)     : 284, 1.5%

#limma
dge <- DGEList(counts=raw_counts)
dge
design <- model.matrix(~ subjects+ tissue, designTable)
design

summary(decideTests(fit))
      (Intercept) patientB patientC tissuetumour
Down            16         0         0            0
NotSig        2244     13903     13903        13903
Up           11643         0         0            0

ADD REPLY • link 23 months ago k___ ▴ 30

2

Entering edit mode

limma generally gives stricter FDR control than does DESeq2 because the DESeq test doesn't take into account uncertainty in estimating the dispersion parameter. So it is perfectly possible to get a moderate number of DE genes from DESeq2 and none from limma. That doesn't indicate any error.

It is very surprising however that you also have no baseline differences between the subjects either. If this is human data, then I have never seen a dataset with no baseline differences between human subjects. This would suggest there is something unusual about your dataset.

There is no rule that says you need to use FDR < 0.05. You could obviously use FDR < 0.1 if you wanted to see more results.

Generally speaking, if you don't get DE results when you think you should, the next step is to troubleshoot the quality of your data with diagnostic plots. Changing the statistical test in hope of getting a different result isn't usually the right approach.

ADD REPLY • link 23 months ago Gordon Smyth 50k