同名公主號：BBio
rstatix包支持管道符吆寨，與tidyverse包設(shè)計(jì)理念相一致伐蒂，用于執(zhí)行基本的統(tǒng)計(jì)檢驗(yàn)遏暴，包括t檢驗(yàn)般卑、Wilcoxon檢驗(yàn)、ANOVA矢门、Kruskal-Wallis和相關(guān)性分析撒桨。
ggpubr包提供了一些易于使用的函數(shù)煤禽，基于ggplot2液样，可繪制用于文章發(fā)表的各種圖形振亮。當(dāng)然，ggpubr也提供了一些統(tǒng)計(jì)函數(shù)鞭莽，但是rstatix在統(tǒng)計(jì)方面更加豐富坊秸。
如何使用兩者繪制出下面好看的圖形呢？以單細(xì)胞數(shù)據(jù)為例澎怒。
//細(xì)胞比例組間差異
-
生成測試數(shù)據(jù)
library(Seurat)
library(tidyverse)
library(rstatix)
library(ggpubr)
library(ggplot2)
#生成測試數(shù)據(jù)
pbmc <- pbmc_small
func <- function(x)str_c('S', x%%12+1)
sample <- pbmc[['groups']] %>% add_column(order=1:dim(pbmc)[2]) %>% mutate(remainder=func(order)) %>% pull(remainder)
df <- table(Idents(pbmc), sample) %>%
prop.table(2) %>%
as.data.frame.array() %>%
rownames_to_column("cluster") %>%
pivot_longer(unique(sample), names_to="sample", values_to="percent") %>%
left_join(data.frame(sample=str_c("S", 1:12), group=rep(c("G1", "G2"), each = 6)))
df
# A tibble: 18 × 4
# cluster sample percent group
# <chr> <chr> <dbl> <chr>
# 1 0 S2 0.571 G1
# 2 0 S3 0.429 G1
# 3 0 S4 0.385 G2
# 4 0 S5 0.385 G2
# 5 0 S6 0.462 G2
# 6 0 S1 0.462 G1
# 7 1 S2 0.214 G1
# 8 1 S3 0.357 G1
# 9 1 S4 0.385 G2
#10 1 S5 0.385 G2
#11 1 S6 0.231 G2
#12 1 S1 0.308 G1
#13 2 S2 0.214 G1
#14 2 S3 0.214 G1
#15 2 S4 0.231 G2
#16 2 S5 0.231 G2
#17 2 S6 0.308 G2
#18 2 S1 0.231 G1
-
細(xì)胞比例組間差異柱狀圖
#只畫cluster 2
data <- df[df$cluster==2, ]
comparisons <- list(c("G1", "G2"))
ggbarplot(data, x = "group", y = "percent",
color = "group",
palette = c("#00AFBB", "#E7B800"),
width = 0.2,
x.text.angle = 90,
add=c('mean_se')) +
labs(x=NULL) + theme(legend.position="none", panel.border=element_rect(color='black', fill=NA, size=0.7, linetype="solid")) +
stat_compare_means(comparisons = comparisons, method = 't.test', label = 'p.signif')
image-20220927212134435.png
-
細(xì)胞比例組間差異箱線圖
#只畫cluster2
data <- df[df$cluster==2, ]
comparisons <- list(c("G1", "G2"))
p1 <- ggboxplot(data, x = "group", y = "percent",
color = "group",
palette =c("#00AFBB", "#E7B800"),
add = "jitter", shape = "group") +
labs(x=NULL) + theme(legend.position="none", panel.border=element_rect(color='black', fill=NA, size=0.7, linetype="solid")) +
stat_compare_means(comparisons = comparisons, method = 't.test', label = 'p.signif', size=4)
#畫所有cluster
p2 <- ggboxplot(df, x = "group", y = "percent",
color = "group",
facet.by = "cluster",
palette =c("#00AFBB", "#E7B800"),
add = "jitter", shape = "group") +
labs(x=NULL) + theme(legend.position="none", panel.border=element_rect(color='black', fill=NA, size=0.7, linetype="solid")) +
stat_compare_means(comparisons = comparisons, method = 't.test', label = 'p.signif', size=4)
p1 + p2
image-20220927213718138.png
#以上圖形褒搔，很坐標(biāo)都是組別，ggpubr中的stat_compare_means可以直接畫圖喷面。但是如果橫坐標(biāo)是cluster呢站超？
stat.test <- df %>%
group_by(cluster) %>%
t_test(percent ~ group) %>%
adjust_pvalue(method = "bonferroni") %>%
add_significance("p")
stat.test <- stat.test %>%
add_xy_position(x = "cluster", dodge = 0.8)
# A tibble: 3 × 16
# cluster .y. group1 group2 n1 n2 statistic df p p.adj
# <chr> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
#1 0 percent G1 G2 6 6 2.03 9.72 0.0706 0.212
#2 1 percent G1 G2 6 6 -2.52 9.02 0.0326 0.0978
#3 2 percent G1 G2 6 6 0.210 9.04 0.838 1
# p.signif y.position groups x xmin xmax
# <chr> <dbl> <named list> <dbl> <dbl> <dbl>
#1 ns 0.707 <chr [2]> 1 0.8 1.2
#2 ns 0.540 <chr [2]> 2 1.8 2.2
#3 ns 0.469 <chr [2]> 3 2.8 3.2
p3 <- ggboxplot(df, x = "cluster", y = "percent",
color = "group",
palette =c("#00AFBB", "#E7B800", "#FC4E07"),
add = "jitter", shape = "group") +
labs(x=NULL) + theme(legend.position="none", panel.border=element_rect(color='black', fill=NA, size=0.7, linetype="solid"))
p4 <- p1 + stat_pvalue_manual(stat.test, label = "p", tip.length = 0)
p5 <- p1 + stat_pvalue_manual(stat.test, label = "p.signif", tip.length = 0)
p4 + p5
image-20220927222432462.png
//組間基因表達(dá)差異
-
生成測試數(shù)據(jù)
library(Seurat)
library(tidyverse)
library(rstatix)
library(ggpubr)
library(ggplot2)
#生成測試數(shù)據(jù)，添加一個(gè)分組
pbmc <- pbmc_small
g3_cells <- sample(colnames(pbmc), 27)
pbmc[['groups']][g3_cells, ] <- "g3"
unique(pbmc$groups)
-
組間基因差異表達(dá)小提琴圖
#繪制組間所有細(xì)胞中基因表達(dá)差異
data <- pbmc
Idents(data) <- data$groups
levels(data) <- c('g1', 'g2', 'g3')
comparisons <- list(c("g2", "g3"))
p <- VlnPlot(data, features='LYZ', cols=c("#00AFBB", "#E7B800", "#FC4E07")) +
theme(legend.position="none", panel.border=element_rect(color='black', fill=NA, size=0.7, linetype="solid"), axis.text.x=element_text(angle=0, hjust=0.5)) +
labs(x=NULL)
p6 <- p + stat_compare_means(comparisons = comparisons, method = 't.test', label = 'p.signif', size=8) + lims(y=c(0, 9))
comparisons <- list(c("g2", "g3"), c('g1', 'g2'), c('g1', 'g3'))
p7 <- p + stat_compare_means(comparisons = comparisons, method = 't.test', label = 'p.signif', size=8) + lims(y=c(0, 11))
p6 + p7
png('test.png', height=200, width=360)
print(p)
dev.off()
image-20220928204123772.png
#如果橫坐標(biāo)是細(xì)胞類型呢乖酬？
data <- pbmc
data$groups <- factor(data$groups, levels=c('g1', 'g2', 'g3'))
comparisons <- list(c("g2", "g3"), c('g1', 'g3'))
p <- VlnPlot(data, features='LYZ', cols=c("#00AFBB", "#E7B800", "#FC4E07"), split.by='groups') +
theme(legend.position="none", panel.border=element_rect(color='black', fill=NA, size=0.7, linetype="solid"), axis.text.x=element_text(angle=0, hjust=0.5)) +
labs(x=NULL)
#執(zhí)行統(tǒng)計(jì)檢驗(yàn)
df <- p$data
head(df)
# LYZ ident split
#ATGCCAGAACGACT 4.968834e+00 0 g3
#CATGGCCTGTGCAT 4.776148e+00 0 g3
#GAACCTGATGAACC 4.753098e+00 0 g2
#TGACTGGATTCTCA 6.328626e-06 0 g2
#AGTCAGACTGCACA 4.042683e-06 0 g2
#TCTGATACACGTGT 4.968820e+00 0 g3
stat.test <- df %>%
group_by(ident) %>%
t_test(LYZ ~ split, comparisons = comparisons) %>%
adjust_pvalue(method = "bonferroni") %>%
add_significance("p")
stat.test <- stat.test %>%
add_xy_position(x = "ident", dodge = 0.8)
print(stat.test, width=Inf)
p8 <- p + stat_pvalue_manual(stat.test, label = "p.signif", tip.length = 0, step.increase = 0)
p8
#Error in `check_aesthetics()`:
#! Aesthetics must be either length 1 or the same as the data (6): fill
#無解報(bào)錯(cuò)，不知所以融求，試著用ggplot2重畫咬像，還是同樣報(bào)錯(cuò)，ggpubr中的ggviolin就不會報(bào)錯(cuò)
p <- ggviolin(df, x = "ident", y = "LYZ", fill = "split",
palette = c("#00AFBB", "#E7B800", "#FC4E07"),
trim = TRUE) +
labs(x=NULL, fill=NULL)
p8 <- p + stat_pvalue_manual(stat.test, label = "{p} {p.signif}", tip.length = 0, step.increase = 0)
p8
image-20220928221841351.png
