load("../1.align/txi.rsem.RData")
metadata <- read.csv("metadata.csv", stringsAsFactors = TRUE)
metadata$Trt <- factor(metadata$Trt, levels=c("SAL","LPS"))
rownames(metadata) <- metadata$ID
all.equal(colnames(txi.rsem$counts), rownames(metadata))
## [1] TRUE
Differential gene expression for design: Group + Drug + Trt
dds <- DESeqDataSetFromTximport(txi.rsem, colData = metadata, design = ~ Group + Drug + Trt)
# filter low expressed genes
keep <- rowSums(counts(dds)) >= 50
table(keep)
## keep
## FALSE TRUE
## 38866 16621
dds <- dds[keep,]
dds <- DESeq(dds)
resultsNames(dds)
## [1] "Intercept" "Group_MIA_vs_CON" "Drug_PLX_vs_CON" "Trt_LPS_vs_SAL"
PCA
vsd <- vst(dds, blind=FALSE)
m <- plotPCA(vsd, intgroup=c("Group", "Drug"), returnData = TRUE)
ggplot(m, aes(PC1, PC2, color = group)) + geom_point(size = 3) + theme_cowplot() + xlab("PC1: 71% variance") + ylab("PC2: 7% variance") + geom_text_repel(aes(label = name)) + scale_color_manual(values = c("black", "turquoise3", "red", "purple2"))
Cluster
tree = hclust(dist(t(assay(vsd))), method = "average")
plot(tree, main = "hclust, average VST", sub="", xlab="", cex.lab = 1, cex.axis = 1, cex.main = 1)
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] ggfortify_0.4.12 ggrepel_0.9.1
## [3] AnnotationHub_2.22.1 BiocFileCache_1.14.0
## [5] dbplyr_2.1.1 org.Mm.eg.db_3.12.0
## [7] AnnotationDbi_1.52.0 DESeq2_1.30.1
## [9] SummarizedExperiment_1.20.0 Biobase_2.50.0
## [11] MatrixGenerics_1.2.1 matrixStats_0.59.0
## [13] GenomicRanges_1.42.0 GenomeInfoDb_1.26.7
## [15] IRanges_2.24.1 S4Vectors_0.28.1
## [17] BiocGenerics_0.36.1 cowplot_1.1.1
## [19] forcats_0.5.1 stringr_1.4.0
## [21] dplyr_1.0.7 purrr_0.3.4
## [23] readr_1.4.0 tidyr_1.1.3
## [25] tibble_3.1.2 ggplot2_3.3.5
## [27] tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-2 ellipsis_0.3.2
## [3] XVector_0.30.0 fs_1.5.0
## [5] rstudioapi_0.13 farver_2.1.0
## [7] bit64_4.0.5 interactiveDisplayBase_1.28.0
## [9] fansi_0.5.0 lubridate_1.7.10
## [11] xml2_1.3.2 splines_4.0.3
## [13] cachem_1.0.5 geneplotter_1.68.0
## [15] knitr_1.33 jsonlite_1.7.2
## [17] broom_0.7.8 annotate_1.68.0
## [19] shiny_1.6.0 BiocManager_1.30.16
## [21] compiler_4.0.3 httr_1.4.2
## [23] backports_1.2.1 assertthat_0.2.1
## [25] Matrix_1.3-4 fastmap_1.1.0
## [27] cli_3.0.0 later_1.2.0
## [29] htmltools_0.5.1.1 tools_4.0.3
## [31] gtable_0.3.0 glue_1.4.2
## [33] GenomeInfoDbData_1.2.4 rappdirs_0.3.3
## [35] Rcpp_1.0.7 cellranger_1.1.0
## [37] jquerylib_0.1.4 vctrs_0.3.8
## [39] xfun_0.24 rvest_1.0.0
## [41] mime_0.11 lifecycle_1.0.0
## [43] XML_3.99-0.6 zlibbioc_1.36.0
## [45] scales_1.1.1 hms_1.1.0
## [47] promises_1.2.0.1 RColorBrewer_1.1-2
## [49] yaml_2.2.1 curl_4.3.2
## [51] gridExtra_2.3 memoise_2.0.0
## [53] sass_0.4.0 stringi_1.6.2
## [55] RSQLite_2.2.7 highr_0.9
## [57] BiocVersion_3.12.0 genefilter_1.72.1
## [59] BiocParallel_1.24.1 rlang_0.4.11
## [61] pkgconfig_2.0.3 bitops_1.0-7
## [63] evaluate_0.14 lattice_0.20-44
## [65] labeling_0.4.2 bit_4.0.4
## [67] tidyselect_1.1.1 magrittr_2.0.1
## [69] R6_2.5.0 generics_0.1.0
## [71] DelayedArray_0.16.3 DBI_1.1.1
## [73] pillar_1.6.1 haven_2.4.1
## [75] withr_2.4.2 survival_3.2-11
## [77] RCurl_1.98-1.3 modelr_0.1.8
## [79] crayon_1.4.1 utf8_1.2.1
## [81] rmarkdown_2.9 locfit_1.5-9.4
## [83] grid_4.0.3 readxl_1.3.1
## [85] blob_1.2.1 reprex_2.0.0
## [87] digest_0.6.27 xtable_1.8-4
## [89] httpuv_1.6.1 munsell_0.5.0
## [91] bslib_0.2.5.1