新闻媒体-太阳成tyc7111cc

发稿时间：2020-07-29来源：天昊生物

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

1. 丰度数据和分组信息读取

未处理的数据读取，以data表示；分组信息以info表示。

in [1]:

data = read.table('phylum.taxon.abundance.xls',header = t,sep = 't',row.names = 1)
data

out[1]:

in [2]:

info = read.table('sample_groups.xls',header = f,sep = 't',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 = 't',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]:

往期相关链接：

1、r基础篇

；；
；；
；
；
；；
；

2、r进阶

；

；

；

；

；

；

3、数据提交

;

；

；

；

4、表达谱分析

；；；

5、医学数据分析

；；；；

天昊客户服务中心

手机/微信号：18964693703