args <- commandArgs(trailingOnly = TRUE)

input = args[1]
summaryinput = args[2]
outputdir = args[3]
setwd(outputdir)

#dat = read.table("NWK276_MID6_25NT/8_V-REGION-nt-mutation-statistics_NWK276_MID6_25NT_051113.txt", header=T, sep="\t", fill=T, stringsAsFactors=F)
dat = read.table(input, header=T, sep="\t", fill=T, stringsAsFactors=F)

datSum = read.table(summaryinput, header=T, sep="\t", fill=T, stringsAsFactors=F)
datSum = datSum[,c("Sequence.ID", "AA.JUNCTION")]

dat = merge(dat, datSum, by="Sequence.ID", all.x=T)

#dat = dat[dat$Functionality == "productive",]

dat$VGene = gsub("^Homsap ", "", dat$V.GENE.and.allele)
dat$VGene = gsub("[*].*", "", dat$VGene) 

dat$past = paste(dat$AA.JUNCTION, dat$VGene)

#dat = dat[!duplicated(dat$past), ]

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")
}

#print(paste("After filtering on 'productive' and deduplicating based on V and AA JUNCTION there are", length(dat$Sequence.ID), "sequences left"))

result = data.frame(x = 1:5, y = 1:5, z = 1:5)
row.names(result) = c("Number of Mutations (%)", "Transition (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C G (%)")

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")

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
}

VRegionMutations =  sum(dat$FR1.IMGT.Nb.of.mutations + 
                        dat$CDR1.IMGT.Nb.of.mutations + 
                        dat$FR2.IMGT.Nb.of.mutations + 
                        dat$CDR2.IMGT.Nb.of.mutations +
                        dat$FR3.IMGT.Nb.of.mutations)

VRegionNucleotides = sum( dat$FR1.IMGT.Nb.of.nucleotides +
                          dat$CDR1.IMGT.Nb.of.nucleotides +
                          dat$FR2.IMGT.Nb.of.nucleotides +
                          dat$CDR2.IMGT.Nb.of.nucleotides +
                          dat$FR3.IMGT.Nb.of.nucleotides)

result[1,"x"] = VRegionMutations
result[1,"y"] = VRegionNucleotides
result[1,"z"] = round(result[1,"x"] / result[1,"y"] * 100, digits=1)

transitionMutations = sum(dat$FR1.IMGT.a.g +
                          dat$FR1.IMGT.g.a +
                          dat$FR1.IMGT.c.t +
                          dat$FR1.IMGT.t.c +
                          dat$CDR1.IMGT.a.g +
                          dat$CDR1.IMGT.g.a +
                          dat$CDR1.IMGT.c.t +
                          dat$CDR1.IMGT.t.c +
                          dat$FR2.IMGT.a.g +
                          dat$FR2.IMGT.g.a +
                          dat$FR2.IMGT.c.t +
                          dat$FR2.IMGT.t.c +
                          dat$CDR2.IMGT.a.g +
                          dat$CDR2.IMGT.g.a +
                          dat$CDR2.IMGT.c.t +
                          dat$CDR2.IMGT.t.c +
                          dat$FR3.IMGT.a.g +
                          dat$FR3.IMGT.g.a +
                          dat$FR3.IMGT.c.t +
                          dat$FR3.IMGT.t.c)

result[2,"x"] = transitionMutations
result[2,"y"] = VRegionMutations
result[2,"z"] = round(result[2,"x"] / result[2,"y"] * 100, digits=1)

transversionMutations = sum(  dat$FR1.IMGT.a.c +
                              dat$FR1.IMGT.c.a +
                              dat$FR1.IMGT.a.t +
                              dat$FR1.IMGT.t.a +
                              dat$FR1.IMGT.g.c +
                              dat$FR1.IMGT.c.g +
                              dat$FR1.IMGT.g.t +
                              dat$FR1.IMGT.t.g +
                              dat$CDR1.IMGT.a.c +
                              dat$CDR1.IMGT.c.a +
                              dat$CDR1.IMGT.a.t +
                              dat$CDR1.IMGT.t.a +
                              dat$CDR1.IMGT.g.c +
                              dat$CDR1.IMGT.c.g +
                              dat$CDR1.IMGT.g.t +
                              dat$CDR1.IMGT.t.g +
                              dat$FR2.IMGT.a.c +
                              dat$FR2.IMGT.c.a +
                              dat$FR2.IMGT.a.t +
                              dat$FR2.IMGT.t.a +
                              dat$FR2.IMGT.g.c +
                              dat$FR2.IMGT.c.g +
                              dat$FR2.IMGT.g.t +
                              dat$FR2.IMGT.t.g +
                              dat$CDR2.IMGT.a.c +
                              dat$CDR2.IMGT.c.a +
                              dat$CDR2.IMGT.a.t +
                              dat$CDR2.IMGT.t.a +
                              dat$CDR2.IMGT.g.c +
                              dat$CDR2.IMGT.c.g +
                              dat$CDR2.IMGT.g.t +
                              dat$CDR2.IMGT.t.g +
                              dat$FR3.IMGT.a.c +
                              dat$FR3.IMGT.c.a +
                              dat$FR3.IMGT.a.t +
                              dat$FR3.IMGT.t.a +
                              dat$FR3.IMGT.g.c +
                              dat$FR3.IMGT.c.g +
                              dat$FR3.IMGT.g.t +
                              dat$FR3.IMGT.t.g)

result[3,"x"] = transversionMutations
result[3,"y"] = VRegionMutations
result[3,"z"] = round(result[3,"x"] / result[3,"y"] * 100, digits=1)

transitionMutationsAtGC = sum(dat$FR1.IMGT.g.a +
                              dat$FR1.IMGT.c.t +
                              dat$CDR1.IMGT.g.a +
                              dat$CDR1.IMGT.c.t +
                              dat$FR2.IMGT.g.a +
                              dat$FR2.IMGT.c.t +
                              dat$CDR2.IMGT.g.a +
                              dat$CDR2.IMGT.c.t +
                              dat$FR3.IMGT.g.a +
                              dat$FR3.IMGT.c.t)

totalMutationsAtGC = sum(dat$FR1.IMGT.g.a +
                         dat$FR1.IMGT.c.t +
                         dat$FR1.IMGT.c.a +
                         dat$FR1.IMGT.g.c +
                         dat$FR1.IMGT.c.g +
                         dat$FR1.IMGT.g.t +
                         dat$CDR1.IMGT.g.a +
                         dat$CDR1.IMGT.c.t +
                         dat$CDR1.IMGT.c.a +
                         dat$CDR1.IMGT.g.c +
                         dat$CDR1.IMGT.c.g +
                         dat$CDR1.IMGT.g.t +
                         dat$FR2.IMGT.g.a +
                         dat$FR2.IMGT.c.t +
                         dat$FR2.IMGT.c.a +
                         dat$FR2.IMGT.g.c +
                         dat$FR2.IMGT.c.g +
                         dat$FR2.IMGT.g.t +
                         dat$CDR2.IMGT.g.a +
                         dat$CDR2.IMGT.c.t +
                         dat$CDR2.IMGT.c.a +
                         dat$CDR2.IMGT.g.c +
                         dat$CDR2.IMGT.c.g +
                         dat$CDR2.IMGT.g.t +
                         dat$FR3.IMGT.g.a +
                         dat$FR3.IMGT.c.t +
                         dat$FR3.IMGT.c.a +
                         dat$FR3.IMGT.g.c +
                         dat$FR3.IMGT.c.g +
                         dat$FR3.IMGT.g.t)

result[4,"x"] = transitionMutationsAtGC
result[4,"y"] = totalMutationsAtGC
result[4,"z"] = round(result[4,"x"] / result[4,"y"] * 100, digits=1)

result[5,"x"] = totalMutationsAtGC
result[5,"y"] = VRegionMutations
result[5,"z"] = round(result[5,"x"] / result[5,"y"] * 100, digits=1)


#transitiontable

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
nts = c("a", "c", "g", "t")


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="")
    transitionTable[NT1,NT2] = sum( dat[,FR1] + 
                                    dat[,CDR1] +
                                    dat[,FR2] +
                                    dat[,CDR2] +
                                    dat[,FR3])
  }
}

setwd(outputdir)
write.table(x=result, file="mutations.txt", 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)
cat(length(dat$Sequence.ID), file="n.txt")