decideTests vs topTable results
1
0
Entering edit mode
Anna ▴ 20
@9d8a27d9
Last seen 3 days ago
Germany

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
decideTests topTable limma • 1.1k views
ADD COMMENT
0
Entering edit mode
@gordon-smyth
Last seen 2 hours ago
WEHI, Melbourne, Australia

a) Looks ok.

b) You've almost answered the question yourself. topTable() is doing an anova F-test on 5 df for all the contrasts at once, whereas decideTests() as you can see is just giving results for each of the contrasts individually. The combined test is the one you want and in this case is more powerful.

c) It is perfectly possible for a gene to be DE for the anova F-test test but not statistically significant for any contrast t-test individually. The anova F-test combines information from all the comparisons.

BTW, you are using limma and Bioconductor software from nearly 3 years ago. I'd recommend upgrading to current versions of R and Bioconductor.

ADD COMMENT
0
Entering edit mode

Thanks Gordon! unfortunately in our company we only have R 4.0.1 so latest bioconductor versions cause some issues :(. Thank you very much!

ADD REPLY

Login before adding your answer.

Traffic: 702 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6