利用phyloseq可视化结果
phyloseq：https://github.com/joey711/phyloseq加载 phyloseq 包library(phyloseq);packageVersion("phyloseq")
library(Biostrings);packageVersion("Biostrings")
library(ggplot2); packageVersion("ggplot2")
theme_set(theme_bw())从 dada2 结果中生成 phyloseq 输入文件
#metadata
samples.out <- rownames(seqtab.nochim)
samples.out
subject <- sapply(strsplit(samples.out, "D"), `[`, 1)
subject
gender <- substr(subject,1,1)
subject <- substr(subject,2,999)
day <- as.integer(sapply(strsplit(samples.out, "D"), `[`, 2))
samdf <- data.frame(Subject=subject, Gender=gender, Day=day)
samdf$When <- "Early"
samdf$When <- "Late"
rownames(samdf) <- samples.out
samdf#读入文件，otu_table,sample_data,tax_table
ps <- phyloseq(otu_table(seqtab.nochim, taxa_are_rows=FALSE),
sample_data(samdf),
tax_table(taxid),
phy_tree(fitGTR$tree))
#方法二：
ps1 <- phyloseq(otu_table(seqtab.nochim, taxa_are_rows=FALSE),
sample_data(samdf),
tax_table(taxa))
library("ape")
random_tree = rtree(ntaxa(ps1), rooted=TRUE, tip.label=taxa_names(ps1))
plot(random_tree)
ps = merge_phyloseq(ps1, random_tree)
#去除Mock
ps <- prune_samples(sample_names(ps) != "Mock", ps) # Remove mock sample
dna <- Biostrings::DNAStringSet(taxa_names(ps))
dna
names(dna) <- taxa_names(ps)
names(dna)
ps <- merge_phyloseq(ps, dna)
ps
taxa_names(ps) <- paste0("ASV", seq(ntaxa(ps)))
psalpha 多样性
plot_richness(ps, x="Day", measures=c("Shannon", "Simpson"), color="When")
# Transform data to proportions as appropriate for Bray-Curtis distances
ps.prop <- transform_sample_counts(ps, function(otu) otu/sum(otu))
ord.nmds.bray <- ordinate(ps.prop, method="NMDS", distance="bray")
plot_ordination(ps.prop, ord.nmds.bray, color="When", title="Bray NMDS")条形图
top20 <- names(sort(taxa_sums(ps), decreasing=TRUE))
ps.top20 <- transform_sample_counts(ps, function(OTU) OTU/sum(OTU))
ps.top20 <- prune_taxa(top20, ps.top20)
plot_bar(ps.top20, x="Day", fill="family") + facet_wrap(~When, scales="free_x")画树
plot_tree(ps)
plot_tree(ps, "treeonly")
plot_tree(ps, "treeonly", nodeplotblank)
plot_tree(ps, "treeonly", nodeplotblank, ladderize="left")
plot_tree(ps, "treeonly", nodeplotblank, ladderize=TRUE)
plot_tree(ps, nodelabf=nodeplotblank, label.tips="taxa_names", ladderize="left")
# plot_tree(ps, "anythingelse")
# plot_tree(ps, nodelabf=nodeplotboot(), ladderize="left", color="SampleType")
# plot_tree(ps, nodelabf=nodeplotboot(), ladderize="left", color="Class")
# plot_tree(ps, nodelabf=nodeplotboot(), ladderize="left", color="SampleType", shape="Class")
# # The default
# plot_tree(ps, color="SampleType", ladderize="left")
页:
[1]