Mercurial > repos > davidvanzessen > mutation_analysis

comparison aa_histogram.r @ 49:5c6b9e99d576 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Wed, 18 Nov 2015 05:55:04 -0500
parents 57d197f149c3
children d4e72eeea640
comparison
equal deleted inserted replaced
48:d09b1bdfd388 49:5c6b9e99d576
3 args <- commandArgs(trailingOnly = TRUE) 3 args <- commandArgs(trailingOnly = TRUE)
4 4
5 input = args[1] 5 input = args[1]
6 outfile = args[2] 6 outfile = args[2]
7 7
8 dat = read.table(input, header=F, sep=",") 8 dat = read.table(input, sep="\t", fill=T, header=T, quote="")
9 dat=as.numeric(dat[1,])
10 dat_sum = sum(dat)
11 9
12 dat_freq = dat / dat_sum * 100 10
11
12 mutations.at.position = as.numeric(dat[1,])
13 aa.at.position = as.numeric(dat[2,])
14
15 dat_freq = mutations.at.position / aa.at.position
13 dat_dt = data.frame(i=1:length(dat_freq), freq=dat_freq) 16 dat_dt = data.frame(i=1:length(dat_freq), freq=dat_freq)
17
18 options(width=220)
19
20 print(dat[,20:40])
14 21
15 m = ggplot(dat_dt, aes(x=i, y=freq)) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) 22 m = ggplot(dat_dt, aes(x=i, y=freq)) + theme(axis.text.x = element_text(angle = 90, hjust = 1))
16 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)) 23 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))
17 m = m + annotate("segment", x = 0.5, y = -0.3, xend=26.5, yend=-0.3, colour="darkgreen", size=1) + annotate("text", x = 13, y = -0.2, label="FR1") 24 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")
18 m = m + annotate("segment", x = 26.5, y = -0.4, xend=38.5, yend=-0.4, colour="darkblue", size=1) + annotate("text", x = 32.5, y = -0.3, label="CDR1") 25 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")
19 m = m + annotate("segment", x = 38.5, y = -0.3, xend=55.5, yend=-0.3, colour="darkgreen", size=1) + annotate("text", x = 47, y = -0.2, label="FR2") 26 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")
20 m = m + annotate("segment", x = 55.5, y = -0.4, xend=65.5, yend=-0.4, colour="darkblue", size=1) + annotate("text", x = 60.5, y = -0.3, label="CDR2") 27 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")
21 m = m + annotate("segment", x = 65.5, y = -0.3, xend=104.5, yend=-0.3, colour="darkgreen", size=1) + annotate("text", x = 85, y = -0.2, label="FR3") 28 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")
29 m = m + expand_limits(y=c(-0.1,1)) + xlab("AA position") + ylab("Frequency") + ggtitle("AA mutation frequency")
22 write.table(dat_dt, paste(dirname(outfile), "/aa_histogram.txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T) 30 write.table(dat_dt, paste(dirname(outfile), "/aa_histogram.txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T)
23 31
24 png(filename=outfile, width=1280, height=720) 32 png(filename=outfile, width=1280, height=720)
25 print(m) 33 print(m)
26 dev.off() 34 dev.off()