interpreting LRT results DESeq2
1
3
Entering edit mode
@solgakarbitechnionacil-6453
Last seen 4.5 years ago
European Union

Hello,

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)

locale:
[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,

Olga.

deseq2 LRT RNA-Seq • 4.5k views
ADD COMMENT
1
Entering edit mode
@mikelove
Last seen 14 hours ago
United States

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 COMMENT
0
Entering edit mode

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.

Login before adding your answer.

Traffic: 226 users visited in the last hour
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