R stackedbarplot fig4

### 16s - frequencies stacked barplots - figure 4
# applied to rarefaction_script output
# paper version



rm(list=ls())


setwd()

freq=read.table("filename_taxotofreq_rarefied.txt", header=T, sep="\t")
colnames(freq)[4]="exp"
for (j in 5:ncol(freq)){
set(freq,which(freq[[j]]<0),j,0)
}


# removal of unclassified sequences
freq$unclassified=NULL

library(data.table, quietly = T)


####################################################################################

library(RColorBrewer)
color_palette=function(freq){
all_table=colnames(freq)[5:ncol(freq)]
color=colorRampPalette(c("darkred","red","orange","yellow","green","blue", "purple","pink"))(length(all_table))
all_color=data.frame(Family=all_table, color)
return(all_color)
}


description_echantillon=function(ech, colors="new", sample_name="none", threshold=0, size=0.7, max=8){
# Sommes des fr?quences des ?chantillons choisis
subfreq=freq[freq$sample.ID%in%ech,]
rownames(subfreq)=subfreq[,1]
data_table=subfreq[,-c(1:4)]

# Proportions
data_prop=apply(data_table,1,function(x){round(100*x/sum(x),2)})
other=data_prop[rowSums(data_prop)<threshold,]
other=colSums(other)
data_prop=data_prop[rowSums(data_prop)>=threshold,]
fam=rownames(data_prop)
data_prop=data_prop[rowSums(data_prop)!=0,]
if(length(colors)==1){
colors=color_palette(subfreq[,c(1:4,match(rownames(data_prop),colnames(subfreq)))])
}
colors=colors[match(rownames(data_prop), colors$Family),]
data_prop=rbind(data_prop,other)
colors$Family=as.character(colors$Family)
colors=rbind(colors, c("other","#838B8B"))
# barplot
par(mfrow=c(1,2),mar=c(10,5,4,5), cex=size)
bp=barplot(data_prop, col=as.character(colors$color), horiz=T, las=1, xlim=c(0,100), width = 0.8, ylim=c(0,max))
if(sample_name=="none"){
sample_title=paste(ech, collapse = " + ")
}else{
sample_title=sample_name
}
title(paste("proportion of each family in : \n", sample_title))
plot.new()
legend("top", legend=colors$Family[order(colors$Family,decreasing=T)], col=as.character(colors$color[order(colors$Family,decreasing=T)]), pch=15, cex=1.5)
}


COLORS=color_palette(freq)
COLORS$color=as.character(COLORS$color)

####################################################################################
## used for paper figures:
#colors adjustment - Orders
COLORS[3,2]="#DC143C"
COLORS[5,2]="#FF8C00"
COLORS[9,2]="#FFD700"
COLORS[11,2]="#ADFF2F"
COLORS[12,2]="#006400"
COLORS[14,2]="#008080"
COLORS[15,2]="#1E90FF"
COLORS[18,2]="#00008B"
COLORS[22,2]="#800080"
COLORS[23,2]="#8A2BE2"
COLORS[24,2]="#9370DB"
COLORS[25,2]="#FF7F50"
COLORS[26,2]="#FF69B4"

#colors adjustment - Families
COLORS[1,2]="#DC143C"
COLORS[5,2]="#F08080"
COLORS[8,2]="#FF8C00"
COLORS[9,2]="#FFD700"
COLORS[16,2]="#ADFF2F"
COLORS[17,2]="#6B8E23"
COLORS[20,2]="#006400"
COLORS[21,2]="#3CB371"
COLORS[24,2]="#20B2AA"
COLORS[27,2]="#008080"
COLORS[30,2]="#00008B"
COLORS[31,2]="#B0E0E6"
COLORS[32,2]="#1E90FF"
COLORS[33,2]="#00CED1"
COLORS[49,2]="#8A2BE2"
COLORS[51,2]="#9370DB"
COLORS[56,2]="#800080"
COLORS[57,2]="#FF69B4"


# graphes
nb_max=8
mix_ech=freq$sample.ID[freq$sample.type=="fly" & freq$level=="wild founders"]
description_echantillon(ech=mix_ech, color=COLORS, sample_name = "wild founder flies samples", threshold = 1, size=0.8, max=nb_max)
mix_ech=freq$sample.ID[freq$sample.type=="food" & freq$level=="wild founders"]
description_echantillon(ech=mix_ech, color=COLORS, sample_name = "wild founder food samples", threshold = 1, size=0.8, max=nb_max)

mix_ech=freq$sample.ID[freq$sample.type=="fly" & freq$level=="yw parents"]
description_echantillon(ech=mix_ech, color=COLORS, sample_name = "yw parents flies samples", threshold = 1, size=0.8, max=nb_max)

mix_ech=freq$sample.ID[freq$sample.type=="fly" & freq$level=="yw progeny"]
description_echantillon(ech=mix_ech, color=COLORS, sample_name = "yw progeny flies samples", threshold = 1, size=0.8, max=nb_max)
mix_ech=freq$sample.ID[freq$sample.type=="food" & freq$level=="yw progeny"]
description_echantillon(ech=mix_ech, color=COLORS, sample_name = "yw progeny food samples", threshold = 1, size=0.8, max=nb_max)