前面学习了一些箱图的画图技巧，这次我们试着重复一个paper中的图。从这个paper中来看，是在普通箱图的基础上又分成了Tumor和Normal2个group，并且每个group进行了统计学的显著性检验。另外按照mean/median进行了降序排序。但是并不是每个project都有2组。利用统计的方法是wilcoxon rank-sum test。
library(ggplot2)
library(ggpubr)
library(tidyverse)
library(latex2exp)
library(rstatix)
data <- read.csv("hsa-let-7a-5p.TCGA_PanCancer_Expression_Data.csv", header = T)
head(data)
data_new <- data %>%
group_by(project) %>%
mutate(median = median(expr), group_max = max(expr)) %>% #按project来计算均值
arrange(desc(median)) #按median降序排序
# 调整因子顺序，因为geom_boxplot在画图的时候，对于字符变量，会按照ASCII排序，所以要指定他们的顺序。
data_new$project <- factor(data_new$project, levels = unique(data_new$project))
ggplot(data_new,aes(project, expr, fill = group))+
geom_boxplot(outlier.shape = 21, outlier.fill = "white")+
scale_fill_manual(values = c("Tumor" = "#d6503a", "Normal" = "#5488ef"))+
theme_bw()+
theme(axis.text = element_text(face = "bold"),
axis.text.x = element_text(angle = 45, hjust = 1))
从图中来看，实现了2个group的划分，也实现了从大到小的排序。
下面我们来添加如何添加统计检验的结果，对于每个project对于tumor和normal进行显著性检验。
如果只针对含有normal和tumor的project的话，用上个帖子的stat_compare_means就可以实现。
data_new1 <- data_new %>% group_by(project) %>% mutate(group_number = length(unique(group))) #这里我们主要计算每个project含有group的数量是1还是2
data_new2 <- data_new1[which(data_new1$group_number == 2),] #挑选group数量是2的project，也就是挑选了含有tumor和normal的project
ggplot(data_new2, aes(x = project, y =expr, fill = group))+
geom_boxplot(outlier.shape = 21, outlier.fill = "white")+
scale_fill_manual(name = NULL, values = c("Tumor" = "#d6503a", "Normal" = "#5488ef"))+
theme_bw()+
theme(axis.text = element_text(face = "bold"),
axis.text.x = element_text(angle = 45, hjust = 1),
axis.line = element_line(colour = "black"),
panel.border = element_blank(),
panel.background = element_blank())+
stat_compare_means(data = data_new2,aes(x=project,y=expr,group=group),label = "p.signif")
但是需要注意的是并不是每一个都含有normal，所以我们可能需要首先计算一下显著性检验，然后手动添加上去这个label。
data_new3 <- unique(data_new2[,c(2, 7,8)])
# 此数据用于添加显著性检验的label：
stat_data <- data_new2 %>%
group_by(project) %>%
wilcox_test(expr ~ group) %>%
adjust_pvalue(method = "fdr") %>%
add_significance('p')%>%
add_xy_position(x = "project")
stat_data2 <- merge(stat_data, data_new3, by = "row.names", all = T) #主要为了加入group的最大值，来控制label的Y轴位置。
ggplot(data_new, aes(x = project, y =expr, fill = group))+
geom_boxplot(outlier.shape = 21, outlier.fill = "white")+
scale_fill_manual(name = NULL, values = c("Tumor" = "#d6503a", "Normal" = "#5488ef"))+
theme_bw()+
theme(axis.text = element_text(face = "bold"),
axis.text.x = element_text(angle = 45, hjust = 1),
axis.line = element_line(colour = "black"),
panel.border = element_blank(),
panel.background = element_blank(),
legend.position = c(0.99, 0.99),
legend.justification = c(1,1))+
annotate(geom = "text", x=stat_data2$x, y = stat_data2$group_max+0.25, label = as.character(stat_data2$p.signif)) #来手动的添加统计检验的文本上去
最后修改一下坐标注释等。
ggplot(data_new, aes(x = project, y =expr, fill = group))+
geom_boxplot(outlier.shape = 21, outlier.fill = "white")+
scale_fill_manual(name = NULL, values = c("Tumor" = "#d6503a", "Normal" = "#5488ef"))+
theme_bw()+
labs(y = TeX("miRNA Level ($log_{2}$CPM)", bold = T),x="")+
theme(axis.text = element_text(face = "bold"),
axis.text.x = element_text(angle = 45, hjust = 1),
axis.line = element_line(colour = "black"),
panel.border = element_blank(),
panel.background = element_blank(),
legend.position = c(0.99, 0.99),
legend.justification = c(1,1))+
annotate(geom = "text", x=stat_data2$x, y = stat_data2$group_max+0.25, label = as.character(stat_data2$p.signif))
但是，还有一个中括号，我不太清楚怎么添加比较好。
