tWT0h=read.table("/data/images/proton/DKlab/orsalia/perExon/wt01.cnt.rpmPerGene.csv",header=T)
xWT0hRNA=log2(0.5*(tWT0h$wt01.bam_rpm+tWT0h$wt02.bam_rpm));yWT0hPoly=log2(0.5*(tWT0h$WT_0h_replicate_I.bam_rpm+tWT0h$WT_0h_replicate_II_.bam_rpm));
xWT0hRNA[mapply(is.infinite, xWT0hRNA)] <- NA;yWT0hPoly[mapply(is.infinite, yWT0hPoly)] <- NA;

tKO0h=read.table("/data/images/proton/DKlab/orsalia/perExon/Ko01.cnt.rpmPerGene.csv",header=T)
xKO0hRNA=log2(0.5*(tKO0h$Ko01.bam_rpm+tKO0h$Ko02.bam_rpm));yKO0hPoly=log2(0.5*(tKO0h$KO_0h_replicate_I.bam_rpm+tKO0h$KO_0h_replicate_II.bam_rpm));
xKO0hRNA[mapply(is.infinite, xKO0hRNA)] <- NA;yKO0hPoly[mapply(is.infinite, yKO0hPoly)] <- NA;

bin<-hexbin(xWT0hRNA,yWT0hPoly, xbins=50,xbnds=c(-10,18),ybnds=c(-8,18));bin@count=log2(bin@count);
png("wt0h-PolysomeVSRNAseq-hist.png",width=1024,height=1024,pointsize = 22)
plot(bin, main="wt0h",xlab="log2_reads_per_exon_per_million_RNAseq",ylab="log2_reads_per_exon_per_million_PolysomesSeq")
dev.off()
bin<-hexbin(xKO0hRNA,yKO0hPoly, xbins=50,xbnds=c(-10,18),ybnds=c(-8,18));bin@count=log2(bin@count);
png("Ko0h-PolysomeVSRNAseq-hist.png",width=1024,height=1024,pointsize = 22)
plot(bin, main="ko0h",xlab="log2_reads_per_exon_per_million_RNAseq",ylab="log2_reads_per_exon_per_million_PolysomesSeq")
dev.off()


dd1=xWT0hRNA-yWT0hPoly;
dd2=xKO0hRNA-yKO0hPoly;

#centroids <- aggregate(cbind(dd1,dd2),mean)


png("0h-wtVSko-PolysomeVSRNAseq_perGene.png",width=2048,height=2048,pointsize = 22)
plot (dd1,dd2,pch=".",xlab="log2_reads_per_gene_per_million_RNA_minus_Polysomes_wt0h",ylab="delta_log2_reads_per_gene_per_million_RNA_minus_Polysomes_ko0h",xlim=c(-20,20),ylim=c(-20,20))
points(mean(na.omit(cbind(dd1))),mean(na.omit(cbind(dd2))),col="green")
abline(lm(dd1~dd2), col="red") # regression line (y~x)
abline(0,1, col="blue") # regression line (y~x)
grid(nx=NULL, ny=NULL)
set=as.logical(abs(dd1-dd2)>=5)
text(dd1[set],dd2[set], labels = tWT0h$Geneid[set], pos = 3,cex=0.5)
dev.off()



tWT2h=read.table("/data/images/proton/DKlab/orsalia/perExon/wt21.cnt.rpmPerGene.csv",header=T)
xWT2hRNA=log2(0.5*(tWT2h$wt21.bam_rpm+tWT2h$wt22.bam_rpm));yWT2hPoly=log2(0.5*(tWT2h$WT_2h_replicate_I.bam_rpm+tWT2h$WT_2h_replicate_II.bam_rpm));
xWT2hRNA[mapply(is.infinite, xWT2hRNA)] <- NA;yWT2hPoly[mapply(is.infinite, yWT2hPoly)] <- NA;

tKO2h=read.table("/data/images/proton/DKlab/orsalia/perExon/Ko21.cnt.rpmPerGene.csv",header=T)
xKO2hRNA=log2(0.5*(tKO2h$Ko21.bam_rpm+tKO2h$Ko22.bam_rpm));yKO2hPoly=log2(0.5*(tKO2h$KO_2h_replicate_I.bam_rpm+tKO2h$KO_2h_replicate_II.bam_rpm));
xKO2hRNA[mapply(is.infinite, xKO2hRNA)] <- NA;yKO2hPoly[mapply(is.infinite, yKO2hPoly)] <- NA;

bin<-hexbin(xWT2hRNA,yWT2hPoly, xbins=50,xbnds=c(-11,18),ybnds=c(-8,18));bin@count=log2(bin@count);
png("wt2h-PolysomeVSRNAseq-hist.png",width=1024,height=1024,pointsize = 22)
plot(bin, main="wt2h",xlab="log2_reads_per_exon_per_million_RNAseq",ylab="log2_reads_per_exon_per_million_PolysomesSeq")
dev.off()
bin<-hexbin(xKO2hRNA,yKO2hPoly, xbins=50,xbnds=c(-11,18),ybnds=c(-8,18));bin@count=log2(bin@count);
png("Ko2h-PolysomeVSRNAseq-hist.png",width=1024,height=1024,pointsize = 22)
plot(bin, main="ko2h",xlab="log2_reads_per_exon_per_million_RNAseq",ylab="log2_reads_per_exon_per_million_PolysomesSeq")
dev.off()

dd1=xWT2hRNA-yWT2hPoly;
dd2=xKO2hRNA-yKO2hPoly;

png("2h-wtVSko-PolysomeVSRNAseq_perGene.png",width=2048,height=2048,pointsize = 22)
plot (dd1,dd2,pch=".",xlab="log2_reads_per_gene_per_million_RNA_minus_Polysomes_wt2h",ylab="delta_log2_reads_per_gene_per_million_RNA_minus_Polysomes_ko2h",xlim=c(-20,20),ylim=c(-20,20))
points(mean(na.omit(cbind(dd1))),mean(na.omit(cbind(dd2))),col="green")
abline(lm(dd1~dd2), col="red") # regression line (y~x)
abline(0,1, col="blue") # regression line (y~x)
grid(nx=NULL, ny=NULL)
set=as.logical(abs(dd1-dd2)>=5)
text(dd1[set],dd2[set], labels = tWT2h$Geneid[set], pos = 3,cex=0.5)
dev.off()
tWT6h=read.table("/data/images/proton/DKlab/orsalia/perExon/wt61.cnt.rpmPerGene.csv",header=T)
xWT6hRNA=log2(0.5*(tWT6h$wt61.bam_rpm+tWT6h$wt62.bam_rpm));yWT6hPoly=log2(0.5*(tWT6h$WT_6h_replicate_I.bam_rpm+tWT6h$WT_6h_replicate_II.bam_rpm));
xWT6hRNA[mapply(is.infinite, xWT6hRNA)] <- NA;yWT6hPoly[mapply(is.infinite, yWT6hPoly)] <- NA;

tKO6h=read.table("/data/images/proton/DKlab/orsalia/perExon/Ko61.cnt.rpmPerGene.csv",header=T)
xKO6hRNA=log2(0.5*(tKO6h$Ko61.bam_rpm+tKO6h$Ko62.bam_rpm));yKO6hPoly=log2(0.5*(tKO6h$KO_6h_replicate_I.bam_rpm+tKO6h$KO_6h_replicate_II.bam_rpm));
xKO6hRNA[mapply(is.infinite, xKO6hRNA)] <- NA;yKO6hPoly[mapply(is.infinite, yKO6hPoly)] <- NA;

bin<-hexbin(xWT6hRNA,yWT6hPoly, xbins=50,xbnds=c(-11,18),ybnds=c(-8,18));bin@count=log2(bin@count);
png("wt6h-PolysomeVSRNAseq-hist.png",width=1024,height=1024,pointsize = 22)
plot(bin, main="wt6h",xlab="log2_reads_per_exon_per_million_RNAseq",ylab="log2_reads_per_exon_per_million_PolysomesSeq")
dev.off()
bin<-hexbin(xKO6hRNA,yKO6hPoly, xbins=50,xbnds=c(-11,18),ybnds=c(-8,18));bin@count=log2(bin@count);
png("Ko6h-PolysomeVSRNAseq-hist.png",width=1024,height=1024,pointsize = 22)
plot(bin, main="ko6h",xlab="log2_reads_per_exon_per_million_RNAseq",ylab="log2_reads_per_exon_per_million_PolysomesSeq")
dev.off()

dd1=xWT6hRNA-yWT6hPoly;
dd2=xKO6hRNA-yKO6hPoly;

png("6h-wtVSko-PolysomeVSRNAseq_perGene.png",width=2048,height=2048,pointsize = 22)
plot (dd1,dd2,pch=".",xlab="log2_reads_per_gene_per_million_RNA_minus_Polysomes_wt6h",ylab="delta_log2_reads_per_gene_per_million_RNA_minus_Polysomes_ko6h",xlim=c(-20,20),ylim=c(-20,20))
points(mean(na.omit(cbind(dd1))),mean(na.omit(cbind(dd2))),col="green")
abline(lm(dd1~dd2), col="red") # regression line (y~x)
abline(0,1, col="blue") # regression line (y~x)
grid(nx=NULL, ny=NULL)
set=as.logical(abs(dd1-dd2)>=5)
text(dd1[set],dd2[set], labels = tWT6h$Geneid[set], pos = 3,cex=0.5)
dev.off()


tWTIFN=read.table("/data/images/proton/DKlab/orsalia/perExon/wtIFN1.cnt.rpmPerGene.csv",header=T)
xWTIFNRNA=log2(0.5*(tWTIFN$wtIFN1.bam_rpm+tWTIFN$wtIFN2.bam_rpm));yWTIFNPoly=log2(0.5*(tWTIFN$WT_IFN.gamma_replicate_I.bam_rpm+tWTIFN$WT_IFN.gamma_replicate_II.bam_rpm));
xWTIFNRNA[mapply(is.infinite, xWTIFNRNA)] <- NA;yWTIFNPoly[mapply(is.infinite, yWTIFNPoly)] <- NA;

tKOIFN=read.table("/data/images/proton/DKlab/orsalia/perExon/KOIFN1.cnt.rpmPerGene.csv",header=T)
xKOIFNRNA=log2(0.5*(tKOIFN$KoIFN1.bam_rpm+tKOIFN$KoIFN2.bam_rpm));yKOIFNPoly=log2(0.5*(tKOIFN$KO_IFN.gamma_replicate_I.bam_rpm+tKOIFN$KO_IFN.gamma_replicate_II.bam_rpm));
xKOIFNRNA[mapply(is.infinite, xKOIFNRNA)] <- NA;yKOIFNPoly[mapply(is.infinite, yKOIFNPoly)] <- NA;

bin<-hexbin(xWTIFNRNA,yWTIFNPoly, xbins=50,xbnds=c(-11,18),ybnds=c(-8,18));bin@count=log2(bin@count);
png("wtIFN-PolysomeVSRNAseq-hist.png",width=1024,height=1024,pointsize = 22)
plot(bin, main="wtIFN",xlab="log2_reads_per_exon_per_million_RNAseq",ylab="log2_reads_per_exon_per_million_PolysomesSeq")
dev.off()
bin<-hexbin(xKOIFNRNA,yKOIFNPoly, xbins=50,xbnds=c(-11,18),ybnds=c(-8,18));bin@count=log2(bin@count);
png("KoIFN-PolysomeVSRNAseq-hist.png",width=1024,height=1024,pointsize = 22)
plot(bin, main="koIFN",xlab="log2_reads_per_exon_per_million_RNAseq",ylab="log2_reads_per_exon_per_million_PolysomesSeq")
dev.off()

dd1=xWTIFNRNA-yWTIFNPoly;
dd2=xKOIFNRNA-yKOIFNPoly;

png("IFN-wtVSko-PolysomeVSRNAseq_perGene.png",width=2048,height=2048,pointsize = 22)
plot (dd1,dd2,pch=".",xlab="log2_reads_per_gene_per_million_RNA_minus_Polysomes_wtIFN",ylab="delta_log2_reads_per_gene_per_million_RNA_minus_Polysomes_koIFN",xlim=c(-20,20),ylim=c(-20,20))
points(mean(na.omit(cbind(dd1))),mean(na.omit(cbind(dd2))),col="green")
abline(lm(dd1~dd2), col="red") # regression line (y~x)
abline(0,1, col="blue") # regression line (y~x)
grid(nx=NULL, ny=NULL)
set=as.logical(abs(dd1-dd2)>=5)
text(dd1[set],dd2[set], labels = tWTIFN$Geneid[set], pos = 3,cex=0.5)
dev.off()

tWTIL4=read.table("/data/images/proton/DKlab/orsalia/perExon/wt01IL4.cnt.rpmPerGene.csv",header=T)
xWTIL4RNA=log2(0.5*(tWTIL4$wt01IL4.bam_rpm+tWTIL4$wt02IL4.bam_rpm));yWTIL4Poly=log2(0.5*(tWTIL4$WT_IL4_replicate_I.bam_rpm+tWTIL4$WT_IL4_replicate_II.bam_rpm));
xWTIL4RNA[mapply(is.infinite, xWTIL4RNA)] <- NA;yWTIL4Poly[mapply(is.infinite, yWTIL4Poly)] <- NA;

tKOIL4=read.table("/data/images/proton/DKlab/orsalia/perExon/Ko1IL4.cnt.rpmPerGene.csv",header=T)
xKOIL4RNA=log2(0.5*(tKOIL4$Ko1IL4.bam_rpm+tKOIL4$Ko2IL4.bam_rpm));yKOIL4Poly=log2(0.5*(tKOIL4$KO_IL4_replicate_I.bam_rpm+tKOIL4$KO_IL4_replicate_II.bam_rpm));
xKOIL4RNA[mapply(is.infinite, xKOIL4RNA)] <- NA;yKOIL4Poly[mapply(is.infinite, yKOIL4Poly)] <- NA;

bin<-hexbin(xWTIL4RNA,yWTIL4Poly, xbins=50,xbnds=c(-11,18),ybnds=c(-8,18));bin@count=log2(bin@count);
png("wtIL4-PolysomeVSRNAseq-hist.png",width=1024,height=1024,pointsize = 22)
plot(bin, main="wtIL4",xlab="log2_reads_per_exon_per_million_RNAseq",ylab="log2_reads_per_exon_per_million_PolysomesSeq")
dev.off()
bin<-hexbin(xKOIL4RNA,yKOIL4Poly, xbins=50,xbnds=c(-11,18),ybnds=c(-8,18));bin@count=log2(bin@count);
png("KoIL4-PolysomeVSRNAseq-hist.png",width=1024,height=1024,pointsize = 22)
plot(bin, main="koIL4",xlab="log2_reads_per_exon_per_million_RNAseq",ylab="log2_reads_per_exon_per_million_PolysomesSeq")
dev.off()

dd1=xWTIL4RNA-yWTIL4Poly;
dd2=xKOIL4RNA-yKOIL4Poly;

png("IL4-wtVSko-PolysomeVSRNAseq_perGene.png",width=2048,height=2048,pointsize = 22)
plot (dd1,dd2,pch=".",xlab="log2_reads_per_gene_per_million_RNA_minus_Polysomes_wtIL4",ylab="delta_log2_reads_per_gene_per_million_RNA_minus_Polysomes_koIL4",xlim=c(-20,20),ylim=c(-20,20))
points(mean(na.omit(cbind(dd1))),mean(na.omit(cbind(dd2))),col="green")
abline(lm(dd1~dd2), col="red") # regression line (y~x)
abline(0,1, col="blue") # regression line (y~x)
grid(nx=NULL, ny=NULL)
set=as.logical(abs(dd1-dd2)>=5)
text(dd1[set],dd2[set], labels = tWTIL4$Geneid[set], pos = 3,cex=0.5)
dev.off()

dd1=xWT0hRNA-yWT0hPoly;
dd2=xWT2hRNA-yWT2hPoly;
png("WT-0hVS2h-PolysomeVSRNAseq_perGene.png",width=2048,height=2048,pointsize = 22)
plot (dd1,dd2,pch=".",xlab="log2_reads_per_gene_per_million_RNA_minus_Polysomes_wt0h",ylab="delta_log2_reads_per_gene_per_million_RNA_minus_Polysomes_wt2h",xlim=c(-20,20),ylim=c(-20,20))
points(mean(na.omit(cbind(dd1))),mean(na.omit(cbind(dd2))),col="green")
abline(lm(dd1~dd2), col="red") # regression line (y~x)
abline(0,1, col="blue") # regression line (y~x)
grid(nx=NULL, ny=NULL)
set=as.logical(abs(dd1-dd2)>=5)
text(dd1[set],dd2[set], labels = tWT0h$Geneid[set], pos = 3,cex=0.5)
dev.off()

dd1=xKO0hRNA-yKO0hPoly;
dd2=xKO2hRNA-yKO2hPoly;
png("KO-0hVS2h-PolysomeVSRNAseq_perGene.png",width=2048,height=2048,pointsize = 22)
plot (dd1,dd2,pch=".",xlab="log2_reads_per_gene_per_million_RNA_minus_Polysomes_ko0h",ylab="delta_log2_reads_per_gene_per_million_RNA_minus_Polysomes_ko2h",xlim=c(-20,20),ylim=c(-20,20))
points(mean(na.omit(cbind(dd1))),mean(na.omit(cbind(dd2))),col="green")
abline(lm(dd1~dd2), col="red") # regression line (y~x)
abline(0,1, col="blue") # regression line (y~x)
grid(nx=NULL, ny=NULL)
set=as.logical(abs(dd1-dd2)>=5)
text(dd1[set],dd2[set], labels = tKO0h$Geneid[set], pos = 3,cex=0.5)
dev.off()

dd1=xWT0hRNA-yWT0hPoly;
dd2=xWT6hRNA-yWT6hPoly;
png("WT-0hVS6h-PolysomeVSRNAseq_perGene.png",width=2048,height=2048,pointsize = 22)
plot (dd1,dd2,pch=".",xlab="log2_reads_per_gene_per_million_RNA_minus_Polysomes_wt0h",ylab="delta_log2_reads_per_gene_per_million_RNA_minus_Polysomes_wt6h",xlim=c(-20,20),ylim=c(-20,20))
points(mean(na.omit(cbind(dd1))),mean(na.omit(cbind(dd2))),col="green")
abline(lm(dd1~dd2), col="red") # regression line (y~x)
abline(0,1, col="blue") # regression line (y~x)
grid(nx=NULL, ny=NULL)
set=as.logical(abs(dd1-dd2)>=5)
text(dd1[set],dd2[set], labels = tWT0h$Geneid[set], pos = 3,cex=0.5)
dev.off()

dd1=xKO0hRNA-yKO0hPoly;
dd2=xKO6hRNA-yKO6hPoly;
png("KO-0hVS6h-PolysomeVSRNAseq_perGene.png",width=2048,height=2048,pointsize = 22)
plot (dd1,dd2,pch=".",xlab="log2_reads_per_gene_per_million_RNA_minus_Polysomes_ko0h",ylab="delta_log2_reads_per_gene_per_million_RNA_minus_Polysomes_ko6h",xlim=c(-20,20),ylim=c(-20,20))
points(mean(na.omit(cbind(dd1))),mean(na.omit(cbind(dd2))),col="green")
abline(lm(dd1~dd2), col="red") # regression line (y~x)
abline(0,1, col="blue") # regression line (y~x)
grid(nx=NULL, ny=NULL)
set=as.logical(abs(dd1-dd2)>=5)
text(dd1[set],dd2[set], labels = tKO0h$Geneid[set], pos = 3,cex=0.5)
dev.off()