1.0.1_gene_recovery_heatmap_ggplot.R

#setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) # only for RStudio

#Here, set the file name of the table generated by get_seq_lengths.py
wd=getwd()
sample.filename = file.path(wd, "1_phylo_reconstruction","1.0_hybpiper","seq_lengths.txt")

#Calculate sizes for gene and sample labels, increase the multiplier to make text bigger
gene.size.multiplier = .03
sample.size.multiplier = .5

#Below is an example of a ggplot save command
#ggsave("plastid_heatmap.png", height = 4, width = 6, units = "in")

#You should not need to change anything below this line!
#-------------------------------------------------------

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

library(ggplot2)
library(reshape2)

sample.data= as.matrix(read.table(sample.filename,header=T,row.names=1,sep="\t"))
sample.len = sample.data[2:nrow(sample.data),]
reference.len = as.numeric(sample.data[1,])

#Calculate the percentage length recovered relative to the reference.
percent.len=sweep(sample.len,2,as.numeric(reference.len),'/')
percent.len = ifelse(percent.len>1,1,percent.len)
percent.long = melt(percent.len)
percent.long$Var1 = as.factor(percent.long$Var1)

gene.size = dim(percent.len)[2] * gene.size.multiplier
sample.size = dim(percent.len)[1] * sample.size.multiplier

pdf(file.path(wd,"plots","hybpiper_heatmap.pdf"),width=12)
ggplot(data = percent.long, aes(x=Var2, y=Var1, fill = value))+
  geom_raster()+
  #guides(fill=FALSE)+ #remove this line if you want the heat scale to appear
  scale_fill_gradient(high = "#132B43", low = "#FFFFFF")+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  ylab(NULL)+xlab(NULL)+
  theme(axis.text.y=element_text(face="italic",size = sample.size),axis.text.x =element_text(size=gene.size))
dev.off()