Hi colin @colinbrislawn,
Im making a heatmap now. Read your reply, I think use plot_heatmap is more easier than @jbisanz poseted, but I dont know, how to use it in the qiime2R. What‘s the "physeq" you said?
Is it this?
physeq<-qza_to_phyloseq( features="table-dada2-3.qza", tree="rooted-tree.qza", taxonomy="taxonomy.qza", metadata = "qiime2/metadata.tsv" )
Then I want to filter my sample,but it failed
physeq %>%
filter(Source=="Plastic") %>%
tax_glom("Phylum") %>%
plot_heatmap() + facet_grid(~Plastic-type)
ggsave("qiime2/image/heatmap.pdf", height=4, width=20, device="pdf") # save a PDF 3 inches by 4 inches
what' wrong?
I already used the command by @jbisanz to make a heat map, and I changed a little,
like this:
library(tidyverse) library(qiime2R)
setwd(dir="D:/Rdata/micro_plastic_2020/16s/try1")
metadata<-read_q2metadata("qiime2/metadata.tsv")
head(metadata)
taxonomy<-read_qza("taxonomy.qza")$data
head(taxonomy)
SVs<-read_qza("table-dada2-3.qza")$d
head(SVs)
SVs<-apply(SVs, 2, function(x) x/sum(x)*100) #convert to percent
head(SVs)
TSVs<-SVs %>%
as.data.frame() %>%
rownames_to_column("Feature.ID") %>%
gather(-Feature.ID, key="SampleID", value="Abundance")
filter the data, extract top 30 specie
Table <-TSVs %>% left_join(metadata) %>% left_join(taxonomy) %>% filter(`Source`=="Plastic")
TableToPlot<- group_by(Table, Feature.ID) %>%
summarize(MeanAbundance=mean(Abundance)) %>%
arrange(desc(MeanAbundance)) %>%
top_n(30, MeanAbundance) %>%
pull(Feature.ID) #extract only the names from the table
SVs %>%
as.data.frame() %>%
rownames_to_column("Feature.ID") %>%
gather(-Feature.ID, key="SampleID", value="Abundance") %>%
mutate(Feature.ID=if_else(Feature.ID %in% TableToPlot, Feature.ID, "Other")) %>% #flag features to be collapsed
group_by(SampleID, Feature.ID) %>%
summarize(Abundance=sum(Abundance)) %>%
left_join(metadata) %>%
mutate(NormAbundance=log10(Abundance+0.01)) %>% # do a log10 transformation after adding a 0.01% pseudocount. Could also add 1 read before transformation to percent
left_join(taxonomy) %>%
filter(`Source`=="Plastic") %>%
mutate(Feature=paste(Feature.ID, Taxon)) %>%
mutate(Feature=gsub("[kpcofgs]__", "", Feature)) %>% # trim out leading text from taxonomy string %>%
mutate(`Fertilizer-type`=factor(`Fertilizer-type`, levels=c("CK","CF","OF"))) %>%
ggplot(aes(x=`Fertilizer-type`, y=Feature, fill=NormAbundance)) +
geom_tile() +
facet_grid(~`Plastic-type`, scales="free_x") +
theme_q2r() +
theme(axis.text.x=element_text(angle=45, hjust=1)) +
scale_fill_viridis_c(name="log10(% Abundance)")
ggsave("qiime2/image/heatmap.pdf", height=4, width=20, device="pdf") # save a PDF 3 inches by 4 inches
Result
But,I want only display the "Phylum" level, So how should I do ? I just want filter my sample and display the "phylum" level ,with the abundance in heatmap.
Yu.