Dear all, I am a bit confused with my results. I read userGuide especially chapter 9 and 12 but I am still not sure I have figured it out. I have an experiment with 8 subjects (used as blocking factor) and 6 timepoints (one is the timepoint before treatment, aka screening) and also I know I have batch effect (so I add my 3 batches in the model) .
Subject Timepoint Batch
1001 Screening 3
1001 2 1
1001 4 1
1001 8 1
1001 43 1
1001 91 1
1002 Screening 3
1002 2 1
....
1008 Screening 2
1008 2 2
1008 4 2
1008 8 2
1008 43 2
1008 91 2
What I would like to see is if I have genes that are changing between timepoints, compared to screening timepoint What confuses me is that when I set my contrasts and run the top table (which is an anova like approach as far as I understood) I get 140 genes DE (not a lot based on what I would expect). But when I run the decideTest I get only a handful of genes. my design is
designr = model.matrix( ~ 0 +Timepoint + Batch , metadata)
And contrasts
contrast.matrix<-makeContrasts(Timepoint2-TimepointScreening, Timepoint4-TimepointScreening , Timepoint8-TimepointScreening , Timepoint43-TimepointScreening, Timepoint91-TimepointScreening,levels = designr)
My question is
a)Am I doing something wrong?
b)Most importantly.. If not how can I interpret this difference?
c) Also what I am still not clear. If I a gene is DE between timepoint 8 vs 2 but not between timepoint 8 vs screening or timepoint 2 vs screening, I dont except to see that popping up in my toptable, right? Is a thought while trying to understand the difference between decideTests and topTable..
# apply duplicateCorrelation is two rounds
set.seed(1234)
designr = model.matrix( ~ 0 +Timepoint + Batch , metadata)
vobj_tmpr = voom( geneExpr, designr, plot=FALSE)
dupcorr <- duplicateCorrelation(vobj_tmpr,designr,block=metadata$Subject)
# run voom considering the duplicateCorrelation results
# in order to compute more accurate precision weights
vobjr = voom( geneExpr, designr, plot=TRUE, block=metadata$Subject, correlation=dupcorr$consensus)
dupcorr <- duplicateCorrelation(vobjr, designr, block=metadata$Subject)
# Fit the model for each gene,
# this step uses only the genome-wide average for the random effect
fitDupCor <- lmFit(vobjr, designr, block=metadata$Subject, correlation=dupcorr$consensus)
fit0r <- eBayes(fitDupCor)
#make contrasts
contrast.matrix<-makeContrasts(Timepoint2-TimepointScreening,
Timepoint4-TimepointScreening ,
Timepoint8-TimepointScreening ,
Timepoint43-TimepointScreening,
Timepoint91-TimepointScreening,levels = designr)
#fit again
fitr<-contrasts.fit(fit0r,contrast.matrix)
fitr<-eBayes(fitr)
#get results from Top table
pr<-topTable(fitr,adjust="BH",p.value=0.1, number = "inf")
dim(pr)
## Result: 145 9
#decide tests using separate method
degr <-decideTests(fitr, p=0.1, adjust="BH")
summary(degr)
# Timepoint2 - TimepointScreening Timepoint4 - TimepointScreening
#Down 0 0
#NotSig 19277 19277
#Up 0 0
# Timepoint8 - TimepointScreening Timepoint43 - TimepointScreening
#Down 4 1
#NotSig 19266 19276
#Up 7 0
# Timepoint91 - TimepointScreening
#Down 0
#NotSig 19277
#Up 0
sessionInfo( )
R version 4.0.1 (2020-06-06)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.3 LTS
Matrix products: default
BLAS: /opt/R/4.0.1/lib/R/lib/libRblas.so
LAPACK: /opt/R/4.0.1/lib/R/lib/libRlapack.so
locale:
[1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8 LC_COLLATE=C.UTF-8
[5] LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8 LC_PAPER=C.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats4 parallel stats graphics grDevices utils datasets methods base
other attached packages:
[1] factoextra_1.0.7 GenomicFeatures_1.42.3 RColorBrewer_1.1-2
[4] igraph_1.4.3 WGCNA_1.72-1 fastcluster_1.2.3
[7] dynamicTreeCut_1.63-1 AnnotationDbi_1.52.0 janitor_2.2.0
[10] ggpubr_0.6.0 nlme_3.1-162 kableExtra_1.3.4
[13] variancePartition_1.20.0 scales_1.2.1 BiocParallel_1.24.1
[16] pheatmap_1.0.12 lubridate_1.9.2 forcats_1.0.0
[19] stringr_1.5.0 dplyr_1.1.2 purrr_1.0.1
[22] readr_2.1.4 tidyr_1.3.0 tibble_3.2.1
[25] ggplot2_3.4.2 tidyverse_2.0.0 edgeR_3.32.1
[28] limma_3.46.0 DESeq2_1.30.1 SummarizedExperiment_1.20.0
[31] Biobase_2.50.0 MatrixGenerics_1.2.1 matrixStats_0.58.0
[34] GenomicRanges_1.42.0 GenomeInfoDb_1.26.7 IRanges_2.24.1
[37] S4Vectors_0.28.1 AnnotationHub_2.22.1 BiocFileCache_1.14.0
[40] dbplyr_2.3.2 BiocGenerics_0.36.0
loaded via a namespace (and not attached):
[1] utf8_1.2.1 tidyselect_1.2.0 lme4_1.1-33
[4] RSQLite_2.2.4 htmlwidgets_1.5.3 grid_4.0.1
[7] munsell_0.5.0 codetools_0.2-18 preprocessCore_1.52.1
[10] statmod_1.4.35 DT_0.17 withr_2.5.0
[13] colorspace_2.0-0 highr_0.8 knitr_1.31
[16] rstudioapi_0.14 ggsignif_0.6.4 labeling_0.4.2
[19] GenomeInfoDbData_1.2.4 farver_2.1.0 bit64_4.0.5
[22] rprojroot_2.0.2 vctrs_0.6.2 generics_0.1.0
[25] xfun_0.22 timechange_0.2.0 R6_2.5.0
[28] doParallel_1.0.16 locfit_1.5-9.4 bitops_1.0-6
[31] cachem_1.0.4 DelayedArray_0.16.2 vroom_1.6.3
[34] promises_1.2.0.1 nnet_7.3-15 gtable_0.3.0
[37] rlang_1.1.1 genefilter_1.72.1 systemfonts_1.0.4
[40] splines_4.0.1 rtracklayer_1.50.0 rstatix_0.7.2
[43] impute_1.64.0 broom_1.0.5 checkmate_2.0.0
[46] BiocManager_1.30.10 yaml_2.3.7 reshape2_1.4.4
[49] abind_1.4-5 crosstalk_1.1.1 backports_1.2.1
[52] httpuv_1.5.5 Hmisc_4.5-0 tools_4.0.1
[55] ellipsis_0.3.2 gplots_3.1.1 jquerylib_0.1.3
[58] Rcpp_1.0.6 plyr_1.8.6 base64enc_0.1-3
[61] progress_1.2.2 zlibbioc_1.36.0 RCurl_1.98-1.12
[64] prettyunits_1.1.1 rpart_4.1-15 openssl_2.0.6
[67] cowplot_1.1.1 ggrepel_0.9.3 cluster_2.1.1
[70] here_1.0.1 fs_1.6.2 colorRamps_2.3.1
[73] tinytex_0.30 magrittr_2.0.3 data.table_1.14.0
[76] stringfish_0.15.8 qs_0.25.5 hms_1.1.3
[79] mime_0.10 evaluate_0.14 xtable_1.8-4
[82] pbkrtest_0.5.1 XML_3.99-0.14 jpeg_0.1-8.1
[85] gridExtra_2.3 compiler_4.0.1 biomaRt_2.46.3
[88] KernSmooth_2.23-18 crayon_1.4.1 minqa_1.2.4
[91] htmltools_0.5.5 later_1.1.0.1 tzdb_0.4.0
[94] Formula_1.2-4 geneplotter_1.68.0 RcppParallel_5.1.7
[97] RApiSerialize_0.1.2 DBI_1.1.1 MASS_7.3-53.1
[100] rappdirs_0.3.3 boot_1.3-27 Matrix_1.3-2
[103] car_3.1-2 cli_3.6.1 pkgconfig_2.0.3
[106] GenomicAlignments_1.26.0 foreign_0.8-81 xml2_1.3.4
[109] foreach_1.5.1 svglite_2.0.0 annotate_1.68.0
[112] bslib_0.4.2 webshot_0.5.2 XVector_0.30.0
[115] rvest_1.0.3 snakecase_0.11.0 digest_0.6.31
[118] Biostrings_2.58.0 rmarkdown_2.7 htmlTable_2.1.0
[121] curl_5.0.1 Rsamtools_2.6.0 shiny_1.6.0
[124] gtools_3.8.2 nloptr_1.2.2.2 jsonlite_1.8.5
[127] lifecycle_1.0.3 carData_3.0-4 viridisLite_0.3.0
[130] askpass_1.1 fansi_0.4.2 pillar_1.9.0
[133] lattice_0.20-41 fastmap_1.1.0 httr_1.4.6
[136] survival_3.2-10 GO.db_3.12.1 interactiveDisplayBase_1.28.0
[139] glue_1.6.2 png_0.1-7 iterators_1.0.13
[142] BiocVersion_3.12.0 bit_4.0.4 sass_0.4.6
[145] stringi_1.5.3 blob_1.2.1 latticeExtra_0.6-29
[148] caTools_1.18.1 memoise_2.0.1
Thanks Gordon! unfortunately in our company we only have R 4.0.1 so latest bioconductor versions cause some issues :(. Thank you very much!