How to extract all comparison in DESeq2 LRT test
1
0
Entering edit mode
Yijing • 0
@0fcbf96d
Last seen 20 months ago
Germany

Hi,

I want to compare the differences among 3 genotypes, namely Hom, Het and WT.

Since there are three groups, I use test = LRT.

Code as following:

# Construct dds 
dds <- DESeqDataSetFromMatrix(countData = cts_br_filter,
                          colData = metadata_br_filter,
                          design = ~ TissueArea + Genotype)

# Run DESeq2 
dds <- DESeq(dds, test = "LRT", reduced = ~ TissueArea)

When I used resultsNames(dds), I do not understand the output: why there is no "Genotype_Hom_vs_Het"? For the TissueArea, what do those mean?

> resultsNames(dds) 
[1] "Intercept"                                  "TissueArea_corpus.callosum_vs_cerebellum"  
[3] "TissueArea_hippocampus_vs_cerebellum"       "TissueArea_prefrontal.cortex_vs_cerebellum"
[5] "Genotype_Het_vs_WT"                         "Genotype_Hom_vs_WT"

When I want to extract the results for "Hom" vs "Het", an error occurs.

res_HomHet <- results(dds, contrast = c("Genotype", "Hom", "Het"))
res_HomHet_shrunken <- lfcShrink(dds, coef = "Genotype_Hom_vs_Het", type = "apeglm")

> res_HomHet_shrunken <- lfcShrink(dds, coef = "Genotype_Hom_vs_Het", type = "apeglm")
    Error in lfcShrink(dds, coef = "Genotype_Hom_vs_Het", type = "apeglm") : 
      coef %in% resultsNamesDDS is not TRUE

If I want to get the results for "Hom" vs "Het", what should I do?

Look forward to your reply!
DESeq2 Likelihood-ratiotest • 2.0k views
ADD COMMENT •
0
Entering edit mode
@james-w-macdonald-5106
Last seen just now
United States

By default you will set one of the genotypes to be the baseline, and the other two coefficients will be comparisons to the baseline. In your case the baseline is WT, and you compare either Het or Hom to that baseline. If you want to compare Het vs Hom, you would conventionally add contrast = c("Genotype", "Het","Hom") as you have done above. But you can't use the contrast argument with lfcShrink

> fake0 <- makeExampleDESeqDataSet()
> colData(fake0) <- DataFrame(condition = factor(rep(c("Control","Treat1","Treat2"), each = 4)))
> dds <- DESeq(fake0, "LRT", reduced = ~1)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
> resultsNames(dds)
[1] "Intercept"                   "condition_Treat1_vs_Control"
[3] "condition_Treat2_vs_Control"

> lfcShrink(dds, contrast = c("condition","Treat1","Treat2"))
Error in lfcShrink(dds, contrast = c("condition", "Treat1", "Treat2")) : 
  type='apeglm' shrinkage only for use with 'coef'

## urf. Change contrast
> colData(fake0)$condition <- relevel(colData(fake0)$condition, "Treat1")
> dds <- DESeq(fake0, "LRT", reduced = ~1)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

> resultsNames(dds)
[1] "Intercept"                   "condition_Control_vs_Treat1"
[3] "condition_Treat2_vs_Treat1" 

> lfcShrink(dds, "condition_Treat2_vs_Treat1")
using 'apeglm' for LFC shrinkage. If used in published research, please cite:
    Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
    sequence count data: removing the noise and preserving large differences.
    Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
log2 fold change (MAP): condition Treat2 vs Treat1 
LRT p-value: '~ condition' vs '~ 1' 
DataFrame with 1000 rows and 5 columns
          baseMean log2FoldChange     lfcSE    pvalue      padj
         <numeric>      <numeric> <numeric> <numeric> <numeric>
gene1     61.93415      0.0302668  0.124168  0.722282  0.993072
gene2      4.44099     -0.0163095  0.127315  0.261012  0.979982
gene3     11.18672      0.0197914  0.127609  0.549889  0.993072
gene4      7.39364      0.0279736  0.130055  0.209022  0.979982
gene5    304.41708      0.0131412  0.120535  0.578371  0.993072
...            ...            ...       ...       ...       ...
gene996    2.05200    -0.00407019  0.126433 0.3975133  0.993072
gene997    6.84232    -0.00416264  0.125284 0.4163598  0.993072
gene998    3.17878    -0.02259350  0.129013 0.2057943  0.979982
gene999  162.11266    -0.02686071  0.122492 0.7490227  0.993072
gene1000  11.14948     0.02935432  0.129245 0.0865619  0.976984
ADD COMMENT •
0
Entering edit mode

Many thanks to your reply!

I have one more question:

I can use res_HomHet <- results(dds, contrast = c("Genotype", "Hom", "Het")) to get the result, but just cannot shrink. Is this result the real Hom vs Het? Can I use this unshrunk value for further analysis like extracting significant DEGs?

Or do I need to relevel just as you did colData(fake0)$condition <- relevel(colData(fake0)$condition, "Treat1") to get the Treat2 vs Treat1?

ADD REPLY •
0
Entering edit mode

Note that the lfc is not part of the LRT testing. LRT os a goodness of fit test. I would use Wald here, and reserve LRT for something like timecourse experiments with many TPs. You can use contrasts to shrink if using the ashr method.

ADD REPLY •
0
Entering edit mode

Thank you for your reply!

But according to DESeq2 tutorial, The LRT is therefore useful for testing multiple terms at once, for example testing 3 or more levels of a factor at once, or all interactions between two variables. I thought the Wald test is just like t-test can only used for two groups comparison. You think it is better to use Wald here?

ADD REPLY •
0
Entering edit mode

Yes, but you only have three groups, wouldn't it (interpret-wise) be much more intuituve to do geno1-geno2, geno2-geno3 and geno1-geno3? The Wald test tests fold changes so you have quantitiative information. It's on you of course but LRT in my mind is when you have for example a time course across a dozen TPs and parwise would give excessively much contrasts.

ADD REPLY •
0
Entering edit mode

Thank you for your reply!

I am not sure if it is the same as the t-test and ANOVA. Because for more than 3 groups, if using a t-test, the p-value is not correct.

ADD REPLY •

Login before adding your answer.

Similar Posts
Loading Similar Posts
Traffic: 954 users visited in the last hour
Content Search
Users
Tags
Badges
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6