--- a
+++ b/b_DownstreamAnalysisScript/bulkRNAana_2a_DESeq.R
@@ -0,0 +1,63 @@
+
+# MESSAGE -----------------------------------------------------------------
+#
+# author: Yulin Lyu
+# email: lvyulin@pku.edu.cn
+#
+# require: R < 4.0 (recommended)
+#
+# note: DESeq is only for visualization usage in my protocols.
+#       It was no longer available since bioconductor 3.11 for R4.
+#       The last version of DESeq was built on R3.
+#       So this package may not work on R4.
+#
+# ---
+
+# 1. Load packages --------------------------------------------------------
+
+setwd("exampleData/RNA")
+
+# grammar
+library(tidyverse)
+library(magrittr)
+library(glue)
+
+# analysis
+library(DESeq)
+
+# for more information, please refer to:
+# http://bioconductor.org/packages/3.10/bioc/vignettes/DESeq/inst/doc/DESeq.pdf
+
+# 2. Load data ------------------------------------------------------------
+
+dataMtx <- readRDS("mid/dataMtx.rds")
+
+colnames(dataMtx)
+usedMtx <- dataMtx[c(
+  c("wanted sample names"),
+  NULL
+)]
+
+# 3. Analyze --------------------------------------------------------------
+
+dir.create("DESeq")
+
+condition <- factor(
+  c("a", "b", "c"),
+  levels = c("a", "b", "c")
+)
+
+cds <- newCountDataSet(usedMtx, condition) %>%
+  estimateSizeFactors() %>%
+  estimateDispersions(method = "blind", sharingMode = "fit-only")
+
+# save vsd profile for visualization later
+vsd <- varianceStabilizingTransformation(cds)
+saveRDS(vsd, "DESeq/vsd.rds")
+
+# perform DEG test (BETTER DO NOT, only if there is no other choice)
+# DEGres <- list()
+# DEGres$a_vs_b <- nbinomTest(cds, "a", "b")
+# DEGres$a_vs_c <- nbinomTest(cds, "a", "c")
+
+# iwalk(DEGres, ~ write_csv(.x, glue("DESeq/{.y}.DEG.csv")))