bed = read.table('L:/members/Schon/Writing/nanoPARE/manuscript_draft/manuscript_data/supplemental_tables/ncRNA_features.tsv',header=T,sep='\t',stringsAsFactors = F)
bed = bed[bed[,'Transcript.type']%in%c('unit','tandem'),]
feature_type = bed[,"End.type"]

plot_colors = c(WT='#1C75BC',
                XRN1='#F7941E',
                XRN4='#00CC66')

comparison_df = data.frame(
  type1=c('tRNA_TSS','tRNA_TSS','snoRNA_TSS','miRNA_TSS','miRNA_TSS'),
  type2=c('tRNA_5P','tRNA_3P','snoRNA_5P','miRNA_3p5P','miRNA_3p3P'),
  xpos1=c(1,1,5,7,7),
  xpos2=c(2,3,6,8,9),
  ypos=c(8,10,8,8,10)
)

sigstars = function(pvalue,thresholds = c(10^-2,10^-3,10^-4)){
  stars = sapply(
    1:length(thresholds),
    function(x)paste(rep('*',x),collapse = '')
  )
  if(any(pvalue < thresholds)){
    return(stars[max(which(pvalue < thresholds))])
  }else{
    return("")
  }
}

pseudo = 0.01

# XRN1 DIGESTION #
nc_lfc = log2((rowMeans(bed[,12:14])+pseudo)/(rowMeans(bed[,9:11])+pseudo))

feature_list = list()
all_feature_types = sort(unique(feature_type),decreasing = T)

for(ft in all_feature_types){
  feature_list[[ft]] = as.numeric(nc_lfc[feature_type == ft])
}

library(beeswarm)
pdf('ncRNA_features_xrn1_fc.pdf',width = 8,height = 4,useDingbats = F,pointsize = 8)
beeswarm(feature_list,labels = gsub('_','\n',names(feature_list)),
         pch=20,spacing = .5,method = 'center',col='white',cex=.5,ylim=c(-12,12),
         main="Detected noncoding RNA 5' ends",xlab="5' end type",ylab='log2( Xrn1+ / Xrn1- ) RPM',las=1)

for(i in 1:length(feature_list)){
  dens_vals = density(feature_list[[i]],na.rm = T,kernel = 'gaussian',bw = 'SJ',n = 1000)
  polygon(x=c((-dens_vals$y)/2+i,(dens_vals$y[length(dens_vals$y):1])/2+i),y=c(dens_vals$x,dens_vals$x[length(dens_vals$x):1]),
          col=adjustcolor(plot_colors["XRN1"],alpha.f = .25),border=NA)
}
beeswarm(feature_list,pch=20,spacing = 1,method = 'center',col=plot_colors['XRN1'],cex=.25,add = T)
beeswarm(lapply(feature_list,function(x)mean(x,na.rm = T)),pch=15,col=plot_colors['XRN1'],cex=2,add = T)
beeswarm(lapply(feature_list,function(x)mean(x,na.rm = T)),pch=15,col='white',cex=1,add = T)
abline(h=0,v=c(3.5,4.5,6.5))

for(i in 1:length(feature_list)){
  p = t.test(feature_list[[i]])$p.value
  text(x=i,y=10,labels = sigstars(p),pos = 3,cex=2)
}

dev.off()






# XRN4 MUTANT #
nc_lfc = log2((rowMeans(bed[,15:17])+pseudo)/(rowMeans(bed[,9:11])+pseudo))
names(nc_lfc) = bed[,4]

feature_list = list()
all_feature_types = sort(unique(feature_type),decreasing = T)

for(ft in all_feature_types){
  feature_list[[ft]] = as.numeric(nc_lfc[feature_type == ft])
}

library(beeswarm)
pdf('ncRNA_features_xrn4_fc.pdf',width = 8,height = 4,useDingbats = F,pointsize = 8)
beeswarm(feature_list,labels = gsub('_','\n',names(feature_list)),
         pch=20,spacing = .5,method = 'center',col='white',cex=.5,ylim=c(-12,12),
         main="Detected noncoding RNA 5' ends",xlab="5' end type",ylab='log2( Xrn1+ / Xrn1- ) RPM',las=1)

for(i in 1:length(feature_list)){
  dens_vals = density(feature_list[[i]],na.rm = T,kernel = 'gaussian',bw = 'SJ',n = 1000)
  polygon(x=c((-dens_vals$y)/2+i,(dens_vals$y[length(dens_vals$y):1])/2+i),y=c(dens_vals$x,dens_vals$x[length(dens_vals$x):1]),
          col=adjustcolor(plot_colors["XRN4"],alpha.f = .25),border=NA)
}
beeswarm(feature_list,pch=20,spacing = 1,method = 'center',col=plot_colors['XRN4'],cex=.25,add = T)
beeswarm(lapply(feature_list,function(x)mean(x,na.rm = T)),pch=15,col=plot_colors['XRN4'],cex=2,add = T)
beeswarm(lapply(feature_list,function(x)mean(x,na.rm = T)),pch=15,col='white',cex=1,add = T)
abline(h=0,v=c(3.5,4.5,6.5))

for(i in 1:length(feature_list)){
  p = t.test(feature_list[[i]])$p.value
  text(x=i,y=10,labels = sigstars(p),pos = 3,cex=2)
}

dev.off()

