Mercurial > repos > davidvanzessen > mutation_analysis

comparison sequence_overview.r.old @ 90:f0e8dac22c6e draft

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