|
26
|
1 library(ggplot2)
|
|
|
2
|
|
|
3 args <- commandArgs(trailingOnly = TRUE)
|
|
|
4
|
|
|
5 input = args[1]
|
|
|
6 outfile = args[2]
|
|
104
|
7 gene = args[3]
|
|
26
|
8
|
|
56
|
9 print("---------------- read input ----------------")
|
|
|
10
|
|
49
|
11 dat = read.table(input, sep="\t", fill=T, header=T, quote="")
|
|
|
12
|
|
56
|
13 print("---------------- as numeric ----------------")
|
|
49
|
14
|
|
|
15 mutations.at.position = as.numeric(dat[1,])
|
|
|
16 aa.at.position = as.numeric(dat[2,])
|
|
26
|
17
|
|
56
|
18 print("---------------- freq data.frame ----------------")
|
|
|
19
|
|
49
|
20 dat_freq = mutations.at.position / aa.at.position
|
|
26
|
21 dat_dt = data.frame(i=1:length(dat_freq), freq=dat_freq)
|
|
|
22
|
|
56
|
23 print("---------------- weird stuff ----------------")
|
|
|
24
|
|
49
|
25 options(width=220)
|
|
|
26
|
|
|
27 print(dat[,20:40])
|
|
|
28
|
|
51
|
29 print(dat_dt) #need this or it will fail????
|
|
|
30
|
|
56
|
31 print("---------------- plot ----------------")
|
|
|
32
|
|
67
|
33 #also need these two, I don't even know...
|
|
|
34 print(paste("dat_freq", length(dat_freq)))
|
|
|
35 print(paste("dat_dt", nrow(dat_dt)))
|
|
|
36
|
|
26
|
37 m = ggplot(dat_dt, aes(x=i, y=freq)) + theme(axis.text.x = element_text(angle = 90, hjust = 1))
|
|
82
|
38 m = m + geom_bar(stat="identity", colour = "black", fill = "darkgrey", alpha=0.8) + scale_x_continuous(breaks=1:length(dat_freq), labels=1:length(dat_freq))
|
|
49
|
39 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")
|
|
|
40 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")
|
|
|
41 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")
|
|
|
42 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")
|
|
|
43 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")
|
|
104
|
44 m = m + expand_limits(y=c(-0.1,1)) + xlab("AA position") + ylab("Frequency") + ggtitle(paste(gene, "AA mutation frequency"))
|
|
57
|
45
|
|
|
46 print("---------------- write/print ----------------")
|
|
|
47
|
|
102
|
48 print("writing dat_dt") #need this
|
|
104
|
49 write.table(dat_dt, paste(dirname(outfile), "/aa_histogram_", gene, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T)
|
|
102
|
50 print("writing png") #also need this, file is haunted
|
|
26
|
51 png(filename=outfile, width=1280, height=720)
|
|
|
52 print(m)
|
|
|
53 dev.off()
|
