R绘图:ggplot2绘制火山图

火山图是测序分析报告中最为核心的图片之一。绘制火山图的方法有许多,Excel和第三方软件等,本文主要运用ggplot2和ggrepel两个R包演示[1][2]
横坐标为Fold change(倍数),越偏离中心差异倍数越大;纵坐标为P value(P值),值越大差异越显著。

(一):绘图分析:

基本要素: 区分上调、下调的散点辅助线坐标轴图例标签
运用参数: geom_point()、geom_vline and geom_hline、labs()、theme()、geom_text_repel()

注:区分上、下调基因,就要对差异基因进行统计学分析,主要包括统计显著性和生物学差异显著性;需要用到ifelse函数 [3]

(二):数据处理

rm(list = ls())
library(ggplot2)
# 数据来自https://www.jianshu.com/p/a973c29c9103 [参考资料2]
dataset <- read.table(file = "C:/Users/Administrator/Documents/dataset_volcano.txt", 
                      header = TRUE, sep = "")
# 设置p_value和logFC的阈值
cut_off_pvalue = 0.0000001  #统计显著性
cut_off_logFC = 1           #差异倍数值

# 根据阈值参数,上调基因设置为‘up’,下调基因设置为‘Down’,无差异设置为‘Stable’,并保存到change列中
dataset$change = ifelse(dataset$P.Value < cut_off_pvalue & abs(dataset$logFC) >= cut_off_logFC, 
                          ifelse(dataset$logFC> cut_off_logFC ,'Up','Down'),
                          'Stable')
head(dataset)
          gene     logFC  P.Value change
1 LOC100909539 -4.514013 4.33e-12   Down
2       Steap2 -3.678175 1.12e-11   Down
3        Trpt1  2.594115 1.21e-11     Up
4        Cers6 -3.630773 1.45e-11   Down
5       Hs3st1 -3.129658 1.74e-11   Down
6        Lama5 -2.772380 2.51e-11   Down

先从整体判断条件的真假:
  1. 不满足,记作"stable";
  2. 满足;
     logFC大于1,记作"up";
     logFC小于1,记作"down";

(三):绘图

p <- ggplot(
  # 数据、映射、颜色
  dataset, aes(x = logFC, y = -log10(P.Value), colour=change)) +
  geom_point(alpha=0.4, size=3.5) +
  scale_color_manual(values=c("#546de5", "#d2dae2","#ff4757"))+
  # 辅助线
  geom_vline(xintercept=c(-1,1),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = -log10(cut_off_pvalue),lty=4,col="black",lwd=0.8) +
  # 坐标轴
  labs(x="log2(fold change)",
       y="-log10 (p-value)")+
  theme_bw()+
  # 图例
  theme(plot.title = element_text(hjust = 0.5), 
        legend.position="right", 
        legend.title = element_blank())
p
  • Rplot01.png

(四):添加标签

# 将需要标记的基因放置在label列(logFC >= 5)
library(ggrepel)
dataset$label <- ifelse(dataset$P.Value > cut_off_pvalue & abs(dataset$logFC) >= 5,
                      as.character(dataset$gene), "")
dataset[2600:2606, ]
     gene      logFC  P.Value change label
2600   Ano5  0.1883082 1.20e-07 Stable      
2601  Rab3a  0.8501624 1.21e-07 Stable      
2602  Pank1  5.2116194 1.21e-07 Stable Pank1
2603  Ghitm  2.4541362 1.21e-07 Stable      
2604 Ndufa8  1.4141038 1.21e-07 Stable      
2605 Zc3h7a -0.9070267 1.21e-07 Stable      
2606  Trpm3  0.1879326 1.21e-07 Stable

p + geom_label_repel(data = dataset, aes(x = dataset$logFC, 
                                     y = -log10(dataset$P.Value), 
                                     label = label),
                 size = 3, box.padding = unit(0.5, "lines"),
                 point.padding = unit(0.8, "lines"), 
                 segment.color = "black", 
                 show.legend = FALSE)

关于ggrepel包,后期会写参数使用相关教程

  • Rplot02.png

(五):Volcano plot[4]

# 火山图应用
genes <- read.table(file = "C:/Users/Administrator/Documents/volcano.txt", header = TRUE)
genes$Significant <- ifelse(genes$padj < 0.05 & abs(genes$log2FoldChange) >= 1, 
                            ifelse(genes$log2FoldChange > 1, "Up", "Down"), "Stable")
ggplot(
  # 数据、映射、颜色
  genes, aes(x = log2FoldChange, y = -log10(pvalue))) +
  geom_point(aes(color = Significant), size=2) +
  scale_color_manual(values = c("blue","grey", "red")) +
  # 注释
  geom_text_repel(
    data = subset(genes, padj < 0.05 & abs(genes$log2FoldChange) >= 1),
    aes(label = Gene),
    size = 5,
    box.padding = unit(0.35, "lines"),
    point.padding = unit(0.3, "lines")) +
  # 辅助线
  geom_vline(xintercept=c(-1,1),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = 1,lty=4,col="black",lwd=0.8) +
  # 坐标轴
  labs(x="log2(fold change)",
       y="-log10 (p-value)") +
  # 图例
  theme(legend.position = "bottom")
Rplot01.png

进一步修饰

# 进一步修饰
ggplot(
  # 数据、映射、颜色
  genes, aes(x = log2FoldChange, y = -log10(pvalue))) +
  geom_point(aes(color = Significant), size=2) +
  scale_color_manual(values = c("blue","grey", "red")) +
  # 注释
  geom_label_repel(
    data = subset(genes, padj < 0.05 & abs(genes$log2FoldChange) >= 1),
    aes(label = Gene),
    size = 5,fill = "darkred", color = "white",
    box.padding = unit(0.35, "lines"),
    point.padding = unit(0.3, "lines")) +
  # 辅助线
  geom_vline(xintercept=c(-0.5, 0.5),lty=4,col="#666666",lwd=0.8) +
  geom_hline(yintercept = 1,lty=4,col="#666666",lwd=0.8) +
  # 坐标轴
  labs(x="log2(fold change)",
       y="-log10 (p-value)") +
  # 图例
  theme(legend.position = "bottom")
Rplot02.png

参考资料:


  1. 使用ggplot2和ggrepel包做火山图

  2. R 数据可视化 02 | 火山图

  3. R语言_ifelse()函数用法

  4. ggplot2 texts : Add text annotations to a graph in R software

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 203,362评论 5 477
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,330评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 150,247评论 0 337
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,560评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,580评论 5 365
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,569评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,929评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,587评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,840评论 1 297
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,596评论 2 321
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,678评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,366评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,945评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,929评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,165评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 43,271评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,403评论 2 342

推荐阅读更多精彩内容