【绘图进阶】之配对箱体图绘制(七)

原创生信团上海天昊生物

前面基础课程中介绍过箱体图的绘制，本期介绍成对样品的箱体图绘制，如果您有成对样品的数据，这方面的学习可千万不要错过哦。本期课程从数据处理和数据绘图两个方面介绍，同时增加更细致的文字描述，如果您仍有不理解的地方，留言给我们，我们会及时的回复。

1. 丰度数据和分组信息读取
未处理的数据读取，以data表示；分组信息以info表示。
In [1]:
data = read.table('phylum.taxon.Abundance.xls',header = T,sep = ' ',row.names = 1)
data

Out[1]:

In [2]:
info = read.table('sample_groups.xls',header = F,sep = ' ',col.names = c('sample','group'))
info

Out[2]:

2. 提取数据框
通过grep获取All和Verrucomicrobia的行号
In [3]:
grep("All|Verrucomicrobia", rownames(data))

Out[3]:
10 12

处理后的数据以数据框df表示。- grep("All|Verrucomicrobia", rownames(data)) 用于提取非All和Verrucomicrobia的行，info$sample提取样品列

In [4]:
df = data[-grep("All|Unassigned", rownames(data)),info$sample]
df

Out[4]:

3. 求两组的均值、方差和标准差；并计算成对样品T检验P值。
In [5]:
df['H_mean'] = apply(df[,info[info$group=='H',1]],1,mean)
df['H_var'] = apply(df[,info[info$group=='H',1]],1,var)
df['H_sd'] = apply(df[,info[info$group=='H',1]],1,sd)

df['L_mean'] = apply(df[,info[info$group=='L',1]],1,mean)
df['L_var'] = apply(df[,info[info$group=='L',1]],1,var)
df['L_sd'] = apply(df[,info[info$group=='L',1]],1,sd)

df['p.value'] = apply(df,1,function(x){ t.test(as.numeric(x[info[info$group=='H',1]]),as.numeric(x[info[info$group=='L',1]]),paired = T)$p.value })

df

									

									Out[5]:

数据保存至文本文件phylum.taxon.Abundance.new.xls 中。

In [6]:
write.table(df,'phylum.taxon.Abundance.new.xls',sep = ' ',row.names = T, col.names = NA,quote = F)

4. 提取P值显著的数据
In [7]:
df_diff = df[df$p.value<0.05,info$sample]
df_diff

Out[7]:

1. 数据合并

In [8]:
df_diff = data.frame(t(df_diff),group = info$group)
df_diff

Out[8]:

2. 加载R包

In [9]:
.libPaths("C:/Program Files/R/R-3.6.1/library")
library(ggpubr)

3. 绘图

method包括：t.test、wilcox.test、anova、kruskal.test；本次分析采用t.test。
ggarrange可以将多个图绘制在一张图上，ncol=2 按照每行2张图绘制，行不限制。

In [10]:
p1= ggpaired(df_diff, x="group", y="Bacteroidetes", color = "group", line.color = "gray",xlab = 'Bacteroidetes',
line.size = 0.4, palette = "npg")+ stat_compare_means(method = "t.test",paired = TRUE)

p2= ggpaired(df_diff, x="group", y="Proteobacteria", color = "group", line.color = "gray", xlab = 'Proteobacteria', line.size = 0.4, palette = "npg")+ stat_compare_means(method = "t.test",paired = TRUE)
ggarrange(p1, p2, ncol = 2,common.legend = T)

Out[10]: