library(reshape2)

args <- commandArgs(trailingOnly = TRUE)

input.file = args[1]
outputdir = args[2]
gene.classes = unlist(strsplit(args[3], ","))
hotspot.analysis.sum.file = args[4]
NToverview.file = paste(outputdir, "ntoverview.txt", sep="/")
NTsum.file = paste(outputdir, "ntsum.txt", sep="/")
main.html = "index.html"

setwd(outputdir)

merged = read.table(input.file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
hotspot.analysis.sum = read.table(hotspot.analysis.sum.file, header=F, sep=",", fill=T, stringsAsFactors=F, quote="")

merged$seq_conc = paste(merged$CDR1.IMGT.seq, merged$FR2.IMGT.seq, merged$CDR2.IMGT.seq, merged$FR3.IMGT.seq, merged$CDR3.IMGT.seq)

IDs = merged[,c("Sequence.ID", "seq_conc", "best_match", "Functionality")]
IDs$best_match = as.character(IDs$best_match)

#dat = data.frame(data.table(dat)[, list(freq=.N), by=c("best_match", "seq_conc")])

dat = data.frame(table(merged$seq_conc))
#dat = data.frame(table(merged$seq_conc, merged$Functionality))

#dat = dat[dat$Freq > 1,]

#names(dat) = c("seq_conc", "Functionality", "Freq")
names(dat) = c("seq_conc", "Freq")

dat$seq_conc = factor(dat$seq_conc)

dat = dat[order(as.character(dat$seq_conc)),]

#writing html from R...
td = function(val) { paste("<td>", val, "</td>", sep="") }
tr = function(val) { capture.output(cat("<tr>", td(val), "</tr>", sep="")) }
make.link = function(id, clss, val) { paste("<a href='", clss, "_", id, ".html'>", val, "</a>", sep="") }
tbl = function(df) { res = "<table border='1'>"; for(i in 1:nrow(df)){ res = paste(res, tr(df[i,]), sep=""); }; res = paste(res, "</table>"); }

cat("<table border='1'>", file=main.html, append=F)
cat("<caption>CDR1+FR2+CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T)
cat("<tr><th>Sequence</th><th>Functionality</th><th>ca1</th><th>ca2</th><th>cg1</th><th>cg2</th><th>cg3</th><th>cg4</th><th>cm</th></tr>", file=main.html, append=T)



single.sequences=0 #sequence only found once, skipped
in.multiple=0 #same sequence across multiple subclasses
multiple.in.one=0 #same sequence multiple times in one subclass
matched=0 #should be the same als matched sequences

for(i in 1:nrow(dat)){
	
	ca1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^ca1", IDs$best_match),]
	ca2 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^ca2", IDs$best_match),]
	
	cg1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^cg1", IDs$best_match),]
	cg2 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^cg2", IDs$best_match),]
	cg3 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^cg3", IDs$best_match),]
	cg4 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^cg4", IDs$best_match),]
	
	cm = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^cm", IDs$best_match),]
	allc = rbind(ca1, ca2, cg1, cg2, cg3, cg4, cm)
	
	classes = c(nrow(ca1), nrow(ca2), nrow(cg1), nrow(cg2), nrow(cg3), nrow(cg4), nrow(cm))
	
	classes.sum = sum(classes)
	
	if(classes.sum == 1){
		single.sequences = single.sequences + 1
		next
	}
	
	matched = matched + sum(classes > 0) #count in how many subclasses the sequence occurs.
	
	if(any(classes  == classes.sum)){
		in.multiple = in.multiple + 1
	} else {
		multiple.in.one = multiple.in.one + 1
	}
	
	id = as.numeric(dat[i,"seq_conc"])
	
	functionality = paste(unique(allc[,"Functionality"]), collapse=",")

	if(nrow(ca1) > 0){
		cat(tbl(ca1), file=paste("ca1_", id, ".html", sep=""))
	}

	if(nrow(ca2) > 0){
		cat(tbl(ca2), file=paste("ca2_", id, ".html", sep=""))
	}

	if(nrow(cg1) > 0){
		cat(tbl(cg1), file=paste("cg1_", id, ".html", sep=""))
	}

	if(nrow(cg2) > 0){
		cat(tbl(cg2), file=paste("cg2_", id, ".html", sep=""))
	}

	if(nrow(cg3) > 0){
		cat(tbl(cg3), file=paste("cg3_", id, ".html", sep=""))
	}

	if(nrow(cg4) > 0){
		cat(tbl(cg4), file=paste("cg4_", id, ".html", sep=""))
	}

	if(nrow(cm) > 0){
		cat(tbl(cm), file=paste("cm_", id, ".html", sep=""))
	}
	
	ca1.html = make.link(id, "ca1", nrow(ca1))
	ca2.html = make.link(id, "ca2", nrow(ca2))
	
	cg1.html = make.link(id, "cg1", nrow(cg1))
	cg2.html = make.link(id, "cg2", nrow(cg2))
	cg3.html = make.link(id, "cg3", nrow(cg3))
	cg4.html = make.link(id, "cg4", nrow(cg4))
	
	cm.html = make.link(id, "cm", nrow(cm))
	
	rw = c(as.character(dat[i,"seq_conc"]), functionality, ca1.html, ca2.html, cg1.html, cg2.html, cg3.html, cg4.html, cm.html)

	cat(tr(rw), file=main.html, append=T)
}

cat("</table>", file=main.html, append=T)

print(paste("Single sequences:", single.sequences))
print(paste("Sequences in multiple subclasses:", in.multiple))
print(paste("Multiple sequences in one subclass:", multiple.in.one))
print(paste("Count that should match 'matched' sequences:", matched))

#ACGT overview

NToverview = merged
NToverview$seq = paste(NToverview$CDR1.IMGT.seq, NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq, sep="_")

NToverview$A = nchar(gsub("[^Aa]", "", NToverview$seq))
NToverview$C = nchar(gsub("[^Cc]", "", NToverview$seq))
NToverview$G = nchar(gsub("[^Gg]", "", NToverview$seq))
NToverview$T = nchar(gsub("[^Tt]", "", NToverview$seq))

#Nsum = data.frame(Sequence.ID="-", best_match="Sum", seq="-", A = sum(NToverview$A), C = sum(NToverview$C), G = sum(NToverview$G), T = sum(NToverview$T))

#NToverview = rbind(NToverview, NTsum)

NTresult = data.frame(nt=c("A", "C", "T", "G"))

for(clazz in gene.classes){
	NToverview.sub = NToverview[grepl(paste("^", clazz, sep=""), NToverview$best_match),]
	new.col.x = c(sum(NToverview.sub$A), sum(NToverview.sub$C), sum(NToverview.sub$T), sum(NToverview.sub$G))
	new.col.y = sum(new.col.x)
	new.col.z = round(new.col.x / new.col.y * 100, 2)
	
	tmp = names(NTresult)
	NTresult = cbind(NTresult, data.frame(new.col.x, new.col.y, new.col.z))
	names(NTresult) = c(tmp, paste(clazz, c("x", "y", "z"), sep=""))
}

write.table(NToverview[,c("Sequence.ID", "best_match", "seq", "A", "C", "G", "T")], NToverview.file, quote=F, sep="\t", row.names=F, col.names=T)

NToverview = NToverview[!grepl("unmatched", NToverview$best_match),]

new.col.x = c(sum(NToverview$A), sum(NToverview$C), sum(NToverview$T), sum(NToverview$G))
new.col.y = sum(new.col.x)
new.col.z = round(new.col.x / new.col.y * 100, 2)

tmp = names(NTresult)
NTresult = cbind(NTresult, data.frame(new.col.x, new.col.y, new.col.z))
names(NTresult) = c(tmp, paste("all", c("x", "y", "z"), sep=""))

names(hotspot.analysis.sum) = names(NTresult)

hotspot.analysis.sum = rbind(hotspot.analysis.sum, NTresult)

write.table(hotspot.analysis.sum, hotspot.analysis.sum.file, quote=F, sep=",", row.names=F, col.names=F, na="0")