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.

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"))