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

0. 准备

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

---

1. 火山图

画不带基因标签的火山图

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）筛选指定行

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

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

for_label <- dat %>% 
  head(10)

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

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. 合并上述基因

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

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对应的列作为标签加到图中，设置数据来源，设置标签颜色（包括文字和方框颜色）

保存火山图

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

---

 2. 差异基因热图

load(file = 'step2output.Rdata')

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

挑选差异基因的两种方法

1. 选择全部差异基因

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

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

2. 差异基因数目

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

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. 差异基因数目

获取差异基因

n = exp[cg,]
dim(n)

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

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

画热图

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. 设置差异基因分组

保存热图

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

---

3. 合并图片

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

3. 合并图片，加标签

---