Question: interpreting LRT results DESeq2
gravatar for
4.2 years ago by
European Union


I am working on RNA-Seq data which consists of 10 samples. My model has two factor: Time and Genotype, having 3 time point: Tbase, T45, Tend and 2 genotypes: WT, mutant (for Tend and Tbase each genotype has 2 replicates, for T45 each genotype has a single replicate).

I used the LRT test of DESeq2 v.1.6.3 since I am interested in the ratio of ratios results:

- (mutant_Tend /mutant_Tbase) / (WT_Tend /WT_Tbase)

- (mutant_T45 /mutant_Tbase) / (WT_T45 /WT_Tbase)

The commands that I used are:

> dds = DESeqDataSetFromMatrix(countData, colData, design = ~ Time + Genotype + Genotype:Time)

> dds = DESeq( dds, test = "LRT", reduced = ~ Time + Genotype)

> resultsNames(dds_LRT)

[1] "Intercept"               "Time_T45_vs_Tbase"       "Time_Tend_vs_Tbase"      "Genotype_mutant_vs_WT"   "TimeT45.Genotypemutant" 
[6] "TimeTend.Genotypemutant"

> res_Tend_vs_Tbase = results(dds_LRT, name = "TimeTend.Genotypemutant")

> res_T45_vs_Tbase = results(dds_LRT, name = "TimeT45.Genotypemutant") 

I have 3 questions:

1. Is this the correct way to receive the comparisons that I am interested in?

2. If so, I am also interested in comparing the two time points: Tend vs T45, what is the correct way to call results function for that comparison?

3. For both results: res_Tend_vs_Tbase, res_T45_vs_Tbase, I received different log2FoldChange values but the pvalue and padj columns are identical in both comparisons. Is there an explanation to that?

In case the sessionInfo() output is neede:

> sessionInfo()
R version 3.1.1 (2014-07-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)

[1] LC_COLLATE=Hebrew_Israel.1255  LC_CTYPE=Hebrew_Israel.1255    LC_MONETARY=Hebrew_Israel.1255 LC_NUMERIC=C                  
[5] LC_TIME=Hebrew_Israel.1255    

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] DESeq2_1.6.3              RcppArmadillo_0.4.650.1.1 Rcpp_0.11.4               GenomicRanges_1.18.4      GenomeInfoDb_1.2.4       
[6] IRanges_2.0.1             S4Vectors_0.4.0           BiocGenerics_0.12.1      

loaded via a namespace (and not attached):
 [1] acepack_1.3-3.3      annotate_1.44.0      AnnotationDbi_1.28.1 base64enc_0.1-2      BatchJobs_1.5        BBmisc_1.9          
 [7] Biobase_2.26.0       BiocParallel_1.0.3   brew_1.0-6           checkmate_1.5.1      cluster_2.0.1        codetools_0.2-10    
[13] colorspace_1.2-4     DBI_0.3.1            digest_0.6.8         fail_1.2             foreach_1.4.2        foreign_0.8-63      
[19] Formula_1.2-0        genefilter_1.48.1    geneplotter_1.44.0   ggplot2_1.0.0        grid_3.1.1           gtable_0.1.2        
[25] Hmisc_3.15-0         iterators_1.0.7      lattice_0.20-30      latticeExtra_0.6-26  locfit_1.5-9.1       MASS_7.3-39         
[31] munsell_0.4.2        nnet_7.3-9           plyr_1.8.1           proto_0.3-10         RColorBrewer_1.1-2   reshape2_1.4.1      
[37] rpart_4.1-9          RSQLite_1.0.0        scales_0.2.4         sendmailR_1.2-1      splines_3.1.1        stringr_0.6.2       
[43] survival_2.38-1      tools_3.1.1          XML_3.98-1.1         xtable_1.7-4         XVector_0.6.0

Thank you very much for the help,


deseq2 rna-seq lrt • 3.5k views
ADD COMMENTlink modified 4.2 years ago by Michael Love24k • written 4.2 years ago by
Answer: interpreting LRT results DESeq2
gravatar for Michael Love
4.2 years ago by
Michael Love24k
United States
Michael Love24k wrote:

regarding why the LRT p-value doesn't change when you change 'name',

from ?results:

"For analyses using the likelihood ratio test (using nbinomLRT), the p-values are determined solely by the difference in deviance between the full and reduced model formula. A log2 fold change is included, which can be controlled using the name argument."

DESeq2 offers tests for specific terms using the Wald test. You don't have to rerun the model though. In order to calculate a p-value for a specific interaction term, you can just add test="Wald" to the results() call.

Yes those are the correct interaction terms for your two ratios of ratios.

To contrast the two time points T45 and Tend you can contrast the two interaction terms, using the list() style of contrasts:

results(dds, contrast=list("TimeTend.Genotypemutant", "TimeT45.Genotypemutant"))


ADD COMMENTlink written 4.2 years ago by Michael Love24k

As always, thank you very much for your quick reply!

I am sorry to bother you with questions that are already answered in the help section, I don't know how I've missed it.

ADD REPLYlink written 4.2 years ago by
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 163 users visited in the last hour