modified: makesrc3.mk
[GalaxyCodeBases.git] / BioInfo / BS-Seq / bwa-meth / compare / src / plot-quals.R
blob7c3005e2e80b8660457a9c36cbca1485edbc338c
1 # make supplemental figures for paper
2 options(stringsAsFactors=FALSE)
3 library(ggplot2)
5 #X11.options(antialias="subpixel", type="cairo")
7 args = commandArgs(TRUE)
9 df = read.delim(args[1])
10 out_eps_or_png = args[2]
12 df[df$method == "bis2", "method"] = "bismark-bt2"
13 df[df$method == "bis1", "method"] = "bismark-bt1"
14 df[df$method == "bwastrand", "method"] = "bwa-strand"
15 df[df$method == "bwa", "method"] = "bwameth"
17 df[grep("sim_R1", df$method, fixed=TRUE), "method"] = "bison"
18 df[grep("real_R1", df$method, fixed=TRUE), "method"] = "bison"
19 df$method = factor(df$method)#, levels=c("last", "bsmap", "gsnap", "bwa", "bismark", "bwa-strand", "bsmooth", "bison"))
20 df$method = factor(as.character(df$method), levels=rev(levels(df$method)))
21 df = df[df$qual > 0,]
23 df = df[order(df$qual),]
25 df = df[(df$on + df$off) != 0,]
26 df$on = df$on * 100
27 df$off = df$off * 100
29 require("grid")
30 print(dim(df))
32 SIM=any(grep("sim", args))
33 #         geom_point(color="grey60", size=ifelse(any(grep("sim", args)), 1.26, 0.2)) +
35 p = ggplot(df, aes(x=off, y=on, by=method)) +
36          geom_line(aes(color=method), linetype="dashed", linewidth=ifelse(SIM, 0.1, 0.25) ) +
37          geom_point(aes(color=method), size=ifelse(SIM, 1.0, 0.55)) +
38          scale_shape(solid = FALSE) 
39 p = p + ylab("% Reads On Target")
40 p = p + xlab("% Reads Off Target")
41 p = p + scale_color_brewer(palette="Set1")
42 if(any(grep("sim", args))){
43     p = p + xlim(xmin=0, xmax=1.5)
44     p = p + ylim(ymin=65, ymax=92)
45     p = p + geom_point(data=df[df$method %in% c('bismark-bt1', 'bismark-bt2', 'gsnap', "bsmap"),],
46                        aes(x=off, y=on, by=method, color=method))
47 } else {
48     # make the points larger, but don't try this for the bwa strand comparison
49     if(!any(grep("strand", args))){
50         p = p + geom_point(data=df[df$method %in% c('bismark-bt1', 'bismark-bt2', 'gsnap', "bsmap"),],
51                        aes(x=off, y=on, by=method, color=method))
52     }
53     p = p + xlim(xmin=4, xmax=14)
54     p = p + ylim(ymin=55, ymax=80)
56 #p = p + theme_bw()
57 p = p + theme(
58              #legend.position = c(0.55, 0.25),
59              #legend.position = c(ifelse(SIM, 0.69, 0.34), 0.25),
60              #legend.position = ifelse(SIM, c(0.69, 0.25), c(0.29 ,0.75)),
61              legend.position = c(0.71 ,0.25),
62               legend.text=element_text(size=5, lineheight=4),
63               axis.text=element_text(size=5),
64               axis.title=element_text(size=7),
65               legend.key.size=unit(5, "mm")
67 #              legend.key.height=5
68               )
69 if(any(grep("strand", args))){
70    p = p + theme(
71              legend.position = c(ifelse(SIM, 0.69, 0.34), 0.25)
72     )
74 p = p + guides(color=guide_legend(ncol=2, title=NULL))
75 #print(p)
76 ggsave(file=out_eps_or_png, units="cm", width=8.6, height=6.3,
77     dpi=1200)