annotate sequence_overview.r @ 81:b6f9a640e098 draft

Uploaded

library(reshape2)

args <- commandArgs(trailingOnly = TRUE)

before.unique.file = args[1]
merged.file = args[2]
outputdir = args[3]
gene.classes = unlist(strsplit(args[4], ","))
hotspot.analysis.sum.file = args[5]
NToverview.file = paste(outputdir, "ntoverview.txt", sep="/")
NTsum.file = paste(outputdir, "ntsum.txt", sep="/")
main.html = "index.html"
empty.region.filter = args[6]


setwd(outputdir)

before.unique = read.table(before.unique.file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
merged = read.table(merged.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="")

#before.unique = before.unique[!grepl("unmatched", before.unique$best_match),]

if(empty.region.filter == "leader"){
	before.unique$seq_conc = paste(before.unique$FR1.IMGT.seq, before.unique$CDR1.IMGT.seq, before.unique$FR2.IMGT.seq, before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq, before.unique$CDR3.IMGT.seq)
} else if(empty.region.filter == "FR1"){
	before.unique$seq_conc = paste(before.unique$CDR1.IMGT.seq, before.unique$FR2.IMGT.seq, before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq, before.unique$CDR3.IMGT.seq)
} else if(empty.region.filter == "CDR1"){
	before.unique$seq_conc = paste(before.unique$FR2.IMGT.seq, before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq, before.unique$CDR3.IMGT.seq)
} else if(empty.region.filter == "FR2"){
	before.unique$seq_conc = paste(before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq, before.unique$CDR3.IMGT.seq)
}

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

dat = data.frame(table(before.unique$seq_conc))

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...
get.bg.color = function(val){
	if(val %in% c("TRUE", "FALSE", "T", "F")){ #if its a logical value, give the background a green/red color
		return(ifelse(val,"#eafaf1","#f9ebea"))
	} else if (!is.na(as.numeric(val))) { #if its a numerical value, give it a grey tint if its >0
		return(ifelse(val > 0,"#eaecee","white"))
	} else {
		return("white")
	}
}
td = function(val) {
  return(paste("<td bgcolor='", get.bg.color(val), "'>", val, "</td>", sep=""))
}
tr = function(val) { 
	return(paste(c("<tr>", sapply(val, td), "</tr>"), collapse="")) 
}

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("<center><img src='data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAA8AAAAPCAYAAAA71pVKAAAAzElEQVQoka2TwQ2CQBBFpwTshw4ImW8ogJMlUIMmhNCDxgasAi50oSXA8XlAjCG7aqKTzGX/vsnM31mzR0gk7tTudO5MEizpzvQ4ryUSe408J3Xn+grE0p1rnpOamVmWsZG4rS+dzzAMsN8Hi9yyjI1JNGtxu4VxBJgLRLpoTKIPiW0LlwtUVRTubW2OBGUJu92cZRmdfbKQMAw8o+vi5v0fLorZ7Y9waGYJjsf38DJz0O1PsEQffOcv4Sa6YYfDDJ5Obzbsp93+5VfdATueO1fdLdI0AAAAAElFTkSuQmCC'> Please note that this tab is based on all sequences before filter unique sequences and the remove duplicates based on filters are applied. In this table only sequences occuring more than once are included. </center>", file=main.html, append=F)
cat("<table border='1' class='pure-table pure-table-striped'>", file=main.html, append=T)

if(empty.region.filter == "leader"){
	cat("<caption>FR1+CDR1+FR2+CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T)
} else if(empty.region.filter == "FR1"){
	cat("<caption>CDR1+FR2+CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T)
} else if(empty.region.filter == "CDR1"){
	cat("<caption>FR2+CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T)
} else if(empty.region.filter == "FR2"){
	cat("<caption>CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T)
}

cat("<tr>", file=main.html, append=T)
cat("<th>Sequence</th><th>Functionality</th><th>IGA1</th><th>IGA2</th><th>IGG1</th><th>IGG2</th><th>IGG3</th><th>IGG4</th><th>IGM</th><th>IGE</th><th>UN</th>", file=main.html, append=T)
cat("<th>total IGA</th><th>total IGG</th><th>total IGM</th><th>total IGE</th><th>number of subclasses</th><th>present in both IGA and IGG</th><th>present in IGA, IGG and IGM</th><th>present in IGA, IGG and IGE</th><th>present in IGA, IGG, IGM and IGE</th><th>IGA1+IGA2</th>", file=main.html, append=T)
cat("<th>IGG1+IGG2</th><th>IGG1+IGG3</th><th>IGG1+IGG4</th><th>IGG2+IGG3</th><th>IGG2+IGG4</th><th>IGG3+IGG4</th>", file=main.html, append=T)
cat("<th>IGG1+IGG2+IGG3</th><th>IGG2+IGG3+IGG4</th><th>IGG1+IGG2+IGG4</th><th>IGG1+IGG3+IGG4</th><th>IGG1+IGG2+IGG3+IGG4</th>", file=main.html, append=T)
cat("</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
unmatched=0 #all of the sequences are unmatched
some.unmatched=0 #one or more sequences in a clone are unmatched
matched=0 #should be the same als matched sequences

sequence.id.page="by_id.html"

for(i in 1:nrow(dat)){
	
	ca1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGA1", IDs$best_match),]
	ca2 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGA2", IDs$best_match),]
	
	cg1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGG1", IDs$best_match),]
	cg2 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGG2", IDs$best_match),]
	cg3 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGG3", IDs$best_match),]
	cg4 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGG4", IDs$best_match),]
	
	cm = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGM", IDs$best_match),]
	
	ce = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGE", IDs$best_match),]
	
	un = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^unmatched", IDs$best_match),]
	
	allc = rbind(ca1, ca2, cg1, cg2, cg3, cg4, cm, ce, un)
	
	ca1.n = nrow(ca1)
	ca2.n = nrow(ca2)
	
	cg1.n = nrow(cg1)
	cg2.n = nrow(cg2)
	cg3.n = nrow(cg3)
	cg4.n = nrow(cg4)
	
	cm.n = nrow(cm)
	
	ce.n = nrow(ce)
	
	un.n = nrow(un)
	
	classes = c(ca1.n, ca2.n, cg1.n, cg2.n, cg3.n, cg4.n, cm.n, ce.n, un.n)
	
	classes.sum = sum(classes)
	
	if(classes.sum == 1){
		single.sequences = single.sequences + 1
		next
	}
	
	if(un.n == classes.sum){
		unmatched = unmatched + 1
		next
	}
	
	classes.no.un = classes[-length(classes)]
	
	in.classes = sum(classes.no.un > 0)
	
	matched = matched + in.classes #count in how many subclasses the sequence occurs.
	
	if(any(classes == classes.sum)){
		multiple.in.one = multiple.in.one + 1
	} else if (un.n > 0) {
		some.unmatched = some.unmatched + 1
	} else {
		in.multiple = in.multiple + 1
	}
	
	id = as.numeric(dat[i,"seq_conc"])
	
	functionality = paste(unique(allc[,"Functionality"]), collapse=",")
	
	by.id.row = c()
	
	if(ca1.n > 0){
		cat(tbl(ca1), file=paste("IGA1_", id, ".html", sep=""))
	}

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

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

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

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

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

	if(cm.n > 0){
		cat(tbl(cm), file=paste("IGM_", id, ".html", sep=""))
	}

	if(ce.n > 0){
		cat(tbl(ce), file=paste("IGE_", id, ".html", sep=""))
	}

	if(un.n > 0){
		cat(tbl(un), file=paste("un_", id, ".html", sep=""))
	}
	
	ca1.html = make.link(id, "IGA1", ca1.n)
	ca2.html = make.link(id, "IGA2", ca2.n)
	
	cg1.html = make.link(id, "IGG1", cg1.n)
	cg2.html = make.link(id, "IGG2", cg2.n)
	cg3.html = make.link(id, "IGG3", cg3.n)
	cg4.html = make.link(id, "IGG4", cg4.n)
	
	cm.html = make.link(id, "IGM", cm.n)
	
	ce.html = make.link(id, "IGE", ce.n)
	
	un.html = make.link(id, "un", un.n)
	
	#extra columns
	ca.n = ca1.n + ca2.n
	
	cg.n = cg1.n + cg2.n + cg3.n + cg4.n
	
	#in.classes
	
	in.ca.cg = (ca.n > 0 & cg.n > 0)
	
	in.ca.cg.cm = (ca.n > 0 & cg.n > 0 & cm.n > 0)
	
	in.ca.cg.ce = (ca.n > 0 & cg.n > 0 & ce.n > 0)
	
	in.ca.cg.cm.ce = (ca.n > 0 & cg.n > 0 & cm.n > 0 & ce.n > 0)
	
	in.ca1.ca2 = (ca1.n > 0 & ca2.n > 0)
	
	in.cg1.cg2 = (cg1.n > 0 & cg2.n > 0)
	in.cg1.cg3 = (cg1.n > 0 & cg3.n > 0)
	in.cg1.cg4 = (cg1.n > 0 & cg4.n > 0)
	in.cg2.cg3 = (cg2.n > 0 & cg3.n > 0)
	in.cg2.cg4 = (cg2.n > 0 & cg4.n > 0)
	in.cg3.cg4 = (cg3.n > 0 & cg4.n > 0)
	
	in.cg1.cg2.cg3 = (cg1.n > 0 & cg2.n > 0 & cg3.n > 0)
	in.cg2.cg3.cg4 = (cg2.n > 0 & cg3.n > 0 & cg4.n > 0)
	in.cg1.cg2.cg4 = (cg1.n > 0 & cg2.n > 0 & cg4.n > 0)
	in.cg1.cg3.cg4 = (cg1.n > 0 & cg3.n > 0 & cg4.n > 0)
	
	in.cg.all = (cg1.n > 0 & cg2.n > 0 & cg3.n > 0 & cg4.n > 0)
	
	#rw = c(as.character(dat[i,"seq_conc"]), functionality, ca1.html, ca2.html, cg1.html, cg2.html, cg3.html, cg4.html, cm.html, un.html)
	rw = c(as.character(dat[i,"seq_conc"]), functionality, ca1.html, ca2.html, cg1.html, cg2.html, cg3.html, cg4.html, cm.html, ce.html, un.html)
	rw = c(rw, ca.n, cg.n, cm.n, ce.n, in.classes, in.ca.cg, in.ca.cg.cm, in.ca.cg.ce, in.ca.cg.cm.ce, in.ca1.ca2, in.cg1.cg2, in.cg1.cg3, in.cg1.cg4, in.cg2.cg3, in.cg2.cg4, in.cg3.cg4, in.cg1.cg2.cg3, in.cg2.cg3.cg4, in.cg1.cg2.cg4, in.cg1.cg3.cg4, in.cg.all)
	
	

	cat(tr(rw), file=main.html, append=T)
	
	
	for(i in 1:nrow(allc)){ #generate html by id
		html = make.link(id, allc[i,"best_match"], allc[i,"Sequence.ID"])
		cat(paste(html, "<br />"), file=sequence.id.page, 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("Matched with unmatched:", some.unmatched))
print(paste("Count that should match 'matched' sequences:", matched))

#ACGT overview

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

if(empty.region.filter == "leader"){
	NToverview$seq = paste(NToverview$FR1.IMGT.seq, NToverview$CDR1.IMGT.seq, NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq)
} else if(empty.region.filter == "FR1"){
	NToverview$seq = paste(NToverview$CDR1.IMGT.seq, NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq)
} else if(empty.region.filter == "CDR1"){
	NToverview$seq = paste(NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq)
} else if(empty.region.filter == "FR2"){
	NToverview$seq = paste(NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq)
}

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){
	print(paste("class:", clazz))
	NToverview.sub = NToverview[grepl(paste("^", clazz, sep=""), NToverview$best_match),]
	print(paste("nrow:", nrow(NToverview.sub)))
	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=""))
}

NToverview.tmp = NToverview[,c("Sequence.ID", "best_match", "seq", "A", "C", "G", "T")]

names(NToverview.tmp) = c("Sequence.ID", "best_match", "Sequence of the analysed region", "A", "C", "G", "T")

write.table(NToverview.tmp, 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")