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

R语言添加p-value和显著性标记,原文链接 https://mp.weixin.qq.com/s/gRw0krS3LY7c0QK9y47EJw 作者简介 Introduction taoyan:伪码农,R语言爱好者,爱开源。 个人博客: https://ytlogos.gi...

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-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两组间和Y40位置
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-value
  stat_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-value
  stat_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 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

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

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

猜你喜欢

写在后面

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

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

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

  • 发表于 2018-02-18 13:38
  • 阅读 ( 12665 )
  • 分类:其他组学

0 条评论

请先 登录 后评论
不写代码的码农
刘永鑫

工程师

64 篇文章

作家榜 »

  1. 祝让飞 118 文章
  2. 柚子 91 文章
  3. 刘永鑫 64 文章
  4. admin 57 文章
  5. 生信分析流 55 文章
  6. SXR 44 文章
  7. 张海伦 31 文章
  8. 爽儿 25 文章