热图

library(dplyr)
library(pheatmap)

绘制热图的数据是做完log10(tpm + 1)后的基因表达矩阵，不能用raw counts。一般基因总数可能有两万左右，不适合全部展示在热图中，需要挑选上调和下调最多的基因。

首先读取基因表达量信息，如果是tpm要计算log10(tpm+1)，如果是raw counts要先计算tpm。一般不用fpkm。

# read expression matrix ----
f <- file.path('_.csv')
gene_counts <- readr::read_csv(f)
glimpse(gene_expr)

## calculate log10(tpm + 1) ----
gene_tpm <- gene_counts |>
  dplyr::select(gene_symbol, starts_with('tpm')) |> 
  tibble::column_to_rownames(var = 'gene_symbol')
log_tpm_1p <- log10(gene_tpm + 1)

再读取差异表达基因，选择上调、下调倍数最高的各15个进行绘图。

## read DEG file ---- 
f <- file.path('_.csv')
gene_diff_sig <- readr::read_csv(f)

## select genes for plot ----
gene_plot <- gene_diff_sig |>
  dplyr::filter(signif != 'ns') |> 
  dplyr::mutate(log2FC_abs = abs(log2FC)) |> 
  dplyr::group_by(signif) |> 
  slice_max(log2FC_abs, n = 15) |> 
  ungroup() |> 
  dplyr::arrange(log2FC)
glimpse(gene_plot)

heatmap_data <- log_tpm_1p[gene_plot$gene_symbol, ]

此时得到的heatmap_data变量是行名为基因名、列名为样本名的data.frame。 接下来可以自行按行做scale运算，也可以用pheatmap()函数自带的scale功能。 另外，pheatmap()生成的图片不属于ggplot2体系，无法用ggsave()保存，用自带的filename参数设置保存路径。

# scale by row
heatmap_data <- heatmap_data |> t() |> scale() |> t()

## set color palette ----
color_palette <- colorRampPalette(c("dodgerblue4", "peachpuff", "deeppink4"))(50)

## plot heatmap ----
pheatmap(
  heatmap_data,
  scale = 'none', # don't scale again
  color = color_palette,
  cluster_rows = FALSE,
  cluster_cols = FALSE,
  # show_rownames = FALSE,
  # labels_row = data[, 1 , drop = T],
  display_numbers = F,
  border_color = 'white',
  number_color = 'white',
  # filename = '_heatmap.pdf'
)