Question

Interpretation of LRT results

0

Entering edit mode

Laia ▴ 10

@239caad3

Last seen 13 months ago

Belgium

Dear Michael,

I am having some doubts when trying to extract relevant info from my gene lists generated by LRT in DEseq2.

Now I have gene lists of differentially expressed genes explained by:

the interaction of the 2 conditions that were applied to my cells,
condition 1
and condition 2

I have been trying to explore this data, by splitting the lists into Up and Down-reg gene lists, and entering these signif genes to GO in R, Panther and GOrilla, for example. Nevertheless, I have several levels per each condition, and I am not convinced that I am interpreting the results properly:

The reference condition in LRT is not the reference condition anymore, right? I am comparing the full design with the reduced design, not experimental conditions vs a reference condition. To do that I should do Wald, right? And from Wald, extract Up/Downs. But then I need to run many, many contrasts (more than 0).. and I don't feel this is how I should be doing this.

Is it correct to explore and visualize UP/DOWNs from the LRT gene lists? Or does it makes more sense to just say "these genes are affected by the interaction of condition1 and condition2 (LRT p.adj)", but then I have to run Wald to explain when are these genes up/down regulated?

I was also trying to plot my results in a heatmap where I could see the top diff. expressed genes for all the different conditions, but then is where I realize that the comparison is not against a reference anymore, but against the full model, so a heatmap would make no sense in this case. May I ask what would be the best way to visualize this result?

Thank you for your time, once again.

Laia

DESeq2 LRT Heatmap • 1.4k views

ADD COMMENT • link updated 13 months ago by James W. MacDonald 67k • written 2.1 years ago by Laia ▴ 10

score 0 · Answer 1 · 2022-11-07

0

Entering edit mode

James W. MacDonald 67k

@james-w-macdonald-5106

Last seen 22 minutes ago

United States

If you want information about individual contrasts, the only way to do that is to make the individual contrasts you want information about. Using the LRT test is an omnibus test that just tells you that one model fits better than a reduced version of the model without providing information about which particular comparisons are different.

ADD COMMENT • link 2.0 years ago James W. MacDonald 67k

0

Entering edit mode

Hello James,

In my case, I have 3 timepoints that I want to compare timepoint 1 vs. timepoint 0 and timepoint 2 vs. timepoint 0.

I tried to extract the results from res_lrt and res_wald, and they are different.

dds <- DESeqDataSetFromMatrix(cluster3_counts,
                              colData = metadata_3,
                              design = ~ orig.ident + timepoint)

Run DESeq2 with LRT

dds_lrt <- DESeq(dds, test = "LRT", reduced = ~ TissueArea)
res_lrt <- results(dds_lrt)
res_lrt_AvsB <- results(dds_lrt, c("timepoint", "2", "0"))

Run DEseq2 with Wald Test

dds <- DESeq(dds)
res <- results(dds)
res_AvsB <- results(dds_lrt, c("timepoint", "2", "0"))

You mentioned LRT test is an omnibus test, but I still can extract the results for one pair of comparison. I do not know which one is correct.

ADD REPLY • link 13 months ago Yijing • 0

1

Entering edit mode

You need to pay attention to the output.

> dds <- makeExampleDESeqDataSet()

> colData(dds)$condition <- factor(rep(LETTERS[1:3], each = 4))
> dds_lrt <- DESeq(dds, "LRT", reduced = ~1)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
> dds_regular <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

> results(dds_lrt, c("condition","A","B"))
log2 fold change (MLE): condition A vs B     <-- what does that line say?
LRT p-value: '~ condition' vs '~ 1'          <--- and what does that line say?

## Wald
> dds_regular <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

## Now compare.
> all.equal(results(dds_lrt,  c("condition","A","C"))$pvalue, results(dds_lrt,  c("condition","A","B"))$pvalue)
[1] TRUE
> all.equal(results(dds_regular,  c("condition","A","C"))$pvalue, results(dds_regular,  c("condition","A","B"))$pvalue)
[1] "Mean relative difference: 0.6234569"

As I noted before, the LRT is an omnibus test where the p-values test that the model fits better with the term versus without. You can still extract logFC values for any comparison, but that is unrelated to the test. The Wald test is a direct test between two groups, and will differ depending on the groups you compare.

ADD REPLY • link 13 months ago James W. MacDonald 67k