--- a/mutation_analysis.r	Fri Jan 29 08:11:31 2016 -0500
+++ b/mutation_analysis.r	Mon Feb 29 10:49:39 2016 -0500
@@ -170,55 +170,53 @@
 
 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]
+
+calculate_result = function(i, gene, dat, matrx, f, fname, name){
   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)
-  
+     
+  if(nrow(tmp) > 0){
+	  
+	  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(matrx[1,x] / matrx[1,y] * 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)
+  }
   
   transitionTable = data.frame(A=zeros,C=zeros,G=zeros,T=zeros)
   row.names(transitionTable) = c("A", "C", "G", "T")
@@ -250,93 +248,40 @@
   }
   
   
-  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)
+  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(gene, "_value.txt" ,sep=""))
-  cat(length(tmp$Sequence.ID), file=paste(gene, "_n.txt" ,sep=""))
+  cat(matrx[1,x], file=paste(name, "_", fname, "_value.txt" ,sep=""))
+  cat(length(tmp$Sequence.ID), file=paste(name, "_", fname, "_n.txt" ,sep=""))
+  
+  matrx
 }
 
-#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)
+nts = c("a", "c", "g", "t")
+zeros=rep(0, 4)
 
-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)
+funcs = c(median, sum, mean)
+fnames = c("median", "sum", "mean")
 
-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)
+for(i in 1:length(funcs)){
+	func = funcs[[i]]
+	fname = fnames[[i]]
+	
+	matrx = matrix(data = 0, ncol=((length(genes) + 1) * 3),nrow=9)
 
-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(i in 1:length(genes)){
+	  matrx = calculate_result(i, genes[i], dat, matrx, func, fname, genes[i])
+	}
 
-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)])
-		}
-	}
+	matrx = calculate_result(i + 1, ".*", dat, matrx, func, fname, name="all")
+	
+	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=paste("mutations_", fname, ".txt", sep=""), sep=",",quote=F,row.names=T,col.names=F)
 }
-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()))) {