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语言开发一个命令行文件管理工具》这篇文章主要为大家详细介绍了如何使用Go语言开发一款命令行文件管理工具,支持批量重命名,删除,创建,移动文件,需要的小伙伴可以了解下... 目录一、工具功能一览二、核心代码解析1. 主程序结构2. 批量重命名3. 批量删除4. 创建文件/目录5. 批量移动三、如何安
阅读更多...
python使用fastapi实现多语言国际化的操作指南

python使用fastapi实现多语言国际化的操作指南

《python使用fastapi实现多语言国际化的操作指南》本文介绍了使用Python和FastAPI实现多语言国际化的操作指南,包括多语言架构技术栈、翻译管理、前端本地化、语言切换机制以及常见陷阱和... 目录多语言国际化实现指南项目多语言架构技术栈目录结构翻译工作流1. 翻译数据存储2. 翻译生成脚本
阅读更多...
Go语言中三种容器类型的数据结构详解

Go语言中三种容器类型的数据结构详解

《Go语言中三种容器类型的数据结构详解》在Go语言中,有三种主要的容器类型用于存储和操作集合数据：本文主要介绍三者的使用与区别,感兴趣的小伙伴可以跟随小编一起学习一下... 目录基本概念1. 数组（Array）2. 切片（Slice）3. 映射（Map）对比总结注意事项基本概念在 Go 语言中，有三种主要
阅读更多...
C语言中自动与强制转换全解析

C语言中自动与强制转换全解析

《C语言中自动与强制转换全解析》在编写C程序时,类型转换是确保数据正确性和一致性的关键环节,无论是隐式转换还是显式转换,都各有特点和应用场景,本文将详细探讨C语言中的类型转换机制,帮助您更好地理解并在... 目录类型转换的重要性自动类型转换（隐式转换）强制类型转换（显式转换）常见错误与注意事项总结与建议类型
阅读更多...
Go语言利用泛型封装常见的Map操作

Go语言利用泛型封装常见的Map操作

《Go语言利用泛型封装常见的Map操作》Go语言在1.18版本中引入了泛型,这是Go语言发展的一个重要里程碑,它极大地增强了语言的表达能力和灵活性,本文将通过泛型实现封装常见的Map操作,感... 目录什么是泛型泛型解决了什么问题Go泛型基于泛型的常见Map操作代码合集总结什么是泛型泛型是一种编程范式，允
阅读更多...
Android kotlin语言实现删除文件的解决方案

Android kotlin语言实现删除文件的解决方案

《Androidkotlin语言实现删除文件的解决方案》：本文主要介绍Androidkotlin语言实现删除文件的解决方案,在项目开发过程中,尤其是需要跨平台协作的项目,那么删除用户指定的文件的... 目录一、前言二、适用环境三、模板内容1.权限申请2.Activity中的模板一、前言在项目开发过程中，尤
阅读更多...
C语言小项目实战之通讯录功能

C语言小项目实战之通讯录功能

《C语言小项目实战之通讯录功能》：本文主要介绍如何设计和实现一个简单的通讯录管理系统,包括联系人信息的存储、增加、删除、查找、修改和排序等功能,文中通过代码介绍的非常详细,需要的朋友可以参考下... 目录功能介绍：添加联系人模块显示联系人模块删除联系人模块查找联系人模块修改联系人模块排序联系人模块源代码如下
阅读更多...
基于Go语言实现一个压测工具

基于Go语言实现一个压测工具

《基于Go语言实现一个压测工具》这篇文章主要为大家详细介绍了基于Go语言实现一个简单的压测工具,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下... 目录整体架构通用数据处理模块Http请求响应数据处理Curl参数解析处理客户端模块Http客户端处理Grpc客户端处理Websocket客户端
阅读更多...
使用SQL语言查询多个Excel表格的操作方法

使用SQL语言查询多个Excel表格的操作方法

《使用SQL语言查询多个Excel表格的操作方法》本文介绍了如何使用SQL语言查询多个Excel表格,通过将所有Excel表格放入一个.xlsx文件中,并使用pandas和pandasql库进行读取和... 目录如何用SQL语言查询多个Excel表格如何使用sql查询excel内容1. 简介2. 实现思路3
阅读更多...
Go语言实现将中文转化为拼音功能

Go语言实现将中文转化为拼音功能

《Go语言实现将中文转化为拼音功能》这篇文章主要为大家详细介绍了Go语言中如何实现将中文转化为拼音功能,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下... 有这么一个需求：新用户入职 创建一系列账号比较麻烦，打算通过接口传入姓名进行初始化。想把姓名转化成拼音。因为有些账号即需要中文也需要英
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2