GEO代码分析流程 - 5. 差异分析 - 可视化（作图） - 火山图、差异基因热图

0. 准备

?

1 2 3 4

setwd("D:/R/CHOL") rm(list = ls()) load(file = "step1output.Rdata") load(file = "step4output.Rdata")

---

1. 火山图

画不带基因标签的火山图

?

1 2 3 4 5 6 7 8 9 10 11

library(dplyr) library(ggplot2) dat = deg p <- ggplot(data = dat,aes(x = logFC,y = -log10(P.Value))) +   geom_point(alpha = 0.4,size = 3.5,aes(color = change)) +   ylab("-log10(Pvalue)")+   scale_color_manual(values = c("blue", "grey","red"))+   geom_vline(xintercept = c(-logFC_t,logFC_t),lty = 4,col = "black",lwd = 0.8) +   geom_hline(yintercept = -log10(P.Value_t),lty = 4,col = "black",lwd = 0.8) +   theme_bw() p

3. 把deg另存为dat对象，防止后续操作对deg中的数据造成修改。

4. 设置横、纵坐标

5. 画点图，设置点的透明度和大小，以及点的颜色分组映射

6. 添加y轴标签

7. 手动设置点的分组颜色

8. 横参考线：x轴截距，形状，颜色，粗细

9. 竖参考线：y轴截距，形状，颜色，粗细

10. 去掉灰色背景。如果要同时去掉灰色背景和网格线，则代码为+ggplot2::theme_classic()

若要固定图形比例，代码为+ggplot2::coord_fixed()

在火山图中添加基因标签

选择基因的三种方法

1. 自选基因：从deg中按基因名称（symbol）筛选指定行

?

1 2	for_label <- dat %>%   filter(symbol %in% c("TRPM3","SFRP1"))

2. 从deg中选择p值最小的10个基因

?

1 2	for_label <- dat %>%   head(10)

3. 从deg中选择p值最小的前3下调和前3上调

?

1 2 3 4 5 6 7

x1 = dat %>%   filter(change == "up") %>%   head(3) x2 = dat %>%   filter(change == "down") %>%   head(3) for_label <- rbind(x1,x2)

2. p值最小的前3个上调基因

5. p值最小的前3个下调基因

7. 合并上述基因

为选择的基因添加基因标签

?

1 2 3 4

volcano_plot <- p +   geom_point(size = 3, shape = 1, data = for_label) +   ggrepel::geom_label_repel(aes(label = symbol),data = for_label,color="black") volcano_plot

2. 给选择的基因点加圈，设置圈的大小，形状

3. 给选择的基因点加标签：把symbol对应的列作为标签加到图中，设置数据来源，设置标签颜色（包括文字和方框颜色）

保存火山图

?

1	ggsave(plot = volcano_plot,filename = paste0(gse,"volcano.png"))

---

 2. 差异基因热图

?

1	load(file = 'step2output.Rdata')

重新加载数据，防止数据被之前的操作改变

挑选差异基因的两种方法

1. 选择全部差异基因

?

1 2	cg = deg$probe_id[deg$change != "stable"] length(cg)

1. 挑选所有差异基因（上调和下调基因）对应的探针id

2. 差异基因数目

2. 取前30上调和前30下调基因

?

1 2 3 4

x = deg$logFC[deg$change != "stable"] names(x) = deg$probe_id[deg$change != "stable"] cg = c(names(head(sort(x),30)),names(tail(sort(x),30))) length(cg)

1. 获取所有差异基因（上调和下调基因）的logFC

2. 给所有差异基因的logFC取名为相应的探针id

3. 取前30上调和前30下调基因的名字（探针id）

4. 差异基因数目

获取差异基因

?

1 2	n = exp[cg,] dim(n)

1. 在deg数据框中取上述差异基因的子集

2. 查看子集的行数和列数

画热图

?

1 2 3 4 5 6 7 8

library(pheatmap) annotation_col = data.frame(group = group_list) rownames(annotation_col) = colnames(n) library(ggplotify) heatmap_plot <- as.ggplot(pheatmap(n,show_colnames = F,show_rownames = F,                           scale = "row",cluster_cols = F,                           annotation_col = annotation_col)) heatmap_plot

2. 按group_list设置差异基因分组

3. 设置差异基因名称

5. 不显示列名；不显示行名（as.ggplot函数用来同化pheatmap为ggplot图的类型）

6. 按行标度；不显示样本聚类（可修改为显示）

7. 设置差异基因分组

保存热图

?

1	ggsave(heatmap_plot,filename = paste0(gse,"heatmap.png"))

---

3. 合并图片

?

1 2 3

load("pca_plot.Rdata") library(patchwork) (pca_plot + volcano_plot + heatmap_plot) + plot_annotation(tag_levels = "A")

3. 合并图片，加标签

---