dds <- DESeqDataSetFromTximport(txi.rsem, colData = metadata, design = ~ Group)
# filter low expressed genes and samples
dds <- dds[,samp]
keep <- rowSums(counts(dds)) >= 50
table(keep)
## keep
## FALSE TRUE
## 38292 17195
dds <- dds[keep,]
# add in gene symbols
genes <- gsub("\\..*","",rownames(dds))
library(AnnotationHub)
ah <- AnnotationHub()
edb <- ah[["AH53222"]]
symbol <- mapIds(edb, keys=genes, column="SYMBOL", keytype="GENEID", multiVals="first")
mcols(dds) <- cbind(mcols(dds), symbol)
dds <- DESeq(dds)
resultsNames(dds)
## [1] "Intercept" "Group_MIA_vs_CON"
# difference between MIA and CON
summary(results(dds, name="Group_MIA_vs_CON", alpha=0.05))
##
## out of 17195 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 1066, 6.2%
## LFC < 0 (down) : 968, 5.6%
## outliers [1] : 27, 0.16%
## low counts [2] : 2325, 14%
## (mean count < 16)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
# up:1066 down: 968
MvC <- results(dds, name="Group_MIA_vs_CON", alpha=0.05)
MvC$gene <- symbol
library(pheatmap)
vsd <- vst(dds, blind=FALSE)
rownames(vsd) <- symbol
df <- as.data.frame(colData(dds)["Group"])
ann_colors = list(Group = c(CON = "black", MIA = "red"))
DEG <- order(MvC$pvalue)[1:1000]
pheatmap(assay(vsd)[DEG,], cluster_rows=TRUE, show_rownames=FALSE, cluster_cols=F, annotation_col=df, scale="row", border_color = NA, col = colorRampPalette(c("navy", "white", "firebrick3"))(50), annotation_colors = ann_colors, main="MIA vs CON")
# scaled transformed counts
Volcano plots of DEGs between Group
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] pheatmap_1.0.12 ensembldb_2.14.1
## [5] AnnotationFilter_1.14.0 GenomicFeatures_1.42.3
## [7] AnnotationHub_2.22.1 BiocFileCache_1.14.0
## [9] dbplyr_2.1.1 forcats_0.5.1
## [11] stringr_1.4.0 dplyr_1.0.7
## [13] purrr_0.3.4 readr_1.4.0
## [15] tidyr_1.1.3 tibble_3.1.2
## [17] tidyverse_1.3.1 org.Mm.eg.db_3.12.0
## [19] AnnotationDbi_1.52.0 DESeq2_1.30.1
## [21] SummarizedExperiment_1.20.0 Biobase_2.50.0
## [23] MatrixGenerics_1.2.1 matrixStats_0.59.0
## [25] GenomicRanges_1.42.0 GenomeInfoDb_1.26.7
## [27] IRanges_2.24.1 S4Vectors_0.28.1
## [29] BiocGenerics_0.36.1 cowplot_1.1.1
## [31] ggplot2_3.3.5
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 backports_1.2.1
## [3] plyr_1.8.6 lazyeval_0.2.2
## [5] splines_4.0.3 BiocParallel_1.24.1
## [7] digest_0.6.27 htmltools_0.5.1.1
## [9] fansi_0.5.0 magrittr_2.0.1
## [11] memoise_2.0.0 Biostrings_2.58.0
## [13] annotate_1.68.0 modelr_0.1.8
## [15] extrafont_0.17 extrafontdb_1.0
## [17] bdsmatrix_1.3-4 askpass_1.1
## [19] prettyunits_1.1.1 colorspace_2.0-2
## [21] apeglm_1.12.0 blob_1.2.1
## [23] rvest_1.0.0 rappdirs_0.3.3
## [25] haven_2.4.1 xfun_0.24
## [27] crayon_1.4.1 RCurl_1.98-1.3
## [29] jsonlite_1.7.2 genefilter_1.72.1
## [31] survival_3.2-11 glue_1.4.2
## [33] gtable_0.3.0 zlibbioc_1.36.0
## [35] XVector_0.30.0 DelayedArray_0.16.3
## [37] proj4_1.0-10.1 Rttf2pt1_1.3.8
## [39] maps_3.3.0 scales_1.1.1
## [41] mvtnorm_1.1-2 DBI_1.1.1
## [43] Rcpp_1.0.7 emdbook_1.3.12
## [45] xtable_1.8-4 progress_1.2.2
## [47] bit_4.0.4 httr_1.4.2
## [49] RColorBrewer_1.1-2 ellipsis_0.3.2
## [51] pkgconfig_2.0.3 XML_3.99-0.6
## [53] farver_2.1.0 sass_0.4.0
## [55] locfit_1.5-9.4 utf8_1.2.1
## [57] labeling_0.4.2 tidyselect_1.1.1
## [59] rlang_0.4.11 later_1.2.0
## [61] munsell_0.5.0 BiocVersion_3.12.0
## [63] cellranger_1.1.0 tools_4.0.3
## [65] cachem_1.0.5 cli_3.0.0
## [67] generics_0.1.0 RSQLite_2.2.7
## [69] broom_0.7.8 evaluate_0.14
## [71] fastmap_1.1.0 yaml_2.2.1
## [73] knitr_1.33 bit64_4.0.5
## [75] fs_1.5.0 mime_0.11
## [77] ash_1.0-15 ggrastr_0.2.3
## [79] xml2_1.3.2 biomaRt_2.46.3
## [81] compiler_4.0.3 rstudioapi_0.13
## [83] beeswarm_0.4.0 curl_4.3.2
## [85] interactiveDisplayBase_1.28.0 reprex_2.0.0
## [87] geneplotter_1.68.0 bslib_0.2.5.1
## [89] stringi_1.6.2 highr_0.9
## [91] ggalt_0.4.0 lattice_0.20-44
## [93] ProtGenerics_1.22.0 Matrix_1.3-4
## [95] vctrs_0.3.8 pillar_1.6.1
## [97] lifecycle_1.0.0 BiocManager_1.30.16
## [99] jquerylib_0.1.4 bitops_1.0-7
## [101] httpuv_1.6.1 rtracklayer_1.50.0
## [103] R6_2.5.0 promises_1.2.0.1
## [105] KernSmooth_2.23-20 vipor_0.4.5
## [107] MASS_7.3-54 assertthat_0.2.1
## [109] openssl_1.4.4 withr_2.4.2
## [111] GenomicAlignments_1.26.0 Rsamtools_2.6.0
## [113] GenomeInfoDbData_1.2.4 hms_1.1.0
## [115] grid_4.0.3 coda_0.19-4
## [117] rmarkdown_2.9 bbmle_1.0.23.1
## [119] numDeriv_2016.8-1.1 shiny_1.6.0
## [121] lubridate_1.7.10 ggbeeswarm_0.6.0