Dear All,
I am running a DESeq2 within-subject treatment response analysis. I have RNA seq data on 75 individuals, for most of these I have data before and after treatment (but not all). I also have a variable indicating successful treatment response across all individuals. I want to model gene expression as dependent on treatment and treatment response (res01) within each subject. counts is an R object with counts across all samples and genes.
The design matrix:
The analysis is ran as follows:
Originally I though that I could extract the effect of res01 on gene expression as:
At this point I realized that the results looked differently than I expected:
I have two questions regarding this:
1. Why do I get two res01 results and two treatment results? This seems to be due to the within-subject design. I don't understand which model that has been tested here explicitly.
2. Should I run the following code to extract the effect that res01 has on gene expression given all other covariates?
Thanks in advance for any help,
Boel
My sessionInfo:
I am running a DESeq2 within-subject treatment response analysis. I have RNA seq data on 75 individuals, for most of these I have data before and after treatment (but not all). I also have a variable indicating successful treatment response across all individuals. I want to model gene expression as dependent on treatment and treatment response (res01) within each subject. counts is an R object with counts across all samples and genes.
The design matrix:
Code:
> head(A_design) sampleID subject treatment res01 A10a A10a A10 0 0 A10b A10b A10 1 1 A11a A11a A11 0 0 A11b A11b A11 1 0 A12a A12a A12 0 0 A12b A12b A12 1 1
Code:
dds <- DESeqDataSetFromMatrix(countData = counts[,A_design[,1]], colData = as.data.frame(A_design), design = ~ subject + treatment + res01) Atreat.dds <- DESeq(dds)
Code:
> res <- results(Atreat.dds,name='res01',pAdjustMethod='BH') Error in results(Atreat.dds, name = "res01", pAdjustMethod = "BH") : cannot find appropriate results in the DESeqDataSet. possibly nbinomWaldTest or nbinomLRT has not yet been run.
Code:
> resultsNames(Atreat.dds) [1] "Intercept" "subjectA01" "subjectA02" "subjectA03" "subjectA04" [6] "subjectA05" "subjectA06" "subjectA07" "subjectA08" "subjectA09" ..... [76] "subjectA92" "treatment0" "treatment1" "res010" "res011"
1. Why do I get two res01 results and two treatment results? This seems to be due to the within-subject design. I don't understand which model that has been tested here explicitly.
2. Should I run the following code to extract the effect that res01 has on gene expression given all other covariates?
Code:
Ares_res01=results(Atreat.dds, contrast=list("res010",'res011'))
Boel
My sessionInfo:
Code:
> sessionInfo() R version 3.2.1 (2015-06-18) Platform: x86_64-redhat-linux-gnu (64-bit) Running under: Scientific Linux release 6.7 (Carbon) locale: [1] LC_CTYPE=sv_SE.UTF-8 LC_NUMERIC=C [3] LC_TIME=sv_SE.UTF-8 LC_COLLATE=sv_SE.UTF-8 [5] LC_MONETARY=sv_SE.UTF-8 LC_MESSAGES=sv_SE.UTF-8 [7] LC_PAPER=sv_SE.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=sv_SE.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] parallel stats4 stats graphics grDevices utils datasets [8] methods base other attached packages: [1] DESeq2_1.8.1 RcppArmadillo_0.5.400.2.0 [3] Rcpp_0.12.0 GenomicRanges_1.20.5 [5] GenomeInfoDb_1.4.2 IRanges_2.2.7 [7] S4Vectors_0.6.3 BiocGenerics_0.14.0 loaded via a namespace (and not attached): [1] RColorBrewer_1.1-2 futile.logger_1.4.1 plyr_1.8.3 [4] XVector_0.8.0 futile.options_1.0.0 tools_3.2.1 [7] rpart_4.1-10 digest_0.6.8 RSQLite_1.0.0 [10] annotate_1.46.1 gtable_0.1.2 lattice_0.20-33 [13] DBI_0.3.1 proto_0.3-10 gridExtra_2.0.0 [16] genefilter_1.50.0 stringr_1.0.0 cluster_2.0.3 [19] locfit_1.5-9.1 nnet_7.3-10 grid_3.2.1 [22] Biobase_2.28.0 AnnotationDbi_1.30.1 XML_3.98-1.3 [25] survival_2.38-3 BiocParallel_1.2.20 foreign_0.8-66 [28] latticeExtra_0.6-26 Formula_1.2-1 geneplotter_1.46.0 [31] ggplot2_1.0.1 reshape2_1.4.1 lambda.r_1.1.7 [34] magrittr_1.5 scales_0.2.5 Hmisc_3.16-0 [37] MASS_7.3-43 splines_3.2.1 xtable_1.7-4 [40] colorspace_1.2-6 stringi_0.5-5 acepack_1.3-3.3 [43] munsell_0.4.2