使用ComplexHeatmap在热图右侧绘制箱型图

wo_qq111

已于 2023-07-26 16:51:42 修改

阅读量235

点赞数 2

文章标签：学习笔记

于 2023-07-26 16:49:27 首次发布

本文链接：https://blog.csdn.net/weixin_52648333/article/details/131905069

版权

引言

对于单细胞数据，热图可以展示细胞层面的变化趋势，但是却不能定量这种趋势，我在数据分析的过程中需要同时展示一个细胞类型实验组和对照组中细胞的基因表达量水平，为了表明这个基因是否存在显著差异，我在每一行热图的右侧同时绘制实验组和对照组基因表达量的箱型图，并标注了p值。

数据

为了方便重复，我尝试使用ifnb的数据集，但是我没有下载ifnb的seurat对象（不敢InstallData("ifnb")，因为会装到家目录，而家目录已经没空间了），用的是之前不知道从哪里搞来的10x数据，细胞注释那块复现不出来。

之后修改安装路径，在项目路径下添加两个文件，一个.Renviron一个.Rprofile，前者写入R_LIBS_SITE=" "，后者写入.libPaths()里添加所有路径（先打.libPaths()查看已有安装路径，然后.libPaths(c(‘your_path’, 'orig_path1', 'orig_path2')），重启环境。

然后本地安装ifnb.SeuratData包（也就是InstallData("ifnb")做的事），参考安装Seurat内置数据集ifnb.SeuratData的方法（Linux上安装）（SeuratData包比较小，可以直接从GitHub安装），之后参考帮助文档Introduction to scRNA-seq integration，做前期处理，这个标题下的代码块都写在一个文件中，代码如下：

首先导入包

library(Seurat)
library(SeuratData)
set.seed(2746)

导入ifnb数据

# InstallData("ifnb")
# 因本地安装好了，不需要运行这条命令

LoadData("ifnb")
# 加载一下数据

merged_list <- SplitObject(
  ifnb,
  split.by = "stim"
)
# 外周血的数据，分为两组，一组对照，一组干扰素刺激的

接着归一化和寻找高变基因

merged_list <- lapply(
  X = merged_list,
  FUN = function(x) {
    x <- NormalizeData(x)
    x <- FindVariableFeatures(
      x,
      selection.method = "vst",
      nfeatures = 2000
    )
  }
)

选择不同样本共有的特征基因

features <- SelectIntegrationFeatures(
  object.list = merged_list
)

整合不同样本

anchors <- FindIntegrationAnchors(
  object.list = merged_list,
  anchor.features = features
)

merged <- IntegrateData(
  anchorset = anchors
)

saveRDS(merged, "output1/seurat/merged.rds")

使用seurat标准流程

merged <- ScaleData(merged, verbose = FALSE)
merged <- RunPCA(merged, npcs = 50, verbose = FALSE)
DimPlot(merged, reduction = "pca")
ElbowPlot(merged, ndims = 50)
dims <- 30

merged <- RunUMAP(
  merged,
  reduction = "pca",
  dims = 1:dims
)
DimPlot(merged, reduction = "umap")

直接使用ifnb注释好的数据

DimPlot(
  merged,
  reduction = 'umap',
  group.by = 'seurat_annotations',
  label = T
)

结果如下：在这里插入图片描述
seurat对象保存，assay切换回RNA，因为integrate数据只是用来降维聚类的，RNA用来数据分析（差异、打分等）

DefaultAssay(merged) <- "RNA"
Idents(merged) <- 'seurat_annotations'
saveRDS(merged, 'output1/seurat/seurat_obj.rds')

热图绘制

接着新建一个文件，这个大标题下所有子标题的代码都在这个文件中，开始热图的代码：

导入包

library(Seurat)
library(ComplexHeatmap)
library(tidyverse)
library(RColorBrewer)
library(ggpubr)
library(circlize)
set.seed(2746)

热图数据

首先先筛选热图制作用的数据，但是因为我对免疫细胞不了解，不知道对这个数据集的关注点在哪里，我计划绘制B细胞差异基因中在干扰素刺激处理上调和下调logfc绝对值大于0.5的基因，并在热图右侧绘制绘制干扰素刺激和对照组基因表达量的箱型图，并标注p值。

加载seurat对象

seurat_obj <- readRDS('output1/seurat/seurat_obj.rds')

table查看细胞数量

> table(seurat_obj$seurat_annotations)

           B  B Activated    CD14 Mono    CD16 Mono CD4 Memory T 
         978          388         4362         1044         1762 
 CD4 Naive T        CD8 T           DC        Eryth           Mk 
        2504          814          472           55          236 
          NK          pDC  T activated 
         619          132          633

提取B细胞子集对象和矩阵，差异分析

sub_seurat <- subset(
  seurat_obj,
  idents = 'B'
)
sub_matrix <- GetAssayData(
  sub_seurat,
  assay = 'RNA'
) |>
  as.matrix()

Idents(sub_seurat) <- 'stim'
table(sub_seurat$stim)
deg <- FindMarkers(
  sub_seurat,
  ident.1 = 'STIM',
  ident.2 = 'CTRL'
)

筛选p值小于等于0.05且avg_log2FC绝对值大于等于0.5的差异基因，按avg_log2FC降序排列，将行名转为单独一列，并添加上下调信息

select_gene_info <- deg |>
  filter(p_val <= 0.05 & abs(avg_log2FC) >= 0.5) |>
  arrange(desc(avg_log2FC)) |>
  rownames_to_column('name') |>
  mutate(
    info = ifelse(
      avg_log2FC >= 0.5,
      'UP', 'DOWN'
    )
  )

添加细胞注释，一列为细胞id，一列为分组信息，并将细胞按照实验组、对照组排序

cell_info <- data.frame(
  group = as.character(sub_seurat$stim),
  name = names(sub_seurat$stim)
) |>
  arrange(
    factor(group, levels = c('STIM', 'CTRL'))
  )

重新排序矩阵，按照基因和细胞注释中的顺序

sub_matrix <- sub_matrix[
  select_gene_info$name,
  cell_info$name
]

热图注释

行注释，格式为键值对+颜色列表。
键值对可以有多个，每个键值对必须在颜色列表里有对应。
值可以是vector、matrix或注释函数，这里是vector，因为vector里的顺序不会根据热图矩阵中行或列元素的顺序自动调整以一一对应，所以在上一步调整矩阵元素顺序以保证和注释元素顺序一致。
然后就是颜色列表，如果值是离散的，那么就把每个离散的值分配一个颜色即可（格式 named vector），如果是连续的，就用colorRamp2函数，第一个参数是值的所有极值，第二个参数是极值的颜色，所有点的颜色会均匀变化。

row_ha <- HeatmapAnnotation(
  avg_log2FC = select_gene_info$avg_log2FC,
  UP_or_Down = select_gene_info$info,
  col = list(
    avg_log2FC = colorRamp2(
      c(min(select_gene_info$avg_log2FC), -0.5, 0.5, max(select_gene_info$avg_log2FC)),
      c('#A9D5A9', 'white', 'white', '#F8BD7D')
    ),
    UP_or_Down = structure(
      names = c('UP', 'DOWN'),
      c('#F8BD7D', '#A9D5A9')
    )
  ),
  which = 'row'
)

列注释同理

col_ha <- HeatmapAnnotation(
  group = cell_info$group,
  col = list(group = structure(
    names = c('STIM', 'CTRL'),
    c('#E98E29', '#9F9EA3')
  ))
)

箱型图

箱型图通过注释函数实现，参考Chapter 3 Heatmap Annotations | ComplexHeatmap Complete Reference (jokergoo.github.io)
首先准备分组信息，用于自定义函数的输入（只要返回对象是AnnotationFunction生成的即可，所以外面可以再包一个函数），格式为named vector，前两个是对照组和实验组在seurat对象注释中的名称，后两个是对照组和实验组箱型图使用的颜色。

group_info <- structure(
  names = c('ctrl_name', 'expr_name', 'ctrl_color', 'expr_color'),
  c('CTRL', 'STIM', '#9F9EA3', '#E98E29')
)

函数由p值计算和注释函数组成。
p值计算apply函数每行切片，拼接分组信息，然后计算组间的pvalue。

注释函数使用的参数有：cell_fun、var_import 、which 、width ，相比于原本的fun参数，cell_fun可以控制每一行（列）里的绘图情况，而cell_fun用到的外部变量必须由var_import传入，而cell_fun本身只是接受index参数。

cell_fun由pushViewport、绘图函数和popViewport组成，pushViewport里设置了一个格子里x，y坐标轴的范围（以左下角为坐标原点，因为是cell_fun，所以x，y描述的是一行（列）对应绘图的范围，且这个范围指的是x、y轴的刻度范围，而不是绘图区域的绝对大小，详细参数可以在命令行打?viewport）。

绘图函数有grid.rect()、grid.boxplot()、grid.text()、grid.xaxis()

grid.rect()绘制外框，为了围住属于同一行的箱型图
grid.boxplot()绘制箱型图，pos设置箱型图的相对位置，我发现pos的值与坐标轴相对应，如值为1表示在y=1的位置绘制一个箱型图
grid.text()绘制p值，参数x、y表示在这个框中的相对位置
grid.xaxis()绘制x坐标轴，但因为在cell_fun中，为了不每一行都画一个x坐标轴，判断只有到最下面一行才绘制

anno_2box <- function(score_matrix, col_anno_df, group_info){
  # score_matrix is a matrix which colnames is cellnames and rownames
  # is genenames or other valuenames like gsva score
  
  # col_anno_df should be a dataframe contain columns of 'name' which 
  # indicate cell name and 'group' indicate two group of expr and ctrl
  
  # group_info is a named charactor vector which contained two group 
  # names of 'expr_name' and 'ctrl_name' and their color names which
  # are 'expr_color' and 'ctrl_color'
  
  ctrl_group = filter(col_anno_df, group == group_info['ctrl_name'])$name
  expr_group = filter(col_anno_df, group == group_info['expr_name'])$name
  expr_color = group_info['expr_color']
  ctrl_color = group_info['ctrl_color']
  
  
  pval = apply(score_matrix, 1, function(x){
    # View(x)
    # 是一个带name的double
    tmp = data.frame(
      value = x,
      name = names(x)
    )
    tmp = merge(
      tmp, col_anno_df,
      by.x = 'name',
      by.y = 'name'
    )
    compare_means(
      value ~ group,
      data = tmp
    )
  })
  
  AnnotationFunction(
    cell_fun = function(index){
      pushViewport(viewport(
        xscale = c(min(score_matrix), max(score_matrix)*1.1),
        yscale = c(0.5, 2.5)
      ))
      grid.rect()
      grid.boxplot(
        score_matrix[index, ctrl_group],
        pos = 2,
        direction = "horizontal",
        box_width = 0.65,
        size = unit(1, "pt"),
        gp = gpar(fill = ctrl_color)
      )
      grid.text(
        pval[[index]]$p.signif,
        x = 0.95,
        y = 0.5,
        just = c('center', 'centre'),
        rot = -90,
        gp = gpar(fontsize = 5, fontface = 'bold')
      )
      grid.boxplot(
        score_matrix[index, expr_group],
        pos = 1,
        direction = "horizontal",
        box_width = 0.65,
        size = unit(1, "pt"),
        gp = gpar(fill = expr_color)
      )
      if(index == nrow(score_matrix)) grid.xaxis()
      popViewport()
    },
    var_import = list(
      score_matrix = score_matrix,
      ctrl_group = ctrl_group,
      expr_group = expr_group,
      pval = pval,
      expr_color = expr_color,
      ctrl_color = ctrl_color
    ),
    which = 'row',
    width = unit(4, 'cm'),
  )
}

热图

热图这里添加箱型图会覆盖掉行的名称，所以把行名以注释的形式展示，此外参数 cluster_columns 和 cluster_columns 会重新聚类，打乱顺序。

p <- Heatmap(
  sub_matrix,
  column_title = 'B deg',
  name = 'normalized count',
  top_annotation = col_ha,
  left_annotation = row_ha,
  cluster_columns = F,
  cluster_rows = F,
  show_row_dend = F,
  show_column_dend = F,
  show_column_names = F,
  
  width = ncol(sub_matrix) * unit(1, "pt"),
  height = nrow(sub_matrix) * unit(16, "pt"),
  use_raster = F,
  
  col = colorRampPalette(rev(brewer.pal(10, "RdBu")))(30)
) +
  rowAnnotation(
    gene = anno_text(rownames(sub_matrix), gp = gpar(fontsize = 8)),
    Young_vs_Old = anno_2box(sub_matrix, cell_info, group_info),
    annotation_name_rot = 0,
    annotation_name_offset = unit(1, 'cm')
  )

保存时需要根据矩阵大小调整保存大小。

width <- ceiling(ncol(sub_matrix)*0.02) + 6
height <- ceiling(nrow(sub_matrix)*0.25) + 1
pdf('output1/heatmap/B.pdf', width = width, height = height)
print(p)
dev.off()

结果如下：
在这里插入图片描述

总结

这个绘图的代码主要涉及复杂注释的绘制，但是对于大数据量尤其是行数较多的矩阵绘制箱型图不能很直观的体现样本之间的差异（因为看不清了）。
我没有统计学知识，p检验那里不知道是否可以使用compare_means。
我参考了很多博客，如数据清理、pdf绘图中print的使用，但已找不到对应的具体博客了，只好感谢博主们的分享。
感谢生信实验室的师兄师姐提供的很多帮助，基本框架都是师兄师姐提供的。
在R语言学习以及计算机使用的过程中，我因开源和分享而获益良多，所以记录我这次的复杂热图注释学习过程，希望发扬分享的精神。
但是写一篇好累，而且我是个菜鸡，网上的友友们不要问我问题，我回答不了，我的分享精神也用光了😢。

wo_qq111

关注

2
点赞
踩
4

收藏

觉得还不错? 一键收藏
1
评论
使用ComplexHeatmap在热图右侧绘制箱型图

这个绘图的代码主要涉及复杂注释的绘制，但是对于大数据量尤其是行数较多的矩阵绘制箱型图不能很直观的体现样本之间的差异（因为看不清了）。我没有统计学知识，p检验那里不知道是否可以使用。在R语言学习以及计算机使用的过程中，我因开源和分享而获益良多，所以记录我这次的复杂热图注释学习过程，希望发扬分享的精神。但是写一篇好累，而且我是个菜鸡，网上的友友们不要问我问题，我回答不了，我的分享精神也用光了😢。
复制链接

扫一扫