Mercurial > repos > davidvanzessen > shm_csr

comparison aa_histogram.r @ 39:a24f8c93583a draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Thu, 22 Dec 2016 09:39:27 -0500
parents 80c4eebf7bc9
children 1cf60ae234b4
comparison
equal deleted inserted replaced
38:05c62efdc393 39:a24f8c93583a
13 13
14 mutations.by.id = read.table(mutations.by.id.file, sep="\t", fill=T, header=T, quote="") 14 mutations.by.id = read.table(mutations.by.id.file, sep="\t", fill=T, header=T, quote="")
15 absent.aa.by.id = read.table(absent.aa.by.id.file, sep="\t", fill=T, header=T, quote="") 15 absent.aa.by.id = read.table(absent.aa.by.id.file, sep="\t", fill=T, header=T, quote="")
16 16
17 for(gene in genes){ 17 for(gene in genes){
18 graph.title = paste(gene, "AA mutation frequency")
18 if(gene == ""){ 19 if(gene == ""){
19 mutations.by.id.gene = mutations.by.id[!grepl("unmatched", mutations.by.id$best_match),] 20 mutations.by.id.gene = mutations.by.id[!grepl("unmatched", mutations.by.id$best_match),]
20 absent.aa.by.id.gene = absent.aa.by.id[!grepl("unmatched", absent.aa.by.id$best_match),] 21 absent.aa.by.id.gene = absent.aa.by.id[!grepl("unmatched", absent.aa.by.id$best_match),]
22
23 graph.title = "AA mutation frequency all"
21 } else { 24 } else {
22 mutations.by.id.gene = mutations.by.id[grepl(paste("^", gene, sep=""), mutations.by.id$best_match),] 25 mutations.by.id.gene = mutations.by.id[grepl(paste("^", gene, sep=""), mutations.by.id$best_match),]
23 absent.aa.by.id.gene = absent.aa.by.id[grepl(paste("^", gene, sep=""), absent.aa.by.id$best_match),] 26 absent.aa.by.id.gene = absent.aa.by.id[grepl(paste("^", gene, sep=""), absent.aa.by.id$best_match),]
24 } 27 }
25 print(paste("nrow", gene, nrow(absent.aa.by.id.gene))) 28 print(paste("nrow", gene, nrow(absent.aa.by.id.gene)))
31 aa.at.position = colSums(absent.aa.by.id.gene[,-c(1,2,3,4)]) 34 aa.at.position = colSums(absent.aa.by.id.gene[,-c(1,2,3,4)])
32 35
33 dat_freq = mutations.at.position / aa.at.position 36 dat_freq = mutations.at.position / aa.at.position
34 dat_freq[is.na(dat_freq)] = 0 37 dat_freq[is.na(dat_freq)] = 0
35 dat_dt = data.frame(i=1:length(dat_freq), freq=dat_freq) 38 dat_dt = data.frame(i=1:length(dat_freq), freq=dat_freq)
39
36 40
37 print("---------------- plot ----------------") 41 print("---------------- plot ----------------")
38 42
39 m = ggplot(dat_dt, aes(x=i, y=freq)) + theme(axis.text.x = element_text(angle = 90, hjust = 1), text = element_text(size=13, colour="black")) 43 m = ggplot(dat_dt, aes(x=i, y=freq)) + theme(axis.text.x = element_text(angle = 90, hjust = 1), text = element_text(size=13, colour="black"))
40 m = m + geom_bar(stat="identity", colour = "black", fill = "darkgrey", alpha=0.8) + scale_x_continuous(breaks=dat_dt$i, labels=dat_dt$i) 44 m = m + geom_bar(stat="identity", colour = "black", fill = "darkgrey", alpha=0.8) + scale_x_continuous(breaks=dat_dt$i, labels=dat_dt$i)
41 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") 45 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")
42 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") 46 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")
43 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") 47 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")
44 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") 48 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")
45 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") 49 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")
46 m = m + expand_limits(y=c(-0.1,1)) + xlab("AA position") + ylab("Frequency") + ggtitle(paste(gene, "AA mutation frequency")) 50 m = m + expand_limits(y=c(-0.1,1)) + xlab("AA position") + ylab("Frequency") + ggtitle(graph.title)
47 m = m + theme(panel.background = element_rect(fill = "white", colour="black"), panel.grid.major.y = element_line(colour = "black"), panel.grid.major.x = element_blank()) 51 m = m + theme(panel.background = element_rect(fill = "white", colour="black"), panel.grid.major.y = element_line(colour = "black"), panel.grid.major.x = element_blank())
48 #m = m + scale_colour_manual(values=c("black")) 52 #m = m + scale_colour_manual(values=c("black"))
49 53
50 print("---------------- write/print ----------------") 54 print("---------------- write/print ----------------")
51 55