|
微生物生態(tài)學(xué)領(lǐng)域經(jīng)常要使用三種差異分析方法,，分別似乎MRPP，anosim,，adonis,。很多人也問，到底使用哪種方法比較好,，當(dāng)然也有一些博文對他們有過簡單的區(qū)別，但是目前沒什么大的不同,。這三種方法是可以一起使用的,。 不僅僅要一起使用，還要再有多個(gè)分組的情況下,，進(jìn)行兩兩比較,。那這便是咱們這次的故事了： pairMicroTest(ps = ps,dist = "bray",Micromet = "MRPP")
運(yùn)行結(jié)果 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)
#-準(zhǔn)備矩陣和分組文件
map = as.data.frame(sample_data(ps1_rela))
aa = levels(map$Group)
aa
aaa = combn(aa,2)
aaa
dim(aaa)[2]
# 構(gòu)建三個(gè)空列
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)
}
我雖然使用了替代的方法可以繼續(xù)按想法實(shí)現(xiàn)這個(gè)功能，但是%in%的問題并沒有解決，不僅是這個(gè)函數(shù),，還有我曾今推送的maptre中的函數(shù)： # 提取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
本來極提取子集使用一下函數(shù)是最簡單的了,，但是還是一樣的問題，這個(gè)函數(shù)再function中不起作用,。兩個(gè)命令相似的特征都是使用了%in%符號,。 ps <- ps %>%
subset_taxa(
Kingdom %in% "Bacteria"
)
所以%in%符號和function之間有什么隱情呢? 希望可以得到大家的幫助。
|
