查看: 642|回复: 15

[R语言] 比原文献的还好看!绘制添加基因标签的火山图

  [复制链接]

管理员

Rank: 15Rank: 15Rank: 15Rank: 15Rank: 15

主题
907
注册时间
2020.6.16
在线时间
485 小时

发表于 2022.11.14 09:59:11 | 显示全部楼层 |阅读模式
之前在《案例复现(一):表达数据的下载与差异分析》一文中已经介绍了如何从GEO数据库下载数据并进行基于limma的差异分析。完成差异分析后,文章作者首先使用火山图对差异分析结果进行可视化,如下图(A),共得到上调基因2496个,下调基因1130个。


那么,如何使用我们自己得到的差异分析数据,复现出原文的火山图呢?甚至绘制出更新颖更好看的火山图?下面继续为大家演示。

1. 读取数据

  1. #读取差异分析表格;
  2. exp <- read.csv("patient-control_diff_new.csv")
  3. #预览数据;
  4. head(exp)
复制代码


2. 数据过滤

  1. #去掉没有symbol号的探针记录;
  2. expr1 <- exp[exp$GENE_SYMBOL!='',]
  3. #查看筛选后的数据量;
  4. dim(expr1)
  5. #[1] 29833 6
  6. #一个基因对应多个探针,继续过滤;
  7. #查看基因个数;
  8. length(unique(expr1$GENE_SYMBOL))
  9. #[1] 21754
  10. library(dplyr)
  11. #根据平均表达量对表格进行降序排列;
  12. expr2 <- arrange(expr1,desc(AveExpr))
  13. #根据Gene symbol列移除包含重复的行;
  14. expr3<- expr2 %>% distinct(GENE_SYMBOL, .keep_all = TRUE)
  15. #查看表格的行列数;
  16. dim(expr3)
  17. #[1] 21754 6
  18. #预览数据;
  19. head(expr3)
复制代码


3. 绘制火山图

  1. #载入绘图相关的R包;
  2. library(ggplot2)
  3. library(ggrepel)
  4. #转成tibble便于后续使用,去掉不需要的列;
  5. dt <- as_tibble(expr3[c(6,4,3)])
  6. #对p值取对数;
  7. dt$log10PValue <- -log10(dt$P.Value)
  8. #生成显著上下调数据的分组标签;
  9. dt$group <- case_when(dt$logFC > 1 & dt$P.Value < 0.05 ~ "Up",
  10. dt$logFC < -1 & dt$P.Value < 0.05 ~ "Down",
  11. abs(dt$logFC) <= 1 ~ "None",
  12. dt$P.Value >= 0.05 ~ "None")
  13. head(dt)
复制代码


  1. #重新设置阈值(logFC=1.5),生成显著上下调数据的分组标签;
  2. dt$group2 <- case_when(dt$logFC > 1.5 & dt$adj.P.Val < 0.05 ~ "Up",
  3. dt$logFC < -1.5 & dt$adj.P.Val < 0.05 ~ "Down",
  4. abs(dt$logFC) <= 1.5 ~ "None",
  5. dt$adj.P.Val >= 0.05 ~ "None")
  6. #提取上下调基因;
  7. Up <- filter(dt,group2=="Up")
  8. up_genes <- Up$GENE_SYMBOL
  9. Down <- filter(dt,group2=="Down")
  10. down_genes <- Down$GENE_SYMBOL
  11. #确定上下调基因数量;
  12. paste0("The number of up gene is ",length(up_genes))
  13. #[1] "The number of up gene is 2483"
  14. paste0("The number of down gene is ",length(down_genes))
  15. #[1] "The number of down gene is 1165"
复制代码

注意,可能由于上游分析过程中过滤条件、差异分析参数等存在差异,我这里的分析结果与文献原文的结果有些许差异,原文的上调基因数为2496,下调基因数为1130。

  1. #获取表达差异最显著的10个基因;
  2. top10sig <- filter(dt,group!="None") %>% distinct(GENE_SYMBOL,.keep_all = T) %>% top_n(10,abs(logFC))
  3. top10sig
复制代码


  1. #将差异表达Top10的基因表格拆分成up和down两部分;
  2. up <- filter(top10sig,group=="Up")
  3. up
  4. down <- filter(top10sig,group=="Down")
  5. down
复制代码


  1. #新增一列,将Top10的差异基因标记为2,其他的标记为1;
  2. dt$size <- case_when(!(dt$GENE_SYMBOL %in% top10sig$GENE_SYMBOL)~ 1,
  3. dt$GENE_SYMBOL %in% top10sig$GENE_SYMBOL ~ 2)
  4. head(dt)
复制代码


  1. #提取非Top10的基因表格;
  2. df <- filter(dt,size==1)
  3. #指定绘图顺序,将group列转成因子型;
  4. df$group <- factor(df$group,
  5. levels = c("Up","Down","None"),
  6. ordered = T)
  7. #开始绘图,建立映射;
  8. p0 <-ggplot(data=df,aes(logFC,log10PValue,color=group))
  9. #添加散点;
  10. p1 <- p0+geom_point(size=1)
  11. p1
复制代码


  1. #自定义半透明颜色(红绿);
  2. mycolor <- c("#FF9999","#99CC00","gray80")
  3. p2 <- p1 + scale_colour_manual(name="",values=alpha(mycolor,0.9))
  4. p2
复制代码


  1. #继续添加Top10基因对应的点;
  2. p3 <- p2+geom_point(data=up,aes(logFC,log10PValue),
  3. color="#FF9999",size=2,alpha=0.9)+
  4. geom_point(data=down,aes(logFC,log10PValue),
  5. color="#7cae00",size=2,alpha=0.9)
  6. p3
复制代码


  1. #expansion函数设置坐标轴范围两端空白区域的大小;mult为“倍数”模式,add为“加性”模式;
  2. p4<-p3+labs(y="-log10 P Value",x="log2FC")+
  3. scale_y_continuous(expand=expansion(add = c(0.5, 0)),
  4. limits = c(0, 12),
  5. breaks = c(0,3,6,9,12),
  6. label = c("0","3","6","9","12"))+
  7. scale_x_continuous(limits = c(-12, 12),
  8. breaks = c(-8,-4,0,4,8),
  9. label = c("-8","-4","0","4","8"))
  10. p4
复制代码


  1. #添加箭头;
  2. set.seed(007)
  3. p5 <- p4+geom_text_repel(data=top10sig,aes(logFC,log10PValue,label=GENE_SYMBOL),
  4. force=20,color="grey20",
  5. size=2,
  6. point.padding = 0.5,hjust = 0.5,
  7. arrow = arrow(length = unit(0.01, "npc"),
  8. type = "open", ends = "last"),
  9. segment.color="grey20",
  10. segment.size=0.2,
  11. segment.alpha=0.8,
  12. nudge_x=0,
  13. nudge_y=0)
  14. p5
复制代码


  1. #自定义图表主题,对图表做精细调整;
  2. top.mar=0.2
  3. right.mar=0.2
  4. bottom.mar=0.2
  5. left.mar=0.2
  6. #隐藏纵轴,并对字体样式、坐标轴的粗细、颜色、刻度长度进行限定;
  7. mytheme<-theme_classic()+
  8. theme(text=element_text(family = "sans",colour ="gray30",size = 10),
  9. axis.line = element_line(size = 0.6,colour = "gray30"),
  10. axis.ticks = element_line(size = 0.6,colour = "gray30"),
  11. axis.ticks.length = unit(1.5,units = "mm"),
  12. plot.margin=unit(x=c(top.mar,right.mar,bottom.mar,left.mar),
  13. units="inches"))
  14. #应用自定义主题;
  15. p5+mytheme
复制代码


  1. #添加辅助线;
  2. p6 <- p4+geom_hline(yintercept = c(-log10(0.05)),
  3. size = 0.5,
  4. color = "orange",
  5. lty = "dashed")+
  6. geom_vline(xintercept = c(-1,1),
  7. size = 0.5,
  8. color = "orange",
  9. lty = "dashed")
  10. p6
复制代码


  1. #添加其他样式的标签;
  2. #为了方便自定义左右区域的标签,这里使用up、down两个独立的子表格;
  3. p7 <- p6+geom_label_repel(
  4. data = up,aes(logFC,log10PValue,label=GENE_SYMBOL),
  5. nudge_x = 3,
  6. nudge_y = -0.8,
  7. color = "white",
  8. alpha = 0.9,
  9. point.padding = 0.5,
  10. size = 2.5,
  11. fill = "#96C93D",
  12. segment.size = 0.5,
  13. segment.color = "grey50",
  14. direction = "y",
  15. hjust = 0.5) +
  16. geom_label_repel(
  17. data = down,aes(logFC,log10PValue,label=GENE_SYMBOL),
  18. nudge_x = -3,
  19. nudge_y = 0.2,
  20. color = "white",
  21. alpha = 0.9,
  22. point.padding = 0.5,
  23. size = 2.5,
  24. fill = "#9881F5",
  25. segment.size = 0.5,
  26. segment.color = "grey50",
  27. direction = "y",
  28. hjust = 0.5)
  29. #应用自定义主题;
  30. p8 <- p7+mytheme
  31. p8
复制代码


好啦,本次的绘图方法分享就到这里啦!

*未经许可,不得以任何方式复制或抄袭本篇文章之部分或全部内容。版权所有,侵权必究。

本文作者:基迪奥-莫北

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
0
注册时间
2020.11.21
在线时间
129 小时

发表于 2022.11.14 11:24:22 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
3
注册时间
2021.6.22
在线时间
63 小时

发表于 2022.11.14 21:20:59 | 显示全部楼层
学习
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.2.12
在线时间
131 小时

发表于 2022.11.15 08:52:01 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
0
注册时间
2020.11.21
在线时间
129 小时

发表于 2022.11.15 13:07:31 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2016.5.10
在线时间
34 小时

发表于 2022.11.15 20:01:24 | 显示全部楼层
11月的第五天
新的一天加油!
回复 支持 反对

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.2.12
在线时间
131 小时

发表于 2022.11.16 08:50:18 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
0
注册时间
2020.11.21
在线时间
129 小时

发表于 2022.11.16 13:58:02 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
0
注册时间
2020.11.21
在线时间
129 小时

发表于 2022.11.17 10:53:53 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
0
注册时间
2020.11.21
在线时间
129 小时

发表于 2022.11.18 11:29:20 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
0
注册时间
2020.11.21
在线时间
129 小时

发表于 2022.11.19 11:07:07 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2017.11.3
在线时间
337 小时

发表于 2022.11.24 08:59:17 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2019.7.28
在线时间
6 小时

发表于 2022.11.27 09:30:28 | 显示全部楼层
11月的第二十六天
新的一天加油!
回复 支持 反对

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2019.7.28
在线时间
6 小时

发表于 2022.11.27 11:07:45 | 显示全部楼层
11月的第二十六天
新的一天加油!
回复 支持 反对

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2017.11.3
在线时间
337 小时

发表于 2022.11.29 08:04:34 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2022.8.22
在线时间
5 小时

发表于 2022.11.29 11:46:27 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表