利用phyloseq可视化结果

生信喵

phyloseq：

1. https://github.com/joey711/phyloseq

复制代码

加载 phyloseq 包

1. library(phyloseq);packageVersion("phyloseq")
   
2. library(Biostrings);packageVersion("Biostrings")
   
3. library(ggplot2); packageVersion("ggplot2")
   
4. theme_set(theme_bw())

复制代码

从 dada2 结果中生成 phyloseq 输入文件

1. #metadata
   
2. samples.out <- rownames(seqtab.nochim)
   
3. samples.out
   
4. 
   
5. subject <- sapply(strsplit(samples.out, "D"), `[`, 1)
   
6. subject
   
7. gender <- substr(subject,1,1)
   
8. subject <- substr(subject,2,999)
   
9. day <- as.integer(sapply(strsplit(samples.out, "D"), `[`, 2))
   
10. samdf <- data.frame(Subject=subject, Gender=gender, Day=day)
    
11. samdf$When <- "Early"
    
12. samdf$When[samdf$Day>100] <- "Late"
    
13. rownames(samdf) <- samples.out
    
14. samdf

复制代码

1. #读入文件，otu_table,sample_data,tax_table
   
2. ps <- phyloseq(otu_table(seqtab.nochim, taxa_are_rows=FALSE),
   
3.                sample_data(samdf),
   
4.                tax_table(taxid),
   
5.               phy_tree(fitGTR$tree))
   
6. #方法二：
   
7. ps1 <- phyloseq(otu_table(seqtab.nochim, taxa_are_rows=FALSE),
   
8.                sample_data(samdf),
   
9.                tax_table(taxa))
   
10. library("ape")
    
11. random_tree = rtree(ntaxa(ps1), rooted=TRUE, tip.label=taxa_names(ps1))
    
12. plot(random_tree)
    
13. 
    
14. ps = merge_phyloseq(ps1, random_tree)
    
15. 
    
16. 
    
17. 
    
18. #去除Mock
    
19. ps <- prune_samples(sample_names(ps) != "Mock", ps) # Remove mock sample
    
20. 
    
21. 
    
22. dna <- Biostrings::DNAStringSet(taxa_names(ps))
    
23. dna
    
24. names(dna) <- taxa_names(ps)
    
25. names(dna)
    
26. ps <- merge_phyloseq(ps, dna)
    
27. ps
    
28. taxa_names(ps) <- paste0("ASV", seq(ntaxa(ps)))
    
29. ps

复制代码

alpha 多样性

1. plot_richness(ps, x="Day", measures=c("Shannon", "Simpson"), color="When")
   
2. 
   
3. # Transform data to proportions as appropriate for Bray-Curtis distances
   
4. ps.prop <- transform_sample_counts(ps, function(otu) otu/sum(otu))
   
5. ord.nmds.bray <- ordinate(ps.prop, method="NMDS", distance="bray")
   
6. plot_ordination(ps.prop, ord.nmds.bray, color="When", title="Bray NMDS")

复制代码

条形图

1. top20 <- names(sort(taxa_sums(ps), decreasing=TRUE))[1:20]
   
2. ps.top20 <- transform_sample_counts(ps, function(OTU) OTU/sum(OTU))
   
3. ps.top20 <- prune_taxa(top20, ps.top20)
   
4. plot_bar(ps.top20, x="Day", fill="family") + facet_wrap(~When, scales="free_x")

复制代码

画树

1. plot_tree(ps)
   
2. plot_tree(ps, "treeonly")
   
3. plot_tree(ps, "treeonly", nodeplotblank)
   
4. plot_tree(ps, "treeonly", nodeplotblank, ladderize="left")
   
5. plot_tree(ps, "treeonly", nodeplotblank, ladderize=TRUE)
   
6. plot_tree(ps, nodelabf=nodeplotblank, label.tips="taxa_names", ladderize="left")
   
7. 
   
8. # plot_tree(ps, "anythingelse")
   
9. # plot_tree(ps, nodelabf=nodeplotboot(), ladderize="left", color="SampleType")
   
10. # plot_tree(ps, nodelabf=nodeplotboot(), ladderize="left", color="Class")
    
11. # plot_tree(ps, nodelabf=nodeplotboot(), ladderize="left", color="SampleType", shape="Class")
    
12. # # The default
    
13. # plot_tree(ps, color="SampleType", ladderize="left")

复制代码