Hi! I'm pretty new to R, and am having trouble editing and visualizing a heatmap I've made. I used the following code to create the heatmap you can see below:
library(tidyverse)
library(qiime2R)
setwd("/Users/jakecohen/Desktop/QiimeOutput/2019_16S")
metadata<-read_q2metadata("Qiime_16S_merged copy.tsv")
SVs<-read_qza("16S_merged2.qza")$data
taxonomy<-read_qza("merged_16S_data.qza")$data
SVs<-apply(SVs, 2, function(x) x/sum(x)*100) #convert to percent
SVsToPlot<-
data.frame(MeanAbundance=rowMeans(SVs)) %>% #find the average abundance of a SV
rownames_to_column("Feature.ID") %>%
arrange(desc(MeanAbundance)) %>%
top_n(10, 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% SVsToPlot, Feature.ID, "Remainder")) %>% #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) %>%
mutate(Feature=paste(Feature.ID, Taxon)) %>%
mutate(Feature=gsub("[kpcofgs]__", "", Feature)) %>% # trim out leading text from taxonomy string
ggplot(aes(x=SampleID, y=Feature, fill=NormAbundance)) +
geom_tile() +
facet_grid(~Depth
, scales="free_x") +
theme_q2r() +
theme(axis.text.x=element_text(angle=45, hjust=1)) +
scale_fill_viridis_c(name="log10(% Abundance)")
ggsave("heatmap.pdf", height=4, width=11, device="pdf") # save a PDF 3 inches by 4 inches
results<-results %>% mutate(Significant=if_else(we.eBH<0.1,"*", ""))
tree<-drop.tip(tree, tree$tip.label[!tree$tip.label %in% results$Feature.ID]) # remove all the features from the tree we do not have data for
ggtree(tree, layout="circular") %<+% results +
geom_tippoint(aes(fill=diff.btw), shape=21, color="grey50") +
geom_tiplab2(aes(label=Significant), size=10) +
scale_fill_gradient2(low="darkblue",high="darkred", midpoint = 0, mid="white", name="log2(fold-change") +
theme(legend.position="right")
ggsave("tree.pdf", height=10, width=10, device="pdf", useDingbats=F)
However, the long names of the x axis titles is making it hard to view the data. Is it possible to write a line of code to shorten the x axis labels? I'd like to remove the OTU labels (the string of characters at the beginning). Let me know if that's possible! The line "mutate(Feature=gsub("[kpcofgs]__", "", Feature)) %>% # trim out leading text from taxonomy string" is supposed to shorten the axis titles but does not seem to filtering everyting.
Thanks!