Mercurial > repos > davidvanzessen > shm_csr

comparison shm_csr.r @ 24:9955b23db68c draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Fri, 02 Dec 2016 05:51:13 -0500
parents 81453585dfc3
children 8d1c4c75f81b
comparison
equal deleted inserted replaced
23:81453585dfc3 24:9955b23db68c
175 175
176 setwd(outputdir) 176 setwd(outputdir)
177 177
178 write.table(dat, input, sep="\t",quote=F,row.names=F,col.names=T) 178 write.table(dat, input, sep="\t",quote=F,row.names=F,col.names=T)
179 179
180 base.order = data.frame(base=c("A", "T", "C", "G"), order=1:4) 180 base.order.x = data.frame(base=c("A", "C", "G", "T"), order.x=1:4)
181 base.order.y = data.frame(base=c("T", "G", "C", "A"), order.y=1:4)
181 182
182 calculate_result = function(i, gene, dat, matrx, f, fname, name){ 183 calculate_result = function(i, gene, dat, matrx, f, fname, name){
183 tmp = dat[grepl(paste("^", gene, ".*", sep=""), dat$best_match),] 184 tmp = dat[grepl(paste("^", gene, ".*", sep=""), dat$best_match),]
184 185
185 j = i - 1 186 j = i - 1
289 transition = transitionTable 290 transition = transitionTable
290 transition$id = names(transition) 291 transition$id = names(transition)
291 292
292 transition2 = melt(transition, id.vars="id") 293 transition2 = melt(transition, id.vars="id")
293 294
294 transition2 = merge(transition2, base.order, by.x="id", by.y="base") 295 transition2 = merge(transition2, base.order.x, by.x="id", by.y="base")
295 296
296 transition2 = merge(transition2, base.order, by.x="variable", by.y="base") 297 transition2 = merge(transition2, base.order.y, by.x="variable", by.y="base")
297 298
298 transition2[is.na(transition2$value),]$value = 0 299 transition2[is.na(transition2$value),]$value = 0
299 300
300 if(any(transition2$value != 0)){ #having rows of data but a transition table filled with 0 is bad 301 if(any(transition2$value != 0)){ #having a transition table filled with 0 is bad
301 print("Plotting heatmap and transition") 302 print("Plotting heatmap and transition")
302 png(filename=paste("transitions_stacked_", name, ".png", sep="")) 303 png(filename=paste("transitions_stacked_", name, ".png", sep=""))
303 p = ggplot(transition2, aes(factor(reorder(id, order.x)), y=value, fill=factor(reorder(variable, order.y)))) + geom_bar(position="fill", stat="identity", colour="black") #stacked bar 304 p = ggplot(transition2, aes(factor(reorder(id, order.x)), y=value, fill=factor(reorder(variable, order.y)))) + geom_bar(position="fill", stat="identity", colour="black") #stacked bar
304 p = p + xlab("From base") + ylab("To base") + ggtitle("Mutations frequency from base to base") + guides(fill=guide_legend(title=NULL)) 305 p = p + xlab("From base") + ylab("") + ggtitle("Mutations frequency from base to base") + guides(fill=guide_legend(title=NULL))
305 p = p + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=16, colour="black")) + scale_fill_manual(values=c("A" = "blue4", "G" = "lightblue1", "C" = "olivedrab3", "T" = "olivedrab4")) 306 p = p + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=16, colour="black")) + scale_fill_manual(values=c("A" = "blue4", "G" = "lightblue1", "C" = "olivedrab3", "T" = "olivedrab4"))
306 #p = p + scale_colour_manual(values=c("A" = "black", "G" = "black", "C" = "black", "T" = "black")) 307 #p = p + scale_colour_manual(values=c("A" = "black", "G" = "black", "C" = "black", "T" = "black"))
307 print(p) 308 print(p)
308 dev.off() 309 dev.off()
309 png(filename=paste("transitions_heatmap_", name, ".png", sep="")) 310 png(filename=paste("transitions_heatmap_", name, ".png", sep=""))
310 p = ggplot(transition2, aes(factor(reorder(id, order.x)), factor(reorder(variable, order.y)))) + geom_tile(aes(fill = value)) + scale_fill_gradient(low="white", high="steelblue") #heatmap 311 p = ggplot(transition2, aes(factor(reorder(variable, -order.y)), factor(reorder(id, -order.x)))) + geom_tile(aes(fill = value)) + scale_fill_gradient(low="white", high="steelblue") #heatmap
311 p = p + xlab("From base") + ylab("To base") + ggtitle("Mutations frequency from base to base") + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=13, colour="black")) 312 p = p + xlab("To base") + ylab("From Base") + ggtitle("Mutations frequency from base to base") + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=13, colour="black"))
312 print(p) 313 print(p)
313 dev.off() 314 dev.off()
314 } else { 315 } else {
315 #print("No data to plot") 316 #print("No data to plot")
316 } 317 }