# plot the results of CpG island prediction

#define some rgb colours
rgb <- c("#000000", "#E69F00", "#56B4E9", "#009E73", 
    "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# read data which is output of Perl script
data <- read.table("cpg_chr4.out", sep = "\t", header = TRUE)

# make an emtpy plot
plot(0, type = "n", xlim = c(72009075, 72009075 + 
    1520980), ylim = c(1.6, 2.3), xlab = "Position", ylab = "", 
    yaxt = "n", main = "CpG island prediction")

# plot the predicted CpG islands
for (i in 1:length(data$pos)) {
    if (data$cpg[i] == 1) {
        # convert the position numbers to chromosomal positions
        data$pos[i] <- data$pos[i] + 72009075
        #print (data$pos[i])
        lines(c(data$pos[i], data$pos[i]), c(1.7, 1.8), col = rgb[7])
    }
}

# Read file with chr4:72,009,075-73,530,562 region
# annotation. This information was obtained with 
# the Table Browser of the UCSC browser

annot <- read.table("chr4_annotation.txt", sep = "\t", 
    header = TRUE)

color <- 1
prevname <- ""
lines <- length(annot$chr)  # number of lines in the annotation file

for (i in (1:lines)) {
    
    if (annot$category[i] == "exon") {
        # if we consider an exon
        # if it a different gene as compared to the previous line,
        # change the colour
        if (annot$gene_id[i] != prevname) {
            color <- color + 1
        }
        prevname <- annot$gene_id[i]
        
        # draw rectangles for the exons
        rect(annot$beg[i], 1.9, annot$end[i], 2.1, border = rgb[color], 
            col = rgb[color])
    }
    
    if (annot$category[i] == "trans") {
        # to identify the end points of transcripts
        lines(c(annot$beg[i], annot$end[i]), c(2, 2), col = "grey", 
            lw = 2)
        direction <- annot$strand[i]
        if (direction == "+") {
            dir <- 2
        }
        if (direction == "-") {
            dir <- 1
        }
        # print arrows to indicate the location 
        # and direction of transcript
        
        arrows(annot$beg[i], 1.85, annot$end[i], 1.85, col = "grey", 
            code = dir, lw = 5, length = 0.1)
    }
}

# print names of genes (this information
# can not be extracted from the
# chr4_annotation.txt file)
text(72009075 + 2e+05, 2.2, "SLC4A4")
text(72009075 + 610000, 2.2, "GC")
text(72009075 + 950000, 2.2, "NPFFR2")
text(72009075 + 1300000, 2.2, "ADAMTS3")