我在phyloseq
工作中使用了很多。我的数据集通常包含多个条件或参数,需要以相同的方式进行分析(例如,夏季或冬季的细菌和 Lake1 或 Lake2 的相同图),所以我想为此使用函数。我写了一个子集函数,它允许我通过循环组合多个参数。输出存储在列表中以供进一步分析。
但是,这似乎很笨拙。所以我的第一个问题是关于功能的改进。
1)具体来说,我想知道是否
a) 使用多个for loops
来生成子集是一个好主意。
b) 此外,可以优化for loops
和的组合。lapply
和
c)也许有更好的方法来防止现有列表再次无法识别地附加相同对象的新迭代?我实现了这一点,因为我在开发代码时有很多很多测试执行
这里讨论了 for 循环是否比 apply 慢:lapply vs for loop - Performance R
我认为phyloseq
内部调用which
,所以它不必是一个phyloseq
特定的解决方案。
2)我的第二个问题是如何处理这种情况,如果不是所有搜索参数都存在于所有子集中?所以在下面的例子中,如果没有丹麦男性,“丹麦”和“M”的组合就会中断。我想避免这种情况,在这个例子中只有 3 个(丹麦 x F,美国 x F,美国 x M)而不是 4 个子集。目前,该函数需要适应每个特殊的子集,这首先破坏了编写它的目的。
library(phyloseq)
data(enterotype)
# reduce the size of the data set
phyloseq <- filter_taxa(enterotype, function (x) {sum(x > 0.001) >= 1}, prune = TRUE)
# arguments for the subsetting function
phyloseq_object <- phyloseq
Nationality <- c("american", "danish")
Gender <- c("F", "M")
# define a function to obtain sample subsets from the phyloseq object
# per combination of parameters
get_sample_subsets <- function(phyloseq_object, nation, gender) {
sample_subset <- sample_data(phyloseq_object)[ which(sample_data(phyloseq_object)$Nationality == nation &
sample_data(phyloseq_object)$Gender == gender),]
phyloseq_subset <- merge_phyloseq(tax_table(phyloseq_object),
otu_table(phyloseq_object),
#refseq(phyloseq_object),
sample_subset)
phyloseq_subset2 <- filter_taxa(phyloseq_subset, function (x) {sum(x > 0) >= 1 }, prune = TRUE)
return(phyloseq_subset2)
}
# here we pass the arguments for subsetting over two for loops
# to create all possible combinations of the subset parameters etc.
# the subsets are stored within a list, which has to be empty before running the loops
sample_subset_list <- list()
if(length(sample_subset_list) == 0) {
for (nations in Nationality) {
for (gender in Gender) {
tmp <- get_sample_subsets(phyloseq_object = phyloseq_object,
nation = nations, gender = gender)
sample_subset_list[[paste(nations, gender, sep = "_")]] <- tmp
}
}
print(sample_subset_list)
} else {
print("list is not empty, abort to prevent appending...")
}
# You could now for example use the output to calculate ordinations for each subset (this data set has too few entries per subset for that)
# create a list where the distance metrics for the sample subsets are stored
ordination_nmds <- list()
ordination_nmds <- lapply(sample_subset_list, ordinate, method = "NMDS",
dist = "bray", try = 100, autotransform = TRUE)