script to generate stacked barplot for community

barplot.R

# cleanup
rm(list=ls())

# load packages
require(vegan)
require(pander)
require(ggplot2)
require(ape)
library(RColorBrewer)

taxa <- read.table("data/community_matrix/hmp1_class2.txt", header=T, sep='\t',skip=8 )
rownames(taxa) <- taxa$SRS
taxa$SRS <- NULL

# transpose that names are row
taxa.t <- as.data.frame(t(taxa))
mapping <- read.csv2("data/patient_metadata/hmp1-II_public_metadata.tsv", sep="\t", header=T)
mapping <- mapping[which(mapping$SRS != "#N/A"),]

grouping_info<-data.frame(row.names=rownames(taxa.t),annot=mapping[match(rownames(taxa.t), mapping$SRS),]$STArea)

# normalize
x <- taxa.t / rowSums(taxa.t)
x<-x[,order(colSums(x),decreasing=TRUE)]

# subset top N taxa
N<-7
taxa_list<-colnames(x)[1:N]
#remove "__Unknown__" and add it to others
taxa_list<-taxa_list[!grepl("Unclassified",taxa_list)]
N<-length(taxa_list)

new_x<-data.frame(x[,colnames(x) %in% taxa_list],Others=rowSums(x[,!colnames(x) %in% taxa_list]))

df<-NULL
for (i in 1:dim(new_x)[2]){
  tmp<-data.frame(row.names=NULL,Sample=rownames(new_x),Taxa=rep(colnames(new_x)[i],dim(new_x)[1]),Value=new_x[,i],Type=grouping_info[,1])
  if(i==1){df<-tmp} else {df<-rbind(df,tmp)}
}

colours <- c("#F0A3FF", "#0075DC", "#993F00","#4C005C","#2BCE48","#FFCC99","#808080","#94FFB5","#8F7C00","#9DCC00","#C20088","#003380","#FFA405","#FFA8BB","#426600","#FF0010","#5EF1F2","#00998F","#740AFF","#990000","#FFFF00");

library(ggplot2)
p<-ggplot(df,aes(Sample,Value,fill=Taxa))+geom_bar(stat="identity")+facet_grid(. ~ Type, drop=TRUE,scale="free",space="free_x")
p<-p+scale_fill_manual(values=colours[1:(N+1)])
p<-p+theme_bw()+ylab("Proportions")
p<-p+ scale_y_continuous(expand = c(0,0))+theme(strip.background = element_rect(fill="gray85"))
p<-p+theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.5))
p <- p + theme(axis.title.x=element_blank(),
               axis.text.x=element_blank(),
               axis.ticks.x=element_blank())
p