Mercurial > repos > davidvanzessen > shm_csr

comparison pattern_plots.r @ 23:81453585dfc3 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Thu, 01 Dec 2016 09:32:06 -0500
parents 012a738edf5a
children 05c62efdc393
comparison
equal deleted inserted replaced
22:0bea8c187a90 23:81453585dfc3
16 16
17 plot3.path = args[4] 17 plot3.path = args[4]
18 plot3.png = paste(plot3.path, ".png", sep="") 18 plot3.png = paste(plot3.path, ".png", sep="")
19 plot3.txt = paste(plot3.path, ".txt", sep="") 19 plot3.txt = paste(plot3.path, ".txt", sep="")
20 20
21 clean.output = args[5]
22
21 dat = read.table(input.file, header=F, sep=",", quote="", stringsAsFactors=F, fill=T, row.names=1) 23 dat = read.table(input.file, header=F, sep=",", quote="", stringsAsFactors=F, fill=T, row.names=1)
22 24
23 25
24 26
25 classes = c("IGA", "IGA1", "IGA2", "IGG", "IGG1", "IGG2", "IGG3", "IGG4", "IGM", "IGE") 27 classes = c("IGA", "IGA1", "IGA2", "IGG", "IGG1", "IGG2", "IGG3", "IGG4", "IGM", "IGE")
26 xyz = c("x", "y", "z") 28 xyz = c("x", "y", "z")
27 new.names = c(paste(rep(classes, each=3), xyz, sep="."), paste("un", xyz, sep="."), paste("all", xyz, sep=".")) 29 new.names = c(paste(rep(classes, each=3), xyz, sep="."), paste("un", xyz, sep="."), paste("all", xyz, sep="."))
28 30
29 names(dat) = new.names 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)
30 37
31 dat["RGYW.WRCY",] = colSums(dat[c(13,14),], na.rm=T) 38 dat["RGYW.WRCY",] = colSums(dat[c(13,14),], na.rm=T)
32 dat["TW.WA",] = colSums(dat[c(15,16),], na.rm=T) 39 dat["TW.WA",] = colSums(dat[c(15,16),], na.rm=T)
33 40
34 data1 = dat[c("RGYW.WRCY", "TW.WA"),] 41 data1 = dat[c("RGYW.WRCY", "TW.WA"),]
49 #p = p + scale_colour_manual(values=c("RGYW.WRCY" = "black", "TW.WA" = "blue4")) 56 #p = p + scale_colour_manual(values=c("RGYW.WRCY" = "black", "TW.WA" = "blue4"))
50 png(filename=plot1.png, width=480, height=300) 57 png(filename=plot1.png, width=480, height=300)
51 print(p) 58 print(p)
52 dev.off() 59 dev.off()
53 60
54 data2 = dat[5:8,] 61 data2 = dat[c(1, 5:8),]
55
56 data2["sum",] = colSums(data2, na.rm=T)
57 62
58 data2 = data2[,names(data2)[grepl("\\.x", names(data2))]] 63 data2 = data2[,names(data2)[grepl("\\.x", names(data2))]]
59 names(data2) = gsub(".x", "", names(data2)) 64 names(data2) = gsub(".x", "", names(data2))
60 65
61 data2["A/T",] = round(colSums(data2[3:4,]) / data2["sum",] * 100, 1) 66 data2["A/T",] = dat["Targeting of A T (%)",names(dat)[grepl("\\.z", names(dat))]]
62 data2["A/T",is.nan(unlist(data2["A/T",]))] = 0
63 67
64 data2["G/C transversions",] = round(data2[2,] / data2["sum",] * 100, 1) 68 data2["G/C transitions",] = round(data2["Transitions at G C (%)",] / data2["Number of Mutations (%)",] * 100, 1)
65 data2["G/C transitions",] = round(data2[1,] / data2["sum",] * 100, 1)
66 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)
67 72
68 data2["G/C transversions",is.nan(unlist(data2["G/C transversions",]))] = 0 73 data2["G/C transversions",is.nan(unlist(data2["G/C transversions",]))] = 0
69 data2["G/C transversions",is.infinite(unlist(data2["G/C transversions",]))] = 0 74 data2["G/C transversions",is.infinite(unlist(data2["G/C transversions",]))] = 0
70 data2["G/C transitions",is.nan(unlist(data2["G/C transitions",]))] = 0 75 data2["G/C transitions",is.nan(unlist(data2["G/C transitions",]))] = 0
71 data2["G/C transitions",is.infinite(unlist(data2["G/C transitions",]))] = 0 76 data2["G/C transitions",is.infinite(unlist(data2["G/C transitions",]))] = 0
72 77
73 data2 = melt(t(data2[6:8,])) 78 data2 = melt(t(data2[c("A/T","G/C transitions","G/C transversions"),]))
74 79
75 names(data2) = c("Class", "Type", "value") 80 names(data2) = c("Class", "Type", "value")
76 81
77 data2 = data2[order(data2$Type),] 82 data2 = data2[order(data2$Type),]
78 83
90 names(data3) = gsub(".x", "", names(data3)) 95 names(data3) = gsub(".x", "", names(data3))
91 96
92 data3[is.na(data3)] = 0 97 data3[is.na(data3)] = 0
93 #data3[is.infinite(data3)] = 0 98 #data3[is.infinite(data3)] = 0
94 99
95 data3["G/C transitions",] = round(data3[1,] / (data3[5,] + data3[7,]) * 100, 1) 100 data3["G/C transitions",] = round(data3["Transitions at G C (%)",] / (data3["C",] + data3["G",]) * 100, 1)
96 101
97 data3["G/C transversions",] = round(data3[2,] / (data3[5,] + data3[7,]) * 100, 1) 102 data3["G/C transversions",] = round((data3["Targeting of G C (%)",] - data3["Transitions at G C (%)",]) / (data3["C",] + data3["G",]) * 100, 1)
98 103
99 data3["A/T",] = round(data3[3,] / (data3[4,] + data3[6,]) * 100, 1) 104 data3["A/T",] = round(data3["Targeting of A T (%)",] / (data3["A",] + data3["T",]) * 100, 1)
100 105
101 data3["G/C transitions",is.nan(unlist(data3["G/C transitions",]))] = 0 106 data3["G/C transitions",is.nan(unlist(data3["G/C transitions",]))] = 0
102 data3["G/C transitions",is.infinite(unlist(data3["G/C transitions",]))] = 0 107 data3["G/C transitions",is.infinite(unlist(data3["G/C transitions",]))] = 0
103 108
104 data3["G/C transversions",is.nan(unlist(data3["G/C transversions",]))] = 0 109 data3["G/C transversions",is.nan(unlist(data3["G/C transversions",]))] = 0