Phyloseq:相对丰富的otu表和元数据不匹配



我对phyloseq相对陌生,我很难获得一个可接受的相对丰度otu表,用于输入siamcat R代码进行荟萃分析。

# this works: from qza to phyloseq object
ps<-qza_to_phyloseq(
features="all-table.qza",
tree="rooted-tree.qza",
taxonomy = "all-taxonomy.qza",
metadata = "metafinal.tsv"
)
# import metadata
metadata <- read_tsv("metafinal.tsv")
# 30 overlap of the metadata-sample_id with ps, 115 only in metadata
gplots::venn(list(metadata=metadata$sample_id, features=sample_names(ps))
# works: from phyloseq object to relative abundance otu table
table(tax_table(ps)[, "Phylum"])
ps_rel_abund <- transform_sample_counts(ps, function(x){x / sum(x)})
ps_phylum_rel <- tax_glom(ps_rel_abund, "Phylum")
taxa_names(ps_phylum_rel) <- tax_table(ps_phylum_rel)[, "Phylum"]
rel_table <- as(otu_table(ps_phylum_rel), "matrix")
# column names and sample_id are 100% the same
colnames(rel_table)
metadata$sample_id
# 100% overlap:
gplots::venn(list(metadata=metadata$sample_id, featuretable=colnames(rel_table)))
# check that metadata and feature agree
stopifnot(all(colnames(rel_table) == metadata$sample_id))

这里我得到一条错误消息:all(colnames(rel_table(==metadata$sample_id(不是TRUE而下面的siamcat代码根本不起作用。

我的metadata[1:5, 1:5]:sample_id absolute_filepath study experiment_acce…study _title

1 SRR8547628$PWD/Chen_20_da…Chen…SRX5349649 c…2 SRR8547629$PWD/Chen _2020_da…Chen…SRX5349648 c…3 SRR8547630$PWD/Chen _2020_da…Chen…SRX5349647 c…4 SRR8547631$PWD/Chen _2020_da…Chen…SRX5349646 c…5 SRR8547632$PWD/Chen _2020_da…Chen…SRX5349645剖析…

我的rel-table[1:5, 1:5]:SRR5092146芦苇0.0000000 0.00000000 0.00000000脊椎动物0.0000000 0.00000000 0.00000000Apicocomplex 0 0.0000000 0.00000000 0.000000000子囊菌群0.0000000 0.00000000 0.000000000Campilobacterota 0 0.2465222 0.01166882 0.004337051SRR5092150芦苇0.00000000脊椎动物0.00000000Apicocomplex 0.00000000子囊菌门0.00000000Campilobacterota 0.02106281

nrow(metadata)=154ncol(rel_table)=154

拜托,为什么它不起作用?我已经试了好几个星期了,但我无法使代码正常运行。。。

感谢您的时间和帮助。

如果我正确理解你的问题,你会想,为什么元数据和特征表中的样本ID完全重叠,但为什么

stopifnot(all(colnames(rel_table) == metadata$sample_id))

返回CCD_ 5。

我想这是因为你们的样品似乎顺序不同。元数据中的前五个示例是:

SRR8547628
SRR8547629
SRR8547630
SRR8547631
SRR8547632

特征表中的前五个样本是:

SRR5092146
SRR5092147
SRR5092148
SRR5092149
SRR5092150

尝试

stopifnot(all(colnames(rel_table) %in% metadata$sample_id))