diff aa_histogram.r @ 57:16c7fc1c4bf8 draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aa_histogram.r	Fri Mar 18 07:50:34 2016 -0400
@@ -0,0 +1,47 @@
+library(ggplot2)
+
+args <- commandArgs(trailingOnly = TRUE)
+
+input = args[1]
+outfile = args[2]
+
+print("---------------- read input ----------------")
+
+dat = read.table(input, sep="\t", fill=T, header=T, quote="")
+
+print("---------------- as numeric ----------------")
+
+mutations.at.position = as.numeric(dat[1,])
+aa.at.position = as.numeric(dat[2,])
+
+print("---------------- freq data.frame ----------------")
+
+dat_freq = mutations.at.position / aa.at.position
+dat_dt = data.frame(i=1:length(dat_freq), freq=dat_freq)
+
+print("---------------- weird stuff ----------------")
+
+options(width=220)
+
+print(dat[,20:40])
+
+print(dat_dt) #need this or it will fail????
+
+print("---------------- plot ----------------")
+
+m = ggplot(dat_dt, aes(x=i, y=freq)) + theme(axis.text.x = element_text(angle = 90, hjust = 1))
+m = m + geom_histogram(stat="identity", colour = "black", fill = "darkgrey", alpha=0.8) + scale_x_continuous(breaks=1:length(dat_freq), labels=1:length(dat_freq))
+m = m + annotate("segment", x = 0.5, y = -0.05, xend=26.5, yend=-0.05, colour="darkgreen", size=1) + annotate("text", x = 13, y = -0.1, label="FR1")
+m = m + annotate("segment", x = 26.5, y = -0.07, xend=38.5, yend=-0.07, colour="darkblue", size=1) + annotate("text", x = 32.5, y = -0.15, label="CDR1")
+m = m + annotate("segment", x = 38.5, y = -0.05, xend=55.5, yend=-0.05, colour="darkgreen", size=1) + annotate("text", x = 47, y = -0.1, label="FR2")
+m = m + annotate("segment", x = 55.5, y = -0.07, xend=65.5, yend=-0.07, colour="darkblue", size=1) + annotate("text", x = 60.5, y = -0.15, label="CDR2")
+m = m + annotate("segment", x = 65.5, y = -0.05, xend=104.5, yend=-0.05, colour="darkgreen", size=1) + annotate("text", x = 85, y = -0.1, label="FR3")
+m = m + expand_limits(y=c(-0.1,1)) + xlab("AA position") + ylab("Frequency") + ggtitle("AA mutation frequency")
+
+print("---------------- write/print ----------------")
+
+write.table(dat_dt, paste(dirname(outfile), "/aa_histogram.txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T)
+
+png(filename=outfile, width=1280, height=720)
+print(m)
+dev.off()