library(ggplot2)
library(plyr)
library(magrittr)
library(dplyr)

args = commandArgs(trailingOnly=TRUE)

lm_eqn = function(plot)
{
    m = lm(count ~ conc, plot);
#    eq = substitute(italic(y) == a + b %.% italic(x)*","~~italic(R)^2~"="~r2,
#    list(a = format(coef(m)[1], digits = 2), b = format(coef(m)[2], digits = 2),r2 = format(summary(m)$r.squared, digits = 3)))
    eq = paste0("r2: ", format(summary(m)$r.squared, digits = 3))
#    as.character(as.expression(eq));
    eq
}

cbPalette <- c("#22635b", "#87bd5b")
plot <- read.table(args[1], header=F, as.is=T, sep="\t")
colnames(plot)<-c("sample","plat", "spikeIn","count", "conc", "type", "lab", "category")

#plot[ , c("lab", "count")] %>% group_by(lab) %>% summarize( zeroes = n())
missing = plot[ plot$count == 0, c("lab", "count")] %>% group_by(lab) %>% summarize( zeroes = n())

plot$count[plot$count == 0] = 0.0001
plot$count=log10(plot$count)
plot$conc=log10(plot$conc)

forE=plot[plot$count > -2.5,]
finalE=forE[forE$conc > -2.5,]

# text to annotate graphs
graphLabels = data.frame( matrix(nrow = length(unique(plot$lab)), ncol = 2) )
colnames(graphLabels) = c("lab", "label")

index = 1
for (i in unique(plot$lab))
{
    graphLabels[index, "lab"] = i
#    graphLabels[index, "label"] = paste0(as.character(as.expression(lm_eqn(finalE[finalE$lab == i, ]))), ", zeroes: ", missing[missing$lab == i, "zeroes"])   #r2 computed without zeroes
    graphLabels[index, "label"] = paste0(as.character(as.expression(lm_eqn(plot[plot$lab == i, ]))), ", zeroes: ", missing[missing$lab == i, "zeroes"])

    index = index + 1
}

pdf(gsub(".tsv",".pdf",args[1]), bg = "white", width = 7, height = 9)
ggplot(plot, aes(x=conc, y=count, color="typo")) + geom_point(size=5, shape=19, alpha=0.6) + scale_color_manual(values = cbPalette) +
       theme_bw() + xlim(-4,5) + ylim(-4,5) + 
       xlab("log10(Spike-In concentration)") + ylab ("log10(#reads + 0.0001)") + 
       theme(axis.text.x = element_text(size = 16, colour = "black"), axis.text.y = element_text(size = 16,  colour = "black"),
       plot.title = element_text(size = 22), axis.title.x = element_text(size = 18), axis.title.y = element_text(size = 18), 
       legend.position = "none", strip.text.x=element_text(size=22), strip.text.y=element_text(size=22)) + ggtitle(paste0(unique(plot$category), " ", unique(plot$sample))) + 
       facet_wrap(~ lab, ncol=1, scales="free") + geom_smooth(method=lm,se=FALSE,fullrange=TRUE, size=1.5) + theme(legend.position = "none") +
       geom_text(data = graphLabels, aes(x = -3, y = 5, label = label))
dev.off()