结合 OTU 和 tax table 并用 OTU ids (Phyloseq/dada2) 替换实际序列
combine OTU and tax table and replace actual sequences with OTU ids (Phyloseq/dada2)
我按照此处描述的工作流程https://f1000research.com/articles/5-1492/v2 使用示例数据和我自己的数据。这工作正常，但现在我无法生成包含 header 的 OTU table，例如 "OTU00004" 或更好的 "kingdom_phylum_..._Pseudomonas_OTU00004"。我想使用这样的 table 来查找和绘制某个 OTU 在多个样本中的丰度。
我创建了一个名为 ps 的 object，似乎没问题：
ps <- phyloseq(tax_table(taxtab), sample_data(samdf),
otu_table(seqtab, taxa_are_rows = FALSE),phy_tree(fitGTR$tree))
> ps
phyloseq-class experiment-level object
otu_table() OTU Table: [ 454 taxa and 360 samples ]
sample_data() Sample Data: [ 360 samples by 14 sample variables ]
tax_table() Taxonomy Table: [ 454 taxa by 6 taxonomic ranks ]
phy_tree() Phylogenetic Tree: [ 454 tips and 452 internal nodes ]
但 OTU table 中的 header 和分类法 table 中的相应行是实际（此处缩短）序列
> head(otu_table(ps)[1])
GCAAGCGTTACTCGGAATCACTGGGCGTAAAGAGCGCGTAGGCGG#shortened
F3D0 0
> head(tax_table(ps)[1])
Taxonomy Table: [1 taxa by 6 taxonomic ranks]:
Kingdom
GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGCAGGCGGA#shortened "Bacteria"
有没有办法结合otu table和分类学table的信息，并用编号的OTU ids替换序列？我检查了几个 phyloseq 资源和常见问题解答，但我找不到答案。
我想要 table 看起来像这样：
taxonomy_OTU00001 taxonomy_OTU00002 taxonomy_OTU00003
F3D0 #counts #counts #counts
F3D1 #counts #counts #counts
F3D11 #counts #counts #counts
F3D125 #counts #counts #counts
由于此步骤之前的工作流程非常耗时，我不确定如何为该问题提供可重现的示例。
编辑：我根据 dww 的建议生成了一个示例子集。
short_otu2 = short_otu = head(otu_table(ps)[,c(1:6)]) # seq as colnames
short_tax2 = short_tax = tax_table(ps)[colnames(short_otu), ] # seq as rownames
# shorten seqs, must still be unique
colnames(short_otu2) <- substr(colnames(short_otu), 0, 50)
rownames(short_tax2) <- substr(rownames(short_tax), 0, 50)
library(phyloseq)
> dput(short_otu2)
new("otu_table", .Data = structure(c(526L, 375L, 2931L, 994L,
2061L, 419L, 319L, 330L, 1737L, 623L, 1868L, 350L, 402L, 207L,
1880L, 577L, 887L, 303L, 413L, 64L, 838L, 698L, 939L, 484L, 146L,
126L, 496L, 440L, 1183L, 184L, 462L, 37L, 26L, 782L, 271L, 310L
), .Dim = c(6L, 6L), .Dimnames = list(c("F3D0", "F3D1", "F3D11",
"F3D125", "F3D13", "F3D141"), c("GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGCAGGCGGAAGAT",
"GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGCAGGCGGACTCT", "GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGGCTGT",
"GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGGCTTT", "CCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGAGCGTAGGTGGATTGT",
"GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGCCTGC"))), taxa_are_rows = FALSE)
> dput(short_tax2)
new("taxonomyTable", .Data = structure(c("Bacteria", "Bacteria",
"Bacteria", "Bacteria", "Bacteria", "Bacteria", "Bacteroidetes",
"Bacteroidetes", "Bacteroidetes", "Bacteroidetes", "Bacteroidetes",
"Bacteroidetes", "Bacteroidia", "Bacteroidia", "Bacteroidia",
"Bacteroidia", "Bacteroidia", "Bacteroidia", "Bacteroidales",
"Bacteroidales", "Bacteroidales", "Bacteroidales", "Bacteroidales",
"Bacteroidales", "Bacteroidales_S24-7_group", "Bacteroidales_S24-7_group",
"Bacteroidales_S24-7_group", "Bacteroidales_S24-7_group", "Bacteroidaceae",
"Bacteroidales_S24-7_group", NA, NA, NA, NA, "Bacteroides", NA
), .Dim = c(6L, 6L), .Dimnames = list(c("GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGCAGGCGGAAGAT",
"GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGCAGGCGGACTCT", "GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGGCTGT",
"GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGGCTTT", "CCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGAGCGTAGGTGGATTGT",
"GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGCCTGC"), c("Kingdom",
"Phylum", "Class", "Order", "Family", "Genus"))))
这是问题的一部分"replace actual sequences with OTU ids (Phyloseq/dada2)?"
我联系了 phyloseq/dada2 开发人员并根据 Susan Holmes 的回复 (https://github.com/joey711/phyloseq/issues/1030) 我想出了这段代码来用编号的 OTU header 替换扩增子序列.
可以在此处找到进一步的讨论：https://github.com/joey711/phyloseq/issues/213
# this changes the header from the actual sequence to Seq_001, Seq_002 etc
taxa_names(ps)
n_seqs <- seq(ntaxa(ps))
len_n_seqs <- nchar(max(n_seqs))
taxa_names(ps) <- paste("Seq", formatC(n_seqs,
width = len_n_seqs,
flag = "0"), sep = "_")
taxa_names(ps)
将分类法纳入 header 的可能方法如下（接上文）：
# generate a vector containing the full taxonomy path for all OTUs
wholetax <- do.call(paste, c(as.data.frame(tax_table(ps))
[c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus")],
sep = "__")) # to distinguish from "_" within tax ranks
# turn the otu_table into a data.frame
otu_export <- as.data.frame(otu_table(ps))
tmp <- names(otu_export)
# paste wholetax and OTU_ids together
for(i in 1:length(tmp)){
names(tmp)[i] = paste(wholetax[i], tmp[i], sep = "__")
}
# overwrite old names
names(otu_export) <- names(tmp)
> head(otu_export)[5]
# output:
Bacteria__Bacteroidetes__Bacteroidia__Bacteroidales__Bacteroidaceae__Bacteroides__Seq_005
F3D0 146
F3D1 126
F3D11 496
F3D125 440
F3D13 1183
F3D141 184
这还不包括表格之间正确排序的测试！所以请确保粘贴和覆盖是正确的。
这样你就有了一个 data.frame 包含每个分类排名的分类 "splittable"，OTU id，样本名称和一个文件中的计数。但是除了导出文件之外，您仍然维护 phyloseq-structure，其中 OTU_ids link 不同的表，例如 otu_table() 和 tax_table()。另一种方法是向 taxa_names() 命令提供 wholetax 向量，我还没有测试过。
非常欢迎提出改进建议！
我按照此处描述的工作流程https://f1000research.com/articles/5-1492/v2 使用示例数据和我自己的数据。这工作正常，但现在我无法生成包含 header 的 OTU table，例如 "OTU00004" 或更好的 "kingdom_phylum_..._Pseudomonas_OTU00004"。我想使用这样的 table 来查找和绘制某个 OTU 在多个样本中的丰度。
我创建了一个名为 ps 的 object，似乎没问题：
ps <- phyloseq(tax_table(taxtab), sample_data(samdf),
otu_table(seqtab, taxa_are_rows = FALSE),phy_tree(fitGTR$tree))
> ps
phyloseq-class experiment-level object
otu_table() OTU Table: [ 454 taxa and 360 samples ]
sample_data() Sample Data: [ 360 samples by 14 sample variables ]
tax_table() Taxonomy Table: [ 454 taxa by 6 taxonomic ranks ]
phy_tree() Phylogenetic Tree: [ 454 tips and 452 internal nodes ]
但 OTU table 中的 header 和分类法 table 中的相应行是实际（此处缩短）序列
> head(otu_table(ps)[1])
GCAAGCGTTACTCGGAATCACTGGGCGTAAAGAGCGCGTAGGCGG#shortened
F3D0 0
> head(tax_table(ps)[1])
Taxonomy Table: [1 taxa by 6 taxonomic ranks]:
Kingdom
GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGCAGGCGGA#shortened "Bacteria"
有没有办法结合otu table和分类学table的信息，并用编号的OTU ids替换序列？我检查了几个 phyloseq 资源和常见问题解答，但我找不到答案。
我想要 table 看起来像这样：
taxonomy_OTU00001 taxonomy_OTU00002 taxonomy_OTU00003
F3D0 #counts #counts #counts
F3D1 #counts #counts #counts
F3D11 #counts #counts #counts
F3D125 #counts #counts #counts
由于此步骤之前的工作流程非常耗时，我不确定如何为该问题提供可重现的示例。
编辑：我根据 dww 的建议生成了一个示例子集。
short_otu2 = short_otu = head(otu_table(ps)[,c(1:6)]) # seq as colnames
short_tax2 = short_tax = tax_table(ps)[colnames(short_otu), ] # seq as rownames
# shorten seqs, must still be unique
colnames(short_otu2) <- substr(colnames(short_otu), 0, 50)
rownames(short_tax2) <- substr(rownames(short_tax), 0, 50)
library(phyloseq)
> dput(short_otu2)
new("otu_table", .Data = structure(c(526L, 375L, 2931L, 994L,
2061L, 419L, 319L, 330L, 1737L, 623L, 1868L, 350L, 402L, 207L,
1880L, 577L, 887L, 303L, 413L, 64L, 838L, 698L, 939L, 484L, 146L,
126L, 496L, 440L, 1183L, 184L, 462L, 37L, 26L, 782L, 271L, 310L
), .Dim = c(6L, 6L), .Dimnames = list(c("F3D0", "F3D1", "F3D11",
"F3D125", "F3D13", "F3D141"), c("GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGCAGGCGGAAGAT",
"GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGCAGGCGGACTCT", "GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGGCTGT",
"GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGGCTTT", "CCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGAGCGTAGGTGGATTGT",
"GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGCCTGC"))), taxa_are_rows = FALSE)
> dput(short_tax2)
new("taxonomyTable", .Data = structure(c("Bacteria", "Bacteria",
"Bacteria", "Bacteria", "Bacteria", "Bacteria", "Bacteroidetes",
"Bacteroidetes", "Bacteroidetes", "Bacteroidetes", "Bacteroidetes",
"Bacteroidetes", "Bacteroidia", "Bacteroidia", "Bacteroidia",
"Bacteroidia", "Bacteroidia", "Bacteroidia", "Bacteroidales",
"Bacteroidales", "Bacteroidales", "Bacteroidales", "Bacteroidales",
"Bacteroidales", "Bacteroidales_S24-7_group", "Bacteroidales_S24-7_group",
"Bacteroidales_S24-7_group", "Bacteroidales_S24-7_group", "Bacteroidaceae",
"Bacteroidales_S24-7_group", NA, NA, NA, NA, "Bacteroides", NA
), .Dim = c(6L, 6L), .Dimnames = list(c("GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGCAGGCGGAAGAT",
"GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGCAGGCGGACTCT", "GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGGCTGT",
"GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGGCTTT", "CCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGAGCGTAGGTGGATTGT",
"GCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGCCTGC"), c("Kingdom",
"Phylum", "Class", "Order", "Family", "Genus"))))
这是问题的一部分"replace actual sequences with OTU ids (Phyloseq/dada2)?"
我联系了 phyloseq/dada2 开发人员并根据 Susan Holmes 的回复 (https://github.com/joey711/phyloseq/issues/1030) 我想出了这段代码来用编号的 OTU header 替换扩增子序列.
可以在此处找到进一步的讨论：https://github.com/joey711/phyloseq/issues/213
# this changes the header from the actual sequence to Seq_001, Seq_002 etc
taxa_names(ps)
n_seqs <- seq(ntaxa(ps))
len_n_seqs <- nchar(max(n_seqs))
taxa_names(ps) <- paste("Seq", formatC(n_seqs,
width = len_n_seqs,
flag = "0"), sep = "_")
taxa_names(ps)
将分类法纳入 header 的可能方法如下（接上文）：
# generate a vector containing the full taxonomy path for all OTUs
wholetax <- do.call(paste, c(as.data.frame(tax_table(ps))
[c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus")],
sep = "__")) # to distinguish from "_" within tax ranks
# turn the otu_table into a data.frame
otu_export <- as.data.frame(otu_table(ps))
tmp <- names(otu_export)
# paste wholetax and OTU_ids together
for(i in 1:length(tmp)){
names(tmp)[i] = paste(wholetax[i], tmp[i], sep = "__")
}
# overwrite old names
names(otu_export) <- names(tmp)
> head(otu_export)[5]
# output:
Bacteria__Bacteroidetes__Bacteroidia__Bacteroidales__Bacteroidaceae__Bacteroides__Seq_005
F3D0 146
F3D1 126
F3D11 496
F3D125 440
F3D13 1183
F3D141 184
这还不包括表格之间正确排序的测试！所以请确保粘贴和覆盖是正确的。
这样你就有了一个 data.frame 包含每个分类排名的分类 "splittable"，OTU id，样本名称和一个文件中的计数。但是除了导出文件之外，您仍然维护 phyloseq-structure，其中 OTU_ids link 不同的表，例如 otu_table() 和 tax_table()。另一种方法是向 taxa_names() 命令提供 wholetax 向量，我还没有测试过。
非常欢迎提出改进建议！