R語言添加p-value和顯著性標記

2022-01-09 宏基因組


taoyan：偽碼農，R語言愛好者，愛開源。

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

往期回顧


上篇文章中提了一下如何通過ggpubr包為ggplot圖添加p-value以及顯著性標記,本文將詳細介紹。利用數據集ToothGrowth進行演示。

# 先加載包

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()函數
該函數主要用用法如下：

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

注釋：

比較獨立的兩組

compare_means(len~supp, data=ToothGrowth)

結果解釋：

   ## 繪製箱線圖

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()映射來更改：

舉個慄子：

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

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

compare_means(len~dose, data=ToothGrowth, method = "anova")

可視化

ggboxplot(ToothGrowth, x="dose", y="len", color = "dose", palette = "jco")+

stat_compare_means()

#使用其他的方法

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)

#可以指定比較哪些組

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.即所有的平均值

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


接下來利用survminer包中的數據集myeloma來講解一下為什麼有時候我們需要將ref.group設置為.all.

library(survminer)#沒安裝的先安裝再加載

data("myeloma")

head(myeloma)

我們將根據患者的分組來繪製DEPDC1基因的表達譜，看不同組之間是否存在顯著性的差異，我們可以在7組之間進行比較，但是這樣的話組間比較的組合就太多了，因此我們可以將7組中每一組與全部平均值進行比較，看看DEPDC1基因在不同的組中是否過表達還是低表達。

compare_means(DEPDC1~molecular_group, data = myeloma, ref.group = ".all.", method = "t.test")

# 可視化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中

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)

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

sessionInfo()

## R version 3.4.0 (2017-04-21)

## Platform: x86_64-w64-mingw32/x64 (64-bit)

## Running under: Windows 8.1 x64 (build 9600)

##

## 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.0 ggpubr_0.1.3 magrittr_1.5 ggplot2_2.2.1

##

## loaded via a namespace (and not attached):

## [1] Rcpp_0.12.11 compiler_3.4.0 plyr_1.8.4

## [4] tools_3.4.0 digest_0.6.12 evaluate_0.10

## [7] tibble_1.3.3 gtable_0.2.0 nlme_3.1-131

## [10] lattice_0.20-35 rlang_0.1.1 Matrix_1.2-10

## [13] psych_1.7.5 ggsci_2.4 DBI_0.6-1

## [16] cmprsk_2.2-7 yaml_2.1.14 parallel_3.4.0

## [19] gridExtra_2.2.1 dplyr_0.5.0 stringr_1.2.0

## [22] knitr_1.16 survMisc_0.5.4 rprojroot_1.2

## [25] grid_3.4.0 data.table_1.10.4 KMsurv_0.1-5

## [28] R6_2.2.1 km.ci_0.5-2 survival_2.41-3

## [31] foreign_0.8-68 rmarkdown_1.5 reshape2_1.4.2

## [34] tidyr_0.6.3 purrr_0.2.2.2 splines_3.4.0

## [37] backports_1.1.0 scales_0.4.1 htmltools_0.3.6

## [40] assertthat_0.2.0 mnormt_1.5-5 xtable_1.8-2

## [43] colorspace_1.3-2 ggsignif_0.2.0 labeling_0.3

## [46] stringi_1.1.5 lazyeval_0.2.0 munsell_0.4.3

## [49] broom_0.4.2 zoo_1.8-0

猜你喜歡寫在後面

為鼓勵讀者交流、快速解決科研困難，我們建立了「宏基因組」專業討論群，目前己有國內外120+ PI，1200+ 一線科研人員加入。參與討論，獲得專業解答，歡迎分享此文至朋友圈，並掃碼加主編好友帶你入群，務必備註「姓名-單位-研究方向-職稱/年級」。技術問題尋求幫助，首先閱讀《如何優雅的提問》學習解決問題思路，仍末解決群內討論，問題不私聊，幫助同行。

學習16S擴增子、宏基因組科研思路和分析實戰，關注「宏基因組」

相關焦點

  • R語言可視化學習筆記之添加p-value和顯著性標記
    上篇文章中提了一下如何通過ggpubr包為ggplot圖添加p-value以及顯著性標記,本文將詳細介紹。利用數據集ToothGrowth進行演示。添加p-value主要利用ggpubr包中的兩個函數：##compare_means()函數該函數主要用用法如下：compare_means(formula, data, method = "wilcox.test", paired = FALSE,group.by =
  • 使用ggpubr包添加p值和顯著性標記
    使用ggpubr包添加p值和顯著性標記今天我們的學習內容主要有兩個方面如何簡單比較兩組或多組的平均值如何自動化為箱線圖添加
  • 差異分析、顯著性標記及統計作圖的自動實現R代碼示例
    本人的很多經驗學自《R語言實戰 第二版》，它的154頁有這一段話。因此，對於非參數的檢驗方法，我也是一直在使用wilcoxon檢驗+p值校正的方法，代替多重比較。然後，還是大神師妹提醒的（就在今天，你們信嗎？
  • R語言 | 差異顯著性檢驗
    我們經常要比較兩組或多組數據是否具有顯著差異，同時我們還會用差異顯著性檢驗識別不同組樣品中具有顯著差異的變量。這篇推文會分別介紹經常使用的差異顯著性檢驗方法在R語言中的實現。方法選擇差異顯著性檢驗具有多種方法，分別針對不同的情況，我們要根據自身情況選擇合適的方法進行分析。
  • 69-R可視化12-用easylabel輕鬆手動添加標記
    Date : [[2021-12-06_Mon]]Tags : #R/index/02 #R/R可視化 #R/R數據科學 #R/R包 #R/偷懶者easylabel (r-project.org)[1]前言先前我介紹過ggrepel 這個包：[[67-R可視化11-用ggrepel更加美觀的添加標記（火山圖的實現）]]
  • R語言添加特殊文本以及annotate()函數的使用小技巧
    p+ annotate('text',x=1.2,y=7.5, #定義添加文本和其位置 label=expression(-log[10]*'(p value)'), size=8,color='red')+ #其他個性化定義 labs(y=expression(-log[10
  • R語言學習筆記之相關性矩陣分析及其可視化
    cor()只能計算出相關係數，無法給出顯著性水平p-value,Hmisc包裡的rcorr()函數能夠同時給出相關係數以及顯著性水平p-value。library(Hmisc)#加載包res2 <- rcorr(as.matrix(mydata))res2#可以用res2$r、res2$P來提取相關係數以及顯著性p-valueres2$r
  • R繪圖---分組樣本見兩兩T-test並添加顯著性標記
    小夥伴們在進行數據分析的時候，一定會遇到分組T測驗的問題，那麼這一類問題怎麼才能在R語言中實現呢，分析的結果如何以可視化的形式給大家展現出來呢？
  • 每天學習一點R:43.差異顯著性檢驗
    我們經常要比較兩組或多組數據是否具有顯著差異，同時我們還會用差異顯著性檢驗識別不同組樣品中具有顯著差異的變量。這篇推文會分別介紹經常使用的差異顯著性檢驗方法在R語言中的實現。方法選擇差異顯著性檢驗具有多種方法，分別針對不同的情況，我們要根據自身情況選擇合適的方法進行分析。
  • 差異表達基因分析:差異倍數(fold change), 差異的顯著性(P-value) | 火山圖
    Furthermore, when the denominator is close to zero, the ratio is not stable, and the fold change value can be disproportionately affected by measurement noise.
  • 顯著性分析後如何標記「abc」?
    科學論文和雜誌中最常見的是用「字母」和「星號」標記差異顯著性，如果所要比較的數據較少，用「*」號標記比較簡單明了，而數據較多時就可能會用到了今天要講的「字母標記法」。今天的內容分兩個部分：「心法」和「招式」。在進行標記之前先把下面的「心法」認真看一遍，內容如下： 1.
  • [小工具] Wilcoxon-tk小工具計算分組數值的秩和檢驗p-value
    Wilcoxon-tk 主要解決的問題是根據不同分組信息計算不同顯著性差異 p-value, 使用wilcoxon-command-line的C程序函數,簡單封裝,完成直接讀數值表,然後批處理計算P values。
  • R語言 | 一個多元線性回歸在R中的實現示例
    R2_adj <- c(R2_adj, fit_stat$adj.r.squared) #提取校正後 R2 p_value <- c(p_value, fit_stat$coefficients[2,4]) #提取顯著性 p 值} env_stat <- data.frame(row.names = env, R2_adj,
  • 技術貼 | R語言:組學關聯分析和pheatmap可視化
    split將ko-metabolite結果列拆成兩列結果保留r_value p_value顯著相關標記「*」library(psych)library(stringr)correlate = function(other, metabo, route){    #  讀取方式：check.name=F, row.names
  • 對於P_value的理解
    通常有t檢驗（用於樣本含量較小，倆樣本均數以及兩本均數與總體均數的之間的比較）、z檢驗（用於樣本含量較大，兩個平均數的差異是否顯著）、卡方檢驗（用於倆定類變量，實際觀測值和理論推導值的偏離程度）等~     通常，我們會設定原假設為H0，指兩樣本無差異，Ha為H0的補，指兩者有差異，而所謂P值：「p-value is the probability we get this sample or
  • 小孩都看得懂的 p-value
    ，那什麼是 p-value 呢？Well Done, Wikipedia, 這下連大人都徹底不懂 p-value 了。但希望下面極簡的講解能讓小孩懂什麼是 p-value。做試驗：p-value牢記：p-value 就是硬幣公平時觀測結果發生的概率。
  • R語言繪製分組柱狀圖示例
    ggplot2繪製帶誤差線以及顯著性標記的分組柱狀圖    簡介了堆疊柱狀圖，順便再簡介一下分組柱狀圖在R中的繪製。
  • r語言一元回歸模型專題及常見問題 - CSDN
    F檢驗法：F檢驗用於對所有的自變量X在整體上看對於Y的線性顯著性，也是用P-value判斷顯著性，小於0.01更小時說明整體上自變量與Y相關關係顯著。Error，為參數的標準差，sd(a), sd(b)t value，為t值，為T檢驗的值Pr(>|t|) ，表示P-value值，用於T檢驗判定，匹配顯著性標記顯著性標記
  • 數據挖掘常見的 p-value 解讀
    p-value，那什麼是 p-value 呢？Well Done, Wikipedia, 這下連大人都徹底不懂 p-value 了。但希望下面極簡的講解能讓小孩懂什麼是 p-value。做試驗：p-value牢記：p-value 就是硬幣公平時觀測結果發生的概率。 第一次硬幣是反面，p-value 是多少？你覺得硬幣不公平嗎？絕壁不會啊！
  • R語言完美重現STAMP結果圖
    統計學檢驗在繪製Extended error bar之前首先要對數據進行差異顯著性檢驗，以選出豐度在不同組間具有顯著差異的物種或功能基因，這裡以兩組數據為例，使用的檢驗方式通常為t-test和Wilcox秩和檢驗。