I need help in making sure my model is correct for my experiment and that I am running the analysis correct using the LRT for timecourse data with Deseq2 in R. I am trying to run a differential gene expression analysis with RNAseq data for seeds from multiple years (3) and imbibed over multiple timepoints (3). I have tried using the LRT example and wanted to check if I am doing it correctly? I also want to do comparisons between years across certain timepoints. When I try to write contrast statements for the later part I get an error on the names I am using. Could you help me identify the correct trems to use if I am trying to run a contrast for example for seeds from 1999 and 1995 time point 0. Error message contained after the problem code. Also, the resultsNames(dds) do not seem to contain all possible combinations, so I am not sure if I am doing my analysis correctly.
Another question, I would like to have my T0 (timepoint 0) be the level that expression is based off for DEG, have I done that correctly with the 'relevel' usage?
Here is my metadata: metadata Year Timepoint Group gm10 H1999 T24 1999_T24 gm11 H1995 T24 1995_T24 gm12 H1995 T24 1995_T24 gm13 H1999 T24 1999_T24 gm14 H1999 T48 1999_T48 gm15 H2019 T48 2019_T48 gm16 H1999 T48 1999_T48 gm17 H2019 T24 2019_T24 gm18 H1995 T24 1995_T24 gm19 H1999 T0 1999_T0 gm1 H2019 T24 2019_T24 gm20 H1995 T0 1995_T0 gm21 H1999 T0 1999_T0 gm22 H1999 T0 1999_T0 gm23 H1995 T0 1995_T0 gm24 H1999 T48 1999_T48 gm25 H2019 T48 2019_T48 gm26 H1999 T48 1999_T48 gm27 H1995 T48 1995_T48 gm28 H2019 T24 2019_T24 gm29 H1999 T0 1999_T0 gm2 H2019 T0 2019_T0 gm30 H1999 T24 1999_T24 gm31 H1995 T0 1995_T0 gm32 H2019 T0 2019_T0 gm33 H1999 T24 1999_T24 gm34 H2019 T48 2019_T48 gm35 H1995 T24 1995_T24 gm36 H2019 T24 2019_T24 gm3 H1995 T0 1995_T0 gm4 H1995 T48 1995_T48 gm5 H2019 T48 2019_T48 gm6 H1995 T48 1995_T48 gm7 H2019 T0 2019_T0 gm8 H1995 T48 1995_T48 gm9 H2019 T0 2019_T0
Code run for Differential expression: ```dds<- DESeqDataSetFromMatrix(countData = mycounts, colData = metadata, design=~Year+Timepoint+Year:Timepoint)
dds$Timepoint <- relevel(dds$Timepoint, ref="T0") dds<- DESeq(dds)
res <- results(dds)
resultsNames(dds)
[1] "Intercept" "Year_H1999_vs_H1995" "Year_H2019_vs_H1995" "Timepoint_T24_vs_T0"
[5] "Timepoint_T48_vs_T0" "YearH1999.TimepointT24" "YearH2019.TimepointT24" "YearH1999.TimepointT48"
[9] "YearH2019.TimepointT48"
ddsLRT <- DESeq(dds, reduced=~Year+Timepoint, test="LRT")
res<-results(ddsLRT)
resultsNames(ddsLRT)
[1] "Intercept" "Year_H1999_vs_H1995" "Year_H2019_vs_H1995" "Timepoint_T24_vs_T0"
[5] "Timepoint_T48_vs_T0" "YearH1999.TimepointT24" "YearH2019.TimepointT24" "YearH1999.TimepointT48"
[9] "YearH2019.TimepointT48"
res_LRT_IHW <- results(ddsLRT, filterFun = ihw )
res_LRT_IHW <- results(ddsLRT, filterFun = ihw ) head(res_LRT_IHW) summary(res_LRT_IHW) LRTpadj0.05 <-subset(res_LRT_IHW, padj<0.05) summary(LRTpadj0.05) write.csv(LRTpadj0.05, file = "LRT_03302022.csv")
For contrasts within/across timepoints and across seed harvest years
res99_24<strong text- results(dds, contrast = c("Year", "H1999_T0", "H1999_T24"))
Error in cleanContrast(object, contrast, expanded = isExpanded, listValues = listValues, : H1999_T0 and H1999_T24 should be levels of Year such that Year_H1999_T0_vs_H1995 and Year_H1999_T24_vs_H1995 are contained in 'resultsNames(object)'
sessionInfo( )
sessionInfo( ) R version 4.0.1 (2020-06-06) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows >= 8 x64 (build 9200)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252
attached base packages: [1] parallel stats4 stats graphics grDevices utils datasets methods base
other attached packages:
[1] IHW_1.16.0 forcats_0.5.1 stringr_1.4.0
[4] purrr_0.3.4 readr_2.1.2 tidyr_1.2.0
[7] tibble_3.1.6 ggplot2_3.3.5 tidyverse_1.3.1
[10] dplyr_1.0.8 BiocManager_1.30.16 DESeq2_1.28.1
[13] SummarizedExperiment_1.18.2 DelayedArray_0.14.1 matrixStats_0.61.0
[16] Biobase_2.48.0 GenomicRanges_1.40.0 GenomeInfoDb_1.24.2
[19] IRanges_2.22.2 S4Vectors_0.26.1 BiocGenerics_0.34.0
loaded via a namespace (and not attached):
[1] bitops_1.0-7 fs_1.5.2 lubridate_1.8.0 bit64_4.0.5
[5] RColorBrewer_1.1-2 httr_1.4.2 tools_4.0.1 backports_1.4.1
[9] utf8_1.2.2 R6_2.5.1 DBI_1.1.2 colorspace_2.0-2
[13] withr_2.5.0 tidyselect_1.1.2 bit_4.0.4 compiler_4.0.1
[17] fdrtool_1.2.17 cli_3.1.1 rvest_1.0.2 xml2_1.3.3
[21] slam_0.1-50 scales_1.1.1 genefilter_1.70.0 XVector_0.28.0
[25] pkgconfig_2.0.3 lpsymphony_1.16.0 dbplyr_2.1.1 fastmap_1.1.0
[29] rlang_1.0.1 readxl_1.3.1 rstudioapi_0.13 RSQLite_2.2.10
[33] generics_0.1.2 jsonlite_1.8.0 BiocParallel_1.22.0 RCurl_1.98-1.6
[37] magrittr_2.0.2 GenomeInfoDbData_1.2.3 Matrix_1.4-0 Rcpp_1.0.8
[41] munsell_0.5.0 fansi_1.0.2 lifecycle_1.0.1 stringi_1.7.6
[45] yaml_2.2.2 zlibbioc_1.34.0 grid_4.0.1 blob_1.2.2
[49] crayon_1.5.1 lattice_0.20-45 haven_2.4.3 splines_4.0.1
[53] annotate_1.66.0 hms_1.1.1 locfit_1.5-9.4 pillar_1.7.0
[57] geneplotter_1.66.0 reprex_2.0.1 XML_3.99-0.8 glue_1.6.1
[61] modelr_0.1.8 vctrs_0.3.8 tzdb_0.2.0 cellranger_1.1.0
[65] gtable_0.3.0 assertthat_0.2.1 cachem_1.0.6 xtable_1.8-4
[69] broom_0.7.12 survival_3.2-13 AnnotationDbi_1.50.3 memoise_2.0.1
[73] ellipsis_0.3.2
``` Thank you for all your help.
