rm(list=ls())
require("shape")
require("RColorBrewer")
setwd("/Users/Helena/Documents/UWOceanography/FISH 554/Class assignment/Fig2/")

d1 = read.table("/Users/claireellis/Desktop/methratio_out_CG5x_M1_fig_test.txt", header=F)  								# Table of CDS and tRNA with start position, end position, GC, and annotation
d2 = read.table("fig2_data2.txt", header=F)									# Same for intergenic regions
head(d1)
d1$annot <- NA
colnames(d1) <- c("type","start","end","meth","annot")
colnames(d2) <- c("type","start","end","GC","annot")
d1$start_rad <- (2*pi*d1$start)/1948860									# 2952962 is the total number of bp in the genome, you need to edit this
d1$end_rad <-(2*pi*d1$end)/1948860  
d2$start_rad <- (2*pi*d2$start)/1948860
d2$end_rad <-(2*pi*d2$end)/1948860

d <- rbind(d1,d2)															# combine the genes and intergenic regions if you want to plot everything (not a great idea, I tried it)
d<- d[order(d$start),]

high_GC = max(d1$meth)							# picks out the indices of genes with especially high or low GC
low_GC = min(d1$meth)

B2R<- colorRampPalette(brewer.pal(9,"Blues"))(100)
d1$color <- B2R[round((d1$meth*50)/mean(d1$meth),0)]							# Pick a color to plot based on how different the GC is from the mean


# plot

par(mar=c(1,1,1,1),oma=c(1,1,1,1))
mat<- matrix(c(1,2),ncol=2,nrow=1)
layout(mat, heights=c(1), widths=c(170,30))


plot(1,xlim=c(-1.5, 1.5),ylim=c(-1.5, 1.5),type='n', axes=F, asp=1)			# the shape package requires an initial blank plot, change the xlim and ylim parameters to zoom in or out


for (row in 1:nrow(d1)) {
  filledcircle(r1 = 1.3, r2=.9, mid = c(0,0), dr = 0.05, col = d1$color[row], 
               from=d1$start_rad[row], to=d1$end_rad[row])
  #   if (row %in% c(high_GC,low_GC)) {										# this commented out section will plot the annotations of genes that have particularly high or low GC
  #     x = 1.5*cos(median(c(d1$start_rad[row],d1$end_rad[row])))
  #     y = 1.5*sin(median(c(d1$start_rad[row],d1$end_rad[row])))
  #     if (row %in% high_GC) {
  #       label = c("dark red")
  #     }
  #     if (row %in% low_GC) {
  #       label = c("dark blue")
  #     }
  #     if (x > 0 & y > 0) {
  #       shift = c(0,0)
  #     }
  #     if (x < 0 & y < 0) {
  #       shift = c(1,1)
  #     }
  #     if (x > 0 & y < 0) {
  #       shift = c(0,1)
  #     }
  #     if (x < 0 & y > 0) {
  #       shift = c(1,0)
  #     }
  #     text(x,y,as.character(d1$annot[row]), cex=.6, col=label,srt=d1$start_rad[row]*180/pi, adj=c(0,0))
  #   }
}
text(0,0,expression(italic("Male 1")),cex=2)
text(0,-.2,"Spermatozoa",cex=2)

# legend

legend_image <- as.raster(matrix(rev(B2R)[2:67], ncol=1))
plot(c(0,2),c(0,1),type = 'n', axes = F,xlab = '', ylab = '')
mtext(text='%meth',side=3,cex=1)
text(x=1.5, y = c(0,.25,1), labels = round(c(0),mean(d1$meth),max(d1$meth)),1),cex=1)
rasterImage(legend_image, 0, 0, 1, 1)