Mercurial > repos > davidvanzessen > mutation_analysis

comparison aa_histogram.r @ 108:6add3e66f4fa draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Wed, 13 Jul 2016 09:04:15 -0400
parents 01c9993865af
children 0096cd454380
comparison
equal deleted inserted replaced
107:01c9993865af 108:6add3e66f4fa
24 absent.aa.by.id.gene = absent.aa.by.id[grepl(paste("^", gene, sep=""), absent.aa.by.id$best_match),] 24 absent.aa.by.id.gene = absent.aa.by.id[grepl(paste("^", gene, sep=""), absent.aa.by.id$best_match),]
25 } 25 }
26 if(nrow(mutations.by.id.gene) == 0){ 26 if(nrow(mutations.by.id.gene) == 0){
27 next 27 next
28 } 28 }
29 29
30 print(paste("nrow", gene, nrow(absent.aa.by.id.gene)))
31
30 mutations.at.position = colSums(mutations.by.id.gene[,-c(1,2)]) 32 mutations.at.position = colSums(mutations.by.id.gene[,-c(1,2)])
31 aa.at.position = colSums(absent.aa.by.id.gene[,-c(1,2,3,4)]) 33 aa.at.position = colSums(absent.aa.by.id.gene[,-c(1,2,3,4)])
32 34
33 dat_freq = mutations.at.position / aa.at.position 35 dat_freq = mutations.at.position / aa.at.position
34 dat_dt = data.frame(i=1:length(dat_freq), freq=dat_freq) 36 dat_dt = data.frame(i=1:length(dat_freq), freq=dat_freq)
43 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 = 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")
44 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 + 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")
45 m = m + expand_limits(y=c(-0.1,1)) + xlab("AA position") + ylab("Frequency") + ggtitle(paste(gene, "AA mutation frequency")) 47 m = m + expand_limits(y=c(-0.1,1)) + xlab("AA position") + ylab("Frequency") + ggtitle(paste(gene, "AA mutation frequency"))
46 48
47 print("---------------- write/print ----------------") 49 print("---------------- write/print ----------------")
48 50
51 dat.sums = data.frame(index=1:length(mutations.at.position), mutations.at.position=mutations.at.position, aa.at.position=aa.at.position)
52
53 write.table(dat.sums, paste(outdir, "/aa_histogram_sum_", gene, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T)
54 write.table(mutations.by.id.gene, paste(outdir, "/aa_histogram_count_", gene, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T)
55 write.table(absent.aa.by.id.gene, paste(outdir, "/aa_histogram_absent_", gene, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T)
49 write.table(dat_dt, paste(outdir, "/aa_histogram_", gene, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T) 56 write.table(dat_dt, paste(outdir, "/aa_histogram_", gene, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T)
50 png(filename=paste(outdir, "/aa_histogram_", gene, ".png", sep=""), width=1280, height=720) 57 png(filename=paste(outdir, "/aa_histogram_", gene, ".png", sep=""), width=1280, height=720)
51 print(m) 58 print(m)
52 dev.off() 59 dev.off()
53 } 60 }