🤔 艾录叔 | 多组比较及绘图

生信漫卷

发布于 2022-10-31 17:08:13

3140

发布于 2022-10-31 17:08:13

文章被收录于专栏：R语言及实用科研软件R语言及实用科研软件

1写在前面

写毕业课题统计时编写的一段代码，大量数据很快就可以统计出结果并作用，方便的很。统计使用的是r基础stat包，绘图使用的ggplot2包。都是很常见的，网上教程也很多。

2示例数据

示例为利用excel随机生成的一列数字

y <- read.table("clipboard", header = F)

该法是直接访问的剪贴板，可以用read.xlsx，read.table，read.csv等函数读取已经整理好的数据。分组信息也可读入。

3输入分组信息

a1 <- factor(c(rep(c('模型组','低剂组','高剂组'),
                  each = 5)))
a2 <- factor(c(rep(c('空白组','模型组','治疗组'),
                  time = 3)))

each和time参数的不同，可以自己运行一下。

4计算各组均值和标准差

mean <- tapply(y$V1,a1,mean)
sd <- tapply(y$V1,a1,sd)
data <- data.frame(mean,sd)
#排序
data$group <- factor(row.names(data), levels = c('模型组','低剂组','高剂组'))

5正态性检验

tapply(y$V1,a1,chisq.test)

6方差齐性检验

bartlett.test(y$V1,a1)

7单因素方差分析

result.aov<-aov (y$V1~a1)
summary(result.aov)
# 两两比较
pairwise.t.test(y$V1,a1,p.adjust.method = 'bonferroni')

8另一种多重比较方法

TukeyHSD(result.aov)
plot(TukeyHSD(result.aov))

9非参数检验

kruskal.test(y$V1,a1)
#两两比较
pairwise.wilcox.test(y$V1,a1,p.adjust.method = "bonferroni")

10使用sink()函数获取结果

sink("calc.txt")
print(data)
print(tapply(y$V1,a1,chisq.test))
print(bartlett.test(y$V1,a1))
#多组数据符合正态分布，且方差齐则使用单因素方差分析。
result.aov<-aov (y$V1~a1)
print(summary(result.aov))
print(pairwise.t.test(y$V1,a1,p.adjust.method = 'bonferroni'))
#多组数据不符合正态分布；或虽符合正态分布，但方差不齐，则使用非参数检验。
print(kruskal.test(y$V1,a1))
print(pairwise.wilcox.test(y$V1,a1,p.adjust.method = 'bonferroni'))
sink()

输出结果

      mean       sd  group
低剂组 14.2 2.167948 低剂组
高剂组 25.0 1.581139 高剂组
模型组  7.2 3.420526 模型组
$低剂组

	Chi-squared test for given probabilities

data:  X[[i]]
X-squared = 1.3239, df = 4, p-value = 0.8573


$高剂组

	Chi-squared test for given probabilities

data:  X[[i]]
X-squared = 0.4, df = 4, p-value = 0.9825


$模型组

	Chi-squared test for given probabilities

data:  X[[i]]
X-squared = 6.5, df = 4, p-value = 0.1648



	Bartlett test of homogeneity of variances

data:  y$V1 and a1
Bartlett's K-squared = 2.1535, df = 2, p-value = 0.3407

            Df Sum Sq Mean Sq F value   Pr(>F)    
f            2  804.1   402.1   63.82 4.03e-07 ***
Residuals   12   75.6     6.3                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

	Pairwise comparisons using t tests with pooled SD 

data:  y$V1 and a1

       低剂组  高剂组 
高剂组 5.7e-05 -      
模型组 0.0026  3.1e-07

P value adjustment method: bonferroni 
	Kruskal-Wallis rank sum test

data:  y$V1 and a1
Kruskal-Wallis chi-squared = 12.545, df = 2, p-value = 0.001888

	Pairwise comparisons using Wilcoxon rank sum test with continuity correction 

data:  y$V1 and a1 

       低剂组 高剂组
高剂组 0.036  -     
模型组 0.035  0.036 

P value adjustment method: bonferroni

11绘制条形图加误差线

library(ggplot2)
p<-ggplot(data,aes(group,mean))+
  geom_col(
    fill = c('white','black','gray80'),
    col ='black',
    width = 0.6,
    position = position_dodge(0.8),
    lwd = 1.5,
  )+
  theme_classic(
    base_size = 32,
    base_family = 'serif'#字体
  )+
  scale_y_continuous(
    expand = c(0,0)
  )+
  scale_x_discrete(
    expand = c(0.2,0.2)
  )+
  geom_blank(aes(y= (mean+sd)*1.2))+
  geom_errorbar(
    aes(ymin = mean, ymax = mean+sd),
    width = 0.3,lwd = 1.5
  )+
  xlab("这是x轴")+
  ylab("这是y轴")
p

12把图保存下来

tiff('barplot.tif')
p
dev.off()

13小提琴图，加误差线，不要图例

小提琴图和箱线图用到的是所有数据，需要构建包含所有数据的表格。

data2 <- as.data.frame(y$V1)
data2$fenzu <- a1
#注意两个表中的同样数据列名一致，要不会报错。
colnames(data2) <-c('mean','group')
#排序
data2$group <- factor(data2$group, levels = c('模型组','低剂组','高剂组'))

p2<-ggplot(data = data2,aes(group,mean,fill=group))+
  geom_violin(trim=FALSE,color="black",
              lwd = 1.5
  )+
  theme_classic(
    base_size = 32,
    base_family = 'serif'
  )+
  geom_blank(data = data, aes(y= (mean+sd)*1.2))+
  geom_point(data = data, aes(group,mean),pch = 18, size = 6)+
  geom_errorbar(data = data,
                aes(ymin = mean-sd, ymax = mean+sd),
    width = 0.1,lwd = 1.5
  )+
  xlab("这里是x轴")+
  ylab("这里是y轴")+
  theme(legend.position="none")
p2

14小提琴图加箱线图

p3<- ggplot(data2, aes(group, mean ,fill = group))+
  geom_violin(trim=FALSE,color="black",
              lwd = 1.5
  )+
  geom_boxplot(width=0.2,position=position_dodge(0.9),lwd= 1.5)+
  theme_classic(
    base_size = 32,
    base_family = 'serif'
  )+
  geom_blank(data = data, aes(y= (mean+sd)*1.2))+
  xlab("这里是x轴")+
  ylab("这里是y轴")+
  theme(legend.position="none")
p3