私は phyloseq に比較的慣れていないため、メタ分析用の siamcat R コードへの入力に適した相対存在量の otu-table を取得するのに苦労しています。
# 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) is not TRUE and the following siamcat code is not working at all.
my metadata[1:5, 1:5]: sample_id absolute_filepath study experiment_acce… study_title
1 SRR8547628 $PWD/Chen_2020_da… Chen… SRX5349649 c の解剖… 2 SRR8547629 $PWD/Chen_2020_da… Chen… SRX5349648 c の解剖… 3 SRR8547630 $PWD/Chen_2020_da… SRX5349648 の解剖… 3 SRX79… Chen 4 SRR8547631 $PWD/Chen_2020_da… Chen… SRX5349646 c の解剖… 5 SRR8547632 $PWD/Chen_2020_da… Chen… SRX5349645 c の解剖…</p>
my rel-table[1:5, 1:5]: SRR5092146 SRR5092147 SRR5092148 SRR5092149 Phragmoplastophyta 0 0.0000000 0.00000000 0.000000000 Vertebrata 0 0.0000000 0.00000000 0.000000000 Apicomplexa 0 0.0000000 0.00000000 0.000000000 Ascomycota 0 0.0000000 0.00000000 0.000000000 Campilobacterota 0 0.2465222 0.01166882 0.004337051 SRR5092150 Phragmoplastophyta 0.00000000 Vertebrata 0.00000000 Apicomplexa 0.00000000 Ascomycota 0.00000000 Campilobacterota 0.02106281
nrow(metadata)= 154
ncol(rel_table)= 154
なぜ機能しないのですか？何週間も試しましたが、コードを正しく実行できません...
お時間をいただき、ありがとうございました。