R语言添加p-value和显著性标记
R语言添加p-value和显著性标记，原文链接 https://mp.weixin.qq.com/s/gRw0krS3LY7c0QK9y47EJw
作者简介 Introduction
taoyan：伪码农，R语言爱好者，爱开源。
个人博客: https://ytlogos.github.io/
往期回顾
上篇文章中提了一下如何通过ggpubr包为ggplot图添加p-value以及显著性标记,本文将详细介绍。利用数据集ToothGrowth进行演示。
ggpubr安装和加载
# 直接从CRAN安装 install.packages("ggpubr", repo="http://cran.us.r-project.org") #先加载包 library(ggpubr) #加载数据集ToothGrowth data("ToothGrowth") head(ToothGrowth)
数据格式如下：
len supp dose 1 4.2 VC 0.5 2 11.5 VC 0.5 3 7.3 VC 0.5 4 5.8 VC 0.5 5 6.4 VC 0.5 6 10.0 VC 0.5
比较方法
R中常用的比较方法主要有下面几种：
| T-test | t.test() | 比较两组(参数) |
| Wilcoxon test | wilcox.test() | 比较两组(非参数) |
| ANOVA | aov()或anova() | 比较多组(参数) |
| Kruskal-Wallis | kruskal.test() | 比较多组(非参数) |
添加p-value
主要利用ggpubr包中的两个函数：
compare_means():可以进行一组或多组间的比较
stat_compare_mean():自动添加p-value、显著性标记到ggplot图中
compare_means()函数
该函数主要用用法如下：
compare_means(formula, data, method = "wilcox.test", paired = FALSE,
group.by = NULL, ref.group = NULL, ...)
注释：
formula:形如x~group，其中x是数值型变量，group是因子，可以是一个或者多个
data：数据集
method:比较的方法，默认为"wilcox.test", 其他可选方法为："t.test"、"anova"、"kruskal.test"
paired:是否要进行paired test(TRUE or FALSE)
group_by: 比较时是否要进行分组
ref.group: 是否需要指定参考组
stat_compare_means()函数
主要用法：
stat_compare_means(mapping = NULL, comparisons = NULL hide.ns = FALSE,
label = NULL, label.x = NULL, label.y = NULL, ...)
注释：
mapping:由aes()创建的一套美学映射
comparisons:指定需要进行比较以及添加p-value、显著性标记的组
hide.ns:是否要显示显著性标记ns
label:显著性标记的类型，可选项为：p.signif(显著性标记)、p.format(显示p-value)
label.x、label.y:显著性标签调整
…:其他参数
比较独立的两组
compare_means(len~supp, data=ToothGrowth)
统计结果如下:
# A tibble: 1 x 8 .y. group1 group2 p p.adj p.format p.signif method <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr> 1 len OJ VC 0.06449067 0.06449067 0.064 ns Wilcoxon
结果解释：
.y:测试中使用的y变量
p:p-value
p.adj:调整后的p-value。默认为p.adjust.method="holm"
p.format:四舍五入后的p-value
p.signif:显著性水平
method:用于统计检验的方法
绘制箱线图
# 绘制箱线图 p <- ggboxplot(ToothGrowth, x="supp", y="len", color = "supp", palette = "jco", add = "jitter") # 添加p-value, 默认是Wilcoxon test p+stat_compare_means()
# 使用t-test统计检验方法 p+stat_compare_means(method = "t.test")
上述显著性标记可以通过label.x、label.y、hjust及vjust来调整 显著性标记可以通过aes()映射来更改：
aes(label=..p.format..)或aes(lebel=paste0("p=",..p.format..))：只显示p-value，不显示统计检验方法
aes(label=..p.signif..)：仅显示显著性水平
aes(label=paste0(..method..,"\n", "p=",..p.format..))：p-value与显著性水平分行显示
举个栗子：
# 显示显著性水平，位置在1.5两组间和Y轴40位置 p+stat_compare_means(aes(label=..p.signif..), label.x = 1.5, label.y = 40) # 也可以将标签指定为字符向量，不要映射，只需将p.signif两端的..去掉即可 p+stat_compare_means(label = "p.signif", label.x = 1.5, label.y = 40)
比较两个paired sample
# 比较两个paired sample compare_means(len~supp, data=ToothGrowth, paired = TRUE) # 利用ggpaired()进行可视化 ggpaired(ToothGrowth, x="supp", y="len", color = "supp", line.color = "gray", line.size = 0.4, palette = "jco")+ stat_compare_means(paired = TRUE)
多组比较 Global test
anova进行多组比较
compare_means(len~dose, data=ToothGrowth, method = "anova")
统计结果如下：
# A tibble: 1 x 6 .y. p p.adj p.format p.signif method <chr> <dbl> <dbl> <chr> <chr> <chr> 1 len 9.532727e-16 9.532727e-16 9.5e-16 **** Anova
可视化, default Kruskal-Wallis
ggboxplot(ToothGrowth, x="dose", y="len", color = "dose", palette = "jco")+ stat_compare_means()
使用其他的方法，anova
#使用其他的方法, anova ggboxplot(ToothGrowth, x="dose", y="len", color = "dose", palette = "jco")+ stat_compare_means(method = "anova")
成对比较
默认进行两两比较
# Pairwise comparisons:如果分组变量中包含两个以上的水平，那么会自动进行pairwise test,默认方法为”wilcox.test” compare_means(len~dose, data=ToothGrowth)
比较结果：
# A tibble: 3 x 8 .y. group1 group2 p p.adj p.format p.signif method <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr> 1 len 0.5 1 7.020855e-06 1.404171e-05 7.0e-06 **** Wilcoxon 2 len 0.5 2 8.406447e-08 2.521934e-07 8.4e-08 **** Wilcoxon 3 len 1 2 1.772382e-04 1.772382e-04 0.00018 *** Wilcoxon
可以指定比较哪些组
# 可以指定比较哪些组 my_comparisons <- list(c("0.5", "1"), c("1", "2"), c("0.5", "2")) ggboxplot(ToothGrowth, x="dose", y="len", color = "dose",palette = "jco")+ stat_compare_means(comparisons=my_comparisons)+ # Add pairwise comparisons p-value stat_compare_means(label.y = 50) # Add global p-value
可以通过修改参数label.y来更改标签的位置
ggboxplot(ToothGrowth, x="dose", y="len", color = "dose",palette = "jco")+ stat_compare_means(comparisons=my_comparisons, label.y = c(29, 35, 40))+ # Add pairwise comparisons p-value stat_compare_means(label.y = 45) # Add global p-value
至于通过添加线条来连接比较的两组，这一功能已由包ggsignif实现
设定参考组设定参考组
compare_means(len~dose, data=ToothGrowth, ref.group = "0.5", #以dose=0.5组为参考组 method = "t.test" ) # 可视化 ggboxplot(ToothGrowth, x="dose", y="len", color = "dose", palette = "jco")+ stat_compare_means(method = "anova", label.y = 40)+ # Add global p-value stat_compare_means(label = "p.signif", method = "t.test", ref.group = "0.5") # Pairwise comparison against reference
参考组也可以设置为.all.即所有的平均值
# 参考组也可以设置为.all.即所有的平均值 compare_means(len~dose, data=ToothGrowth, ref.group = ".all.", method = "t.test") #可视化 ggboxplot(ToothGrowth, x="dose", y="len", color = "dose", palette = "jco")+ stat_compare_means(method = "anova", label.y = 40)+# Add global p-value stat_compare_means(label = "p.signif", method = "t.test", ref.group = ".all.")#Pairwise comparison against all
为什么有时需要将ref.group设置为.all
接下来利用survminer包中的数据集myeloma来讲解一下为什么有时候我们需要将ref.group设置为.all.
# 利用survminer包中的数据集myeloma来讲解一下为什么有时候我们需要将ref.group设置为.all. install.packages("survminer") library(survminer) #没安装的先安装再加载 data("myeloma") head(myeloma)
我们将根据患者的分组来绘制DEPDC1基因的表达谱，看不同组之间是否存在显著性的差异，我们可以在7组之间进行比较，但是这样的话组间比较的组合就太多了，因此我们可以将7组中每一组与全部平均值进行比较，看看DEPDC1基因在不同的组中是否过表达还是低表达。
compare_means(DEPDC1~molecular_group, data = myeloma, ref.group = ".all.", method = "t.test")
比较结果如下：
# A tibble: 7 x 8 .y. group1 group2 p p.adj p.format p.signif <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> 1 DEPDC1 .all. Cyclin D-1 2.877529e-01 1.000000e+00 0.29 ns 2 DEPDC1 .all. Cyclin D-2 4.244240e-01 1.000000e+00 0.42 ns 3 DEPDC1 .all. Hyperdiploid 2.725486e-08 1.907840e-07 2.7e-08 **** 4 DEPDC1 .all. Low bone disease 5.258400e-06 3.155040e-05 5.3e-06 **** 5 DEPDC1 .all. MAF 2.538126e-01 1.000000e+00 0.25 ns 6 DEPDC1 .all. MMSET 5.784193e-01 1.000000e+00 0.58 ns 7 DEPDC1 .all. Proliferation 2.393921e-05 1.196961e-04 2.4e-05 **** # ... with 1 more variables: method <chr>
可视化DEPDC1基因表达谱
ggboxplot(myeloma, x="molecular_group", y="DEPDC1", color = "molecular_group", add = "jitter", legend="none")+ rotate_x_text(angle = 45)+ geom_hline(yintercept = mean(myeloma$DEPDC1), linetype=2)+# Add horizontal line at base mean stat_compare_means(method = "anova", label.y = 1600)+ # Add global annova p-value stat_compare_means(label = "p.signif", method = "t.test", ref.group = ".all.")# Pairwise comparison against all
从图中可以看出，DEPDC1基因在Proliferation组中显著性地过表达，而在Hyperdiploid和Low bone disease显著性地低表达
我们也可以将非显著性标记ns去掉，只需要将参数hide.ns=TRUE
ggboxplot(myeloma, x="molecular_group", y="DEPDC1", color = "molecular_group", add = "jitter", legend="none")+ rotate_x_text(angle = 45)+ geom_hline(yintercept = mean(myeloma$DEPDC1), linetype=2)+# Add horizontal line at base mean stat_compare_means(method = "anova", label.y = 1600)+ # Add global annova p-value stat_compare_means(label = "p.signif", method = "t.test", ref.group = ".all.", hide.ns = TRUE)# Pairwise comparison against all
分面多组同时比较
按另一个变量进行分组之后进行统计检验，比如按变量dose进行分组：
compare_means(len~supp, data=ToothGrowth, group.by = "dose") #可视化 p <- ggboxplot(ToothGrowth, x="supp", y="len", color = "supp", palette = "jco", add = "jitter", facet.by = "dose", short.panel.labs = FALSE)#按dose进行分面 #label只绘制p-value p+stat_compare_means(label = "p.format")
# label绘制显著性水平 p+stat_compare_means(label = "p.signif", label.x = 1.5)
将所有箱线图绘制在一个panel中
# 将所有箱线图绘制在一个panel中 p <- ggboxplot(ToothGrowth, x="dose", y="len", color = "supp", palette = "jco", add = "jitter") p+stat_compare_means(aes(group=supp)) # 只显示p-value p+stat_compare_means(aes(group=supp), label = "p.format") # 显示显著性水平 p+stat_compare_means(aes(group=supp), label = "p.signif")
进行paired sample检验
compare_means(len~supp, data=ToothGrowth, group.by = "dose", paired = TRUE) # 可视化 p <- ggpaired(ToothGrowth, x="supp", y="len", color = "supp", palette = "jco", line.color="gray", line.size=0.4, facet.by = "dose", short.panel.labs = FALSE) # 按dose分面 # 只显示p-value p+stat_compare_means(label = "p.format", paired = TRUE)
其他图形
有误差棒的条形图
实际上我以前的文章里有纯粹用ggplot2实现
ggbarplot(ToothGrowth, x="dose", y="len", add = "mean_se")+ stat_compare_means()+ stat_compare_means(ref.group = "0.5", label = "p.signif", label.y = c(22, 29))
有误差棒的线图
ggline(ToothGrowth, x="dose", y="len", add = "mean_se")+ stat_compare_means()+ stat_compare_means(ref.group = "0.5", label = "p.signif", label.y = c(22, 29))
条形图(两个分组变量)
ggbarplot(ToothGrowth, x="dose", y="len", add = "mean_se", color = "supp", palette = "jco", position = position_dodge(0.8))+ stat_compare_means(aes(group=supp), label = "p.signif", label.y = 29)
线图-两分组
# line, multiply pair group ggline(ToothGrowth, x="dose", y="len", add = "mean_se", color = "supp", palette = "jco")+ stat_compare_means(aes(group=supp), label = "p.signif", label.y = c(16, 25, 29))
环境信息
sessionInfo()
信息如下：
R version 3.4.3 (2017-11-30) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows >= 8 x64 (build 9200) Matrix products: default locale: [1] LC_COLLATE=Chinese (Simplified)_China.936 [2] LC_CTYPE=Chinese (Simplified)_China.936 [3] LC_MONETARY=Chinese (Simplified)_China.936 [4] LC_NUMERIC=C [5] LC_TIME=Chinese (Simplified)_China.936 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] survminer_0.4.2 bindrcpp_0.2 ggpubr_0.1.6.999 magrittr_1.5 [5] ggplot2_2.2.1 loaded via a namespace (and not attached): [1] Rcpp_0.12.13 compiler_3.4.3 plyr_1.8.4 [4] bindr_0.1 tools_3.4.3 digest_0.6.15 [7] tibble_1.3.4 gtable_0.2.0 nlme_3.1-131 [10] lattice_0.20-35 pkgconfig_2.0.1 rlang_0.1.4 [13] Matrix_1.2-12 psych_1.7.8 ggsci_2.8 [16] cmprsk_2.2-7 yaml_2.1.15 parallel_3.4.3 [19] gridExtra_2.3 knitr_1.19 dplyr_0.7.4 [22] stringr_1.2.0 survMisc_0.5.4 grid_3.4.3 [25] tidyselect_0.2.3 data.table_1.10.4-3 glue_1.2.0 [28] KMsurv_0.1-5 R6_2.2.2 km.ci_0.5-2 [31] survival_2.41-3 foreign_0.8-69 tidyr_0.8.0 [34] purrr_0.2.4 reshape2_1.4.2 splines_3.4.3 [37] scales_0.5.0 assertthat_0.2.0 mnormt_1.5-5 [40] xtable_1.8-2 colorspace_1.3-2 ggsignif_0.4.0 [43] labeling_0.3 stringi_1.1.5 lazyeval_0.2.1 [46] munsell_0.4.3 broom_0.4.3 zoo_1.8-1
猜你喜欢
- 热文：1高分文章 2不可或缺的人 3图表规范
- 一文读懂：1微生物组 2寄生虫益处 3进化树
- 必备技能：1提问 2搜索 3Endnote
- 文献阅读 1热心肠 2SemanticScholar3geenmedical
- 扩增子分析：1图表解读 2分析流程 3统计绘图 4功能预测
- 科研经验：1云笔记 2云协作 3公众号
- 系列教程：1Biostar 2微生物组 3宏基因组
- 生物科普 1肠道细菌 2人体上的生命 3生命大跃进4细胞的暗战 5人体奥秘
写在后面
为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外1200+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。 
学习扩增子、宏基因组科研思路和分析实战，关注“宏基因组” 
点击阅读原文，跳转最新文章目录阅读https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA
- 发表于 2018-02-18 13:38
- 阅读 ( 14226 )
- 分类：其他组学
