Mercurial > repos > davidvanzessen > mutation_analysis
diff mutation_analysis.r @ 114:e7b550d52eb7 draft
Uploaded
| author | davidvanzessen |
|---|---|
| date | Tue, 09 Aug 2016 07:20:41 -0400 |
| parents | ade5cf6fd2dc |
| children | ede6c4ee5196 |
line wrap: on
line diff
--- a/mutation_analysis.r Thu Aug 04 04:52:51 2016 -0400 +++ b/mutation_analysis.r Tue Aug 09 07:20:41 2016 -0400 @@ -169,6 +169,8 @@ setwd(outputdir) +base.order = data.frame(base=c("A", "T", "C", "G"), order=1:4) + calculate_result = function(i, gene, dat, matrx, f, fname, name){ tmp = dat[grepl(paste("^", gene, ".*", sep=""), dat$best_match),] @@ -179,67 +181,67 @@ if(nrow(tmp) > 0){ - if(fname == "sum"){ + if(fname == "sum"){ matrx[1,x] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) matrx[1,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1) matrx[1,z] = round(f(matrx[1,x] / matrx[1,y]) * 100, digits=1) - } else { + } else { matrx[1,x] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) matrx[1,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1) matrx[1,z] = round(f(tmp$VRegionMutations / tmp$VRegionNucleotides) * 100, digits=1) - } - - matrx[2,x] = round(f(tmp$transitionMutations, na.rm=T), digits=1) - matrx[2,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) - matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1) - - matrx[3,x] = round(f(tmp$transversionMutations, na.rm=T), digits=1) - matrx[3,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) - matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1) - - matrx[4,x] = round(f(tmp$transitionMutationsAtGC, na.rm=T), digits=1) - matrx[4,y] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1) - matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1) - - matrx[5,x] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1) - matrx[5,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) - matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1) - - matrx[6,x] = round(f(tmp$transitionMutationsAtAT, na.rm=T), digits=1) - matrx[6,y] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1) - matrx[6,z] = round(matrx[6,x] / matrx[6,y] * 100, digits=1) - - matrx[7,x] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1) - matrx[7,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) - matrx[7,z] = round(matrx[7,x] / matrx[7,y] * 100, digits=1) - - matrx[8,x] = round(f(tmp$nonSilentMutationsFR, na.rm=T), digits=1) - matrx[8,y] = round(f(tmp$silentMutationsFR, na.rm=T), digits=1) - matrx[8,z] = round(matrx[8,x] / matrx[8,y], digits=1) - - matrx[9,x] = round(f(tmp$nonSilentMutationsCDR, na.rm=T), digits=1) - matrx[9,y] = round(f(tmp$silentMutationsCDR, na.rm=T), digits=1) - matrx[9,z] = round(matrx[9,x] / matrx[9,y], digits=1) - - if(fname == "sum"){ - matrx[10,x] = round(f(rowSums(tmp[,c("FR2.IMGT.Nb.of.nucleotides", "FR3.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1) - matrx[10,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1) - matrx[10,z] = round(matrx[10,x] / matrx[10,y], digits=1) - - matrx[11,x] = round(f(rowSums(tmp[,c("CDR1.IMGT.Nb.of.nucleotides", "CDR2.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1) - matrx[11,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1) - matrx[11,z] = round(matrx[11,x] / matrx[11,y], digits=1) - } - } + } + + matrx[2,x] = round(f(tmp$transitionMutations, na.rm=T), digits=1) + matrx[2,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) + matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1) + + matrx[3,x] = round(f(tmp$transversionMutations, na.rm=T), digits=1) + matrx[3,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) + matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1) + + matrx[4,x] = round(f(tmp$transitionMutationsAtGC, na.rm=T), digits=1) + matrx[4,y] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1) + matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1) + + matrx[5,x] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1) + matrx[5,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) + matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1) + + matrx[6,x] = round(f(tmp$transitionMutationsAtAT, na.rm=T), digits=1) + matrx[6,y] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1) + matrx[6,z] = round(matrx[6,x] / matrx[6,y] * 100, digits=1) + + matrx[7,x] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1) + matrx[7,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) + matrx[7,z] = round(matrx[7,x] / matrx[7,y] * 100, digits=1) + + matrx[8,x] = round(f(tmp$nonSilentMutationsFR, na.rm=T), digits=1) + matrx[8,y] = round(f(tmp$silentMutationsFR, na.rm=T), digits=1) + matrx[8,z] = round(matrx[8,x] / matrx[8,y], digits=1) + + matrx[9,x] = round(f(tmp$nonSilentMutationsCDR, na.rm=T), digits=1) + matrx[9,y] = round(f(tmp$silentMutationsCDR, na.rm=T), digits=1) + matrx[9,z] = round(matrx[9,x] / matrx[9,y], digits=1) + + if(fname == "sum"){ + matrx[10,x] = round(f(rowSums(tmp[,c("FR2.IMGT.Nb.of.nucleotides", "FR3.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1) + matrx[10,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1) + matrx[10,z] = round(matrx[10,x] / matrx[10,y], digits=1) + + matrx[11,x] = round(f(rowSums(tmp[,c("CDR1.IMGT.Nb.of.nucleotides", "CDR2.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1) + matrx[11,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1) + matrx[11,z] = round(matrx[11,x] / matrx[11,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 + 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){ + if(nrow(tmp) > 0){ for(nt1 in nts){ for(nt2 in nts){ if(nt1 == nt2){ @@ -259,20 +261,40 @@ } } } - } - - - print(paste("writing value file: ", name, "_", fname, "_value.txt" ,sep="")) - - write.table(x=transitionTable, file=paste("transitions_", name ,"_", fname, ".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_", name , "_", fname, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T) - - cat(matrx[1,x], file=paste(name, "_", fname, "_value.txt" ,sep="")) - cat(nrow(tmp), file=paste(name, "_", fname, "_n.txt" ,sep="")) - - print(paste(fname, name, nrow(tmp))) - - matrx + transition = transitionTable + transition$id = names(transition) + + transition2 = melt(transition, id.vars="id") + + transition2 = merge(transition2, base.order, by.x="id", by.y="base") + transition2 = merge(transition2, base.order, by.x="variable", by.y="base") + + transition2[is.na(transition2$value),]$value = 0 + + png(filename=paste("transitions_stacked_", name, ".png", sep="")) + p = ggplot(transition2, aes(factor(reorder(id, order.x)), y=value, fill=factor(reorder(variable, order.y)))) + geom_bar(position="fill", stat="identity") #stacked bar + p = p + xlab("From base") + ylab("To base") + ggtitle("Mutations frequency from base to base") + guides(fill=guide_legend(title=NULL)) + print(p) + dev.off() + + png(filename=paste("transitions_heatmap_", name, ".png", sep="")) + p = ggplot(transition2, aes(factor(reorder(id, order.x)), factor(reorder(variable, order.y)))) + geom_tile(aes(fill = value), colour="white") + scale_fill_gradient(low="white", high="steelblue") #heatmap + p = p + xlab("From base") + ylab("To base") + ggtitle("Mutations frequency from base to base") + print(p) + dev.off() + } + + #print(paste("writing value file: ", name, "_", fname, "_value.txt" ,sep="")) + + write.table(x=transitionTable, file=paste("transitions_", name ,"_", fname, ".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_", name , "_", fname, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T) + + cat(matrx[1,x], file=paste(name, "_", fname, "_value.txt" ,sep="")) + cat(nrow(tmp), file=paste(name, "_", fname, "_n.txt" ,sep="")) + + #print(paste(fname, name, nrow(tmp))) + + matrx } nts = c("a", "c", "g", "t") @@ -322,12 +344,6 @@ write.table(x=new.table, file="mutations_sum.txt", sep=",",quote=F,row.names=F,col.names=F) - - -if (!("ggplot2" %in% rownames(installed.packages()))) { - install.packages("ggplot2", repos="http://cran.xl-mirror.nl/") -} - dat = dat[!grepl("^unmatched", dat$best_match),] #blegh