Mercurial > repos > davidvanzessen > shm_csr
comparison aa_histogram.r @ 81:b6f9a640e098 draft
Uploaded
| author | davidvanzessen |
|---|---|
| date | Fri, 19 Feb 2021 15:10:54 +0000 |
| parents | |
| children | 729738462297 |
comparison
equal
deleted
inserted
replaced
| 80:a4617f1d1d89 | 81:b6f9a640e098 |
|---|---|
| 1 library(ggplot2) | |
| 2 | |
| 3 args <- commandArgs(trailingOnly = TRUE) | |
| 4 | |
| 5 mutations.by.id.file = args[1] | |
| 6 absent.aa.by.id.file = args[2] | |
| 7 genes = strsplit(args[3], ",")[[1]] | |
| 8 genes = c(genes, "") | |
| 9 outdir = args[4] | |
| 10 | |
| 11 | |
| 12 print("---------------- read input ----------------") | |
| 13 | |
| 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="") | |
| 16 | |
| 17 for(gene in genes){ | |
| 18 graph.title = paste(gene, "AA mutation frequency") | |
| 19 if(gene == ""){ | |
| 20 mutations.by.id.gene = mutations.by.id[!grepl("unmatched", mutations.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" | |
| 24 } else { | |
| 25 mutations.by.id.gene = mutations.by.id[grepl(paste("^", gene, sep=""), mutations.by.id$best_match),] | |
| 26 absent.aa.by.id.gene = absent.aa.by.id[grepl(paste("^", gene, sep=""), absent.aa.by.id$best_match),] | |
| 27 } | |
| 28 print(paste("nrow", gene, nrow(absent.aa.by.id.gene))) | |
| 29 if(nrow(mutations.by.id.gene) == 0){ | |
| 30 next | |
| 31 } | |
| 32 | |
| 33 mutations.at.position = colSums(mutations.by.id.gene[,-c(1,2)]) | |
| 34 aa.at.position = colSums(absent.aa.by.id.gene[,-c(1,2,3,4)]) | |
| 35 | |
| 36 dat_freq = mutations.at.position / aa.at.position | |
| 37 dat_freq[is.na(dat_freq)] = 0 | |
| 38 dat_dt = data.frame(i=1:length(dat_freq), freq=dat_freq) | |
| 39 | |
| 40 | |
| 41 print("---------------- plot ----------------") | |
| 42 | |
| 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")) | |
| 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) | |
| 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") | |
| 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") | |
| 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") | |
| 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") | |
| 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") | |
| 50 m = m + expand_limits(y=c(-0.1,1)) + xlab("AA position") + ylab("Frequency") + ggtitle(graph.title) | |
| 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()) | |
| 52 #m = m + scale_colour_manual(values=c("black")) | |
| 53 | |
| 54 print("---------------- write/print ----------------") | |
| 55 | |
| 56 | |
| 57 dat.sums = data.frame(index=1:length(mutations.at.position), mutations.at.position=mutations.at.position, aa.at.position=aa.at.position) | |
| 58 | |
| 59 write.table(dat.sums, paste(outdir, "/aa_histogram_sum_", gene, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T) | |
| 60 write.table(mutations.by.id.gene, paste(outdir, "/aa_histogram_count_", gene, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T) | |
| 61 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) | |
| 62 write.table(dat_dt, paste(outdir, "/aa_histogram_", gene, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T) | |
| 63 | |
| 64 png(filename=paste(outdir, "/aa_histogram_", gene, ".png", sep=""), width=1280, height=720) | |
| 65 print(m) | |
| 66 dev.off() | |
| 67 | |
| 68 ggsave(paste(outdir, "/aa_histogram_", gene, ".pdf", sep=""), m, width=14, height=7) | |
| 69 } |