Main effects
dds <- DESeqDataSetFromTximport(txi.rsem, colData = metadata, design = ~ Group + Trt + Region)## using counts and average transcript lengths from tximportdds <- dds[,-21] # exclude MSL3 outlier by cluster
# filter low expressed genes
keep <- rowSums(counts(dds)) >= 50
table(keep)## keep
## FALSE TRUE
## 37491 17996dds <- dds[keep,]
# add in gene symbols
genes <- gsub("\\..*","",rownames(dds))
symbol <- mapIds(org.Mm.eg.db, keys=genes, column="SYMBOL", keytype="ENSEMBL", multiVals="first")## 'select()' returned 1:many mapping between keys and columnsmcols(dds) <- cbind(mcols(dds), symbol)
dds <- DESeq(dds)## estimating size factors## using 'avgTxLength' from assays(dds), correcting for library size## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingresultsNames(dds)## [1] "Intercept" "Group_MIA_vs_CON" "Trt_LPS_vs_SAL"
## [4] "Region_Str_vs_Ctx"group = results(dds, contrast=c("Group","MIA","CON"), alpha=0.05)
summary(group) # up: 35 down: 128##
## out of 17996 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 35, 0.19%
## LFC < 0 (down) : 128, 0.71%
## outliers [1] : 71, 0.39%
## low counts [2] : 4495, 25%
## (mean count < 17)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?resultstrt = results(dds, contrast=c("Trt","LPS","SAL"), alpha=0.05)
summary(trt) # up: 3480 down: 3359##
## out of 17996 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 3480, 19%
## LFC < 0 (down) : 3359, 19%
## outliers [1] : 71, 0.39%
## low counts [2] : 683, 3.8%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?resultsregion = results(dds, contrast=c("Region","Str","Ctx"), alpha=0.05)
summary(region) # up: 461 down: 190##
## out of 17996 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 461, 2.6%
## LFC < 0 (down) : 190, 1.1%
## outliers [1] : 71, 0.39%
## low counts [2] : 0, 0%
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?resultsVolcano plots of DEGs between Group, Treatment, and Brain Region
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/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] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] EnhancedVolcano_1.9.13 ggrepel_0.9.1
## [3] org.Mm.eg.db_3.12.0 AnnotationDbi_1.52.0
## [5] DESeq2_1.30.1 SummarizedExperiment_1.20.0
## [7] Biobase_2.50.0 MatrixGenerics_1.2.1
## [9] matrixStats_0.59.0 GenomicRanges_1.42.0
## [11] GenomeInfoDb_1.26.7 IRanges_2.24.1
## [13] S4Vectors_0.28.1 BiocGenerics_0.36.1
## [15] cowplot_1.1.1 forcats_0.5.1
## [17] stringr_1.4.0 dplyr_1.0.7
## [19] purrr_0.3.4 readr_1.4.0
## [21] tidyr_1.1.3 tibble_3.1.2
## [23] ggplot2_3.3.5 tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
## [1] ggbeeswarm_0.6.0 colorspace_2.0-2 ellipsis_0.3.2
## [4] XVector_0.30.0 fs_1.5.0 rstudioapi_0.13
## [7] farver_2.1.0 bit64_4.0.5 mvtnorm_1.1-2
## [10] apeglm_1.12.0 fansi_0.5.0 lubridate_1.7.10
## [13] xml2_1.3.2 splines_4.0.3 extrafont_0.17
## [16] cachem_1.0.5 geneplotter_1.68.0 knitr_1.33
## [19] jsonlite_1.7.2 broom_0.7.8 Rttf2pt1_1.3.8
## [22] annotate_1.68.0 dbplyr_2.1.1 compiler_4.0.3
## [25] httr_1.4.2 backports_1.2.1 assertthat_0.2.1
## [28] Matrix_1.3-4 fastmap_1.1.0 cli_2.5.0
## [31] htmltools_0.5.1.1 tools_4.0.3 coda_0.19-4
## [34] gtable_0.3.0 glue_1.4.2 GenomeInfoDbData_1.2.4
## [37] maps_3.3.0 Rcpp_1.0.6 bbmle_1.0.23.1
## [40] cellranger_1.1.0 jquerylib_0.1.4 vctrs_0.3.8
## [43] ggalt_0.4.0 extrafontdb_1.0 xfun_0.24
## [46] rvest_1.0.0 lifecycle_1.0.0 XML_3.99-0.6
## [49] zlibbioc_1.36.0 MASS_7.3-54 scales_1.1.1
## [52] hms_1.1.0 proj4_1.0-10.1 RColorBrewer_1.1-2
## [55] yaml_2.2.1 memoise_2.0.0 ggrastr_0.2.3
## [58] emdbook_1.3.12 sass_0.4.0 bdsmatrix_1.3-4
## [61] stringi_1.6.2 RSQLite_2.2.7 highr_0.9
## [64] genefilter_1.72.1 BiocParallel_1.24.1 rlang_0.4.11
## [67] pkgconfig_2.0.3 bitops_1.0-7 evaluate_0.14
## [70] lattice_0.20-44 labeling_0.4.2 bit_4.0.4
## [73] tidyselect_1.1.1 plyr_1.8.6 magrittr_2.0.1
## [76] R6_2.5.0 generics_0.1.0 DelayedArray_0.16.3
## [79] DBI_1.1.1 pillar_1.6.1 haven_2.4.1
## [82] withr_2.4.2 survival_3.2-11 RCurl_1.98-1.3
## [85] ash_1.0-15 modelr_0.1.8 crayon_1.4.1
## [88] KernSmooth_2.23-20 utf8_1.2.1 rmarkdown_2.9
## [91] locfit_1.5-9.4 grid_4.0.3 readxl_1.3.1
## [94] blob_1.2.1 reprex_2.0.0 digest_0.6.27
## [97] xtable_1.8-4 numDeriv_2016.8-1.1 munsell_0.5.0
## [100] beeswarm_0.4.0 vipor_0.4.5 bslib_0.2.5.1