annotate sequence_overview.r @ 22:0bea8c187a90 draft

Uploaded

1 library(reshape2)
2 
3 args <- commandArgs(trailingOnly = TRUE)
4 
5 before.unique.file = args[1]
6 merged.file = args[2]
7 outputdir = args[3]
8 gene.classes = unlist(strsplit(args[4], ","))
9 hotspot.analysis.sum.file = args[5]
10 NToverview.file = paste(outputdir, "ntoverview.txt", sep="/")
11 NTsum.file = paste(outputdir, "ntsum.txt", sep="/")
12 main.html = "index.html"
13 empty.region.filter = args[6]
14 
15 
16 setwd(outputdir)
17 
18 before.unique = read.table(before.unique.file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
19 merged = read.table(merged.file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
20 hotspot.analysis.sum = read.table(hotspot.analysis.sum.file, header=F, sep=",", fill=T, stringsAsFactors=F, quote="")
21 
22 #before.unique = before.unique[!grepl("unmatched", before.unique$best_match),]
23 
24 if(empty.region.filter == "leader"){
25 	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)
26 } else if(empty.region.filter == "FR1"){
27 	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)
28 } else if(empty.region.filter == "CDR1"){
29 	before.unique$seq_conc = paste(before.unique$FR2.IMGT.seq, before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq)
30 } else if(empty.region.filter == "FR2"){
31 	before.unique$seq_conc = paste(before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq)
32 }
33 
34 IDs = before.unique[,c("Sequence.ID", "seq_conc", "best_match", "Functionality")]
35 IDs$best_match = as.character(IDs$best_match)
36 
37 dat = data.frame(table(before.unique$seq_conc))
38 
39 names(dat) = c("seq_conc", "Freq")
40 
41 dat$seq_conc = factor(dat$seq_conc)
42 
43 dat = dat[order(as.character(dat$seq_conc)),]
44 
45 #writing html from R...
46 get.bg.color = function(val){
47 	if(val %in% c("TRUE", "FALSE", "T", "F")){ #if its a logical value, give the background a green/red color
48 		return(ifelse(val,"#eafaf1","#f9ebea"))
49 	} else if (!is.na(as.numeric(val))) { #if its a numerical value, give it a grey tint if its >0
50 		return(ifelse(val > 0,"#eaecee","white"))
51 	} else {
52 		return("white")
53 	}
54 }
55 td = function(val) {
56   return(paste("<td bgcolor='", get.bg.color(val), "'>", val, "</td>", sep=""))
57 }
58 tr = function(val) { 
59 	return(paste(c("<tr>", sapply(val, td), "</tr>"), collapse="")) 
60 }
61 
62 make.link = function(id, clss, val) { 
63 	paste("<a href='", clss, "_", id, ".html'>", val, "</a>", sep="") 
64 }
65 tbl = function(df) {
66 	res = "<table border='1'>"
67 	for(i in 1:nrow(df)){ 
68 		res = paste(res, tr(df[i,]), sep="")
69 	}
70 	res = paste(res, "</table>")
71 }
72 
73 cat("<table border='1' class='pure-table pure-table-striped'>", file=main.html, append=F)
74 
75 if(empty.region.filter == "leader"){
76 	cat("<caption>FR1+CDR1+FR2+CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T)
77 } else if(empty.region.filter == "FR1"){
78 	cat("<caption>CDR1+FR2+CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T)
79 } else if(empty.region.filter == "CDR1"){
80 	cat("<caption>FR2+CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T)
81 } else if(empty.region.filter == "FR2"){
82 	cat("<caption>CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T)
83 }
84 
85 cat("<tr>", file=main.html, append=T)
86 cat("<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><th>un</th>", file=main.html, append=T)
87 cat("<th>total CA</th><th>total CG</th><th>number of subclasses</th><th>present in both Ca and Cg</th><th>Ca1+Ca2</th>", file=main.html, append=T)
88 cat("<th>Cg1+Cg2</th><th>Cg1+Cg3</th><th>Cg1+Cg4</th><th>Cg2+Cg3</th><th>Cg2+Cg4</th><th>Cg3+Cg4</th>", file=main.html, append=T)
89 cat("<th>Cg1+Cg2+Cg3</th><th>Cg2+Cg3+Cg4</th><th>Cg1+Cg2+Cg4</th><th>Cg1+Cg3+Cg4</th><th>Cg1+Cg2+Cg3+Cg4</th>", file=main.html, append=T)
90 cat("</tr>", file=main.html, append=T)
91 
92 
93 
94 single.sequences=0 #sequence only found once, skipped
95 in.multiple=0 #same sequence across multiple subclasses
96 multiple.in.one=0 #same sequence multiple times in one subclass
97 unmatched=0 #all of the sequences are unmatched
98 some.unmatched=0 #one or more sequences in a clone are unmatched
99 matched=0 #should be the same als matched sequences
00 
01 sequence.id.page="by_id.html"
02 
03 for(i in 1:nrow(dat)){
04 	
05 	ca1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGA1", IDs$best_match),]
06 	ca2 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGA2", IDs$best_match),]
07 	
08 	cg1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGG1", IDs$best_match),]
09 	cg2 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGG2", IDs$best_match),]
10 	cg3 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGG3", IDs$best_match),]
11 	cg4 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGG4", IDs$best_match),]
12 	
13 	cm = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^IGM", IDs$best_match),]
14 	
15 	un = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^unmatched", IDs$best_match),]
16 	allc = rbind(ca1, ca2, cg1, cg2, cg3, cg4, cm, un)
17 	
18 	ca1.n = nrow(ca1)
19 	ca2.n = nrow(ca2)
20 	
21 	cg1.n = nrow(cg1)
22 	cg2.n = nrow(cg2)
23 	cg3.n = nrow(cg3)
24 	cg4.n = nrow(cg4)
25 	
26 	cm.n = nrow(cm)
27 	
28 	un.n = nrow(un)
29 	
30 	classes = c(ca1.n, ca2.n, cg1.n, cg2.n, cg3.n, cg4.n, cm.n, un.n)
31 	
32 	classes.sum = sum(classes)
33 	
34 	if(classes.sum == 1){
35 		single.sequences = single.sequences + 1
36 		next
37 	}
38 	
39 	if(un.n == classes.sum){
40 		unmatched = unmatched + 1
41 		next
42 	}
43 	
44 	in.classes = sum(classes > 0)
45 	
46 	matched = matched + in.classes #count in how many subclasses the sequence occurs.
47 	
48 	if(any(classes  == classes.sum)){
49 		multiple.in.one = multiple.in.one + 1
50 	} else if (un.n > 0) {
51 		some.unmatched = some.unmatched + 1
52 	} else {
53 		in.multiple = in.multiple + 1
54 	}
55 	
56 	id = as.numeric(dat[i,"seq_conc"])
57 	
58 	functionality = paste(unique(allc[,"Functionality"]), collapse=",")
59 	
60 	by.id.row = c()
61 	
62 	if(ca1.n > 0){
63 		cat(tbl(ca1), file=paste("IGA1_", id, ".html", sep=""))
64 	}
65 
66 	if(ca2.n > 0){
67 		cat(tbl(ca2), file=paste("IGA2_", id, ".html", sep=""))
68 	}
69 
70 	if(cg1.n > 0){
71 		cat(tbl(cg1), file=paste("IGG1_", id, ".html", sep=""))
72 	}
73 
74 	if(cg2.n > 0){
75 		cat(tbl(cg2), file=paste("IGG2_", id, ".html", sep=""))
76 	}
77 
78 	if(cg3.n > 0){
79 		cat(tbl(cg3), file=paste("IGG3_", id, ".html", sep=""))
80 	}
81 
82 	if(cg4.n > 0){
83 		cat(tbl(cg4), file=paste("IGG4_", id, ".html", sep=""))
84 	}
85 
86 	if(cm.n > 0){
87 		cat(tbl(cm), file=paste("IGM_", id, ".html", sep=""))
88 	}
89 
90 	if(un.n > 0){
91 		cat(tbl(un), file=paste("un_", id, ".html", sep=""))
92 	}
93 	
94 	ca1.html = make.link(id, "IGA1", ca1.n)
95 	ca2.html = make.link(id, "IGA2", ca2.n)
96 	
97 	cg1.html = make.link(id, "IGG1", cg1.n)
98 	cg2.html = make.link(id, "IGG2", cg2.n)
99 	cg3.html = make.link(id, "IGG3", cg3.n)
00 	cg4.html = make.link(id, "IGG4", cg4.n)
01 	
02 	cm.html = make.link(id, "IGM", cm.n)
03 	
04 	un.html = make.link(id, "un", un.n)
05 	
06 	#extra columns
07 	ca.n = ca1.n + ca2.n
08 	
09 	cg.n = cg1.n + cg2.n + cg3.n + cg4.n
10 	
11 	#in.classes
12 	
13 	in.ca.cg = (ca.n > 0 & cg.n > 0)
14 	
15 	in.ca1.ca2 = (ca1.n > 0 & ca2.n > 0)
16 	
17 	in.cg1.cg2 = (cg1.n > 0 & cg2.n > 0)
18 	in.cg1.cg3 = (cg1.n > 0 & cg3.n > 0)
19 	in.cg1.cg4 = (cg1.n > 0 & cg4.n > 0)
20 	in.cg2.cg3 = (cg2.n > 0 & cg3.n > 0)
21 	in.cg2.cg4 = (cg2.n > 0 & cg4.n > 0)
22 	in.cg3.cg4 = (cg3.n > 0 & cg4.n > 0)
23 	
24 	in.cg1.cg2.cg3 = (cg1.n > 0 & cg2.n > 0 & cg3.n > 0)
25 	in.cg2.cg3.cg4 = (cg2.n > 0 & cg3.n > 0 & cg4.n > 0)
26 	in.cg1.cg2.cg4 = (cg1.n > 0 & cg2.n > 0 & cg4.n > 0)
27 	in.cg1.cg3.cg4 = (cg1.n > 0 & cg3.n > 0 & cg4.n > 0)
28 	
29 	in.cg.all = (cg1.n > 0 & cg2.n > 0 & cg3.n > 0 & cg4.n > 0)
30 	
31 	
32 	
33 	
34 	#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)
35 	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)
36 	rw = c(rw, ca.n, cg.n, in.classes, in.ca.cg, 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)
37 
38 	cat(tr(rw), file=main.html, append=T)
39 	
40 	
41 	for(i in 1:nrow(allc)){ #generate html by id
42 		html = make.link(id, allc[i,"best_match"], allc[i,"Sequence.ID"])
43 		cat(paste(html, "<br />"), file=sequence.id.page, append=T)
44 	}
45 }
46 
47 cat("</table>", file=main.html, append=T)
48 
49 print(paste("Single sequences:", single.sequences))
50 print(paste("Sequences in multiple subclasses:", in.multiple))
51 print(paste("Multiple sequences in one subclass:", multiple.in.one))
52 print(paste("Matched with unmatched:", some.unmatched))
53 print(paste("Count that should match 'matched' sequences:", matched))
54 
55 #ACGT overview
56 
57 #NToverview = merged[!grepl("^unmatched", merged$best_match),]
58 NToverview = merged
59 
60 if(empty.region.filter == "leader"){
61 	NToverview$seq = paste(NToverview$FR1.IMGT.seq, NToverview$CDR1.IMGT.seq, NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq)
62 } else if(empty.region.filter == "FR1"){
63 	NToverview$seq = paste(NToverview$CDR1.IMGT.seq, NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq)
64 } else if(empty.region.filter == "CDR1"){
65 	NToverview$seq = paste(NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq)
66 } else if(empty.region.filter == "FR2"){
67 	NToverview$seq = paste(NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq)
68 }
69 
70 NToverview$A = nchar(gsub("[^Aa]", "", NToverview$seq))
71 NToverview$C = nchar(gsub("[^Cc]", "", NToverview$seq))
72 NToverview$G = nchar(gsub("[^Gg]", "", NToverview$seq))
73 NToverview$T = nchar(gsub("[^Tt]", "", NToverview$seq))
74 
75 #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))
76 
77 #NToverview = rbind(NToverview, NTsum)
78 
79 NTresult = data.frame(nt=c("A", "C", "T", "G"))
80 
81 for(clazz in gene.classes){
82 	print(paste("class:", clazz))
83 	NToverview.sub = NToverview[grepl(paste("^", clazz, sep=""), NToverview$best_match),]
84 	print(paste("nrow:", nrow(NToverview.sub)))
85 	new.col.x = c(sum(NToverview.sub$A), sum(NToverview.sub$C), sum(NToverview.sub$T), sum(NToverview.sub$G))
86 	new.col.y = sum(new.col.x)
87 	new.col.z = round(new.col.x / new.col.y * 100, 2)
88 	
89 	tmp = names(NTresult)
90 	NTresult = cbind(NTresult, data.frame(new.col.x, new.col.y, new.col.z))
91 	names(NTresult) = c(tmp, paste(clazz, c("x", "y", "z"), sep=""))
92 }
93 
94 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)
95 
96 NToverview = NToverview[!grepl("unmatched", NToverview$best_match),]
97 
98 new.col.x = c(sum(NToverview$A), sum(NToverview$C), sum(NToverview$T), sum(NToverview$G))
99 new.col.y = sum(new.col.x)
00 new.col.z = round(new.col.x / new.col.y * 100, 2)
01 
02 tmp = names(NTresult)
03 NTresult = cbind(NTresult, data.frame(new.col.x, new.col.y, new.col.z))
04 names(NTresult) = c(tmp, paste("all", c("x", "y", "z"), sep=""))
05 
06 names(hotspot.analysis.sum) = names(NTresult)
07 
08 hotspot.analysis.sum = rbind(hotspot.analysis.sum, NTresult)
09 
10 write.table(hotspot.analysis.sum, hotspot.analysis.sum.file, quote=F, sep=",", row.names=F, col.names=F, na="0")
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
24 
25 
26 
27 
28 
29 
30 
31 
32 
33 
34 
35 
36 
37 
38 
39 
40