Hello,

I am working on a RNA-Seq experiment in which there are 3 different cell types analyzed. For each cell type I have 2 biological replicates - each represent a different pool of mice.

I first ran DESeq2 without setting the 'Mice' factor, meaning without pairing the 3 cell types of the same mice pool.

The colData in this case is:

I have used the following design: ~CellType and obtained the 3 pairwise results between each pair of cell types.

An example for a results command I have used is: 

res = results(dds, contrast = c("CellType", "type1", "type2"), alpha = 0.05)

Then, I have created a data set with the 'Mice' factor, creating the following colData:

with the design: ~Mice + CellType   and the same commands for results.

As expected I receive different results, but for some of them I can't seem to understand the differences.

For example, in a certain gene I get the following expression vector:

Meaning this gene is only expressed in the type1 cells.

When requesting the comparison: results(dds,  contrast = c("CellType", "type1", "type2"), alpha = 0.05) I receive the following statistics:

- dds without setting the 'Mice' factor: log2FoldChange = 3.77127, padj = 1.38E-07

- dds with setting the 'Mice' factor: log2FoldChange = 0.313526, padj = 0.87749

Interestingly there are several such examples in which one cell type does not show any expression of the gene and the other shows different expression levels and I receive similar results when comparing them, having the data set with the pairing factor showing non-significant results with a strange log2FoldChange of close to 0.313.

Another example is the following gene:

Where only one of the pools supports any differences between the cell types.

Here, without pairing the Mice pool I get the non-significant (expected) result: log2FoldChange = 1.02, padj = 0.44

and with the setting: log2FoldChange = 1.25, padj = 0.000169

Am I doing something wrong? Is there an explanation for this that I am missing?

Thank you very much,

Olga

> sessionInfo()
R version 3.3.3 (2017-03-06)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

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] grid      splines   parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] matrixStats_0.51.0         ggrepel_0.6.5              PoiClaClu_1.0.2            gridExtra_2.2.1           
 [5] reshape_0.8.6              VennDiagram_1.6.17         futile.logger_1.4.3        colorRamps_2.3            
 [9] colorspace_1.3-2           randomcoloR_1.0.0          cowplot_0.7.0              Hmisc_4.0-2               
[13] Formula_1.2-1              survival_2.40-1            lattice_0.20-34            NOISeq_2.18.0             
[17] Matrix_1.2-8               pheatmap_1.0.8             gplots_3.0.1               RColorBrewer_1.1-2        
[21] DESeq2_1.14.1              SummarizedExperiment_1.4.0 Biobase_2.34.0             GenomicRanges_1.26.4      
[25] GenomeInfoDb_1.10.3        IRanges_2.8.2              S4Vectors_0.12.2           BiocGenerics_0.20.0       
[29] ggplot2_2.2.1              XLConnect_0.2-12           XLConnectJars_0.2-12      

loaded via a namespace (and not attached):
 [1] jsonlite_1.3         gtools_3.5.0         assertthat_0.1       latticeExtra_0.6-28  RSQLite_1.1-2        backports_1.0.5     
 [7] digest_0.6.12        XVector_0.14.1       checkmate_1.8.2      htmltools_0.3.5      plyr_1.8.4           XML_3.98-1.5        
[13] genefilter_1.56.0    zlibbioc_1.20.0      xtable_1.8-2         scales_0.4.1         gdata_2.17.0         BiocParallel_1.8.1  
[19] htmlTable_1.9        tibble_1.2           annotate_1.52.1      nnet_7.3-12          lazyeval_0.2.0       magrittr_1.5        
[25] memoise_1.0.0        foreign_0.8-67       tools_3.3.3          data.table_1.10.4    stringr_1.2.0        V8_1.3              
[31] munsell_0.4.3        locfit_1.5-9.1       cluster_2.0.5        AnnotationDbi_1.36.2 lambda.r_1.1.9       caTools_1.17.1      
[37] RCurl_1.95-4.8       htmlwidgets_0.8      labeling_0.3         bitops_1.0-6         base64enc_0.1-3      gtable_0.2.0        
[43] DBI_0.6              curl_2.4             knitr_1.15.1         futile.options_1.0.0 KernSmooth_2.23-15   rJava_0.9-8         
[49] stringi_1.1.3        Rcpp_0.12.10         geneplotter_1.52.0   rpart_4.1-10         acepack_1.4.1 

hi Olga,

I don't see anything surprising about the results here, from my interpretation of what the model would do when you change the design. Of course, it's important to keep in mind that you have very few samples in this dataset for use to estimate the four coefficients, and so the model has to do the best it can to interpret changes as either true differences or variability, in the absence of information from the counts for individual genes. In combination with each gene's information from the likelihood, there is sharing of information across genes. In situations like you show, where you have one pool of mice or one cell type with expression and the other with 0's, the model is leaning heavily (this happens automatically through the use in DESeq2 of the Bayes formula) on the estimation of inherent biological variability from other genes with more information in the counts.

Update: I just noticed you are showing the adjusted p-values. These depend on the p-values of other genes, which change when you change the design. So it's not enough to just look at the counts for one gene to try to figure out changes in the adjusted p-value (also there is the information sharing explanation above). A simple example:

> p.adjust(c(0,0,.05),method="BH")[3]
[1] 0.05
> p.adjust(c(1,1,.05),method="BH")[3]
[1] 0.15

The second case makes perfect sense to me, what seems like biological variability when you don't inform the model about mice pools is explained, and then the type 1 / type 2 ratio looks more significant afterwards.

The most important thing to appreciate is that there is the bare minimum of sample size to perform statistical testing, and so the model does the best it can by learning about inherent biological variability from other genes with more information.