library(dplyr)
library(pheatmap)
热图
绘制热图的数据是做完log10(tpm + 1)
后的基因表达矩阵,不能用raw counts。一般基因总数可能有两万左右,不适合全部展示在热图中,需要挑选上调和下调最多的基因。
首先读取基因表达量信息,如果是tpm要计算log10(tpm+1),如果是raw counts要先计算tpm。一般不用fpkm。
# read expression matrix ----
<- file.path('_.csv')
f <- readr::read_csv(f)
gene_counts glimpse(gene_expr)
## calculate log10(tpm + 1) ----
<- gene_counts |>
gene_tpm ::select(gene_symbol, starts_with('tpm')) |>
dplyr::column_to_rownames(var = 'gene_symbol')
tibble<- log10(gene_tpm + 1) log_tpm_1p
再读取差异表达基因,选择上调、下调倍数最高的各15个进行绘图。
## read DEG file ----
<- file.path('_.csv')
f <- readr::read_csv(f)
gene_diff_sig
## select genes for plot ----
<- gene_diff_sig |>
gene_plot ::filter(signif != 'ns') |>
dplyr::mutate(log2FC_abs = abs(log2FC)) |>
dplyr::group_by(signif) |>
dplyrslice_max(log2FC_abs, n = 15) |>
ungroup() |>
::arrange(log2FC)
dplyrglimpse(gene_plot)
<- log_tpm_1p[gene_plot$gene_symbol, ] heatmap_data
此时得到的heatmap_data
变量是行名为基因名、列名为样本名的data.frame
。 接下来可以自行按行做scale运算,也可以用pheatmap()
函数自带的scale功能。 另外,pheatmap()
生成的图片不属于ggplot2
体系,无法用ggsave()
保存,用自带的filename
参数设置保存路径。
# scale by row
<- heatmap_data |> t() |> scale() |> t()
heatmap_data
## set color palette ----
<- colorRampPalette(c("dodgerblue4", "peachpuff", "deeppink4"))(50)
color_palette
## 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'
)