数据可视化——R语言为ggplot图形添加P值和显著性水平

• 准备
• • 安装和加载R包
  • 示例数据
• 均值比较的方法
• 用于添加P值的R函数
• • compare_means()
  • stat_compare_means()
• 独立双样本组间比较
• 配对双样本组间比较
• 多组样本的组间比较
• 多个分组变量
• 其他绘图方式

本文对一篇英文博客进行翻译，博客原文链接： Add P-values and Significance Levels to ggplots

概述：本文介绍如何轻松地为ggplot图形添加P值和显著性水平：

• 比较两组或多组的均值
• 自动地将P值和显著性水平添加到ggplot图形中，如箱形图，点图，条形图和折线图等

使用工具： R语言中的ggplot2包和ggpubr包

安装和加载R包

本文使用 ggpubr 包，要求版本高于0.1.3。 ggpubr 是一个基于ggplot2的计算工具包。

• 直接输入以下命令从 CARN 中下载安装

install.packages("ggpubr")
也可以从GitHub中下载安装最新版本
 
if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")
加载ggpubr包
 
library(ggpubr)
ggpubr的官方文档可在以下位置获得：http://www.sthda.com/english/rpkgs/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中用于两组或多组间均值比较的标准统计方法在之前的文章也有描述：comparing means in R
 
均值比较的常见方法：
 
方法
R实现函数
描述
T-test
t.test()
比较两组(参数检验)
Wilcoxon test
wilcox.test()
比较两组(非参数检验)
ANOVA
aov()或anova()
比较多组(参数检验)
Kruskal-Wallis
kruskal.test()
比较多组(非参数检验)
 
以下链接提供了各种方法的详细介绍：
 
 
单样本均值与已知均值的比较
 
  
 
单样本T检验（参数检验）
 
 
单样本Wilcoxon检验（非参数检验）
 
 
 
独立双样本均值比较
 
  
非配对双样本T检验（参数检验）
非配对双样本Wilcoxon检验（非参数检验）
 
 
配对双样本均值比较
 
  
配对双样本T检验（参数检验）
配对双样本Wilcoxon检验（非参数检验）
 
 
两组以上的均值比较
 
  
方差分析（ANOVA, 参数检验） 
    
单因素方差分析
双因素方差分析
 
Kruskal-Wallis检验（替代单因素方差分析的非参数检验）
 
 
用于添加P值的R函数
 
介绍两个ggpubr包中的函数
 
compare_means()：用于执行均值比较
stat_compare_means()：用于在ggplot图形中自动添加P值和显著性水平
 
compare_means()
 
该函数用于执行均值比较。该函数与标准的R函数相比，灵活性更强。
 
简化形式如下：
 
compare_means(formula, data, method = "wilcox.test", paired = FALSE,
  group.by = NULL, ref.group = NULL, ...)
formula：指定一个公式，公式形式为 x ~ group，其中，x 表示一个数值型变量，group 表示一个因子型变量，包含一个或多个水平。例如，一个示例公式为 formula = TP53 ~ cancer_group，表示在 cancer_group 对应的各水平间比较TP53的表达水平；也可以同时指定多个响应变量，如 formula = c(TP53, PTEN) ~ cancer_group。
data：指定一个数据框(data.frame)，数据框需包含formula中的变量。
method：指定统计检验的方法。默认为“wilcox.test”，即Wilcoxon检验（非参数检验）；也可指定其他统计方法： 
  
“t.test”，即T检验（参数检验）。“t.test”和“wilcox.test”用于两组样本间的比较。当超过两组时，将会执行两两比较(pairwise comparison)。
“anova”（参数检验）或 “kruskal.test”（非参数检验），用于执行多组间的单因素方差分析。
 
paired：指定一个逻辑变量，表示是否需要执行配对检验，仅适用于t.test 和wilcox.test。
group.by：指定一个分组变量的字符名，用于在统计检验之前对数据进行分组。当存在group.by指定的变量时，均值比较将在不同水平的各个子集数据中执行。
ref.group：指定一个组别的字符名，作为对照组(reference group)。如果指定，各个分组水平将与对照组水平进行比较。也可指定ref.group为“.all.”，表示每个分组水平将于所有分组水平（如base-mean）进行比较。
 
stat_compare_means()
 
该函数是对ggplot2的扩展，可将均值比较后的P值添加到ggplot图形中，如箱形图、点图、条形图和折线图等。
 
简化形式如下：
 
stat_compare_means(mapping = NULL, comparisons = NULL hide.ns = FALSE,
                   label = NULL,  label.x = NULL, label.y = NULL,  ...)
 
mapping：通过 aes() 设置绘图时的aesthetic
 
 
comparisons：指定一个列表(list)，每个列表元素需为长度等于2的向量。向量的内容可以为X轴的两个组别名（字符型），也可以是两个感兴趣组的组别索引（整数值），表示采用指定的两个组别进行比较。
 
 
hide.ns：逻辑变量，如果设为TRUE，显示显著性水平时将隐藏 ns 字样，即组间差异不显著时不显示 ns 字样。
 
 
label：指定一个字符串，表示标签类型。可为：“p.signif”（显示显著性水平），“p.format”（显示格式化的P值）。
 
 
label.x, label.y：指定一个数值，表示显示标签的绝对坐标位置。
 
 
…：传递给函数compare_means()的参数，如method、paired、ref.group。
 
 
独立双样本组间比较
 
执行统计检验
 
compare_means(len ~ supp, data = ToothGrowth)
示例结果如下：
 
## # A tibble: 1 x 8
##     .y. group1 group2      p  p.adj p.format p.signif   method
## 1   len     OJ     VC 0.0645 0.0645    0.064       ns Wilcoxon
method默认为“wilcox.test”（非参数检验），可指定method = “t.test”，表示T检验（参数检验）
 
返回值为具有以下列的数据框：
 
 
.y.：用于统计检验的数值变量
 
 
p：P值
 
 
p.adj：调整后的P值，调整P值的默认方法为p.adjust.method = “holm”
 
 
p.format: 格式化的P值
 
 
p.signif：显著性水平，即用不同数量的 * 表示显著性水平
 
 
method：用于组间比较的统计方法
 
 
创建添加P值的箱形图
 
p <- ggboxplot(ToothGrowth, x = "supp", y = "len",
          color = "supp", palette = "jco",
          add = "jitter")
#  Add p-value
p + stat_compare_means()
# Change method
p + stat_compare_means(method = "t.test")

 注意：显示P值的标签位置可以通过如下参数来调整：label.x, label.y, hjust 和vjust。
 
显示P值的标签默认为 compare_means() 返回值中的 method 和 p 的组合。也可以通过 aes() 函数指定为其他显示形式。例如：
 
aes(label = ..p.format..) 或 aes(label = paste0(“p =”, ..p.format..))
表示只显示格式化的P值，而不显示method
aes(label = ..p.signif..)
表示仅显示显著性水平
aes(label = paste0(..method.., “\n”, “p =”, ..p.format..))
表示在method名和P值之间添加换行符(“\n”)
示例如下：
 
p + stat_compare_means( aes(label = ..p.signif..), 
                        label.x = 1.5, label.y = 40)

 另外，也可以将参数label指定为字符向量：
 
p + stat_compare_means( label = "p.signif", label.x = 1.5, label.y = 40)
配对双样本组间比较
 
执行统计检验
 
compare_means(len ~ supp, data = ToothGrowth, paired = TRUE)
示例结果如下：
 
## # A tibble: 1 x 8
##     .y. group1 group2       p   p.adj p.format p.signif   method
## 1   len     OJ     VC 0.00431 0.00431   0.0043       ** Wilcoxon
使用函数 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
compare_means(len ~ dose,  data = ToothGrowth, method = "anova")
示例结果如下：




    
 
## # A tibble: 1 x 6
##     .y.        p    p.adj p.format p.signif method
## 1   len 9.53e-16 9.53e-16  9.5e-16     ****  Anova
添加全局检验的P值（所有组比较总的P值）
 
# Default method = "kruskal.test" for multiple groups
ggboxplot(ToothGrowth, x = "dose", y = "len",
          color = "dose", palette = "jco")+
  stat_compare_means()
# Change method to anova
ggboxplot(ToothGrowth, x = "dose", y = "len",
          color = "dose", palette = "jco")+
  stat_compare_means(method = "anova")
两两比较(Pairwise comparisons)
 
如果分组变量包含两个以上的水平，两两比较的检验(pairwise test)将自动执行。默认方法为“wilcox.test”，也可设置为“t.test”。
 
# Perorm pairwise comparisons
compare_means(len ~ dose,  data = ToothGrowth)
示例结果如下：
 
## # A tibble: 3 x 8
##     .y. group1 group2        p    p.adj p.format p.signif   method
## 1   len    0.5      1 7.02e-06 1.40e-05  7.0e-06     **** Wilcoxon
## 2   len    0.5      2 8.41e-08 2.52e-07  8.4e-08     **** Wilcoxon
## 3   len      1      2 1.77e-04 1.77e-04  0.00018      *** Wilcoxon
# Visualize: Specify the comparisons you want
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
如果需要指定标签显示的Y轴位置，可使用参数label.y
 
ggboxplot(ToothGrowth, x = "dose", y = "len",
          color = "dose", palette = "jco")+ 
  stat_compare_means(comparisons = my_comparisons, label.y = c(29, 35, 40))+
  stat_compare_means(label.y = 45)

 注：ggsignif 包也可以很方便的为条形图添加组间比较的P值
 
相对于对照组的多重两两比较的检验(Multiple pairwise tests)
 
# Pairwise comparison against reference
compare_means(len ~ dose,  data = ToothGrowth, ref.group = "0.5",
              method = "t.test")
示例结果如下：
 
## # A tibble: 2 x 8
##     .y. group1 group2        p    p.adj p.format p.signif method
## 1   len    0.5      1 6.70e-09 6.70e-09  6.7e-09     **** T-test
## 2   len    0.5      2 1.47e-16 2.94e-16  < 2e-16     **** T-test
# Visualize
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
相对于所有组(base-mean)的多重两两比较的检验
 
# Comparison of each group against base-mean
compare_means(len ~ dose,  data = ToothGrowth, ref.group = ".all.",
              method = "t.test")
示例结果如下：
 
## # A tibble: 3 x 8
##     .y. group1 group2        p    p.adj p.format p.signif method
## 1   len  .all.    0.5 1.24e-06 3.73e-06  1.2e-06     **** T-test
## 2   len  .all.      1 5.67e-01 5.67e-01     0.57       ns T-test
## 3   len  .all.      2 1.37e-05 2.74e-05  1.4e-05     **** T-test
# Visualize
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
下面使用Github中可用的骨髓瘤数据集展示一些典型情况，其中，与“.all.”的比较将很有用。
 
将依据患者molecular进行分组，绘制各个组别DEPDC1基因的表达水平。目的是比较各个组别之间是否存在差别，如果有差别，差别又在哪里。
 
要回答以上问题，可以在所有7个组之间进行两两比较(pairwise comparison)。 由于组别较多，将导致很多种组别组合，这将很难解释。
 
一个简单的解决办法是将7组中的每一组与“.all.”（如base-mean）比较。当检验结果显著时，可以得出结论：与所有组相比，xxx组的DEPDC1的表达水平显著降低或显著升高。
 
# Load myeloma data from GitHub
myeloma <- read.delim("https://raw.githubusercontent.com/kassambara/data/master/myeloma.txt")
# Perform the test
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 method
## 1 DEPDC1  .all.       Cyclin D-1 0.149690 0.44907  0.14969       ns T-test
## 2 DEPDC1  .all.       Cyclin D-2 0.523143 1.00000  0.52314       ns T-test
## 3 DEPDC1  .all.     Hyperdiploid 0.000282 0.00169  0.00028      *** T-test
## 4 DEPDC1  .all. Low bone disease 0.005084 0.02542  0.00508       ** T-test
## 5 DEPDC1  .all.              MAF 0.086107 0.34443  0.08611       ns T-test
## 6 DEPDC1  .all.            MMSET 0.576291 1.00000  0.57629       ns T-test
## # ... with 1 more rows
注意：上面的R代码中myeloma数据的下载地址存在错误，可从以下链接找到myeloma数据：https://github.com/kassambara/data/blob/master/myeloma.txt
 然后自己将全部数据复制并保存为myeloma.txt，从本地读取myeloma.txt。
 
# Visualize the expression profile
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
根据上图的结果，可以得出proliferation组的DEPDC1表达水平显著升高；Hyperdiploid组和Low bone disease组中DEPDC1表达水平显著降低。




    
 
注意：想要隐藏 ns 标志，可以设置参数 hide.ns = TRUE
 
# Visualize the expression profile
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
多个分组变量
 
使用另一个变量进行分组后再执行独立双样本比较
 
执行统计检验：
 
compare_means(len ~ supp, data = ToothGrowth, 
              group.by = "dose")
示例结果如下：
 
## # A tibble: 3 x 9
##    dose   .y. group1 group2       p  p.adj p.format p.signif   method
## 1   0.5   len     OJ     VC 0.02319 0.0464    0.023        * Wilcoxon
## 2   1.0   len     OJ     VC 0.00403 0.0121    0.004       ** Wilcoxon
## 3   2.0   len     OJ     VC 1.00000 1.0000    1.000       ns Wilcoxon
在上面的示例中，对于分类变量dose的每一个水平，分类变量supp又将数据分为两个水平：OJ和VC，然后在这两个水平上对数值变量len进行均值比较。
 
可视化：创建一个按组划分的多面板框图(此处为“dose”)
 
# Box plot facetted by "dose"
p <- ggboxplot(ToothGrowth, x = "supp", y = "len",
          color = "supp", palette = "jco",
          add = "jitter",
          facet.by = "dose", short.panel.labs = FALSE)
# Use only p.format as label. Remove method name.
p + stat_compare_means(label = "p.format")
# Or use significance symbol as label
p + stat_compare_means(label =  "p.signif", label.x = 1.5)
注意：想要隐藏 ns 标志，可以设置参数hide.ns = TRUE
 
可视化：创建一个包含所有箱形图的单一面板。X表示dose，Y表示len，颜色表示supp
 
p <- ggboxplot(ToothGrowth, x = "dose", y = "len",
          color = "supp", palette = "jco",
          add = "jitter")
p + stat_compare_means(aes(group = supp))
# Show only p-value
p + stat_compare_means(aes(group = supp), label = "p.format")
# Use significance symbol as label
p + stat_compare_means(aes(group = supp), label = "p.signif")
使用另一个变量进行分组后再执行配对双样本比较
 
执行统计检验：
 
compare_means(len ~ supp, data = ToothGrowth, 
              group.by = "dose", paired = TRUE)
实例结果如下：
 
## # A tibble: 3 x 9
##    dose   .y. group1 group2      p  p.adj p.format p.signif   method
## 1   0.5   len     OJ     VC 0.0330 0.0659    0.033        * Wilcoxon
## 2   1.0   len     OJ     VC 0.0191 0.0572    0.019        * Wilcoxon
## 3   2.0   len     OJ     VC 1.0000 1.0000    1.000       ns Wilcoxon
可视化：创建一个按组划分的多面板框图(此处为“dose”)
 
# Box plot facetted by "dose"
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)
# Use only p.format as label. Remove method name.
p + stat_compare_means(label = "p.format", paired = TRUE)
其他绘图方式
 
条形图和折线图（一个分组变量）
 
# Bar plot of mean +/-se
ggbarplot(ToothGrowth, x = "dose", y = "len", add = "mean_se")+
  stat_compare_means() +                                         # Global p-value
  stat_compare_means(ref.group = "0.5", label = "p.signif",
                     label.y = c(22, 29))                   # compare to ref.group
# Line plot of mean +/-se
ggline(ToothGrowth, x = "dose", y = "len", add = "mean_se")+
  stat_compare_means() +                                         # Global p-value
  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)
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))
注意：经过实际测试，笔者发现R语言中的统计方法计算结果的P值与SPSS中的P值存在差异。如，常规的方差分析(ANOVA) + 事后两两组间比较（如Bonferroni校正）使用上述R函数就很难得出与SPSS中一致的结果。如果需要使用SPSS的统计P值，建议对生成的图形进行后期修改。
 
References
 
原文链接：Add P-values and Significance Levels to ggplots
                        原文链接：http://www.sthda.com/english/articles/24-ggpubr-publication-ready-plots/76-add-p-values-and-significance-levels-to-ggplots/
p <- iris %>% 
  ggplot(aes(Species, Sepal.Length)) + 
  geom_violin(aes(fill = Species),trim = F) + 
  stat_boxplot(geom = 'errorbar', width = 0.1) + 
                                    我们的回答是你需要提取相关组的表达量，进行检验后再用ggplot函数添加即可；或者直接提取数据用ggplot作图那么显著添加也就不成问题了。时隔3月，我们这里提供 了一种函数，可以进行基因在两组之间的显著性分析。同时可进行批量的基因分析。这就是所有内容了，其实这样检验你用不用得到倒是其次，主要是这里面包含一些小的细节知识点，学会了就能和其他内容融汇贯通了，自己感悟吧！更多精彩内容请至KS科研分享与服务公众号。显著性检验函数，有点长，可自行保存成R文件，然后每次使用的时候source一下就可以了。
                                    “ggplot2中的次级坐标轴”生活科学哥-R语言科学 2020-06-12 8：35在平时作图中，我们有时希望在一个坐标中进行二个坐标轴的设定，也是为了方便数据的显示。这个过程在EXCEL等当中比较容易实现，但是，如何在R中实现呢？今天我们就来讲一讲操作的过程。数据准备先准备如下数据：library(ggplot2)library(scales)library(magrittr)...
它包含一项评估维生素C对豚鼠牙齿生长的影响的研究数据。
实验在60只豚鼠上进行，其中每只豚鼠通过两种递送方法（橙汁，OJ，或抗坏血酸，VC）分别接受三种剂量水平的维生素C量（0.5、1和2 mg /天， VC）。
                                    上图这样的单细胞StackedVlnPlot在高分文章中出现比较多，比较适合美观的展示多个marker gene的表达分布，而目前Seurat画小提琴图的函数VlnPlot是不能实现这样堆叠效果的。我们在《Seurat结果转为scanpy可处理对象》中介绍过scanpy中的sc.pl.stacked_violin函数可以实现StackedVlnPlot的功能：我们在《生信工程师的自我修养...
data <- data.frame(
  individual=paste( "Mister ", seq(1,60), sep=""),
  group=c( rep('A', 10), rep('B', 30), rep('C', 14), rep('D', 6)) ,
  value1=sample( seq(
这是一个库和一对命令行程序，用于计算p值和临界值，以测试小数单位根和协整。
 该程序是基于James G. MacKinnon和MortenØrregaardNielsen，“分数单位根和协整检验的数值分布函数”，《应用计量经济学》 ，第1卷，第1期中的描述进行建模的。 2014年1月29日，第161-171页。 编译此软件时，会自动下载与该论文相关的数据档案库中的数据 。
 该程序由两个命令行二进制文件组成： fdpval和fdcrit 。 前者计算测试统计的p值，后者根据p值计算关键的测试统计值。 在不带参数的情况下运行每个程序会提供程序的使用信息。
 可从获得该库的最新版本的源代码。
Debian，Ubuntu和其他Debian衍生Linux发行版
二进制版本可从。 通过在/etc/apt/sources.list中添加以下行之一，可以获得包含Debian衍生发行版
                                    该函数中n阶循环是指可能你的数据集存在多种处理，也就是说存在不同处理变量间组合成新的子数据集。该函数的输出结果返回一个数据框，包含不同处理及bootstrap重采样的模型斜率的结果，最后可以利用统计学检验不同处理变量间的显著性差异。第一，可以实现不同类型模型的斜率比较；第二，可以解决不同变量处理间的多阶循环。主要参数的介绍data : 数据集R : bootstrap的次数sub : 需要对哪些变量进行划分成子数据集的循环equa : 对应构建模型的fomula。