1.編示例數(shù)據(jù)
set.seed(10086)
exp = matrix(rnorm(480,mean = 6),ncol = 120)
exp = round(exp,2)
rownames(exp) = paste0("gene",1:4)
colnames(exp) = paste0("test",1:120)
exp[,50:80] = exp[,50:80]+2
exp[,90:120] = exp[,90:120]+4
exp[1:4,1:4]
## test1 test2 test3 test4
## gene1 6.55 5.37 6.32 4.18
## gene2 3.26 6.25 5.63 5.78
## gene3 6.57 6.80 8.11 6.58
## gene4 6.49 7.08 8.49 4.38
到這一步，拿到了一個編的表達矩陣洁灵。
2.格式轉(zhuǎn)換
表達矩陣要畫ggplot2箱線圖饱岸，肯定是要先轉(zhuǎn)變格式的，否則ggplot2代碼無法下手寫徽千。
library(tidyr)
library(tibble)
library(dplyr)
dat = t(exp) %>%
as.data.frame() %>%
rownames_to_column() %>%
mutate(group = rep(c("control","treat1","treat2"),each = 40))
pdat = dat%>%
pivot_longer(cols = starts_with("gene"),
names_to = "gene",
values_to = "value")
head(pdat)
## # A tibble: 6 × 4
## rowname group gene value
## <chr> <chr> <chr> <dbl>
## 1 test1 control gene1 6.55
## 2 test1 control gene2 3.26
## 3 test1 control gene3 6.57
## 4 test1 control gene4 6.49
## 5 test2 control gene1 5.37
## 6 test2 control gene2 6.25
3.畫箱線圖
library(ggplot2)
p1 = ggplot(pdat, aes(y=value, x=gene))+
geom_boxplot(aes(fill=group),outlier.shape = NA)+
theme_bw()+
scale_fill_manual(values = c("#2fa1dd", "#f87669", "#e6b707"))
p1
[圖片上傳失敗...(image-382a55-1692258856792)]
p2 = ggplot(pdat, aes(y=value, x=gene))+
geom_bar(aes(fill=group),position = position_dodge(0.9),stat = "summary")+
theme_bw()+
scale_fill_manual(values = c("#2fa1dd", "#f87669", "#e6b707"))+
stat_summary(aes(group = group), fun.data = 'mean_se',
geom = "errorbar", colour = "black",
width = 0.2,position = position_dodge(0.9))
p2
[圖片上傳失敗...(image-297a28-1692258856792)]
4.加顯著性標記
4.1習以為常版
library(ggpubr)
p1 +
stat_compare_means(aes(group = group, label = after_stat(p.signif)),method = "wilcox.test")
[圖片上傳失敗...(image-ee95bd-1692258856792)]
p2 + stat_compare_means(aes(group = group, label = after_stat(p.signif)),method = "wilcox.test",label.y = 10.5)
[圖片上傳失敗...(image-2a7112-1692258856792)]
這個超級簡單苫费，但顯著性統(tǒng)一在頭頂上，然后也不支持每個基因在兩兩分組中的顯著性單獨計算双抽。
4.2 錯落有致版
用rstatix里的函數(shù)計算出顯著性和添加顯著性的橫縱坐標位置百框，再用stat_pvalue_manual添加上去。
library(rstatix)
stat.test <- pdat %>%
group_by(gene) %>%
wilcox_test(value ~ group) %>%
adjust_pvalue() %>%
add_significance("p.adj") %>%
add_xy_position(x="gene")
p1 + stat_pvalue_manual(stat.test,label = "p.adj.signif",hide.ns = F,
tip.length = 0.01)
[圖片上傳失敗...(image-241141-1692258856792)]
p2 + stat_pvalue_manual(stat.test,label = "p.adj.signif",hide.ns = F,
tip.length = 0.01)
[圖片上傳失敗...(image-16df19-1692258856792)]
4.3 手動添加版
例如上次畫的AUC值柱狀圖牍汹，它的縱坐標和誤差棒不是畫圖時計算出來的铐维，而是拿了現(xiàn)成的數(shù)值來畫柬泽，就無法用上面的方法計算。如果能自行用其他方法計算出組間比較的顯著性如何嫁蛇，可以手動添加锨并。
rm(list = ls())
load("dats.Rdata")
p = ggplot(dats, aes(x = category, y = y, fill = x)) +
geom_bar(stat = "identity", position = "dodge") +
geom_text(aes(label = round(y,2)), vjust = 4, size = 3,
position = position_dodge(width = 0.9)) +
geom_errorbar(aes(ymin = min, ymax = max),
width = 0.2,
position = position_dodge(width = 0.9)) +
scale_fill_manual(values = c("#2fa1dd", "#f87669", "#e6b707"))+
theme_bw()
p
[圖片上傳失敗...(image-c65f9d-1692258856792)]
p + geom_signif(data = dats,y_position=c(1,0.97,0.83,0.89,0.89,0.92),
xmin=c(0.7,1,1.7,2,2.7,3),
xmax=c(1,1.3,2,2.3,3,3.3),
annotation=c("*","**","***","ns","*","**"),
tip_length=0.01,
vjust = 0.3)
[圖片上傳失敗...(image-75aa38-1692258856791)]
