Mercurial > repos > davidvanzessen > mutation_analysis
diff mutation_analysis.r.bak @ 53:7290a88ea202 draft
Uploaded
| author | davidvanzessen |
|---|---|
| date | Mon, 29 Feb 2016 10:49:39 -0500 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mutation_analysis.r.bak Mon Feb 29 10:49:39 2016 -0500 @@ -0,0 +1,469 @@ +library(data.table) +library(ggplot2) + +args <- commandArgs(trailingOnly = TRUE) + +input = args[1] +genes = unlist(strsplit(args[2], ",")) +outputdir = args[3] +print(args[4]) +include_fr1 = ifelse(args[4] == "yes", T, F) +setwd(outputdir) + +dat = read.table(input, header=T, sep="\t", fill=T, stringsAsFactors=F) + +if(length(dat$Sequence.ID) == 0){ + setwd(outputdir) + result = data.frame(x = rep(0, 5), y = rep(0, 5), z = rep(NA, 5)) + row.names(result) = c("Number of Mutations (%)", "Transition (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C G (%)") + write.table(x=result, file="mutations.txt", sep=",",quote=F,row.names=T,col.names=F) + transitionTable = data.frame(A=rep(0, 4),C=rep(0, 4),G=rep(0, 4),T=rep(0, 4)) + row.names(transitionTable) = c("A", "C", "G", "T") + transitionTable["A","A"] = NA + transitionTable["C","C"] = NA + transitionTable["G","G"] = NA + transitionTable["T","T"] = NA + write.table(x=transitionTable, file="transitions.txt", sep=",",quote=F,row.names=T,col.names=NA) + cat("0", file="n.txt") + stop("No data") +} + + + +cleanup_columns = c("FR1.IMGT.c.a", + "FR2.IMGT.g.t", + "CDR1.IMGT.Nb.of.nucleotides", + "CDR2.IMGT.t.a", + "FR1.IMGT.c.g", + "CDR1.IMGT.c.t", + "FR2.IMGT.a.c", + "FR2.IMGT.Nb.of.mutations", + "FR2.IMGT.g.c", + "FR2.IMGT.a.g", + "FR3.IMGT.t.a", + "FR3.IMGT.t.c", + "FR2.IMGT.g.a", + "FR3.IMGT.c.g", + "FR1.IMGT.Nb.of.mutations", + "CDR1.IMGT.g.a", + "CDR1.IMGT.t.g", + "CDR1.IMGT.g.c", + "CDR2.IMGT.Nb.of.nucleotides", + "FR2.IMGT.a.t", + "CDR1.IMGT.Nb.of.mutations", + "CDR1.IMGT.a.g", + "FR3.IMGT.a.c", + "FR1.IMGT.g.a", + "FR3.IMGT.a.g", + "FR1.IMGT.a.t", + "CDR2.IMGT.a.g", + "CDR2.IMGT.Nb.of.mutations", + "CDR2.IMGT.g.t", + "CDR2.IMGT.a.c", + "CDR1.IMGT.t.c", + "FR3.IMGT.g.c", + "FR1.IMGT.g.t", + "FR3.IMGT.g.t", + "CDR1.IMGT.a.t", + "FR1.IMGT.a.g", + "FR3.IMGT.a.t", + "FR3.IMGT.Nb.of.nucleotides", + "FR2.IMGT.t.c", + "CDR2.IMGT.g.a", + "FR2.IMGT.t.a", + "CDR1.IMGT.t.a", + "FR2.IMGT.t.g", + "FR3.IMGT.t.g", + "FR2.IMGT.Nb.of.nucleotides", + "FR1.IMGT.t.a", + "FR1.IMGT.t.g", + "FR3.IMGT.c.t", + "FR1.IMGT.t.c", + "CDR2.IMGT.a.t", + "FR2.IMGT.c.t", + "CDR1.IMGT.g.t", + "CDR2.IMGT.t.g", + "FR1.IMGT.Nb.of.nucleotides", + "CDR1.IMGT.c.g", + "CDR2.IMGT.t.c", + "FR3.IMGT.g.a", + "CDR1.IMGT.a.c", + "FR2.IMGT.c.a", + "FR3.IMGT.Nb.of.mutations", + "FR2.IMGT.c.g", + "CDR2.IMGT.g.c", + "FR1.IMGT.g.c", + "CDR2.IMGT.c.t", + "FR3.IMGT.c.a", + "CDR1.IMGT.c.a", + "CDR2.IMGT.c.g", + "CDR2.IMGT.c.a", + "FR1.IMGT.c.t", + "FR1.IMGT.Nb.of.silent.mutations", + "FR2.IMGT.Nb.of.silent.mutations", + "FR3.IMGT.Nb.of.silent.mutations", + "FR1.IMGT.Nb.of.nonsilent.mutations", + "FR2.IMGT.Nb.of.nonsilent.mutations", + "FR3.IMGT.Nb.of.nonsilent.mutations") + +for(col in cleanup_columns){ + dat[,col] = gsub("\\(.*\\)", "", dat[,col]) + #dat[dat[,col] == "",] = "0" + dat[,col] = as.numeric(dat[,col]) + dat[is.na(dat[,col]),] = 0 +} + +regions = c("FR1", "CDR1", "FR2", "CDR2", "FR3") +if(!include_fr1){ + regions = c("CDR1", "FR2", "CDR2", "FR3") +} + +sum_by_row = function(x, columns) { sum(as.numeric(x[columns]), na.rm=T) } + +VRegionMutations_columns = paste(regions, ".IMGT.Nb.of.mutations", sep="") +dat$VRegionMutations = apply(dat, FUN=sum_by_row, 1, columns=VRegionMutations_columns) + +VRegionNucleotides_columns = paste(regions, ".IMGT.Nb.of.nucleotides", sep="") +dat$VRegionNucleotides = apply(dat, FUN=sum_by_row, 1, columns=VRegionNucleotides_columns) + +transitionMutations_columns = paste(rep(regions, each=4), c(".IMGT.a.g", ".IMGT.g.a", ".IMGT.c.t", ".IMGT.t.c"), sep="") +dat$transitionMutations = apply(dat, FUN=sum_by_row, 1, columns=transitionMutations_columns) + +transversionMutations_columns = paste(rep(regions, each=8), c(".IMGT.a.c",".IMGT.c.a",".IMGT.a.t",".IMGT.t.a",".IMGT.g.c",".IMGT.c.g",".IMGT.g.t",".IMGT.t.g"), sep="") +dat$transversionMutations = apply(dat, FUN=sum_by_row, 1, columns=transversionMutations_columns) + + +transitionMutationsAtGC_columns = paste(rep(regions, each=2), c(".IMGT.g.a",".IMGT.c.t"), sep="") +dat$transitionMutationsAtGC = apply(dat, FUN=sum_by_row, 1, columns=transitionMutationsAtGC_columns) + + +totalMutationsAtGC_columns = paste(rep(regions, each=6), c(".IMGT.c.g",".IMGT.c.t",".IMGT.c.a",".IMGT.g.c",".IMGT.g.a",".IMGT.g.t"), sep="") +#totalMutationsAtGC_columns = paste(rep(regions, each=6), c(".IMGT.g.a",".IMGT.c.t",".IMGT.c.a",".IMGT.c.g",".IMGT.g.t"), sep="") +dat$totalMutationsAtGC = apply(dat, FUN=sum_by_row, 1, columns=totalMutationsAtGC_columns) + +transitionMutationsAtAT_columns = paste(rep(regions, each=2), c(".IMGT.a.g",".IMGT.t.c"), sep="") +dat$transitionMutationsAtAT = apply(dat, FUN=sum_by_row, 1, columns=transitionMutationsAtAT_columns) + +totalMutationsAtAT_columns = paste(rep(regions, each=6), c(".IMGT.a.g",".IMGT.a.c",".IMGT.a.t",".IMGT.t.g",".IMGT.t.c",".IMGT.t.a"), sep="") +#totalMutationsAtAT_columns = paste(rep(regions, each=5), c(".IMGT.a.g",".IMGT.t.c",".IMGT.a.c",".IMGT.g.c",".IMGT.t.g"), sep="") +dat$totalMutationsAtAT = apply(dat, FUN=sum_by_row, 1, columns=totalMutationsAtAT_columns) + + +FRRegions = regions[grepl("FR", regions)] +CDRRegions = regions[grepl("CDR", regions)] + +FR_silentMutations_columns = paste(FRRegions, ".IMGT.Nb.of.silent.mutations", sep="") +dat$silentMutationsFR = apply(dat, FUN=sum_by_row, 1, columns=FR_silentMutations_columns) + +CDR_silentMutations_columns = paste(CDRRegions, ".IMGT.Nb.of.silent.mutations", sep="") +dat$silentMutationsCDR = apply(dat, FUN=sum_by_row, 1, columns=CDR_silentMutations_columns) + +FR_nonSilentMutations_columns = paste(FRRegions, ".IMGT.Nb.of.nonsilent.mutations", sep="") +dat$nonSilentMutationsFR = apply(dat, FUN=sum_by_row, 1, columns=FR_nonSilentMutations_columns) + +CDR_nonSilentMutations_columns = paste(CDRRegions, ".IMGT.Nb.of.nonsilent.mutations", sep="") +dat$nonSilentMutationsCDR = apply(dat, FUN=sum_by_row, 1, columns=CDR_nonSilentMutations_columns) + +mutation.sum.columns = c("Sequence.ID", "VRegionMutations", "VRegionNucleotides", "transitionMutations", "transversionMutations", "transitionMutationsAtGC", "transitionMutationsAtAT", "silentMutationsFR", "nonSilentMutationsFR", "silentMutationsCDR", "nonSilentMutationsCDR") + +write.table(dat[,mutation.sum.columns], "mutation_by_id.txt", sep="\t",quote=F,row.names=F,col.names=T) + +setwd(outputdir) + +nts = c("a", "c", "g", "t") +zeros=rep(0, 4) +matrx = matrix(data = 0, ncol=((length(genes) + 1) * 3),nrow=9) +for(i in 1:length(genes)){ + gene = genes[i] + tmp = dat[grepl(paste(".*", gene, ".*", sep=""), dat$best_match),] + if(gene == "."){ + tmp = dat + } + j = i - 1 + x = (j * 3) + 1 + y = (j * 3) + 2 + z = (j * 3) + 3 + matrx[1,x] = sum(tmp$VRegionMutations) + matrx[1,y] = sum(tmp$VRegionNucleotides) + matrx[1,z] = round(matrx[1,x] / matrx[1,y] * 100, digits=1) + + matrx[2,x] = sum(tmp$transitionMutations) + matrx[2,y] = sum(tmp$VRegionMutations) + matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1) + + matrx[3,x] = sum(tmp$transversionMutations) + matrx[3,y] = sum(tmp$VRegionMutations) + matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1) + + matrx[4,x] = sum(tmp$transitionMutationsAtGC) + matrx[4,y] = sum(tmp$totalMutationsAtGC) + matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1) + + matrx[5,x] = sum(tmp$totalMutationsAtGC) + matrx[5,y] = sum(tmp$VRegionMutations) + matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1) + + matrx[6,x] = sum(tmp$transitionMutationsAtAT) + matrx[6,y] = sum(tmp$totalMutationsAtAT) + matrx[6,z] = round(matrx[6,x] / matrx[6,y] * 100, digits=1) + + matrx[7,x] = sum(tmp$totalMutationsAtAT) + matrx[7,y] = sum(tmp$VRegionMutations) + matrx[7,z] = round(matrx[7,x] / matrx[7,y] * 100, digits=1) + + matrx[8,x] = sum(tmp$nonSilentMutationsFR) + matrx[8,y] = sum(tmp$silentMutationsFR) + matrx[8,z] = round(matrx[8,x] / matrx[8,y], digits=1) + + matrx[9,x] = sum(tmp$nonSilentMutationsCDR) + matrx[9,y] = sum(tmp$silentMutationsCDR) + matrx[9,z] = round(matrx[9,x] / matrx[9,y], digits=1) + + + transitionTable = data.frame(A=zeros,C=zeros,G=zeros,T=zeros) + row.names(transitionTable) = c("A", "C", "G", "T") + transitionTable["A","A"] = NA + transitionTable["C","C"] = NA + transitionTable["G","G"] = NA + transitionTable["T","T"] = NA + + if(nrow(tmp) > 0){ + for(nt1 in nts){ + for(nt2 in nts){ + if(nt1 == nt2){ + next + } + NT1 = LETTERS[letters == nt1] + NT2 = LETTERS[letters == nt2] + FR1 = paste("FR1.IMGT.", nt1, ".", nt2, sep="") + CDR1 = paste("CDR1.IMGT.", nt1, ".", nt2, sep="") + FR2 = paste("FR2.IMGT.", nt1, ".", nt2, sep="") + CDR2 = paste("CDR2.IMGT.", nt1, ".", nt2, sep="") + FR3 = paste("FR3.IMGT.", nt1, ".", nt2, sep="") + if(include_fr1){ + transitionTable[NT1,NT2] = sum(tmp[,c(FR1, CDR1, FR2, CDR2, FR3)]) + } else { + transitionTable[NT1,NT2] = sum(tmp[,c(CDR1, FR2, CDR2, FR3)]) + } + } + } + } + + + write.table(x=transitionTable, file=paste("transitions_", gene ,".txt", sep=""), sep=",",quote=F,row.names=T,col.names=NA) + write.table(x=tmp[,c("Sequence.ID", "best_match", "chunk_hit_percentage", "nt_hit_percentage", "start_locations")], file=paste("matched_", gene ,".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T) + + cat(matrx[1,x], file=paste(gene, "_value.txt" ,sep="")) + cat(length(tmp$Sequence.ID), file=paste(gene, "_n.txt" ,sep="")) +} + +#again for all of the data +tmp = dat +j = i +x = (j * 3) + 1 +y = (j * 3) + 2 +z = (j * 3) + 3 +matrx[1,x] = sum(tmp$VRegionMutations) +matrx[1,y] = sum(tmp$VRegionNucleotides) +matrx[1,z] = round(matrx[1,x] / matrx[1,y] * 100, digits=1) + +matrx[2,x] = sum(tmp$transitionMutations) +matrx[2,y] = sum(tmp$VRegionMutations) +matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1) + +matrx[3,x] = sum(tmp$transversionMutations) +matrx[3,y] = sum(tmp$VRegionMutations) +matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1) + +matrx[4,x] = sum(tmp$transitionMutationsAtGC) +matrx[4,y] = sum(tmp$totalMutationsAtGC) +matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1) + +matrx[5,x] = sum(tmp$totalMutationsAtGC) +matrx[5,y] = sum(tmp$VRegionMutations) +matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1) + +matrx[6,x] = sum(tmp$transitionMutationsAtAT) +matrx[6,y] = sum(tmp$totalMutationsAtAT) +matrx[6,z] = round(matrx[6,x] / matrx[6,y] * 100, digits=1) + +matrx[7,x] = sum(tmp$totalMutationsAtAT) +matrx[7,y] = sum(tmp$VRegionMutations) +matrx[7,z] = round(matrx[7,x] / matrx[7,y] * 100, digits=1) + +matrx[8,x] = sum(tmp$nonSilentMutationsFR) +matrx[8,y] = sum(tmp$silentMutationsFR) +matrx[8,z] = round(matrx[8,x] / matrx[8,y], digits=1) + +matrx[9,x] = sum(tmp$nonSilentMutationsCDR) +matrx[9,y] = sum(tmp$silentMutationsCDR) +matrx[9,z] = round(matrx[9,x] / matrx[9,y], digits=1) + +transitionTable = data.frame(A=1:4,C=1:4,G=1:4,T=1:4) +row.names(transitionTable) = c("A", "C", "G", "T") +transitionTable["A","A"] = NA +transitionTable["C","C"] = NA +transitionTable["G","G"] = NA +transitionTable["T","T"] = NA + + +for(nt1 in nts){ + for(nt2 in nts){ + if(nt1 == nt2){ + next + } + NT1 = LETTERS[letters == nt1] + NT2 = LETTERS[letters == nt2] + FR1 = paste("FR1.IMGT.", nt1, ".", nt2, sep="") + CDR1 = paste("CDR1.IMGT.", nt1, ".", nt2, sep="") + FR2 = paste("FR2.IMGT.", nt1, ".", nt2, sep="") + CDR2 = paste("CDR2.IMGT.", nt1, ".", nt2, sep="") + FR3 = paste("FR3.IMGT.", nt1, ".", nt2, sep="") + if(include_fr1){ + transitionTable[NT1,NT2] = sum(tmp[,c(FR1, CDR1, FR2, CDR2, FR3)]) + } else { + transitionTable[NT1,NT2] = sum(tmp[,c(CDR1, FR2, CDR2, FR3)]) + } + } +} +write.table(x=transitionTable, file="transitions.txt", sep=",",quote=F,row.names=T,col.names=NA) +write.table(x=tmp[,c("Sequence.ID", "best_match", "chunk_hit_percentage", "nt_hit_percentage", "start_locations")], file="matched_all.txt", sep="\t",quote=F,row.names=F,col.names=T) +cat(matrx[1,x], file="total_value.txt") +cat(length(tmp$Sequence.ID), file="total_n.txt") + + + +result = data.frame(matrx) +row.names(result) = c("Number of Mutations (%)", "Transition (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C G (%)", "Transitions at A T (%)", "Targeting of A T (%)", "FR R/S (ratio)", "CDR R/S (ratio)") + +write.table(x=result, file="mutations.txt", sep=",",quote=F,row.names=T,col.names=F) + + +if (!("ggplot2" %in% rownames(installed.packages()))) { + install.packages("ggplot2", repos="http://cran.xl-mirror.nl/") +} + + +genesForPlot = gsub("[0-9]", "", dat$best_match) +genesForPlot = data.frame(table(genesForPlot)) +colnames(genesForPlot) = c("Gene","Freq") +genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq) +write.table(genesForPlot, "all.txt", sep="\t",quote=F,row.names=F,col.names=T) + + +pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=label)) +pc = pc + geom_bar(width = 1, stat = "identity") +pc = pc + coord_polar(theta="y") +pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("Classes", "( n =", sum(genesForPlot$Freq), ")")) + +png(filename="all.png") +pc +dev.off() + + +#blegh +genesForPlot = dat[grepl("ca", dat$best_match),]$best_match +if(length(genesForPlot) > 0){ + genesForPlot = data.frame(table(genesForPlot)) + colnames(genesForPlot) = c("Gene","Freq") + genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq) + + pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=label)) + pc = pc + geom_bar(width = 1, stat = "identity") + pc = pc + coord_polar(theta="y") + pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IgA subclasses", "( n =", sum(genesForPlot$Freq), ")")) + write.table(genesForPlot, "ca.txt", sep="\t",quote=F,row.names=F,col.names=T) + + png(filename="ca.png") + print(pc) + dev.off() +} + +genesForPlot = dat[grepl("cg", dat$best_match),]$best_match +if(length(genesForPlot) > 0){ + genesForPlot = data.frame(table(genesForPlot)) + colnames(genesForPlot) = c("Gene","Freq") + genesForPlot$label = paste(genesForPlot$Gene, "-", genesForPlot$Freq) + + pc = ggplot(genesForPlot, aes(x = factor(1), y=Freq, fill=label)) + pc = pc + geom_bar(width = 1, stat = "identity") + pc = pc + coord_polar(theta="y") + pc = pc + xlab(" ") + ylab(" ") + ggtitle(paste("IgG subclasses", "( n =", sum(genesForPlot$Freq), ")")) + write.table(genesForPlot, "cg.txt", sep="\t",quote=F,row.names=F,col.names=T) + + png(filename="cg.png") + print(pc) + dev.off() +} + +dat$percentage_mutations = round(dat$VRegionMutations / dat$VRegionNucleotides * 100, 2) + +p = ggplot(dat, aes(best_match, percentage_mutations)) +p = p + geom_point(aes(colour=best_match), position="jitter") + geom_boxplot(aes(middle=mean(percentage_mutations)), alpha=0.1, outlier.shape = NA) +p = p + xlab("Subclass") + ylab("Frequency") + ggtitle("Frequency scatter plot") + +png(filename="scatter.png") +print(p) +dev.off() + +write.table(dat[,c("Sequence.ID", "best_match", "VRegionMutations", "VRegionNucleotides", "percentage_mutations")], "scatter.txt", sep="\t",quote=F,row.names=F,col.names=T) + +write.table(dat, input, sep="\t",quote=F,row.names=F,col.names=T) + + + + + + +dat$best_match_class = substr(dat$best_match, 0, 2) +freq_labels = c("0", "0-2", "2-5", "5-10", "10-15", "15-20", "20") +dat$frequency_bins = cut(dat$percentage_mutations, breaks=c(-Inf, 0, 2,5,10,15,20, Inf), labels=freq_labels) + +frequency_bins_data = data.frame(data.table(dat)[, list(frequency_count=.N), by=c("best_match_class", "frequency_bins")]) + +p = ggplot(frequency_bins_data, aes(frequency_bins, frequency_count)) +p = p + geom_bar(aes(fill=best_match_class), stat="identity", position="dodge") +p = p + xlab("Frequency ranges") + ylab("Frequency") + ggtitle("Mutation Frequencies by class") + +png(filename="frequency_ranges.png") +print(p) +dev.off() + +frequency_bins_data_by_class = frequency_bins_data + +write.table(frequency_bins_data_by_class, "frequency_ranges_classes.txt", sep="\t",quote=F,row.names=F,col.names=T) + +frequency_bins_data = data.frame(data.table(dat)[, list(frequency_count=.N), by=c("best_match", "frequency_bins")]) + +write.table(frequency_bins_data, "frequency_ranges_subclasses.txt", sep="\t",quote=F,row.names=F,col.names=T) + + +#frequency_bins_data_by_class +#frequency_ranges_subclasses.txt + + + + + + + + + + + + + + + + + + + + + + + + + + +