Effect size using likelihood ratio test in DESeq2
1
0
Entering edit mode
jd348 • 0
@user-24862
Last seen 1 day ago
Durham

Hello,

I'm hoping to receive some guidance regarding an RNA-sequencing project I'm conducing right now. For this project, I have 9 time points, 3 biological replicates per time point, and 2 genotypes: wildtype vs mutant. I am using a likelihood ratio test in DESeq2 using the following code.

From a DESeq2 training page, the log2 fold change is printed in the results table for consistency with other results table outputs, but is not associated with the actual test. However, I am interested in detecting upregulated and downregulated genes. Seeing as I cannot use the log2 fold change from the DESeq results to measure effect size, is there any other tool I can use in DESeq2 to quantify effect size of these differentially expressed genes?

EDIT: to clarify, I want to find the genes that are consistently upregulated/downregulated across all 9 time points.


dds <- DESeqDataSetFromMatrix(ctdata, coldata, design = ~ genotype + age
dds <- DESeq(dds, test = 'LRT', reduced = ~ age)

sessionInfo( )
R version 4.1.0 (2021-05-18)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 11.4

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
 [1] grid      splines   parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] gprofiler2_0.2.1            ggformula_0.10.1            ggridges_0.5.3              scales_1.1.1               
 [5] ggstance_0.3.5              ggbeeswarm_0.6.0            org.Hs.eg.db_3.13.0         AnnotationDbi_1.54.1       
 [9] biomaRt_2.48.3              EnhancedVolcano_1.10.0      DEGreport_1.28.0            gplots_3.1.1               
[13] VennDiagram_1.6.20          futile.logger_1.4.3         WGCNA_1.70-3                fastcluster_1.2.3          
[17] dynamicTreeCut_1.63-1       WebGestaltR_0.4.4           pheatmap_1.0.12             ggpubr_0.4.0               
[21] openxlsx_4.2.4              BiocParallel_1.26.2         ggrepel_0.9.1               forcats_0.5.1              
[25] stringr_1.4.0               dplyr_1.0.7                 purrr_0.3.4                 readr_2.0.1                
[29] tidyr_1.1.3                 tibble_3.1.4                ggplot2_3.3.5               tidyverse_1.3.1            
[33] DESeq2_1.32.0               SummarizedExperiment_1.22.0 Biobase_2.52.0              MatrixGenerics_1.4.3       
[37] matrixStats_0.60.1          GenomicRanges_1.44.0        GenomeInfoDb_1.28.4         IRanges_2.26.0             
[41] S4Vectors_0.30.0            BiocGenerics_0.38.0        

loaded via a namespace (and not attached):
  [1] rappdirs_0.3.3              bit64_4.0.5                 knitr_1.34                  DelayedArray_0.18.0        
  [5] data.table_1.14.0           rpart_4.1-15                KEGGREST_1.32.0             RCurl_1.98-1.4             
  [9] doParallel_1.0.16           generics_0.1.0              preprocessCore_1.54.0       cowplot_1.1.1              
 [13] lambda.r_1.2.4              RSQLite_2.2.8               bit_4.0.4                   tzdb_0.1.2                 
 [17] xml2_1.3.2                  lubridate_1.7.10            assertthat_0.2.1            xfun_0.26                  
 [21] hms_1.1.0                   fansi_0.5.0                 progress_1.2.2              caTools_1.18.2             
 [25] dbplyr_2.1.1                readxl_1.3.1                igraph_1.2.6                DBI_1.1.1                  
 [29] geneplotter_1.70.0          tmvnsim_1.0-2               htmlwidgets_1.5.4           reshape_0.8.8              
 [33] apcluster_1.4.8             ellipsis_0.3.2              backports_1.2.1             annotate_1.70.0            
 [37] vctrs_0.3.8                 Cairo_1.5-12.2              abind_1.4-5                 cachem_1.0.6               
 [41] withr_2.4.2                 ggforce_0.3.3               lasso2_1.2-21.1             checkmate_2.0.0            
 [45] prettyunits_1.1.1           mnormt_2.0.2                svglite_2.0.0               cluster_2.1.2              
 [49] lazyeval_0.2.2              crayon_1.4.1                genefilter_1.74.0           edgeR_3.34.1               
 [53] pkgconfig_2.0.3             tweenr_1.0.2                nlme_3.1-153                vipor_0.4.5                
 [57] nnet_7.3-16                 rlang_0.4.11                lifecycle_1.0.0             filelock_1.0.2             
 [61] extrafontdb_1.0             BiocFileCache_2.0.0         modelr_0.1.8                polyclip_1.10-0            
 [65] ggrastr_0.2.3               cellranger_1.1.0            rngtools_1.5                Matrix_1.3-4               
 [69] carData_3.0-4               reprex_2.0.1                base64enc_0.1-3             beeswarm_0.4.0             
 [73] whisker_0.4                 GlobalOptions_0.1.2         viridisLite_0.4.0           png_0.1-7                  
 [77] rjson_0.2.20                bitops_1.0-7                ConsensusClusterPlus_1.56.0 KernSmooth_2.23-20         
 [81] Biostrings_2.60.2           blob_1.2.2                  doRNG_1.8.2                 shape_1.4.6                
 [85] jpeg_0.1-9                  rstatix_0.7.0               ggsignif_0.6.3              memoise_2.0.0              
 [89] magrittr_2.0.1              plyr_1.8.6                  zlibbioc_1.38.0             compiler_4.1.0             
 [93] RColorBrewer_1.1-2          ash_1.0-15                  clue_0.3-59                 cli_3.0.1                  
 [97] XVector_0.32.0              htmlTable_2.2.1             formatR_1.11                Formula_1.2-4              
[101] MASS_7.3-54                 tidyselect_1.1.1            stringi_1.7.4               proj4_1.0-10.1             
[105] locfit_1.5-9.4              latticeExtra_0.6-29         tools_4.1.0                 rio_0.5.27                 
[109] circlize_0.4.13             rstudioapi_0.13             foreach_1.5.1               foreign_0.8-81             
[113] logging_0.10-108            gridExtra_2.3               farver_2.1.0                digest_0.6.27              
[117] Rcpp_1.0.7                  car_3.0-11                  broom_0.7.9                 ggalt_0.4.0                
[121] httr_1.4.2                  Nozzle.R1_1.1-1             ggdendro_0.1.22             ComplexHeatmap_2.8.0       
[125] psych_2.1.6                 colorspace_2.0-2            rvest_1.0.1                 XML_3.99-0.7               
[129] fs_1.5.0                    plotly_4.9.4.1              systemfonts_1.0.2           xtable_1.8-4               
[133] jsonlite_1.7.2              futile.options_1.0.1        R6_2.5.1                    Hmisc_4.5-0                
[137] pillar_1.6.2                htmltools_0.5.2             glue_1.4.2                  fastmap_1.1.0              
[141] codetools_0.2-18            maps_3.3.0                  utf8_1.2.2                  lattice_0.20-44            
[145] curl_4.3.2                  gtools_3.9.2                zip_2.2.0                   GO.db_3.13.0               
[149] Rttf2pt1_1.3.9              survival_3.2-13             limma_3.48.3                munsell_0.5.0              
[153] GetoptLong_1.0.5            GenomeInfoDbData_1.2.6      iterators_1.0.13            labelled_2.8.0             
[157] impute_1.66.0               mosaicCore_0.9.0            haven_2.4.3                 gtable_0.3.0               
[161] extrafont_0.17
DESeq2 • 128 views
ADD COMMENT
1
Entering edit mode
ATpoint ▴ 860
@atpoint-13662
Last seen 2 hours ago
Germany

The thing with time-courses is that you do not necessarily have binary patterns of gene expression changes being only up-or downregulated. Over time you could have genes going up initially and then down again. Or down initially and then come back to the baseline, or any combination you could think of. Hence, a single logFC is not really informative. I would take the DEGs (by FDR) and make a heatmap (e.g. ComplexHeatmap package) with the scaled counts to visualize the patterns. I guess these patterns (which you can separate by e.g. hclust) are then what you should use to define the expression status of your genes over time.

ADD COMMENT

Login before adding your answer.

Traffic: 281 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