Mercurial > repos > davidvanzessen > argalaxy_tools

comparison pattern_plots.r @ 4:5ffd52fc35c4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Mon, 12 Dec 2016 05:22:37 -0500
parents
children
comparison
equal deleted inserted replaced
3:beaa487ecf43 4:5ffd52fc35c4
1 library(ggplot2)
2 library(reshape2)
3 library(scales)
4
5 args <- commandArgs(trailingOnly = TRUE)
6
7 input.file = args[1] #the data that's get turned into the "SHM overview" table in the html report "data_sum.txt"
8
9 plot1.path = args[2]
10 plot1.png = paste(plot1.path, ".png", sep="")
11 plot1.txt = paste(plot1.path, ".txt", sep="")
12
13 plot2.path = args[3]
14 plot2.png = paste(plot2.path, ".png", sep="")
15 plot2.txt = paste(plot2.path, ".txt", sep="")
16
17 plot3.path = args[4]
18 plot3.png = paste(plot3.path, ".png", sep="")
19 plot3.txt = paste(plot3.path, ".txt", sep="")
20
21 clean.output = args[5]
22
23 dat = read.table(input.file, header=F, sep=",", quote="", stringsAsFactors=F, fill=T, row.names=1)
24
25
26
27 classes = c("IGA", "IGA1", "IGA2", "IGG", "IGG1", "IGG2", "IGG3", "IGG4", "IGM", "IGE")
28 xyz = c("x", "y", "z")
29 new.names = c(paste(rep(classes, each=3), xyz, sep="."), paste("un", xyz, sep="."), paste("all", xyz, sep="."))
30
31 names(dat) = new.names
32
33 clean.dat = dat
34 clean.dat = clean.dat[,c(paste(rep(classes, each=3), xyz, sep="."), paste("all", xyz, sep="."), paste("un", xyz, sep="."))]
35
36 write.table(clean.dat, clean.output, quote=F, sep="\t", na="", row.names=T, col.names=NA)
37
38 dat["RGYW.WRCY",] = colSums(dat[c(13,14),], na.rm=T)
39 dat["TW.WA",] = colSums(dat[c(15,16),], na.rm=T)
40
41 data1 = dat[c("RGYW.WRCY", "TW.WA"),]
42
43 data1 = data1[,names(data1)[grepl(".z", names(data1))]]
44 names(data1) = gsub("\\..*", "", names(data1))
45
46 data1 = melt(t(data1))
47
48 names(data1) = c("Class", "Type", "value")
49
50 data1 = data1[order(data1$Type),]
51
52 write.table(data1, plot1.txt, quote=F, sep="\t", na="", row.names=F, col.names=T)
53
54 p = ggplot(data1, aes(Class, value)) + geom_bar(aes(fill=Type), stat="identity", position="dodge", colour = "black") + ylab("% of mutations") + guides(fill=guide_legend(title=NULL))
55 p = p + theme(panel.background = element_rect(fill = "white", colour="black"),text = element_text(size=15, colour="black"), axis.text.x = element_text(angle = 45, hjust = 1)) + scale_fill_manual(values=c("RGYW.WRCY" = "white", "TW.WA" = "blue4"))
56 #p = p + scale_colour_manual(values=c("RGYW.WRCY" = "black", "TW.WA" = "blue4"))
57 png(filename=plot1.png, width=480, height=300)
58 print(p)
59 dev.off()
60
61 data2 = dat[c(1, 5:8),]
62
63 data2 = data2[,names(data2)[grepl("\\.x", names(data2))]]
64 names(data2) = gsub(".x", "", names(data2))
65
66 data2["A/T",] = dat["Targeting of A T (%)",names(dat)[grepl("\\.z", names(dat))]]
67
68 data2["G/C transitions",] = round(data2["Transitions at G C (%)",] / data2["Number of Mutations (%)",] * 100, 1)
69
70 data2["mutation.at.gc",] = dat["Transitions at G C (%)",names(dat)[grepl("\\.y", names(dat))]]
71 data2["G/C transversions",] = round((data2["mutation.at.gc",] - data2["Transitions at G C (%)",]) / data2["Number of Mutations (%)",] * 100, 1)
72
73 data2["G/C transversions",is.nan(unlist(data2["G/C transversions",]))] = 0
74 data2["G/C transversions",is.infinite(unlist(data2["G/C transversions",]))] = 0
75 data2["G/C transitions",is.nan(unlist(data2["G/C transitions",]))] = 0
76 data2["G/C transitions",is.infinite(unlist(data2["G/C transitions",]))] = 0
77
78 data2 = melt(t(data2[c("A/T","G/C transitions","G/C transversions"),]))
79
80 names(data2) = c("Class", "Type", "value")
81
82 data2 = data2[order(data2$Type),]
83
84 write.table(data2, plot2.txt, quote=F, sep="\t", na="", row.names=F, col.names=T)
85
86 p = ggplot(data2, aes(x=Class, y=value, fill=Type)) + geom_bar(position="fill", stat="identity", colour = "black") + scale_y_continuous(labels=percent_format()) + guides(fill=guide_legend(title=NULL)) + ylab("% of mutations")
87 p = p + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=15, colour="black"), axis.text.x = element_text(angle = 45, hjust = 1)) + scale_fill_manual(values=c("A/T" = "blue4", "G/C transversions" = "gray74", "G/C transitions" = "white"))
88 #p = p + scale_colour_manual(values=c("A/T" = "blue4", "G/C transversions" = "gray74", "G/C transitions" = "black"))
89 png(filename=plot2.png, width=480, height=300)
90 print(p)
91 dev.off()
92
93 data3 = dat[c(5, 6, 8, 17:20),]
94 data3 = data3[,names(data3)[grepl("\\.x", names(data3))]]
95 names(data3) = gsub(".x", "", names(data3))
96
97 data3[is.na(data3)] = 0
98 #data3[is.infinite(data3)] = 0
99
100 data3["G/C transitions",] = round(data3["Transitions at G C (%)",] / (data3["C",] + data3["G",]) * 100, 1)
101
102 data3["G/C transversions",] = round((data3["Targeting of G C (%)",] - data3["Transitions at G C (%)",]) / (data3["C",] + data3["G",]) * 100, 1)
103
104 data3["A/T",] = round(data3["Targeting of A T (%)",] / (data3["A",] + data3["T",]) * 100, 1)
105
106 data3["G/C transitions",is.nan(unlist(data3["G/C transitions",]))] = 0
107 data3["G/C transitions",is.infinite(unlist(data3["G/C transitions",]))] = 0
108
109 data3["G/C transversions",is.nan(unlist(data3["G/C transversions",]))] = 0
110 data3["G/C transversions",is.infinite(unlist(data3["G/C transversions",]))] = 0
111
112 data3["A/T",is.nan(unlist(data3["A/T",]))] = 0
113 data3["A/T",is.infinite(unlist(data3["A/T",]))] = 0
114
115 data3 = melt(t(data3[8:10,]))
116 names(data3) = c("Class", "Type", "value")
117
118 data3 = data3[order(data3$Type),]
119
120 write.table(data3, plot3.txt, quote=F, sep="\t", na="", row.names=F, col.names=T)
121
122 p = ggplot(data3, aes(Class, value)) + geom_bar(aes(fill=Type), stat="identity", position="dodge", colour = "black") + ylab("% of nucleotides") + guides(fill=guide_legend(title=NULL))
123 p = p + theme(panel.background = element_rect(fill = "white", colour="black"), text = element_text(size=15, colour="black"), axis.text.x = element_text(angle = 45, hjust = 1)) + scale_fill_manual(values=c("A/T" = "blue4", "G/C transversions" = "gray74", "G/C transitions" = "white"))
124 #p = p + scale_colour_manual(values=c("A/T" = "blue4", "G/C transversions" = "gray74", "G/C transitions" = "black"))
125 png(filename=plot3.png, width=480, height=300)
126 print(p)
127 dev.off()
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159