library(data.table)
library(ggplot2)

args <- commandArgs(trailingOnly = TRUE)

input = args[1]
genes = unlist(strsplit(args[2], ","))
outputdir = args[3]
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",
                    "CDR3.IMGT.Nb.of.nucleotides",
                    "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]),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)

calculate_result = function(i, gene, dat, matrx, f, fname, name){
	tmp = dat[grepl(paste(".*", gene, ".*", sep=""), dat$best_match),]

	j = i - 1
	x = (j * 3) + 1
	y = (j * 3) + 2
	z = (j * 3) + 3
	 
	if(nrow(tmp) > 0){
	  
	  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 {
		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)
	  }
  }
  
  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)])
				}
			}
		}
  }
  
  
  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(length(tmp$Sequence.ID), file=paste(name, "_", fname, "_n.txt" ,sep=""))
  
  matrx
}

nts = c("a", "c", "g", "t")
zeros=rep(0, 4)

funcs = c(median, sum, mean)
fnames = c("median", "sum", "mean")

for(i in 1:length(funcs)){
	func = funcs[[i]]
	fname = fnames[[i]]
	
	rows = 9
	if(fname == "sum"){
		rows = 11
	}
	matrx = matrix(data = 0, ncol=((length(genes) + 1) * 3),nrow=rows)

	for(i in 1:length(genes)){
	  matrx = calculate_result(i, genes[i], dat, matrx, func, fname, genes[i])
	}

	matrx = calculate_result(i + 1, ".*", dat, matrx, func, fname, name="all")
	
	result = data.frame(matrx)
	if(fname == "sum"){
		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)", "nt in FR", "nt in CDR")
	} else {
		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)
}


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