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之间有什么隐情呢? 希望可以得到大家的帮助。

赞 (0)

相关推荐