I have previously been using DESeq2 as part of HOMER, but I think I need more control over the analysis of this experiment, and would love some input on how to set up the design for DESeq2 (plus I would learn!).
Experimental setup: RNA-seq (PE) of mice. I have wildtype (WT) and transgenic (TG) mice. Half of my TG mice develop periods of behavioral abnormalities. At four different time points 1 WT mouse, 1 "normal TG" (TGN) and 1 "disease TG" (TGD) has been sampled. This is summarized below:
Where "sampleName" is the uniqe name of each mouse, "Geno" is the genotype, "State" is normal or disease at the time of sampling, "Mate" is which litter the animal belongs to, "Run" is file name, "Batch" is sampling batch, and "Sum" is a collective of "Geno" and "State".
Now, I want to find genes that are differentially expressed in TG compared to WT as well as D compared to N mice. I would like to take into consideration the different sampling batches. After reading around at forums and the manual for DESeq2 I am sure I am not understanding how to properly decide on a design...
FYI, looking at PCA plot and heatmap there is no good separation of my samples based on anything. From the PCA plot I see that "Batch A" clusters somewhat far from the other samples and from the heatmap I see that "Batch D" cluster together while everything else is a hot mess. However, I am expecting only small differences between my conditions.
My first thought was to include an interaction effect between Geno and State, but that gave me an error (because all N animals are TG?)
>library(DESeq2) > dds = DESeqDataSet(se,design= ~Batch + Geno:State) Error in checkFullRank(modelMatrix) : the model matrix is not full rank, so the model cannot be fit as specified. Levels or combinations of levels without any samples have resulted in column(s) of zeros in the model matrix. Please read the vignette section 'Model matrix not full rank': vignette('DESeq2')
So instead I used following:
library(DESeq2) dds = DESeqDataSet(se,design= ~Batch+Sum)
Which has been producing OK results (as is I am getting a few differentially expressed genes). But then I can only get info on the "Sum" and not compare e.g. all WT to all TG using:
results(dds, contrast = c("Geno","TG","WT"), alpha = 0.05)
Which I would like to do.. (Maybe I am doing this wrong?)
I hope this has been somewhat clear. I hope for any help to guide a newbie in the right direction :)
Links to similar questions that I may have missed are also highly appreciated! Trying to learn how properly do this, but am finding it a bit challenging...
> sessionInfo() R version 3.4.2 (2017-09-28) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Ubuntu 16.04.3 LTS Matrix products: default BLAS: /usr/lib/libblas/libblas.so.3.6.0 LAPACK: /usr/lib/lapack/liblapack.so.3.6.0 locale:  LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C  LC_TIME=ja_JP.UTF-8 LC_COLLATE=en_US.UTF-8  LC_MONETARY=ja_JP.UTF-8 LC_MESSAGES=en_US.UTF-8  LC_PAPER=ja_JP.UTF-8 LC_NAME=C  LC_ADDRESS=C LC_TELEPHONE=C  LC_MEASUREMENT=ja_JP.UTF-8 LC_IDENTIFICATION=C attached base packages:  stats4 parallel stats graphics grDevices utils datasets  methods base other attached packages:  org.Mm.eg.db_3.5.0 genefilter_1.60.0  ggbeeswarm_0.6.0 PoiClaClu_1.0.2  RColorBrewer_1.1-2 pheatmap_1.0.8  hexbin_1.27.1 ggplot2_2.2.1  dplyr_0.7.4 DESeq2_1.18.1  BiocParallel_1.12.0 BiocInstaller_1.28.0  GenomicFeatures_1.30.0 AnnotationDbi_1.40.0  GenomicAlignments_1.14.1 Rsamtools_1.30.0  Biostrings_2.46.0 XVector_0.18.0  SummarizedExperiment_1.8.0 DelayedArray_0.4.1  matrixStats_0.52.2 Biobase_2.38.0  GenomicRanges_1.30.0 GenomeInfoDb_1.14.0  IRanges_2.12.0 S4Vectors_0.16.0  BiocGenerics_0.24.0 loaded via a namespace (and not attached):  RMySQL_0.10.13 bit64_0.9-7 splines_3.4.2  Formula_1.2-2 assertthat_0.2.0 latticeExtra_0.6-28  blob_1.1.0 vipor_0.4.5 GenomeInfoDbData_0.99.1  progress_1.1.2 RSQLite_2.0 backports_1.1.1  lattice_0.20-35 glue_1.2.0 digest_0.6.12  checkmate_1.8.5 colorspace_1.3-2 htmltools_0.3.6  Matrix_1.2-11 plyr_1.8.4 XML_3.98-1.9  pkgconfig_2.0.1 biomaRt_2.34.0 zlibbioc_1.24.0  xtable_1.8-2 scales_0.5.0 tibble_1.3.4  htmlTable_1.9 annotate_1.56.1 nnet_7.3-12  lazyeval_0.2.1 survival_2.41-3 magrittr_1.5  memoise_1.1.0 foreign_0.8-69 beeswarm_0.2.3  tools_3.4.2 data.table_1.10.4-3 prettyunits_1.0.2  stringr_1.2.0 munsell_0.4.3 locfit_1.5-9.1  cluster_2.0.6 bindrcpp_0.2 compiler_3.4.2  rlang_0.1.4 grid_3.4.2 RCurl_1.95-4.8  htmlwidgets_0.9 labeling_0.3 bitops_1.0-6  base64enc_0.1-3 gtable_0.2.0 DBI_0.7  reshape2_1.4.2 R6_2.2.2 gridExtra_2.3  knitr_1.17 rtracklayer_1.38.0 bit_1.1-12  bindr_0.1 Hmisc_4.0-3 stringi_1.1.6  Rcpp_0.12.14 geneplotter_1.56.0 rpart_4.1-11  acepack_1.4.1 >