RNA-seq-analysis-using-R/07_plot_key_gene_heatmap.R at main · Bioinformaticslave/RNA-seq-analysis-using-R · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
###############################
# 07_plot_key_gene_heatmap.R
# Heatmap of key immunological DEGs
###############################

source("R/00_config.R")

suppressPackageStartupMessages({
  library(gplots)
  library(RColorBrewer)
  library(DESeq2)
})

vsd <- readRDS(file.path(object_dir, "vsd.rds"))
result_all <- read.csv(file.path(table_dir, paste0(comparison_name, "_all_genes.csv")))

key_genes <- read.table(key_gene_file, sep = ",", stringsAsFactors = FALSE)
key_genes <- unique(as.character(key_genes[, 1]))

comparison_genes <- result_all$gene[
  result_all$padj < padj_cutoff &
    abs(result_all$log2FoldChange) > lfc_cutoff
]

key_genes <- unique(intersect(comparison_genes, key_genes))

exprs_data <- assays(vsd)[[1]]

selected_groups <- c(reference_group, target_group)

selected_samples <- rownames(colData(vsd))[
  colData(vsd)$group %in% selected_groups
]

exprs_data <- exprs_data[, selected_samples]

col_group <- brewer.pal(10, "Set3")[
  as.factor(colData(vsd)[selected_samples, ]$group)
]

mat <- exprs_data[
  intersect(key_genes, rownames(exprs_data)),
  ,
  drop = FALSE
]

pdf(file.path(plot_dir, paste0("05_", comparison_name, "_key_gene_heatmap.pdf")),
    width = 9,
    height = 10)

heatmap.2(
  as.matrix(mat),
  col = colorRampPalette(rev(brewer.pal(10, "RdYlBu")))(256),
  trace = "none",
  srtCol = 60,
  cexCol = 0.7,
  cexRow = 0.6,
  margins = c(18, 12),
  density.info = "none",
  lhei = c(2, 8),
  keysize = 0.9,
  main = "Key immunological DEGs",
  scale = "row",
  ColSideColors = col_group
)

dev.off()

cat("Key gene heatmap saved\\n")