nanoPARE_colors = c(BODY='#808285',CAP='#1C75BC',UNCAP="#F7941E",DIGEST='',XRN4='')

chromosome_lengths = lapply(FASTA,nchar)
binsize = 50000
smallbinsize = 5000

sampletype = 'fb'
folder = 'L:/members/Schon/Writing/nanoPARE/data_generated/'
capped_file = paste(folder,'capped_features/',sampletype,'.W.capped.bed',sep='')
noncapped_file = paste(folder,'noncapped_features/',sampletype,'.W.noncapped.bed',sep='')
capped_features = read.table(capped_file,stringsAsFactors = F)
noncapped_features = read.table(noncapped_file,stringsAsFactors = F)
# 
# capped_features = read.table('L:/members/Schon/Writing/nanoPARE/manuscript_draft/manuscript_data/fb.all_features.capped.bed',stringsAsFactors = F)
# noncapped_features = read.table('L:/members/Schon/Writing/nanoPARE/manuscript_draft/manuscript_data/fb.all_features.noncapped.bed',stringsAsFactors = F)

c_chroms = toupper(gsub('Ath_chr','',capped_features[,1]))
c_positions = capped_features[,2] + capped_features[,7]

n_chroms = toupper(gsub('Ath_chr','',noncapped_features[,1]))
n_positions = noncapped_features[,2] + noncapped_features[,7]

capped_chroms = list()
noncapped_chroms = list()

for(chrom in names(chromosome_lengths)){
  if(chrom%in%c('M','C')){
    b = smallbinsize
  }else{
    b = binsize
  }
  cpos = c_positions[c_chroms == chrom]
  npos = n_positions[n_chroms == chrom]
  numbins = ceiling(chromosome_lengths[[chrom]]/b)
  cbins = vector(length=numbins)
  nbins = vector(length=numbins)
  startpos = 1
  for(i in 1:numbins){
    posrange = c(startpos,startpos + b)
    startpos = posrange[length(posrange)]
    cbins[i] = sum(cpos>=posrange[1] & cpos<=posrange[2])
    nbins[i] = sum(npos>=posrange[1] & npos<=posrange[2])
  }
  capped_chroms[[chrom]] = cbins
  noncapped_chroms[[chrom]] = nbins
}

pdf(
  paste("chromosome_metaplots_",binsize,'.pdf',sep=''),
  useDingbats = F,
  width = 4,
  height = 2,
  pointsize=8
)
par(mfrow=c(1,7),mar=c(4,0,2,0),lend='square',ljoin='mitre',lwd=1,oma=c(0,4,0,4))
xmax = max(max(unlist(capped_chroms)),max(unlist(noncapped_chroms)))
xmax = 50
ymax = max(unlist(lapply(capped_chroms,length)))
lenmax = max(unlist(lapply(capped_chroms,length)))
for(chrom in names(chromosome_lengths)){
  barplot(capped_chroms[[chrom]],horiz = T,xlim=c(-xmax,xmax),ylim=c(0,ymax),space=0,border = NA,col=nanoPARE_colors['CAP'],axes = F)
  barplot(-noncapped_chroms[[chrom]],add=T,horiz = T,xlim=c(-125,125),space=0,border = NA,col=nanoPARE_colors['UNCAP'],axes = F)
  if(chrom%in%c('M','C')){
    segments(x0 = 0,x1 = 0,y0 = 0,y1 = chromosome_lengths[[chrom]]/smallbinsize)
  }else{
    segments(x0 = 0,x1 = 0,y0 = 0,y1 = chromosome_lengths[[chrom]]/binsize)
  }
  title(main = chrom)
  if(chrom == '1'){
    axis(1,at=c(-50,0,50))
    axis(2,at=seq(0,lenmax,by = 50),labels = round(seq(0,lenmax,by = 50)*binsize/10^6,0),las=1)
    abline(h=seq(0,lenmax,by = 50),col='gray')
  }
  if(chrom == '5'){
    abline(v=xmax)
  }
  if(chrom == 'M'){
    axis(1,at=c(-50,0,50))
    axis(4,at=seq(0,lenmax,by = 50),labels = round(seq(0,lenmax,by = 50)*smallbinsize/10^3,0),las=1)
    abline(h=seq(0,lenmax,by = 50),col='gray')
  }
}
dev.off()