library(cobs)
library(quantreg)
library(parallel)
library(gplots)

interpolate<-function(dat,curve){

interpolate_points<-function(row,dat,curve){
    MY_X=dat[row,1]
    MY_Y=dat[row,2]
    VAL1=tail(which(curve[,1]<MY_X),1)
    VAL2=VAL1+1
    X <- curve[c(VAL1,VAL2),1]
    Y <- curve[c(VAL1,VAL2),2]
    INTERP_Y=approx(X,Y,xout=MY_X)$y
    INTERP_Y
}
res<-unlist(mclapply(seq(nrow(dat)),interpolate_points,dat=dat,curve=curve,mc.cores=8))
res
}

################################################
# 1Mbp bins
################################################
x<-read.table("SRP199233.1e6_fmt.tsv",header=T,row.names=1)
x<-x[which(rowSums(x)>=10),]
x<-sweep(x, 2, colSums(x), FUN="/")*1000000
mysd<-apply(x,1,sd)
mean<-apply(x,1,mean)
y<-data.frame(log10(mean),mysd/mean)
colnames(y)=c("logMean","cv")
Rbs.9 <- cobs(y$logMean,y$cv, nknots=20,constraint="none",tau=0.99)
Rbs.median <- cobs(y$logMean,y$cv,nknots=20,constraint="none",tau=0.5)
pred<-data.frame(predict(Rbs.9))

y$interpolated<-interpolate(y,pred)
y$diff=y$cv-y$interpolated
yy <- y[which(y$diff>0.05),]
yy <- yy[order(-yy$diff),]

write.table(yy,file="SRP199233.1e6_regions.tsv")

png("plots/SRP199233.1e6_cv.png",height=960,width=960)
plot(y$logMean,y$cv,pch=18,cex=0.5,xlab="log10(mean)",ylab="CV")
lines(predict(Rbs.9), col = "red", lwd = 1.0)
lines(predict(Rbs.median), col = "blue", lwd = 1.0)
points(yy$logMean,yy$cv)
text(yy$logMean,yy$cv+0.02,labels=rownames(yy),cex=0.8)
dev.off()

png("plots/SRP199233.1e6_hm1.png",height=960,width=960)
my_palette <- colorRampPalette(c("blue", "white", "red"))(n = 25)
zz<-x[which(rownames(x) %in% rownames(yy)),]
heatmap.2(as.matrix(zz),margin=c(8, 22),cexRow=0.65,trace="none",
  cexCol=0.8,col=my_palette,scale="row")
dev.off()

png("plots/SRP199233.1e6_hm2.png",height=960,width=960)
heatmap.2(cor(t(zz)),trace="none",scale="none",margins=c(12,12),
  cexRow=0.8, cexCol=0.8)
dev.off()