# 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)