# plot the results of cpg island prediction # define some colours rgb <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7") # read file being output from Perl script data <- read.table("cpg_short.out", sep = "\t", header = TRUE) # Specify the number of lines of margin to the four sides # of the plot. In this case we want to make room for text # at the right axis. par(mar = c(5, 4, 4, 5) + 0.1) # make first plot, which is empty plot(data$pos, data$cg_obs_exp, type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "") # plot lines to indicate where CpG islands are predicted for (i in 1:length(data$pos)) { if (data$cpg[i] == 1) { lines(c(data$pos[i], data$pos[i]), c(0, 1), col = rgb[2]) } } # before another plot, prevent R from clearing # the graphics device par(new = TRUE) # make second plot with the cg_obs_exp data plot(data$pos, data$cg_obs_exp, type = "l", main = "CpG island prediction", xlab = "Position", ylab = "CpG obs/exp", col = rgb[7]) # print legend legend(8000, 0.9, c("CG ratio", "CpG obs/exp"), col = c(rgb[6], rgb[7]), lwd = 2) par(new = TRUE) # make third plot plot(data$pos, data$cg_ratio, type = "l", xaxt = "n", yaxt = "n", xlab = "", ylab = "", col = rgb[6]) # print ticks for the 2nd y axis axis(4) # print text to 2nd y axis mtext("CpG ratio", side = 4, line = 3)