library(ggplot2) library(plyr) library(magrittr) library(dplyr) args = commandArgs(trailingOnly=TRUE) lm_eqn = function(plot) { m = lm(count ~ conc, plot); eq = paste0("r2: ", format(summary(m)$r.squared, digits = 3)) 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") missing = plot[ plot$count == 0, ] %>% group_by(lab, sample) %>% summarize( zeroes = n()) dim(plot) 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,] dim(finalE) # text to annotate graphs graphLabels = data.frame( matrix(nrow = length(unique(plot$lab)), ncol = 3) ) colnames(graphLabels) = c("lab", "sample", "label") index = 1 for (i in unique(plot$lab)) { for (s in unique(plot$sample)) { graphLabels[index, "lab"] = i graphLabels[index, "sample"] = s definitive = finalE[finalE$lab == i & finalE$sample == s, ] if ( nrow(definitive) != 0 ) { # graphLabels[index, "label"] = paste0(as.character(as.expression(lm_eqn(definitive))), ", zeroes: ", missing[missing$lab == i & missing$sample == s, "zeroes"]) #r2 computed without zeroes graphLabels[index, "label"] = paste0(as.character(as.expression(lm_eqn(plot[plot$lab == i & plot$sample == s,]))), ", zeroes: ", missing[missing$lab == i & missing$sample == s, "zeroes"]) index = index + 1 } } } pdf(gsub(".tsv", "_zeroes.pdf", args[1]), bg = "white", width = 70, height = 45) #pdf(gsub(".tsv", ".pdf", args[1]), bg = "white", width = 70, height = 45) ggplot(plot, aes(x=conc, y=count, color="type")) + 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 = 30, colour = "black"), axis.text.y = element_text(size = 30, colour = "black"), plot.title = element_text(size = 30), axis.title.x = element_text(size = 18), axis.title.y = element_text(size = 18), legend.position = "none", strip.text.x=element_text(size=40), strip.text.y=element_text(size=40)) + ggtitle(unique(plot$category)) + facet_grid( sample ~ lab, scales = "free") + # geom_smooth(data=finalE,method=lm,se=FALSE,fullrange=TRUE, size=1.5) + geom_smooth(method=lm,se=FALSE,fullrange=TRUE, size=1.5) + theme(legend.position = "none") + geom_text(data = graphLabels, aes(x = 1, y = 5, label = label), size = 12) dev.off()