Question

EdgeR glmLRT vs glmQLFTest

0

Entering edit mode

grastalt27 • 0

@grastalt27-10859

Last seen 7.8 years ago

Hi Everyone,

I'm doing some DGE analysis with 4 groups! I have tried to use the glmQLFTest as suggested in the edgeR users manual, but there are barely any differently expressed genes with FDR < 0.99. On the other hand for the same comparisons, glmLRT gives me a few hundred genes.

Is glmQLFTest not valid for multiple contrasts?

If I move forward with my analysis using the glmLRT function, will my results be valid?

Thank you all!

bioconductor edger • 11k views

ADD COMMENT • link updated 7.8 years ago by Steve Lianoglou ★ 13k • written 7.8 years ago by grastalt27 • 0

score 3 · Answer 1 · 2016-06-27

3

Entering edit mode

Steve Lianoglou ★ 13k

@steve-lianoglou-2771

Last seen 14 months ago

United States

What do you mean when you ask if glmQLFTest is "valid for multiple contrasts"?

First thing I'd do would double check that the code is valid in both scenarios -- perhaps you can share a clean up version of your code for us to help diagnose?

Assuming the code is correct, I would strongly urge you to plot the expression values (log2(cpm)) of a subset of the genes being called by glmLRT and not QLF across your comparisons to give you confidence in these calls.

In the general case, I really, really, really, really, like the output from the edgeR's QLF "framework" so I'd be extra careful when working in the corner you seem to be painted in.

Aaron Lun recently outlined a few scenarios where the QLF framework would not be a good choice, so take a minute to read this:

A: EDGE-R exact test vs QL F-test

And while you're off reading, this recent f1000 paper by the edgeR/QLF folks is good reading:

From reads to genes to pathways: differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline

ADD COMMENT • link 7.8 years ago Steve Lianoglou ★ 13k

2

Entering edit mode

I'll second the vote in favor of the QL method over the LRT method. In every dataset I've seen where I did a null comparison between replicates in which no genes should be differentially expressed, glmQLFTest never gave me more than a few significant genes, while glmLRT sometimes gave me many genes. So I trust glmQLFTest much more to control the false positive rate.

ADD REPLY • link 7.8 years ago Ryan C. Thompson ★ 7.9k

0

Entering edit mode

Thanks Steve,

I was wondering if glmQLFTest only works for one comparison like a case v control. I have a case1 v case2, case1 v case 3, case 1 v control, ...

Here is the code:

y <- DGEList(dat, group=group)

batch <- c(0,0,0,0,0,0,1,1,1,1,1,1)
fgroup <- factor(group)
design <- model.matrix(~ 0 + fgroup + batch)
colnames(design) <- levels(c(fgroup, batch))

y <- calcNormFactors(y)
y <- estimateDisp(y, design, robust=TRUE)

# LRT method
fit <- glmFit(y, design, robust=TRUE)
lrt <- glmLRT(fit, contrast=c(-1,1,0,0,0)) # has smallest FDR = 2.15e-39

#QLFTest method
fit <- glmQLFit(y, design, robust=TRUE)
qlf <- glmQLFTest(fit, contrast=c(-1,1,0,0,0)) # top 100 FDR values are exactly 0.9998381

ADD REPLY • link 7.8 years ago grastalt27 • 0

3

Entering edit mode

That's quite an extreme difference. One possible cause is that your genes have highly variable dispersions. If edgeR sees this in the data, it will (correctly) reduce the estimate of the prior d.f. in both estimateDisp and glmQLFit, thus reducing the amount of shrinkage that can be performed. In general, I would consider a prior.d.f below 5 to be small, while anything above 50 is pretty high, i.e., strong shrinkage. For small prior d.f., the power of glmQLFTest will (appropriately) decrease, as inferences are less reliable from weakly-shrunk dispersions that are more unstable. However, glmLRT is blind to this effect, potentially resulting in spuriously low p-values due to imprecisely estimated dispersions with small values.

It's worth having a look at the output of plotBCV and plotQLDisp as diagnostics. The former will reveal if there's some strong trend that hasn't been properly modelled and is subsequently inflating the variability of the dispersions - more aggressive filtering should take care of this, in cases where the strong trend is generated by technical noise at low counts. The QL dispersion plots will allow you to see if there's a mass of outlier genes, which could also inflate the variability of the dispersions (especially if it overwhelms the default winsorizing thresholds). For example, Ig chains are often highly variable if you're working with a data set containing different B-cell clones, so it's worth getting rid of them to avoid large outliers. Some types of batch effects can also increase the variability of the dispersions if they're left unmodelled, as the gene-specific batch effect gets incorporated into the dispersion estimate - check your MDS plot for any patterns in your replicates.

As for your actual question, there are no restrictions on the types of contrasts you can give to glmQLFTest; it'll handle any number of case vs. case/control comparisons. Further, on an unrelated note, robust doesn't do anything in glmFit.

ADD REPLY • link 7.8 years ago Aaron Lun ★ 28k

0

Entering edit mode

As Aaron says, this sounds extreme. If you want to to email it to me, I'd be interested to see test your data to see why the LRT and QL pipelines give such different FDR values for your data.

Gordon

ADD REPLY • link 7.8 years ago Gordon Smyth 50k