Hi,
Within my experiment, I have two age groups (6 months and 14 months) and one continuous variable (cell abundance). The bulk RNA sequencing samples and samples used to quantify cell abundance are biological replicates acquired from mouse BXD stains. To compare these samples, I independently averaged the count data and cell abundance data for these samples by group (strain and age combinations).
My goal is to determine which changes in gene expression are driven by changes in regional cell abundance (ex: Hippo_GFAP, hippocampal astrocyte abundance) vs gene expression changes that are related to aging alone.
I have set my DESeq2 design as ~ Age + Hippo_GFAP + Age:Hippo_GFAP and have run the following code looking at two different contrasts pulled from the results of resultsNames(dds_int)
Depending on the terms I include in the results function I get different values output as "low counts" (305 for name = "Age_14_vs_6" results, and 1218 for name = "Hippo_GFAP" results). Why would different values be output here when I am using the same dds object in the function? Of note, prior to running DESeq, I did perform preliminary filtering parameters to remove genes that had less than 10 counts in half of my samples.
Thanks!
> dds_int = DESeqDataSetFromMatrix(countData = group_count,
+ colData = group_key,
+ design = ~ Age + Hippo_GFAP + Age:Hippo_GFAP)
converting counts to integer mode
> #run_deseq2:
> dds_int = DESeq(dds_int)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
> # view_results_names:
> resultsNames(dds_int)
[1] "Intercept" "Age_14_vs_6" "Hippo_GFAP" "Age14.Hippo_GFAP"
> ######Interaction: Age contrast
> res_int_age <-results(dds_int, name = "Age_14_vs_6")
> summary(res_int_age)
out of 15703 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up) : 117, 0.75%
LFC < 0 (down) : 235, 1.5%
outliers [1] : 0, 0%
low counts [2] : 305, 1.9%
(mean count < 13)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
> ######Interaction: Cell abundance
> res_int_load <-results(dds_int, name = "Hippo_GFAP")
> summary(res_int_load)
out of 15703 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up) : 461, 2.9%
LFC < 0 (down) : 290, 1.8%
outliers [1] : 0, 0%
low counts [2] : 1218, 7.8%
(mean count < 22)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
session info:
sessionInfo( )
R version 4.0.0 (2020-04-24)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 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] DESeq2_1.28.1 SummarizedExperiment_1.18.2 DelayedArray_0.14.1 matrixStats_0.61.0 Biobase_2.48.0
[6] GenomicRanges_1.40.0 GenomeInfoDb_1.24.2 IRanges_2.22.2 S4Vectors_0.26.1 BiocGenerics_0.34.0
[11] sva_3.35.2 BiocParallel_1.22.0 genefilter_1.70.0 mgcv_1.8-40 nlme_3.1-157
[16] ggfortify_0.4.14 ggplot2_3.3.5 dplyr_1.0.7
loaded via a namespace (and not attached):
[1] tidyr_1.0.2 edgeR_3.30.3 bit64_4.0.5 splines_4.0.0 assertthat_0.2.1 blob_1.2.3 GenomeInfoDbData_1.2.3
[8] pillar_1.7.0 RSQLite_2.2.12 lattice_0.20-45 glue_1.4.0 limma_3.44.3 digest_0.6.29 RColorBrewer_1.1-3
[15] XVector_0.28.0 colorspace_1.4-1 cowplot_1.1.1 Matrix_1.4-1 XML_3.99-0.9 pkgconfig_2.0.3 zlibbioc_1.34.0
[22] purrr_0.3.4 xtable_1.8-4 scales_1.2.0 tibble_3.0.1 annotate_1.66.0 farver_2.1.0 generics_0.1.2
[29] ellipsis_0.3.2 cachem_1.0.6 withr_2.5.0 cli_3.1.1 survival_3.3-1 magrittr_2.0.2 crayon_1.5.1
[36] memoise_2.0.1 fansi_0.4.1 tools_4.0.0 lifecycle_1.0.1 stringr_1.4.0 munsell_0.5.0 locfit_1.5-9.4
[43] AnnotationDbi_1.50.3 compiler_4.0.0 tinytex_0.38 rlang_1.0.2 grid_4.0.0 RCurl_1.98-1.6 rstudioapi_0.13
[50] labeling_0.4.2 bitops_1.0-7 gtable_0.3.0 DBI_1.1.2 R6_2.5.1 gridExtra_2.3 fastmap_1.1.0
[57] bit_4.0.4 utf8_1.1.4 stringi_1.4.6 Rcpp_1.0.8.3 vctrs_0.3.8 geneplotter_1.66.0 tidyselect_1.1.2
[64] xfun_0.30
Michael Love