Mercurial > repos > davidvanzessen > mutation_analysis

comparison mutation_analysis.r @ 114:e7b550d52eb7 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Tue, 09 Aug 2016 07:20:41 -0400
parents ade5cf6fd2dc
children ede6c4ee5196
comparison
equal deleted inserted replaced
113:b84477f57318 114:e7b550d52eb7
167 167
168 write.table(dat[,mutation.sum.columns], "mutation_by_id.txt", sep="\t",quote=F,row.names=F,col.names=T) 168 write.table(dat[,mutation.sum.columns], "mutation_by_id.txt", sep="\t",quote=F,row.names=F,col.names=T)
169 169
170 setwd(outputdir) 170 setwd(outputdir)
171 171
172 base.order = data.frame(base=c("A", "T", "C", "G"), order=1:4)
173
172 calculate_result = function(i, gene, dat, matrx, f, fname, name){ 174 calculate_result = function(i, gene, dat, matrx, f, fname, name){
173 tmp = dat[grepl(paste("^", gene, ".*", sep=""), dat$best_match),] 175 tmp = dat[grepl(paste("^", gene, ".*", sep=""), dat$best_match),]
174 176
175 j = i - 1 177 j = i - 1
176 x = (j * 3) + 1 178 x = (j * 3) + 1
177 y = (j * 3) + 2 179 y = (j * 3) + 2
178 z = (j * 3) + 3 180 z = (j * 3) + 3
179 181
180 if(nrow(tmp) > 0){ 182 if(nrow(tmp) > 0){
181 183
182 if(fname == "sum"){ 184 if(fname == "sum"){
183 matrx[1,x] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) 185 matrx[1,x] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
184 matrx[1,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1) 186 matrx[1,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
185 matrx[1,z] = round(f(matrx[1,x] / matrx[1,y]) * 100, digits=1) 187 matrx[1,z] = round(f(matrx[1,x] / matrx[1,y]) * 100, digits=1)
186 } else { 188 } else {
187 matrx[1,x] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) 189 matrx[1,x] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
188 matrx[1,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1) 190 matrx[1,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
189 matrx[1,z] = round(f(tmp$VRegionMutations / tmp$VRegionNucleotides) * 100, digits=1) 191 matrx[1,z] = round(f(tmp$VRegionMutations / tmp$VRegionNucleotides) * 100, digits=1)
190 } 192 }
191 193
192 matrx[2,x] = round(f(tmp$transitionMutations, na.rm=T), digits=1) 194 matrx[2,x] = round(f(tmp$transitionMutations, na.rm=T), digits=1)
193 matrx[2,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) 195 matrx[2,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
194 matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1) 196 matrx[2,z] = round(matrx[2,x] / matrx[2,y] * 100, digits=1)
195 197
196 matrx[3,x] = round(f(tmp$transversionMutations, na.rm=T), digits=1) 198 matrx[3,x] = round(f(tmp$transversionMutations, na.rm=T), digits=1)
197 matrx[3,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) 199 matrx[3,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
198 matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1) 200 matrx[3,z] = round(matrx[3,x] / matrx[3,y] * 100, digits=1)
199 201
200 matrx[4,x] = round(f(tmp$transitionMutationsAtGC, na.rm=T), digits=1) 202 matrx[4,x] = round(f(tmp$transitionMutationsAtGC, na.rm=T), digits=1)
201 matrx[4,y] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1) 203 matrx[4,y] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1)
202 matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1) 204 matrx[4,z] = round(matrx[4,x] / matrx[4,y] * 100, digits=1)
203 205
204 matrx[5,x] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1) 206 matrx[5,x] = round(f(tmp$totalMutationsAtGC, na.rm=T), digits=1)
205 matrx[5,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) 207 matrx[5,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
206 matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1) 208 matrx[5,z] = round(matrx[5,x] / matrx[5,y] * 100, digits=1)
207 209
208 matrx[6,x] = round(f(tmp$transitionMutationsAtAT, na.rm=T), digits=1) 210 matrx[6,x] = round(f(tmp$transitionMutationsAtAT, na.rm=T), digits=1)
209 matrx[6,y] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1) 211 matrx[6,y] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1)
210 matrx[6,z] = round(matrx[6,x] / matrx[6,y] * 100, digits=1) 212 matrx[6,z] = round(matrx[6,x] / matrx[6,y] * 100, digits=1)
211 213
212 matrx[7,x] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1) 214 matrx[7,x] = round(f(tmp$totalMutationsAtAT, na.rm=T), digits=1)
213 matrx[7,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1) 215 matrx[7,y] = round(f(tmp$VRegionMutations, na.rm=T), digits=1)
214 matrx[7,z] = round(matrx[7,x] / matrx[7,y] * 100, digits=1) 216 matrx[7,z] = round(matrx[7,x] / matrx[7,y] * 100, digits=1)
215 217
216 matrx[8,x] = round(f(tmp$nonSilentMutationsFR, na.rm=T), digits=1) 218 matrx[8,x] = round(f(tmp$nonSilentMutationsFR, na.rm=T), digits=1)
217 matrx[8,y] = round(f(tmp$silentMutationsFR, na.rm=T), digits=1) 219 matrx[8,y] = round(f(tmp$silentMutationsFR, na.rm=T), digits=1)
218 matrx[8,z] = round(matrx[8,x] / matrx[8,y], digits=1) 220 matrx[8,z] = round(matrx[8,x] / matrx[8,y], digits=1)
219 221
220 matrx[9,x] = round(f(tmp$nonSilentMutationsCDR, na.rm=T), digits=1) 222 matrx[9,x] = round(f(tmp$nonSilentMutationsCDR, na.rm=T), digits=1)
221 matrx[9,y] = round(f(tmp$silentMutationsCDR, na.rm=T), digits=1) 223 matrx[9,y] = round(f(tmp$silentMutationsCDR, na.rm=T), digits=1)
222 matrx[9,z] = round(matrx[9,x] / matrx[9,y], digits=1) 224 matrx[9,z] = round(matrx[9,x] / matrx[9,y], digits=1)
223 225
224 if(fname == "sum"){ 226 if(fname == "sum"){
225 matrx[10,x] = round(f(rowSums(tmp[,c("FR2.IMGT.Nb.of.nucleotides", "FR3.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1) 227 matrx[10,x] = round(f(rowSums(tmp[,c("FR2.IMGT.Nb.of.nucleotides", "FR3.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1)
226 matrx[10,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1) 228 matrx[10,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
227 matrx[10,z] = round(matrx[10,x] / matrx[10,y], digits=1) 229 matrx[10,z] = round(matrx[10,x] / matrx[10,y], digits=1)
228 230
229 matrx[11,x] = round(f(rowSums(tmp[,c("CDR1.IMGT.Nb.of.nucleotides", "CDR2.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1) 231 matrx[11,x] = round(f(rowSums(tmp[,c("CDR1.IMGT.Nb.of.nucleotides", "CDR2.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1)
230 matrx[11,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1) 232 matrx[11,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
231 matrx[11,z] = round(matrx[11,x] / matrx[11,y], digits=1) 233 matrx[11,z] = round(matrx[11,x] / matrx[11,y], digits=1)
232 } 234 }
233 } 235 }
234 236
235 transitionTable = data.frame(A=zeros,C=zeros,G=zeros,T=zeros) 237 transitionTable = data.frame(A=zeros,C=zeros,G=zeros,T=zeros)
236 row.names(transitionTable) = c("A", "C", "G", "T") 238 row.names(transitionTable) = c("A", "C", "G", "T")
237 transitionTable["A","A"] = NA 239 transitionTable["A","A"] = NA
238 transitionTable["C","C"] = NA 240 transitionTable["C","C"] = NA
239 transitionTable["G","G"] = NA 241 transitionTable["G","G"] = NA
240 transitionTable["T","T"] = NA 242 transitionTable["T","T"] = NA
241 243
242 if(nrow(tmp) > 0){ 244 if(nrow(tmp) > 0){
243 for(nt1 in nts){ 245 for(nt1 in nts){
244 for(nt2 in nts){ 246 for(nt2 in nts){
245 if(nt1 == nt2){ 247 if(nt1 == nt2){
246 next 248 next
247 } 249 }
257 } else { 259 } else {
258 transitionTable[NT1,NT2] = sum(tmp[,c(CDR1, FR2, CDR2, FR3)]) 260 transitionTable[NT1,NT2] = sum(tmp[,c(CDR1, FR2, CDR2, FR3)])
259 } 261 }
260 } 262 }
261 } 263 }
262 } 264 transition = transitionTable
263 265 transition$id = names(transition)
264 266
265 print(paste("writing value file: ", name, "_", fname, "_value.txt" ,sep="")) 267 transition2 = melt(transition, id.vars="id")
266 268
267 write.table(x=transitionTable, file=paste("transitions_", name ,"_", fname, ".txt", sep=""), sep=",",quote=F,row.names=T,col.names=NA) 269 transition2 = merge(transition2, base.order, by.x="id", by.y="base")
268 write.table(x=tmp[,c("Sequence.ID", "best_match", "chunk_hit_percentage", "nt_hit_percentage", "start_locations")], file=paste("matched_", name , "_", fname, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T) 270 transition2 = merge(transition2, base.order, by.x="variable", by.y="base")
269 271
270 cat(matrx[1,x], file=paste(name, "_", fname, "_value.txt" ,sep="")) 272 transition2[is.na(transition2$value),]$value = 0
271 cat(nrow(tmp), file=paste(name, "_", fname, "_n.txt" ,sep="")) 273
272 274 png(filename=paste("transitions_stacked_", name, ".png", sep=""))
273 print(paste(fname, name, nrow(tmp))) 275 p = ggplot(transition2, aes(factor(reorder(id, order.x)), y=value, fill=factor(reorder(variable, order.y)))) + geom_bar(position="fill", stat="identity") #stacked bar
274 276 p = p + xlab("From base") + ylab("To base") + ggtitle("Mutations frequency from base to base") + guides(fill=guide_legend(title=NULL))
275 matrx 277 print(p)
278 dev.off()
279
280 png(filename=paste("transitions_heatmap_", name, ".png", sep=""))
281 p = ggplot(transition2, aes(factor(reorder(id, order.x)), factor(reorder(variable, order.y)))) + geom_tile(aes(fill = value), colour="white") + scale_fill_gradient(low="white", high="steelblue") #heatmap
282 p = p + xlab("From base") + ylab("To base") + ggtitle("Mutations frequency from base to base")
283 print(p)
284 dev.off()
285 }
286
287 #print(paste("writing value file: ", name, "_", fname, "_value.txt" ,sep=""))
288
289 write.table(x=transitionTable, file=paste("transitions_", name ,"_", fname, ".txt", sep=""), sep=",",quote=F,row.names=T,col.names=NA)
290 write.table(x=tmp[,c("Sequence.ID", "best_match", "chunk_hit_percentage", "nt_hit_percentage", "start_locations")], file=paste("matched_", name , "_", fname, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T)
291
292 cat(matrx[1,x], file=paste(name, "_", fname, "_value.txt" ,sep=""))
293 cat(nrow(tmp), file=paste(name, "_", fname, "_n.txt" ,sep=""))
294
295 #print(paste(fname, name, nrow(tmp)))
296
297 matrx
276 } 298 }
277 299
278 nts = c("a", "c", "g", "t") 300 nts = c("a", "c", "g", "t")
279 zeros=rep(0, 4) 301 zeros=rep(0, 4)
280 302
320 342
321 #sum.table = sum.table[c("Number of Mutations (%)", "Median of Number of Mutations (%)", "Transition (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C G (%)", "Transitions at A T (%)", "Targeting of A T (%)", "FR R/S (ratio)", "CDR R/S (ratio)", "nt in FR", "nt in CDR"),] 343 #sum.table = sum.table[c("Number of Mutations (%)", "Median of Number of Mutations (%)", "Transition (%)", "Transversions (%)", "Transitions at G C (%)", "Targeting of C G (%)", "Transitions at A T (%)", "Targeting of A T (%)", "FR R/S (ratio)", "CDR R/S (ratio)", "nt in FR", "nt in CDR"),]
322 344
323 write.table(x=new.table, file="mutations_sum.txt", sep=",",quote=F,row.names=F,col.names=F) 345 write.table(x=new.table, file="mutations_sum.txt", sep=",",quote=F,row.names=F,col.names=F)
324 346
325
326
327 if (!("ggplot2" %in% rownames(installed.packages()))) {
328 install.packages("ggplot2", repos="http://cran.xl-mirror.nl/")
329 }
330
331 dat = dat[!grepl("^unmatched", dat$best_match),] 347 dat = dat[!grepl("^unmatched", dat$best_match),]
332 348
333 #blegh 349 #blegh
334 genesForPlot = dat[grepl("ca", dat$best_match),]$best_match 350 genesForPlot = dat[grepl("ca", dat$best_match),]$best_match
335 if(length(genesForPlot) > 0){ 351 if(length(genesForPlot) > 0){