diff --git a/bulkRNAseq/Fig1F.R b/bulkRNAseq/Fig1F.R
new file mode 100644
index 0000000000000000000000000000000000000000..a3a6f78aa26c792999c847a6537bb41f8c926fb3
--- /dev/null
+++ b/bulkRNAseq/Fig1F.R
@@ -0,0 +1,63 @@
+library(DESeq2)
+library(magrittr)
+library(tibble)
+library(dplyr)
+library(tidyr)
+library(knitr)
+library(ggplot2)
+library(SummarizedExperiment)
+library(DEGreport)
+library(ggrepel)
+library(EnhancedVolcano)
+library(openxlsx)
+library(RColorBrewer)
+library(org.Hs.eg.db)
+library(pheatmap)
+library(clusterProfiler)
+library(ComplexHeatmap)
+library(DOSE)
+
+## using FDR 0.05 ##
+logfc.thr <- 0
+adj.pval.thr <- 0.05
+
+# Volcano plot
+obj <- readRDS(file = "../20191109_120832_DGE.rds")
+DGE.results <- as.data.frame(obj$deseq$dge_res$`HSCMPP_THAL-HSCMPP_ND`)
+colnames(DGE.results) <- c("log2FoldChange","lfcSE","baseMean","PValue","padj")
+
+DGE.results <- readRDS("Fig1F_input.rds") # DGE output with padj < 0.05
+data <- data.frame(gene = row.names(DGE.results), pval = -log10(DGE.results$padj), lfc = DGE.results$log2FoldChange)
+data <- na.omit(data)
+data <- mutate(data, color = case_when(data$lfc > logfc.thr & data$pval > -log10(adj.pval.thr) ~ "Overexpressed",
+                                       data$lfc < -logfc.thr & data$pval > -log10(adj.pval.thr) ~ "Underexpressed",
+                                       abs(data$lfc) < logfc.thr | data$pval < -log10(adj.pval.thr) ~ "NonSignificant"))
+data <- data[order(data$pval, decreasing = TRUE),]
+
+selected_genes <-read.table(file = "Selected_genes", header = TRUE)
+rownames(selected_genes) <- selected_genes$GENE
+
+highl2 <- subset(data, subset = gene %in% selected_genes$GENE, color != "NonSignificant")
+
+vol <- ggplot(data, aes(x = lfc, y = pval, colour = color)) +
+  theme_bw(base_size = 12) +
+  theme(legend.position = "right") +
+  xlab("log2 Fold Change") + 
+  ylab(expression(-log[10]("adjusted p-value"))) +
+  geom_hline(yintercept = -log10(adj.pval.thr), colour = "darkgray")
+if (logfc.thr > 0) {
+  vol <- vol +
+    geom_vline(xintercept = logfc.thr, colour = "darkgray") +
+    geom_vline(xintercept = -logfc.thr, colour = "darkgray")
+}
+vol <- vol +
+  geom_point(size = 2, alpha = 0.8, na.rm = T) +
+  scale_color_manual(name = "Expression",
+                     values = c(Overexpressed = "#CD4F39",
+                                Underexpressed = "#0B5394",
+                                NonSignificant = "darkgray")) +
+  geom_label_repel(data = highl2, aes(label = gene), show.legend = FALSE, max.overlaps = 20)
+ggsave(filename = "Fig1F.png",
+       plot = vol,
+       width = 9, height = 7, dpi = 600)
+
diff --git a/bulkRNAseq/Fig1F_input.rds b/bulkRNAseq/Fig1F_input.rds
new file mode 100644
index 0000000000000000000000000000000000000000..8a4f0c2504996f2d683ca561d34bc7ad5fffeed2
Binary files /dev/null and b/bulkRNAseq/Fig1F_input.rds differ
diff --git a/bulkRNAseq/Selected_genes b/bulkRNAseq/Selected_genes
new file mode 100644
index 0000000000000000000000000000000000000000..a637385cb22f55c38f5af6d05eb4d51d245df81f
--- /dev/null
+++ b/bulkRNAseq/Selected_genes
@@ -0,0 +1,20 @@
+GENE
+NFIB
+HBD
+GYPE
+DLK1
+APOE
+CDKN1A
+MAP1LC3A
+SKIL
+E2F2
+ATG4D
+SQSTM1
+NFKBIE
+RELA
+NFKB2
+TGFB1
+JUN
+ATG14
+TFRC
+NFKB1
\ No newline at end of file