|
76
|
1 library(reshape2)
|
|
|
2
|
|
|
3 args <- commandArgs(trailingOnly = TRUE)
|
|
|
4
|
|
90
|
5 input.file = args[1]
|
|
|
6 outputdir = args[2]
|
|
|
7 gene.classes = unlist(strsplit(args[3], ","))
|
|
|
8 hotspot.analysis.sum.file = args[4]
|
|
81
|
9 NToverview.file = paste(outputdir, "ntoverview.txt", sep="/")
|
|
|
10 NTsum.file = paste(outputdir, "ntsum.txt", sep="/")
|
|
76
|
11 main.html = "index.html"
|
|
|
12
|
|
|
13 setwd(outputdir)
|
|
|
14
|
|
90
|
15 merged = read.table(input.file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
|
|
82
|
16 hotspot.analysis.sum = read.table(hotspot.analysis.sum.file, header=F, sep=",", fill=T, stringsAsFactors=F, quote="")
|
|
76
|
17
|
|
90
|
18 merged$seq_conc = paste(merged$CDR1.IMGT.seq, merged$FR2.IMGT.seq, merged$CDR2.IMGT.seq, merged$FR3.IMGT.seq, merged$CDR3.IMGT.seq)
|
|
82
|
19
|
|
90
|
20 IDs = merged[,c("Sequence.ID", "seq_conc", "best_match", "Functionality")]
|
|
76
|
21 IDs$best_match = as.character(IDs$best_match)
|
|
|
22
|
|
|
23 #dat = data.frame(data.table(dat)[, list(freq=.N), by=c("best_match", "seq_conc")])
|
|
|
24
|
|
90
|
25 dat = data.frame(table(merged$seq_conc, merged$Functionality))
|
|
76
|
26
|
|
87
|
27 #dat = dat[dat$Freq > 1,]
|
|
76
|
28
|
|
82
|
29 names(dat) = c("seq_conc", "Functionality", "Freq")
|
|
76
|
30
|
|
|
31 dat$seq_conc = factor(dat$seq_conc)
|
|
|
32
|
|
82
|
33 dat = dat[order(as.character(dat$seq_conc)),]
|
|
76
|
34
|
|
|
35 #writing html from R...
|
|
|
36 td = function(val) { paste("<td>", val, "</td>", sep="") }
|
|
|
37 tr = function(val) { capture.output(cat("<tr>", td(val), "</tr>", sep="")) }
|
|
|
38 make.link = function(id, clss, val) { paste("<a href='", clss, "_", id, ".html'>", val, "</a>", sep="") }
|
|
|
39 tbl = function(df) { res = "<table border='1'>"; for(i in 1:nrow(df)){ res = paste(res, tr(df[i,]), sep=""); }; res = paste(res, "</table>"); }
|
|
|
40
|
|
|
41 cat("<table border='1'>", file=main.html, append=F)
|
|
88
|
42 cat("<caption>CDR1+FR2+CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T)
|
|
79
|
43 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)
|
|
76
|
44
|
|
|
45 for(i in 1:nrow(dat)){
|
|
|
46 ca1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "ca1",]
|
|
|
47 ca2 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "ca2",]
|
|
|
48
|
|
|
49 cg1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cg1",]
|
|
|
50 cg2 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cg2",]
|
|
|
51 cg3 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cg3",]
|
|
|
52 cg4 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cg4",]
|
|
|
53
|
|
|
54 cm = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cm",]
|
|
|
55
|
|
82
|
56 classes = c(nrow(ca1), nrow(ca2), nrow(cg1), nrow(cg2), nrow(cg3), nrow(cg4), nrow(cm))
|
|
|
57
|
|
|
58 classes.sum = sum(classes)
|
|
|
59
|
|
87
|
60 if(classes.sum == 1){
|
|
82
|
61 next
|
|
|
62 }
|
|
|
63
|
|
76
|
64 id = as.numeric(dat[i,"seq_conc"])
|
|
79
|
65
|
|
|
66 functionality = dat[i,"Functionality"]
|
|
76
|
67
|
|
|
68 if(nrow(ca1) > 0){
|
|
|
69 cat(tbl(ca1), file=paste("ca1_", id, ".html", sep=""))
|
|
|
70 }
|
|
|
71
|
|
|
72 if(nrow(ca2) > 0){
|
|
|
73 cat(tbl(ca2), file=paste("ca2_", id, ".html", sep=""))
|
|
|
74 }
|
|
|
75
|
|
|
76 if(nrow(cg1) > 0){
|
|
|
77 cat(tbl(cg1), file=paste("cg1_", id, ".html", sep=""))
|
|
|
78 }
|
|
|
79
|
|
|
80 if(nrow(cg2) > 0){
|
|
|
81 cat(tbl(cg2), file=paste("cg2_", id, ".html", sep=""))
|
|
|
82 }
|
|
|
83
|
|
|
84 if(nrow(cg3) > 0){
|
|
|
85 cat(tbl(cg3), file=paste("cg3_", id, ".html", sep=""))
|
|
|
86 }
|
|
|
87
|
|
|
88 if(nrow(cg4) > 0){
|
|
|
89 cat(tbl(cg4), file=paste("cg4_", id, ".html", sep=""))
|
|
|
90 }
|
|
|
91
|
|
|
92 if(nrow(cm) > 0){
|
|
|
93 cat(tbl(cm), file=paste("cm_", id, ".html", sep=""))
|
|
|
94 }
|
|
|
95
|
|
|
96 ca1.html = make.link(id, "ca1", nrow(ca1))
|
|
|
97 ca2.html = make.link(id, "ca2", nrow(ca2))
|
|
|
98
|
|
|
99 cg1.html = make.link(id, "cg1", nrow(cg1))
|
|
|
100 cg2.html = make.link(id, "cg2", nrow(cg2))
|
|
|
101 cg3.html = make.link(id, "cg3", nrow(cg3))
|
|
|
102 cg4.html = make.link(id, "cg4", nrow(cg4))
|
|
|
103
|
|
|
104 cm.html = make.link(id, "cm", nrow(cm))
|
|
|
105
|
|
79
|
106 rw = c(as.character(dat[i,"seq_conc"]), as.character(functionality), ca1.html, ca2.html, cg1.html, cg2.html, cg3.html, cg4.html, cm.html)
|
|
76
|
107
|
|
|
108 cat(tr(rw), file=main.html, append=T)
|
|
|
109 }
|
|
|
110
|
|
|
111 cat("</table>", file=main.html, append=T)
|
|
81
|
112
|
|
|
113
|
|
|
114 #ACGT overview
|
|
|
115
|
|
|
116
|
|
|
117
|
|
90
|
118 NToverview = merged
|
|
|
119 NToverview$seq = paste(NToverview$CDR1.IMGT.seq, NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq, sep="_")
|
|
82
|
120
|
|
81
|
121 NToverview$A = nchar(gsub("[^Aa]", "", NToverview$seq))
|
|
|
122 NToverview$C = nchar(gsub("[^Cc]", "", NToverview$seq))
|
|
|
123 NToverview$G = nchar(gsub("[^Gg]", "", NToverview$seq))
|
|
|
124 NToverview$T = nchar(gsub("[^Tt]", "", NToverview$seq))
|
|
|
125
|
|
90
|
126 print(sum(colSums(NToverview[,c("A", "C", "T", "G")])))
|
|
|
127
|
|
82
|
128 #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))
|
|
|
129
|
|
|
130 #NToverview = rbind(NToverview, NTsum)
|
|
|
131
|
|
|
132 NTresult = data.frame(nt=c("A", "C", "T", "G"))
|
|
|
133
|
|
|
134 for(clazz in gene.classes){
|
|
89
|
135 NToverview.sub = NToverview[grepl(clazz, paste("^", NToverview$best_match, sep="")),]
|
|
82
|
136 new.col.x = c(sum(NToverview.sub$A), sum(NToverview.sub$C), sum(NToverview.sub$T), sum(NToverview.sub$G))
|
|
|
137 new.col.y = sum(new.col.x)
|
|
|
138 new.col.z = round(new.col.x / new.col.y * 100, 2)
|
|
|
139
|
|
|
140 tmp = names(NTresult)
|
|
|
141 NTresult = cbind(NTresult, data.frame(new.col.x, new.col.y, new.col.z))
|
|
|
142 names(NTresult) = c(tmp, paste(clazz, c("x", "y", "z"), sep=""))
|
|
|
143 }
|
|
81
|
144
|
|
82
|
145 new.col.x = c(sum(NToverview$A), sum(NToverview$C), sum(NToverview$T), sum(NToverview$G))
|
|
|
146 new.col.y = sum(new.col.x)
|
|
|
147 new.col.z = round(new.col.x / new.col.y * 100, 2)
|
|
|
148
|
|
|
149 tmp = names(NTresult)
|
|
|
150 NTresult = cbind(NTresult, data.frame(new.col.x, new.col.y, new.col.z))
|
|
|
151 names(NTresult) = c(tmp, paste("all", c("x", "y", "z"), sep=""))
|
|
81
|
152
|
|
82
|
153 names(hotspot.analysis.sum) = names(NTresult)
|
|
|
154
|
|
|
155 hotspot.analysis.sum = rbind(hotspot.analysis.sum, NTresult)
|
|
|
156
|
|
|
157 print(hotspot.analysis.sum)
|
|
|
158
|
|
|
159 write.table(hotspot.analysis.sum, hotspot.analysis.sum.file, quote=F, sep=",", row.names=F, col.names=F, na="0")
|
|
|
160
|
|
81
|
161
|
|
|
162 write.table(NToverview, NToverview.file, quote=F, sep="\t", row.names=F, col.names=T)
|
|
|
163
|
|
|
164
|
|
|
165
|
|
|
166
|
|
|
167
|
|
|
168
|
|
|
169
|
|
|
170
|
|
|
171
|
|
|
172
|
|
|
173
|
|
|
174
|
|
|
175
|
|
|
176
|
|
|
177
|
|
|
178
|
|
|
179
|
|
|
180
|
|
|
181
|
|
|
182
|
|
|
183
|
|
|
184
|
|
|
185
|
|
|
186
|
|
|
187
|
|
|
188
|
|
|
189
|
|
|
190
|
|
|
191
|
|
|
192
|
