R语言种%in%的错误到底是什么原因？ / 开普饭

写在前面

微生物生态学领域经常要使用三种差异分析方法，分别似乎MRPP，anosim，adoni。很多人也问，到底使用哪种方法比较好，当然也有一些博文对他们有过简单的区别，但是目前没什么大的不同。这三种方法是可以一起使用的。

不仅仅要一起使用，还要再有多个分组的情况下，进行两两比较。那这便是咱们这次的故事了：

两两比较最重要的是待比较两组的样本otu子表格的选择，再phyloseq中使用：

ps_sub <- phyloseq::subset_samples(ps1_rela,Group %in% c("KO"，"WT");ps_sub

以上代码非常方便，直接指定分组标签就可以嫁给你otu表格在内的全部分析文件取子集和对齐。

因为有多个组，所以我添加了for循环，这也是没有问题的：

for (i in 1:dim(aaa)[2]) { print(i) Desep_group = aaa[,i] map = as.data.frame(sample_data(ps1_rela)) head(map) ps_sub <- phyloseq::subset_samples(ps1_rela,Group %in% Desep_group);ps_sub

ps_sub = phyloseq::filter_taxa(ps_sub, function(x) sum(x ) > 0 , TRUE);ps_sub map = as.data.frame(sample_data(ps_sub)) unif <- phyloseq::distance(ps_sub , method=dist, type="samples") }

但是当我封装为功能的时候却出现了错误：我添加print(ps_sub)参数查看，原来是没有对phyloseq对象取子集，这条命令没有起作用：ps_sub <- phyloseq::subset_samples(ps1_rela,Group %in% Desep_group);ps_sub

但是如果不在function中使用，这是没问题的。

function(ps = ps){ ps1_rela = transform_sample_counts(ps, function(x) x / sum(x) );ps1_rela library(vegan) #-准备矩阵和分组文件 map = as.data.frame(sample_data(ps1_rela))


aa = levels(map$Group)

aa

aaa = combn(aa,2)

aaa

dim(aaa)[2]
# 构建三个空列

ID = rep("a",dim(aaa)[2])

R = rep("a",dim(aaa)[2])

P = rep("a",dim(aaa)[2])

for (i in 1:dim(aaa)[2]) {

print(i)

Desep_group = aaa[,i]

map = as.data.frame(sample_data(ps1_rela))

head(map)

ps_sub <- phyloseq::subset_samples(ps1_rela,Group %in% Desep_group);ps_sub
ps_sub = phyloseq::filter_taxa(ps_sub, function(x) sum(x ) > 0 , TRUE);ps_sub

map = as.data.frame(sample_data(ps_sub))

unif <- phyloseq::distance(ps_sub , method=dist, type="samples")

print(ps_sub)

}

}

所以我将这条代码:

ps_sub <- phyloseq::subset_samples(ps1_rela,Group %in% Desep_group);ps_sub

写了个个相同功能的作为替换：

map$ID = row.names(map) maps<- dplyr::filter(map,Group %in% Desep_group) row.names(maps) = maps$ID ps_sub = ps1_rela sample_data( ps_sub ) = maps

这次就不会出现问题了，也成功运行了这个函数。

三种方法两两比较函数使用以上策略

试运行：

pairMicroTest(ps = ps,dist = "bray",Micromet = "MRPP")

运行结果

ID stat p 1 KO_VS_OE MRPP.delta 0.3 p: 0.002 2 KO_VS_WT MRPP.delta 0.295 p: 0.002 3 OE_VS_WT MRPP.delta 0.275 p: 0.008

主体函数

pairMicroTest = function(ps = ps,dist = "bray",Micromet = "anosim"){


ps1_rela = transform_sample_counts(ps, function(x) x / sum(x) );ps1_rela

library(vegan)

#-准备矩阵和分组文件

map = as.data.frame(sample_data(ps1_rela))
aa = levels(map$Group)

aa

aaa = combn(aa,2)

aaa

dim(aaa)[2]
# 构建三个空列

ID = rep("a",dim(aaa)[2])

R = rep("a",dim(aaa)[2])

P = rep("a",dim(aaa)[2])

# i = 2

for (i in 1:dim(aaa)[2]) {

print(i)

Desep_group = aaa[,i]

map = as.data.frame(sample_data(ps1_rela))

head(map)

map$ID = row.names(map)

maps<- dplyr::filter(map,Group %in% Desep_group)

row.names(maps) = maps$ID

ps_sub = ps1_rela

sample_data( ps_sub ) = maps

# ps_sub <- phyloseq::subset_samples(ps1_rela,Group %in% Desep_group);ps_sub
ps_sub = phyloseq::filter_taxa(ps_sub, function(x) sum(x ) > 0 , TRUE);ps_sub

map = as.data.frame(sample_data(ps_sub))

unif <- phyloseq::distance(ps_sub , method=dist, type="samples")

print(ps_sub)
if (Micromet == "MRPP") {

mrpp = vegan::mrpp(unif, map$Group)

as = round(mrpp$delta,3)

R2 <- paste("MRPP.delta ",as, sep = "")

R2

p_v = paste("p: ",round(mrpp$Pvalue,3), sep = "")

p_v
}
if (Micromet == "anosim") {

dat.ano = anosim(unif, map$Group)

a = round(dat.ano$statistic,3)

R2 <- paste("ANOSIM.r ",a, sep = "")

p_v = paste("p: ",round(dat.ano$signif,3), sep = "")
}

if (Micromet == "adonis") {

ado = vegan::adonis(unif ~ map$Group,permutations = 999)

a = round(as.data.frame(ado$aov.tab[5])[1,1],3)

R2 <- paste("adonis:R ",a, sep = "")

b = as.data.frame(ado$aov.tab[6])[1,1]

p_v = paste("p: ",b, sep = "")

}

ID[i] = paste(Desep_group[1],Desep_group[2],sep = "_VS_")

P[i] = p_v

R[i] = R2

}
result = data.frame(ID = ID,stat = R,p = P)

head(result)
return(result)

}

写在后面

我虽然使用了替代的方法可以继续按想法实现这个功能，但是%in%的问题并没有解决，不仅是这个函数，还有我曾今推送的maptre中的函数：

# 提取otu表格 otu_table = as.data.frame(t(vegan_otu(ps1_rela))) otu_table$mean = rowMeans(otu_table) otu_table$ID = row.names(otu_table) head(otu_table) #按照从大到小排序 otu_table<- arrange(otu_table, desc(mean)) subtab = head(otu_table,N) head(subtab) row.names(subtab) =subtab$ID subtab = subtab[,1:(dim(subtab)[2]-2)] subtab = as.matrix(subtab)

#对phyloseq取子集 ps_sub <- phyloseq(otu_table(subtab, taxa_are_rows=TRUE), tax_table(tax)) ps_sub

本来极提取子集使用一下函数是最简单的了，但是还是一样的问题，这个函数再function中不起作用。两个命令相似的特征都是使用了%in%符号。

ps <- ps %>% subset_taxa( Kingdom %in% "Bacteria" )

所以%in%符号和function之间有什么隐情呢? 希望可以得到大家的帮助。

R语言种%in%的错误到底是什么原因？

相关推荐