# Take the results from "all_features_deseq.R" and try to validate them
# using the public PARE data in wt vs. xrn4
nanoPARE_colors = c(
  BODY='#808285',
  CAP='#1C75BC',
  UNCAP="#F7941E",
  DIGEST='',
  XRN4='#00CC66',
  D234='#CC00FF'
)

german_cap = read.table('../data_tables/german.capped_exon.counts.tsv',stringsAsFactors = F,header = F)

german_noncap = read.table('../data_tables/german.noncapped_exon.counts.tsv',stringsAsFactors = F,header = F)


shared_agis = intersect(german_cap[,1],german_noncap[,1])

cap_plus = german_cap[german_cap[,3]=='+',4:5]
colnames(cap_plus) = c('wt','xrn4')
rownames(cap_plus) = german_cap[german_cap[,3]=='+',1]
cap_minus = german_cap[german_cap[,3]=='-',6:7]
colnames(cap_minus) = c('wt','xrn4')
rownames(cap_minus) = german_cap[german_cap[,3]=='-',1]
german_cap = rbind(cap_plus, cap_minus)
rm(cap_plus,cap_minus)

noncap_plus = german_noncap[german_noncap[,3]=='+',4:5]
colnames(noncap_plus) = c('wt','xrn4')
rownames(noncap_plus) = german_noncap[german_noncap[,3]=='+',1]
noncap_minus = german_noncap[german_noncap[,3]=='-',6:7]
colnames(noncap_minus) = c('wt','xrn4')
rownames(noncap_minus) = german_noncap[german_noncap[,3]=='-',1]
german_noncap = rbind(noncap_plus, noncap_minus)
rm(noncap_plus,noncap_minus)

percent_capped_wt = german_cap[shared_agis,'wt']/(german_cap[shared_agis,'wt']+german_noncap[shared_agis,'wt'])
names(percent_capped_wt) = shared_agis
percent_capped_xrn4 = german_cap[shared_agis,'xrn4']/(german_cap[shared_agis,'xrn4']+german_noncap[shared_agis,'xrn4'])
names(percent_capped_xrn4) = shared_agis

percent_capped_wt_clean = percent_capped_wt[!is.na(percent_capped_wt)]
percent_capped_xrn4_clean = percent_capped_xrn4[!is.na(percent_capped_xrn4)]

percent_change_german = percent_capped_xrn4_clean[shared_agis] - percent_capped_wt_clean[shared_agis]

denslist = list()
meanpoints = vector(length=2)
meanpoints[1] = mean(percent_capped_wt_clean)
meanpoints[2] = mean(percent_capped_xrn4_clean)
denslist[['wt']] = density(percent_capped_wt_clean,bw = .005,n = 500)
denslist[['wt']]$y = denslist[['wt']]$y/max(denslist[['wt']]$y)/2

denslist[['xrn4']] = density(percent_capped_xrn4_clean,bw = .005,n = 500)
denslist[['xrn4']]$y = denslist[['xrn4']]$y/max(denslist[['xrn4']]$y)/2


pdf('german_pare_comparison.pdf',useDingbats = F,width = 3,height = 3,pointsize = 8)
par(lend='square',ljoin='mitre',oma=c(2,0,0,0),mar=c(4,4,2,2))
plot(-10,-10,ylim=c(0,1),xlim=c(0,3),axes=F,xlab=NA,ylab='Percent full-length')
axis(1,c(1,2),labels=c('Col-0','xrn4-5'),tick = F,line = NA,outer = F)
axis(2,seq(0,1,.25),seq(0,100,25),las=1,line = NA,tick = F)
segments(x0 = c(0,0,0),x1 = c(3,3,3),y0 = c(.25,.50,.75),y1 = c(.25,.50,.75),col='gray')
# axis(1,1:9,rep(1:3,3),tick = F,line = NA)
colorvector = c(rep(nanoPARE_colors[["CAP"]],1),rep(nanoPARE_colors[["XRN4"]],1))
for(i in 1:length(denslist)){
  m = max(which(denslist[[i]]$x<=1))
  mn = min(which(denslist[[i]]$x>=0))
  polygon(
    x=c(i+denslist[[i]]$y[mn:m],i-denslist[[i]]$y[m:mn]),
    y=c(denslist[[i]]$x[mn:m],denslist[[i]]$x[m:mn]),
    col=adjustcolor(colorvector[i],alpha.f = .25),border=NA
  )
  points(meanpoints,pch=15,cex=1,col=colorvector)
  # points(meanpoints,pch=0,cex=1,col='black')
}
title(main="German et al. 2008 PARE data")
segments(x0 = c(0,0,0,3),x1 = c(0,3,3,3),y0 = c(0,0,1,0),y1 = c(1,0,1,1))
dev.off()

