R语言添加p-value和显著性标记

2024-02-01 01:40
文章标签 语言 标记 value 显著性

本文主要是介绍R语言添加p-value和显著性标记，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

R语言添加p-value和显著性标记，原文链接 https://mp.weixin.qq.com/s/gRw0krS3LY7c0QK9y47EJw

作者简介 Introduction

taoyan：伪码农，R语言爱好者，爱开源。

个人博客: https://ytlogos.github.io/

往期回顾

  • 一条命令轻松绘制CNS顶级配图-ggpubr
  • R语言聚类分析–cluster, factoextra

上篇文章中提了一下如何通过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中常用的比较方法主要有下面几种：

方法R函数描述
T-testt.test()比较两组(参数)
Wilcoxon testwilcox.test()比较两组(非参数)
ANOVAaov()或anova()比较多组(参数)
Kruskal-Walliskruskal.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()

image

# 使用t-test统计检验方法
p+stat_compare_means(method = "t.test")

image

上述显著性标记可以通过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)

image

多组比较 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()

image

使用其他的方法，anova

#使用其他的方法, anova
ggboxplot(ToothGrowth, x="dose", y="len", color = "dose", palette = "jco")+stat_compare_means(method = "anova")

image

成对比较

默认进行两两比较

# 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

image

可以通过修改参数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-valuestat_compare_means(label = "p.signif", method = "t.test", ref.group = "0.5") # Pairwise comparison against reference

image

参考组也可以设置为.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-valuestat_compare_means(label = "p.signif", method = "t.test",ref.group = ".all.")#Pairwise comparison against all

image

为什么有时需要将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 meanstat_compare_means(method = "anova", label.y = 1600)+ # Add global annova p-valuestat_compare_means(label = "p.signif", method = "t.test", ref.group = ".all.")# Pairwise comparison against all

image

从图中可以看出，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 meanstat_compare_means(method = "anova", label.y = 1600)+ # Add global annova p-valuestat_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")

image

# label绘制显著性水平
p+stat_compare_means(label = "p.signif", label.x = 1.5)

image

将所有箱线图绘制在一个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")

image

进行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)

image

其他图形

有误差棒的条形图

实际上我以前的文章里有纯粹用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))

image

有误差棒的线图

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))

image

条形图(两个分组变量)

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)

image

线图-两分组

# 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))

image

环境信息

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: defaultlocale:
[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热心肠 2SemanticScholar 3geenmedical
  • 扩增子分析：1图表解读 2分析流程 3统计绘图 4功能预测
  • 科研经验：1云笔记 2云协作 3公众号
  • 系列教程：1Biostar 2微生物组 3宏基因组
  • 生物科普 1肠道细菌 2人体上的生命 3生命大跃进 4细胞的暗战 5人体奥秘

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外1200+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。
image

学习扩增子、宏基因组科研思路和分析实战，关注“宏基因组”
image

点击阅读原文，跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA

这篇关于R语言添加p-value和显著性标记的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/665700。

相关文章

Go语言中泄漏缓冲区的问题解决

Go语言中泄漏缓冲区的问题解决

《Go语言中泄漏缓冲区的问题解决》缓冲区是一种常见的数据结构,常被用于在不同的并发单元之间传递数据,然而,若缓冲区使用不当,就可能引发泄漏缓冲区问题,本文就来介绍一下问题的解决,感兴趣的可以了解一下... 目录引言泄漏缓冲区的基本概念代码示例：泄漏缓冲区的产生项目场景：Web 服务器中的请求缓冲场景描述代码
阅读更多...
Go语言如何判断两张图片的相似度

Go语言如何判断两张图片的相似度

《Go语言如何判断两张图片的相似度》这篇文章主要为大家详细介绍了Go语言如何中实现判断两张图片的相似度的两种方法,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下... 在介绍技术细节前，我们先来看看图片对比在哪些场景下可以用得到：图片去重：自动删除重复图片，为存储空间"瘦身"。想象你是一个
阅读更多...
Go语言中Recover机制的使用

Go语言中Recover机制的使用

《Go语言中Recover机制的使用》Go语言的recover机制通过defer函数捕获panic,实现异常恢复与程序稳定性,具有一定的参考价值,感兴趣的可以了解一下... 目录引言Recover 的基本概念基本代码示例简单的 Recover 示例嵌套函数中的 Recover项目场景中的应用Web 服务器中
阅读更多...
Go语言中使用JWT进行身份验证的几种方式

Go语言中使用JWT进行身份验证的几种方式

《Go语言中使用JWT进行身份验证的几种方式》本文主要介绍了Go语言中使用JWT进行身份验证的几种方式,包括dgrijalva/jwt-go、golang-jwt/jwt、lestrrat-go/jw... 目录简介1. github.com/dgrijalva/jwt-go安装：使用示例：解释：2. gi
阅读更多...
Go 语言中的 Struct Tag 的用法详解

Go 语言中的 Struct Tag 的用法详解

《Go语言中的StructTag的用法详解》在Go语言中,结构体字段标签（StructTag）是一种用于给字段添加元信息（metadata）的机制,常用于序列化（如JSON、XML）、ORM映... 目录一、结构体标签的基本语法二、json:"token"的具体含义三、常见的标签格式变体四、使用示例五、使用
阅读更多...
解读@ConfigurationProperties和@value的区别

解读@ConfigurationProperties和@value的区别

《解读@ConfigurationProperties和@value的区别》：本文主要介绍@ConfigurationProperties和@value的区别及说明,具有很好的参考价值,希望对大家... 目录1. 功能对比2. 使用场景对比@ConfigurationProperties@Value3. 核
阅读更多...
Go语言使用slices包轻松实现排序功能

Go语言使用slices包轻松实现排序功能

《Go语言使用slices包轻松实现排序功能》在Go语言开发中,对数据进行排序是常见的需求,Go1.18版本引入的slices包提供了简洁高效的排序解决方案,支持内置类型和用户自定义类型的排序操作,本... 目录一、内置类型排序：字符串与整数的应用1. 字符串切片排序2. 整数切片排序二、检查切片排序状态：
阅读更多...
基于Go语言实现Base62编码的三种方式以及对比分析

基于Go语言实现Base62编码的三种方式以及对比分析

《基于Go语言实现Base62编码的三种方式以及对比分析》Base62编码是一种在字符编码中使用62个字符的编码方式,在计算机科学中,,Go语言是一种静态类型、编译型语言,它由Google开发并开源,... 目录一、标准库现状与解决方案1. 标准库对比表2. 解决方案完整实现代码（含边界处理）二、关键实现细
阅读更多...
如何合理管控Java语言的异常

如何合理管控Java语言的异常

《如何合理管控Java语言的异常》：本文主要介绍如何合理管控Java语言的异常问题,具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地方,望不吝赐教... 目录1、介绍2、Thorwable类3、Error4、Exception类4.1、检查异常4.2、运行时异常5、处理方式5.1. 捕获异常
阅读更多...
C语言中的常见进制转换详解(从二进制到十六进制)

C语言中的常见进制转换详解(从二进制到十六进制)

《C语言中的常见进制转换详解(从二进制到十六进制)》进制转换是计算机编程中的一个常见任务,特别是在处理低级别的数据操作时,C语言作为一门底层编程语言,在进制转换方面提供了灵活的操作方式,今天,我们将深... 目录1、进制基础2、C语言中的进制转换2.1 从十进制转换为其他进制十进制转二进制十进制转八进制十进
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2