日韩黑丝制服一区视频播放|日韩欧美人妻丝袜视频在线观看|九九影院一级蜜桃|亚洲中文在线导航|青草草视频在线观看|婷婷五月色伊人网站|日本一区二区在线|国产AV一二三四区毛片|正在播放久草视频|亚洲色图精品一区

搜索
分享

用ggplot2畫火山圖

 dachuanz 2018-05-16

用ggplot2畫火山圖看代碼似乎還是很簡單，但是真的來畫了，細節(jié)還是蠻多的。

把代碼寫在這里，以便以后用。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
library(ggplot2)
library（Cairo）
data <- read.delim("./data/g1_VS_g3.txt",header = T,sep="\t",na.strings = "")
# 設(shè)置顏色域
data$threshold = as.factor(data$pvalues < 0.05 & abs(log2(data$foldchange)) >=1.5)
Cairo(file="volcan_PNG_300_dpi.png",
      type="png",
      units="in",
      bg="white",
      width=8,
      height=6,
      pointsize=12,
      dpi=300)
ggplot(data=data,
            aes(x=log2(foldchange), y =-log10(pvalues),
                colour=threshold)) +
  geom_point(alpha=0.4, size=1.75) +
  xlim(c(-4, 4)) +
  xlab("log2 fold change") + ylab("-log10 p-value") +
  ggtitle("Volcano picture of DEG")+
  theme_gray() +
  geom_vline(xintercept=c(-1.5,1.5),lty=4,col="grey",lwd=0.5)+ # 在x軸-1.5與1.5的位置畫兩根豎線
  geom_hline(yintercept = -log10(0.05),lty=4,col="grey",lwd=0.5)+ 在p value 0.05的位置畫一根橫線
  theme(legend.position="none")
dev.off()

圖1. volcano of DEGs

這只是一個粗略的圖，要細致，還需要調(diào)整網(wǎng)格，坐標軸，顏色，字體大小，等等……

我們能不能將下調(diào)，上調(diào)的分別改變顏色體現(xiàn)呢？

那就來試試吧，思路其實就是顏色的映射到不同的 見前代碼colour=threshold ，之前的設(shè)置只設(shè)定了2個水平，顯著or不顯著。那其實只要把顯著與不顯著改成上調(diào)、下調(diào)、無顯著改變就可以了。見代碼：

1
data$threshold <- as.factor(ifelse(data$pvalues < 0.05 & abs(log2(data$foldchange)) >=1.5,ifelse(log2(data$foldchange) > 1.5 ,'Up','Down'),'Not'))

建立好了分層，我們還需要做的一點是，分配顏色。

1
scale_color_manual(values=c("blue", "grey","red"))

這樣就搞定了。但是，但是，強迫癥估計還要加上圖歷，映射不同顏色是吧，呵呵。那就看下面的代碼吧

1
theme(legend.position="right")

這樣通過legend.position 指定顯示的位置，或許你還想更改圖例的標題，%……%……這需要你自己動手了，google or bing去尋求學習。

下面放出完整代碼：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
library(ggplot2)
library(ggthemes)
library(Cairo)
data <- read.delim("./data/g1_VS_g3.txt",header = T,sep="\t",na.strings = "")
data$threshold <- as.factor(ifelse(data$pvalues < 0.05 & abs(log2(data$foldchange)) >=1.5,ifelse(log2(data$foldchange) > 1.5 ,'Up','Down'),'Not'))
##Construct the plot object
# with legend
Cairo(file="volcan_PNG_300_lengend_dpi.png",
      type="png",
      units="in",
      bg="white",
      width=5.5,
      height=5,
      pointsize=12,
      dpi=300)
ggplot(data=data,
            aes(x=log2(foldchange), y =-log10(pvalues),
                colour=threshold,fill=threshold)) +
  scale_color_manual(values=c("blue", "grey","red"))+
  geom_point(alpha=0.4, size=1.2) +
  xlim(c(-4, 4)) +
  theme_bw(base_size = 12, base_family = "Times") +
  geom_vline(xintercept=c(-1.5,1.5),lty=4,col="grey",lwd=0.6)+
  geom_hline(yintercept = -log10(0.05),lty=4,col="grey",lwd=0.6)+
  theme(legend.position="right",
        panel.grid=element_blank(),
        legend.title = element_blank(),
        legend.text= element_text(face="bold", color="black",family = "Times", size=8),
        plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(face="bold", color="black", size=12),
        axis.text.y = element_text(face="bold",  color="black", size=12),
        axis.title.x = element_text(face="bold", color="black", size=12),
        axis.title.y = element_text(face="bold",color="black", size=12))+
  labs(x="log2 (fold change)",y="-log10 (p-value)",title="Volcano picture of DEG")
dev.off()
# without legend
Cairo(file="volcan_PNG_300_dpi.png",
      type="png",
      units="in",
      bg="white",
      width=6,
      height=5,
      pointsize=12,
      dpi=300)
ggplot(data=data,
       aes(x=log2(foldchange), y =-log10(pvalues),
           colour=threshold,fill=threshold)) +
  scale_color_manual(values=c("blue", "grey","red"))+
  geom_point(alpha=0.4, size=1.2) +
  xlim(c(-4, 4)) +
  theme_bw(base_size = 12, base_family = "Times") +
  geom_vline(xintercept=c(-1.5,1.5),lty=4,col="grey",lwd=0.6)+
  geom_hline(yintercept = -log10(0.05),lty=4,col="grey",lwd=0.6)+
  theme(legend.position="none",
        panel.grid=element_blank(),
        # legend.title = element_blank(),
        # legend.text= element_text(face="bold", color="black",family = "Times", size=8),
        plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(face="bold", color="black", size=12),
        axis.text.y = element_text(face="bold",  color="black", size=12),
        axis.title.x = element_text(face="bold", color="black", size=12),
        axis.title.y = element_text(face="bold",color="black", size=12))+
  labs(x="log2 (fold change)",y="-log10 (p-value)",title="Volcano picture of DEG")
dev.off()

 

圖2. volcano of DEGs without legend

 

圖3. volcano of DEGs with legend

    本站是提供個人知識管理的網(wǎng)絡(luò)存儲空間，所有內(nèi)容均由用戶發(fā)布，不代表本站觀點。請注意甄別內(nèi)容中的聯(lián)系方式、誘導購買等信息，謹防詐騙。如發(fā)現(xiàn)有害或侵權(quán)內(nèi)容，請點擊一鍵舉報。
    轉(zhuǎn)藏 分享 獻花（0）

    0條評論

    發(fā)表

    請遵守用戶 評論公約

    類似文章 更多