def.par <- par(no.readonly = TRUE) # save default, for resetting...
xhist <- hist(d1, plot=FALSE)
yhist <- hist(d2, breaks=800, freq=F, plot=FALSE)
top <- max(c(xhist$counts, yhist$counts))
nf <- layout(matrix(c(2,0,1,3),2,2,byrow=TRUE), c(5,1), c(1,5), TRUE)
#layout.show(nf)
par(mar=c(4,4,1,1))
plot(d1, d2, xlab="d1 = log2(prediction_m1) - log2(observation)", ylab="d2 = log2(prediction_m2) - log2(observation)", col="#0000ff22", pch=16, cex=1.5, type='n')
abline(h=0,v=0, lty=2)
points(d1, d2, col="#0000ff22", pch=16, cex=1.5)
abline(a=DD,b=1, lty=2, col='red')
abline(a=-DD,b=1, lty=2, col='red')
library(calibrate)
textxy(d1[abs(d1-d2)>DD],d2[abs(d1-d2)>DD], as.character(m$gene_name[abs(d1-d2)>DD]), dcol='red')
a=round(lm(d2[abs(d1-d2)<DD]~d1[abs(d1-d2)<DD])$coefficients[1],3)
b=round(lm(d2[abs(d1-d2)<DD]~d1[abs(d1-d2)<DD])$coefficients[2],3)
abline(a=a,b=b, lty=1, col='red')
legend("topleft", c(paste("d2 = ",b,"*d1 ",ifelse(a>0,'+',''),a,sep=""),"m1: full model","m2: model without repressive marks"), bty='n')
par(mar=c(0,4,1,1))
hist(d1, breaks=800, freq=F, col='#0000ff22', border='#0000ff22', axes=FALSE, ylab="", xlab="",main="H1 cell line")
lines(density(d1, bw=1), col="black")
#barplot(xhist$counts, axes=FALSE, ylim=c(0, top), space=0, main="H1 cell line")
par(mar=c(4,0,1,1))
barplot(yhist$density, axes=FALSE, space=0, horiz=TRUE, col='#0000ff22', border='#0000ff22',)
#hist(d2, breaks=800, freq=F, col='#0000ff22', border='#0000ff22', axes=FALSE, ylab="", xlab="",main="")
lines(density(d2, bw=1)$y, y=seq(1,length(yhist$density), length.out=length(density(d2, bw=1)$y)),col="black")
par(def.par)
