#read gene biotype annotation
setwd("~/Documents/HuR_PARclip_withHeaders/reference_gtf/")
library(GenomicRanges)
library(rtracklayer)
GTFfile = "Mus_musculus.NCBIM37.64-toMM9.gtf"
GTF <- import.gff(GTFfile, format="gtf")
gna=elementMetadata(GTF)$gene_name
gbt=elementMetadata(GTF)$gene_biotype
z=cbind(as.data.frame(gna),as.data.frame(gbt))
z = z[!duplicated(z$gna),]
colnames(z)=c("GeneName", "Biotype")

##### IL4  rpkm- Keep only RNAseq protein Coding RPKM>=1 
setwd("~/Documents/norm_Polys_to_RNAseq_DEGs/")
trl=read.table("Mus_musculus.NCBIM37.64-toMM9.gene-lenghts.txt")
colnames(trl)= c("GeneID", "GeneName", "Length")
trl=trl[!duplicated(trl["GeneName"]),]

setwd("~/giagkas/RNAseq_POLYseq/RNAseq/DEG/")
rn=read.table("edgeRgenes_featurecountExpr-IL4.csv", header = T)
rn$GeneName=rownames(rn)
rn=subset(rn, select=c(9,1:8))
rownames(rn)=c()
colnames(rn)=c("GeneName", "logFC", "logCPM", "PValue", "FDR", "KO1", "KO2", "WT1", "WT2")
setwd("~/giagkas/RNAseq_POLYseq/Polysomes/DEG/")
po=read.table("edgeRgenes_featurecountExpr-IL4.csv", header = T)
colnames(po)=c("GeneName", "logFC", "logCPM", "PValue", "FDR", "KO1", "KO2", "WT1", "WT2")
setwd("~/giagkas/RNAseq_POLYseq/")
m= merge(rn, po, by.x = "GeneName", by.y = "GeneName", all = T )

m=merge(m, trl, by= "GeneName", all.x = T)

#m[2:5]=NULL
#m[6:9]=NULL
#m[is.na(m)] <- 0

# RPKM=CPM/(transcript_length/1000)
rpkm=m
rpkm$KO1.x=m$KO1.x/(m$Length/1000)
rpkm$KO2.x=m$KO2.x/(m$Length/1000)
rpkm$WT1.x=m$WT1.x/(m$Length/1000)
rpkm$WT2.x=m$WT2.x/(m$Length/1000)
rpkm$KO1.y=m$KO1.y/(m$Length/1000)
rpkm$KO2.y=m$KO2.y/(m$Length/1000)
rpkm$WT1.y=m$WT1.y/(m$Length/1000)
rpkm$WT2.y=m$WT2.y/(m$Length/1000)

m=rpkm
m$KOrna=(m$KO1.x+m$KO2.x)/2
m$WTrna=(m$WT1.x+m$WT2.x)/2
m$KOpoly=(m$KO1.y+m$KO2.y)/2
m$WTpoly=(m$WT1.y+m$WT2.y)/2

am=m
#add gene biotype annotation
am=merge(am,z, by="GeneName", all.x=T)
colnames(am)=c("GeneName", "logFC.RNA",  "logCPM.RNA", "PValue.RNA", "FDR.RNA",    "KO1CPM.RNA",
               "KO2CPM.RNA",    "WT1CPM.RNA",    "WT2CPM.RNA",    "logFC.POLY",  "logCPM.POLY", "PValue.POLY",
               "FDR.POLY",    "KO1CPM.POLY",    "KO2CPM.POLY",    "WT1CPM.POLY",    "WT2CPM.POLY",    "GeneID","Length",
               "koRPKM.RNA", "wtRPKM.RNA", "koRPKM.POLY", "wtRPKM.POLY", "Biotype")

write.table(am, file = "RNA-Poly_IL4total.txt", sep = "\t", quote = F, row.names = F)


##protein coding
m=merge(m,z, by="GeneName", all.x=T)
m=subset(m, (m$Biotype=="protein_coding"))
pm=m
colnames(pm)=c("GeneName", "logFC.RNA",  "logCPM.RNA", "PValue.RNA", "FDR.RNA",    "KO1CPM.RNA",
               "KO2CPM.RNA",    "WT1CPM.RNA",    "WT2CPM.RNA",    "logFC.POLY",  "logCPM.POLY", "PValue.POLY",
               "FDR.POLY",    "KO1CPM.POLY",    "KO2CPM.POLY",    "WT1CPM.POLY",    "WT2CPM.POLY",    "GeneID","Length",
               "koRPKM.RNA", "wtRPKM.RNA", "koRPKM.POLY", "wtRPKM.POLY", "Biotype")


setwd("~/giagkas/RNAseq_POLYseq/")
write.table(pm, file = "RNA-Poly_IL4pc.txt", sep = "\t", quote = F, row.names = F)

######
setwd("~/giagkas/RNAseq_POLYseq/plots/")
rpkm=m
rpkm$logKOrna= log2(rpkm$KOrna)
rpkm$logWTrna= log2(rpkm$WTrna)
rpkm$logKOpoly= log2(rpkm$KOpoly)
rpkm$logWTpoly= log2(rpkm$WTpoly)

m2=rpkm
fcru=subset(m2, (m2$WTrna>=1)&(m2$PValue.x<0.05)&(m2$logFC.x>0.584963))
fcrd=subset(m2, (m2$WTrna>=1)&(m2$PValue.x<0.05)&(m2$logFC.x<(-0.584963)))
fcpu=subset(m2, (m2$WTrna>=1)&(m2$WTpoly>=1)&(m2$PValue.y<0.05)&(m2$logFC.y>0.584963))
fcpd=subset(m2, (m2$WTrna>=1)&(m2$WTpoly>=1)&(m2$PValue.y<0.05)&(m2$logFC.y<(-0.584963)))


png(filename = "Abundance_Translation_IL4_PC_wtThreshold.png", width =1200, height = 600)
par(mfrow=c(1,2),font.axis=2, cex.axis=1.5, cex.lab=1.5, lwd=2)
plot(m2$logWTrna,m2$logKOrna, pch=".", ylim = c(0,15),xlim = c(0,15), main= "Abundance M2", xlab = "log2RPKM WT",
     ylab = "log2RPKM KO")
points(fcru$logWTrna, fcru$logKOrna, pch="*", col="red")
points(fcrd$logWTrna, fcrd$logKOrna, pch="*", col="green")
legend("topleft", c(paste("FC>1.5","(n=",nrow(fcru),")", sep = " "),
                    paste("FC<-1.5","(n=",nrow(fcrd),")", sep = " ")), bty = "n", text.font=2, cex=1.5)
plot(m2$logWTpoly,m2$logKOpoly, pch=".", ylim = c(0,15), xlim = c(0,15),main= "Translation M2", xlab = "log2RPKM WT",
     ylab = "log2RPKM KO")
points(fcpu$logWTpoly, fcpu$logKOpoly, pch="*", col="red")
points(fcpd$logWTpoly, fcpd$logKOpoly, pch="*", col="green")
legend("topleft", c(paste("FC>1.5","(n=",nrow(fcpu),")", sep = " "),
                    paste("FC<-1.5","(n=",nrow(fcpd),")", sep = " ")), bty = "n", text.font=2, cex=1.5)

dev.off()




setwd("~/giagkas/RNAseq_POLYseq/plots/")
#m=subset(m, (m$WTrna>=1)&(m$WTpoly>=1))
m=subset(m, (m$WTrna>=1)&(m$WTpoly>=1)&(m$KOrna>=1)&(m$KOpoly>=1))
rpkm=m

rpkm$logKOrna= log2(rpkm$KOrna)
rpkm$logWTrna= log2(rpkm$WTrna)
rpkm$logKOpoly= log2(rpkm$KOpoly)
rpkm$logWTpoly= log2(rpkm$WTpoly)
m2=rpkm

fcru=subset(m2, (m2$PValue.x<0.05)&m2$logFC.x>0.584963)
fcrd=subset(m2, (m2$PValue.x<0.05)&m2$logFC.x<(-0.584963))
fcpu=subset(m2, (m2$PValue.y<0.05)&m2$logFC.y>0.584963)
fcpd=subset(m2, (m2$PValue.y<0.05)&m2$logFC.y<(-0.584963))

png(filename = "Abundance_Translation_IL4_PC.png", width =1200, height = 600)
par(mfrow=c(1,2),font.axis=2, cex.axis=1.5, cex.lab=1.5, lwd=2)
plot(m2$logWTrna,m2$logKOrna, pch=".", ylim = c(0,15),xlim = c(0,15), main= "Abundance M2", xlab = "log2RPKM WT",
     ylab = "log2RPKM KO")
points(fcru$logWTrna, fcru$logKOrna, pch="*", col="red")
points(fcrd$logWTrna, fcrd$logKOrna, pch="*", col="green")
legend("topleft", c(paste("FC>1.5","(n=",nrow(fcru),")", sep = " "),
                    paste("FC<-1.5","(n=",nrow(fcrd),")", sep = " ")), bty = "n", text.font=2, cex=1.5)
plot(m2$logWTpoly,m2$logKOpoly, pch=".", ylim = c(0,15), xlim = c(0,15),main= "Translation M2", xlab = "log2RPKM WT",
     ylab = "log2RPKM KO")
points(fcpu$logWTpoly, fcpu$logKOpoly, pch="*", col="red")
points(fcpd$logWTpoly, fcpd$logKOpoly, pch="*", col="green")
legend("topleft", c(paste("FC>1.5","(n=",nrow(fcpu),")", sep = " "),
                    paste("FC<-1.5","(n=",nrow(fcpd),")", sep = " ")), bty = "n", text.font=2, cex=1.5)

dev.off()



##### IFN  rpkm - Keep only RNAseq protein Coding RPKM>=1 
setwd("~/Documents/norm_Polys_to_RNAseq_DEGs/")
trl=read.table("Mus_musculus.NCBIM37.64-toMM9.gene-lenghts.txt")
colnames(trl)= c("GeneID", "GeneName", "Length")
trl=trl[!duplicated(trl["GeneName"]),]

setwd("~/giagkas/RNAseq_POLYseq/RNAseq/DEG/")
rn=read.table("edgeRgenes_featurecountExpr-IFN.csv", header = T)
rn$GeneName=rownames(rn)
rn=subset(rn, select=c(9,1:8))
rownames(rn)=c()
colnames(rn)=c("GeneName", "logFC", "logCPM", "PValue", "FDR", "KO1", "KO2", "WT1", "WT2")
setwd("~/giagkas/RNAseq_POLYseq/Polysomes/DEG/")
po=read.table("edgeRgenes_featurecountExpr-IFN.csv", header = T)
colnames(po)=c("GeneName", "logFC", "logCPM", "PValue", "FDR", "KO1", "KO2", "WT1", "WT2")
setwd("~/giagkas/RNAseq_POLYseq/")
m= merge(rn, po, by.x = "GeneName", by.y = "GeneName", all = T )

m=merge(m, trl, by= "GeneName", all.x = T)

#m[2:5]=NULL
#m[6:9]=NULL
#m[is.na(m)] <- 0

# RPKM=CPM/(transcript_length/1000)
rpkm=m
rpkm$KO1.x=m$KO1.x/(m$Length/1000)
rpkm$KO2.x=m$KO2.x/(m$Length/1000)
rpkm$WT1.x=m$WT1.x/(m$Length/1000)
rpkm$WT2.x=m$WT2.x/(m$Length/1000)
rpkm$KO1.y=m$KO1.y/(m$Length/1000)
rpkm$KO2.y=m$KO2.y/(m$Length/1000)
rpkm$WT1.y=m$WT1.y/(m$Length/1000)
rpkm$WT2.y=m$WT2.y/(m$Length/1000)

m=rpkm
m$KOrna=(m$KO1.x+m$KO2.x)/2
m$WTrna=(m$WT1.x+m$WT2.x)/2
m$KOpoly=(m$KO1.y+m$KO2.y)/2
m$WTpoly=(m$WT1.y+m$WT2.y)/2

am=m
#add gene biotype annotation
am=merge(am,z, by="GeneName", all.x=T)
colnames(am)=c("GeneName", "logFC.RNA",  "logCPM.RNA", "PValue.RNA", "FDR.RNA",    "KO1CPM.RNA",
               "KO2CPM.RNA",    "WT1CPM.RNA",    "WT2CPM.RNA",    "logFC.POLY",  "logCPM.POLY", "PValue.POLY",
               "FDR.POLY",    "KO1CPM.POLY",    "KO2CPM.POLY",    "WT1CPM.POLY",    "WT2CPM.POLY",    "GeneID","Length",
               "koRPKM.RNA", "wtRPKM.RNA", "koRPKM.POLY", "wtRPKM.POLY", "Biotype")

write.table(am, file = "RNA-Poly_IFNtotal.txt", sep = "\t", quote = F, row.names = F)


##protein coding
m=merge(m,z, by="GeneName", all.x=T)
m=subset(m, (m$Biotype=="protein_coding"))
pm=m
colnames(pm)=c("GeneName", "logFC.RNA",  "logCPM.RNA", "PValue.RNA", "FDR.RNA",    "KO1CPM.RNA",
               "KO2CPM.RNA",    "WT1CPM.RNA",    "WT2CPM.RNA",    "logFC.POLY",  "logCPM.POLY", "PValue.POLY",
               "FDR.POLY",    "KO1CPM.POLY",    "KO2CPM.POLY",    "WT1CPM.POLY",    "WT2CPM.POLY",    "GeneID","Length",
               "koRPKM.RNA", "wtRPKM.RNA", "koRPKM.POLY", "wtRPKM.POLY", "Biotype")


setwd("~/giagkas/RNAseq_POLYseq/")
write.table(pm, file = "RNA-Poly_IFNpc.txt", sep = "\t", quote = F, row.names = F)

#########
setwd("~/giagkas/RNAseq_POLYseq/plots/")
rpkm=m

rpkm$logKOrna= log2(rpkm$KOrna)
rpkm$logWTrna= log2(rpkm$WTrna)
rpkm$logKOpoly= log2(rpkm$KOpoly)
rpkm$logWTpoly= log2(rpkm$WTpoly)
m2=rpkm

fcru=subset(m2, (m2$WTrna>=1)&(m2$PValue.x<0.05)&(m2$logFC.x>0.584963))
fcrd=subset(m2, (m2$WTrna>=1)&(m2$PValue.x<0.05)&(m2$logFC.x<(-0.584963)))
fcpu=subset(m2, (m2$WTrna>=1)&(m2$WTpoly>=1)&(m2$PValue.y<0.05)&(m2$logFC.y>0.584963))
fcpd=subset(m2, (m2$WTrna>=1)&(m2$WTpoly>=1)&(m2$PValue.y<0.05)&(m2$logFC.y<(-0.584963)))


png(filename = "Abundance_Translation_IFN_PC_wtThreshold.png", width =1200, height = 600)
par(mfrow=c(1,2),font.axis=2, cex.axis=1.5, cex.lab=1.5, lwd=2)
plot(m2$logWTrna,m2$logKOrna, pch=".", ylim = c(0,15),xlim = c(0,15), main= "Abundance M1-IFN", xlab = "log2RPKM WT",
     ylab = "log2RPKM KO")
points(fcru$logWTrna, fcru$logKOrna, pch="*", col="red")
points(fcrd$logWTrna, fcrd$logKOrna, pch="*", col="green")
legend("topleft", c(paste("FC>1.5","(n=",nrow(fcru),")", sep = " "),
                    paste("FC<-1.5","(n=",nrow(fcrd),")", sep = " ")), bty = "n", text.font=2, cex=1.5)
plot(m2$logWTpoly,m2$logKOpoly, pch=".", ylim = c(0,15), xlim = c(0,15),main= "Translation M1-IFN", xlab = "log2RPKM WT",
     ylab = "log2RPKM KO")
points(fcpu$logWTpoly, fcpu$logKOpoly, pch="*", col="red")
points(fcpd$logWTpoly, fcpd$logKOpoly, pch="*", col="green")
legend("topleft", c(paste("FC>1.5","(n=",nrow(fcpu),")", sep = " "),
                    paste("FC<-1.5","(n=",nrow(fcpd),")", sep = " ")), bty = "n", text.font=2, cex=1.5)
dev.off()


##

setwd("~/giagkas/RNAseq_POLYseq/plots/")
#m=subset(m, (m$WTrna>=1)&(m$WTpoly>=1))
m=subset(m, (m$WTrna>=1)&(m$WTpoly>=1)&(m$KOrna>=1)&(m$KOpoly>=1))
rpkm=m

rpkm$logKOrna= log2(rpkm$KOrna)
rpkm$logWTrna= log2(rpkm$WTrna)
rpkm$logKOpoly= log2(rpkm$KOpoly)
rpkm$logWTpoly= log2(rpkm$WTpoly)
m2=rpkm

fcru=subset(m2, (m2$PValue.x<0.05)&m2$logFC.x>0.584963)
fcrd=subset(m2, (m2$PValue.x<0.05)&m2$logFC.x<(-0.584963))
fcpu=subset(m2, (m2$PValue.y<0.05)&m2$logFC.y>0.584963)
fcpd=subset(m2, (m2$PValue.y<0.05)&m2$logFC.y<(-0.584963))

png(filename = "Abundance_Translation_IFN_PC.png", width =1200, height = 600)
par(mfrow=c(1,2),font.axis=2, cex.axis=1.5, cex.lab=1.5, lwd=2)
plot(m2$logWTrna,m2$logKOrna, pch=".", ylim = c(0,15),xlim = c(0,15), main= "Abundance M1-IFN", xlab = "log2RPKM WT",
     ylab = "log2RPKM KO")
points(fcru$logWTrna, fcru$logKOrna, pch="*", col="red")
points(fcrd$logWTrna, fcrd$logKOrna, pch="*", col="green")
legend("topleft", c(paste("FC>1.5","(n=",nrow(fcru),")", sep = " "),
                    paste("FC<-1.5","(n=",nrow(fcrd),")", sep = " ")), bty = "n", text.font=2, cex=1.5)
plot(m2$logWTpoly,m2$logKOpoly, pch=".", ylim = c(0,15), xlim = c(0,15),main= "Translation M1-IFN", xlab = "log2RPKM WT",
     ylab = "log2RPKM KO")
points(fcpu$logWTpoly, fcpu$logKOpoly, pch="*", col="red")
points(fcpd$logWTpoly, fcpd$logKOpoly, pch="*", col="green")
legend("topleft", c(paste("FC>1.5","(n=",nrow(fcpu),")", sep = " "),
                    paste("FC<-1.5","(n=",nrow(fcpd),")", sep = " ")), bty = "n", text.font=2, cex=1.5)
dev.off()




##### M0  rpkm - Keep only RNAseq protein Coding RPKM>=1 
setwd("~/Documents/norm_Polys_to_RNAseq_DEGs/")
trl=read.table("Mus_musculus.NCBIM37.64-toMM9.gene-lenghts.txt")
colnames(trl)= c("GeneID", "GeneName", "Length")
trl=trl[!duplicated(trl["GeneName"]),]

setwd("~/giagkas/RNAseq_POLYseq/RNAseq/DEG/")
rn=read.table("edgeRgenes_featurecountExpr-0h.csv", header = T)
rn$GeneName=rownames(rn)
rn=subset(rn, select=c(9,1:8))
rownames(rn)=c()
colnames(rn)=c("GeneName", "logFC", "logCPM", "PValue", "FDR", "KO1", "KO2", "WT1", "WT2")
setwd("~/giagkas/RNAseq_POLYseq/Polysomes/DEG/")
po=read.table("edgeRgenes_featurecountExpr-0h.csv", header = T)
colnames(po)=c("GeneName", "logFC", "logCPM", "PValue", "FDR", "KO1", "KO2", "WT1", "WT2")
setwd("~/giagkas/RNAseq_POLYseq/")
m= merge(rn, po, by.x = "GeneName", by.y = "GeneName", all = T )

m=merge(m, trl, by= "GeneName", all.x = T)

#m[2:5]=NULL
#m[6:9]=NULL
#m[is.na(m)] <- 0

# RPKM=CPM/(transcript_length/1000)
rpkm=m
rpkm$KO1.x=m$KO1.x/(m$Length/1000)
rpkm$KO2.x=m$KO2.x/(m$Length/1000)
rpkm$WT1.x=m$WT1.x/(m$Length/1000)
rpkm$WT2.x=m$WT2.x/(m$Length/1000)
rpkm$KO1.y=m$KO1.y/(m$Length/1000)
rpkm$KO2.y=m$KO2.y/(m$Length/1000)
rpkm$WT1.y=m$WT1.y/(m$Length/1000)
rpkm$WT2.y=m$WT2.y/(m$Length/1000)

m=rpkm
m$KOrna=(m$KO1.x+m$KO2.x)/2
m$WTrna=(m$WT1.x+m$WT2.x)/2
m$KOpoly=(m$KO1.y+m$KO2.y)/2
m$WTpoly=(m$WT1.y+m$WT2.y)/2

am=m
#add gene biotype annotation
am=merge(am,z, by="GeneName", all.x=T)
colnames(am)=c("GeneName", "logFC.RNA",  "logCPM.RNA", "PValue.RNA", "FDR.RNA",    "KO1CPM.RNA",
               "KO2CPM.RNA",    "WT1CPM.RNA",    "WT2CPM.RNA",    "logFC.POLY",  "logCPM.POLY", "PValue.POLY",
               "FDR.POLY",    "KO1CPM.POLY",    "KO2CPM.POLY",    "WT1CPM.POLY",    "WT2CPM.POLY",    "GeneID","Length",
               "koRPKM.RNA", "wtRPKM.RNA", "koRPKM.POLY", "wtRPKM.POLY", "Biotype")

write.table(am, file = "RNA-Poly_0htotal.txt", sep = "\t", quote = F, row.names = F)


##protein coding
m=merge(m,z, by="GeneName", all.x=T)
m=subset(m, (m$Biotype=="protein_coding"))
pm=m
colnames(pm)=c("GeneName", "logFC.RNA",  "logCPM.RNA", "PValue.RNA", "FDR.RNA",    "KO1CPM.RNA",
               "KO2CPM.RNA",    "WT1CPM.RNA",    "WT2CPM.RNA",    "logFC.POLY",  "logCPM.POLY", "PValue.POLY",
               "FDR.POLY",    "KO1CPM.POLY",    "KO2CPM.POLY",    "WT1CPM.POLY",    "WT2CPM.POLY",    "GeneID","Length",
               "koRPKM.RNA", "wtRPKM.RNA", "koRPKM.POLY", "wtRPKM.POLY", "Biotype")


setwd("~/giagkas/RNAseq_POLYseq/")
write.table(pm, file = "RNA-Poly_0hpc.txt", sep = "\t", quote = F, row.names = F)




##
setwd("~/giagkas/RNAseq_POLYseq/plots/")
rpkm=m

rpkm$logKOrna= log2(rpkm$KOrna)
rpkm$logWTrna= log2(rpkm$WTrna)
rpkm$logKOpoly= log2(rpkm$KOpoly)
rpkm$logWTpoly= log2(rpkm$WTpoly)
m2=rpkm
fcru=subset(m2, (m2$WTrna>=1)&(m2$PValue.x<0.05)&(m2$logFC.x>0.584963))
fcrd=subset(m2, (m2$WTrna>=1)&(m2$PValue.x<0.05)&(m2$logFC.x<(-0.584963)))
fcpu=subset(m2, (m2$WTrna>=1)&(m2$WTpoly>=1)&(m2$PValue.y<0.05)&(m2$logFC.y>0.584963))
fcpd=subset(m2, (m2$WTrna>=1)&(m2$WTpoly>=1)&(m2$PValue.y<0.05)&(m2$logFC.y<(-0.584963)))


png(filename = "Abundance_Translation_0h_PC_wtThreshold.png", width =1200, height = 600)
par(mfrow=c(1,2),font.axis=2, cex.axis=1.5, cex.lab=1.5, lwd=2)
plot(m2$logWTrna,m2$logKOrna, pch=".", ylim = c(0,15),xlim = c(0,15), main= "Abundance M0", xlab = "log2RPKM WT",
     ylab = "log2RPKM KO")
points(fcru$logWTrna, fcru$logKOrna, pch="*", col="red")
points(fcrd$logWTrna, fcrd$logKOrna, pch="*", col="green")
legend("topleft", c(paste("FC>1.5","(n=",nrow(fcru),")", sep = " "),
                    paste("FC<-1.5","(n=",nrow(fcrd),")", sep = " ")), bty = "n", text.font=2, cex=1.5)
plot(m2$logWTpoly,m2$logKOpoly, pch=".", ylim = c(0,15), xlim = c(0,15),main= "Translation M0", xlab = "log2RPKM WT",
     ylab = "log2RPKM KO")
points(fcpu$logWTpoly, fcpu$logKOpoly, pch="*", col="red")
points(fcpd$logWTpoly, fcpd$logKOpoly, pch="*", col="green")
legend("topleft", c(paste("FC>1.5","(n=",nrow(fcpu),")", sep = " "),
                    paste("FC<-1.5","(n=",nrow(fcpd),")", sep = " ")), bty = "n", text.font=2, cex=1.5)

dev.off()



##
setwd("~/giagkas/RNAseq_POLYseq/plots/")
#m=subset(m, (m$WTrna>=1)&(m$WTpoly>=1))
m=subset(m, (m$WTrna>=1)&(m$WTpoly>=1)&(m$KOrna>=1)&(m$KOpoly>=1))
rpkm=m

rpkm$logKOrna= log2(rpkm$KOrna)
rpkm$logWTrna= log2(rpkm$WTrna)
rpkm$logKOpoly= log2(rpkm$KOpoly)
rpkm$logWTpoly= log2(rpkm$WTpoly)
m2=rpkm

fcru=subset(m2, (m2$PValue.x<0.05)&m2$logFC.x>0.584963)
fcrd=subset(m2, (m2$PValue.x<0.05)&m2$logFC.x<(-0.584963))
fcpu=subset(m2, (m2$PValue.y<0.05)&m2$logFC.y>0.584963)
fcpd=subset(m2, (m2$PValue.y<0.05)&m2$logFC.y<(-0.584963))

png(filename = "Abundance_Translation_0h_PC.png", width =1200, height = 600)
par(mfrow=c(1,2),font.axis=2, cex.axis=1.5, cex.lab=1.5, lwd=2)
plot(m2$logWTrna,m2$logKOrna, pch=".", ylim = c(0,15),xlim = c(0,15), main= "Abundance M0", xlab = "log2RPKM WT",
     ylab = "log2RPKM KO")
points(fcru$logWTrna, fcru$logKOrna, pch="*", col="red")
points(fcrd$logWTrna, fcrd$logKOrna, pch="*", col="green")
legend("topleft", c(paste("FC>1.5","(n=",nrow(fcru),")", sep = " "),
                    paste("FC<-1.5","(n=",nrow(fcrd),")", sep = " ")), bty = "n", text.font=2, cex=1.5)
plot(m2$logWTpoly,m2$logKOpoly, pch=".", ylim = c(0,15), xlim = c(0,15),main= "Translation M0", xlab = "log2RPKM WT",
     ylab = "log2RPKM KO")
points(fcpu$logWTpoly, fcpu$logKOpoly, pch="*", col="red")
points(fcpd$logWTpoly, fcpd$logKOpoly, pch="*", col="green")
legend("topleft", c(paste("FC>1.5","(n=",nrow(fcpu),")", sep = " "),
                    paste("FC<-1.5","(n=",nrow(fcpd),")", sep = " ")), bty = "n", text.font=2, cex=1.5)

dev.off()









##### M1-2hs  rpkm - Keep only RNAseq protein Coding RPKM>=1 
setwd("~/Documents/norm_Polys_to_RNAseq_DEGs/")
trl=read.table("Mus_musculus.NCBIM37.64-toMM9.gene-lenghts.txt")
colnames(trl)= c("GeneID", "GeneName", "Length")
trl=trl[!duplicated(trl["GeneName"]),]

setwd("~/giagkas/RNAseq_POLYseq/RNAseq/DEG/")
rn=read.table("edgeRgenes_featurecountExpr-2h.csv", header = T)
rn$GeneName=rownames(rn)
rn=subset(rn, select=c(9,1:8))
rownames(rn)=c()
colnames(rn)=c("GeneName", "logFC", "logCPM", "PValue", "FDR", "KO1", "KO2", "WT1", "WT2")
setwd("~/giagkas/RNAseq_POLYseq/Polysomes/DEG/")
po=read.table("edgeRgenes_featurecountExpr-2h.csv", header = T)
colnames(po)=c("GeneName", "logFC", "logCPM", "PValue", "FDR", "KO1", "KO2", "WT1", "WT2")
setwd("~/giagkas/RNAseq_POLYseq/")
m= merge(rn, po, by.x = "GeneName", by.y = "GeneName", all = T )

m=merge(m, trl, by= "GeneName", all.x = T)

#m[2:5]=NULL
#m[6:9]=NULL
#m[is.na(m)] <- 0

# RPKM=CPM/(transcript_length/1000)
rpkm=m
rpkm$KO1.x=m$KO1.x/(m$Length/1000)
rpkm$KO2.x=m$KO2.x/(m$Length/1000)
rpkm$WT1.x=m$WT1.x/(m$Length/1000)
rpkm$WT2.x=m$WT2.x/(m$Length/1000)
rpkm$KO1.y=m$KO1.y/(m$Length/1000)
rpkm$KO2.y=m$KO2.y/(m$Length/1000)
rpkm$WT1.y=m$WT1.y/(m$Length/1000)
rpkm$WT2.y=m$WT2.y/(m$Length/1000)

m=rpkm
m$KOrna=(m$KO1.x+m$KO2.x)/2
m$WTrna=(m$WT1.x+m$WT2.x)/2
m$KOpoly=(m$KO1.y+m$KO2.y)/2
m$WTpoly=(m$WT1.y+m$WT2.y)/2

am=m
#add gene biotype annotation
am=merge(am,z, by="GeneName", all.x=T)
colnames(am)=c("GeneName", "logFC.RNA",  "logCPM.RNA", "PValue.RNA", "FDR.RNA",    "KO1CPM.RNA",
               "KO2CPM.RNA",    "WT1CPM.RNA",    "WT2CPM.RNA",    "logFC.POLY",  "logCPM.POLY", "PValue.POLY",
               "FDR.POLY",    "KO1CPM.POLY",    "KO2CPM.POLY",    "WT1CPM.POLY",    "WT2CPM.POLY",    "GeneID","Length",
               "koRPKM.RNA", "wtRPKM.RNA", "koRPKM.POLY", "wtRPKM.POLY", "Biotype")

write.table(am, file = "RNA-Poly_2htotal.txt", sep = "\t", quote = F, row.names = F)


##protein coding
m=merge(m,z, by="GeneName", all.x=T)
m=subset(m, (m$Biotype=="protein_coding"))
pm=m
colnames(pm)=c("GeneName", "logFC.RNA",  "logCPM.RNA", "PValue.RNA", "FDR.RNA",    "KO1CPM.RNA",
               "KO2CPM.RNA",    "WT1CPM.RNA",    "WT2CPM.RNA",    "logFC.POLY",  "logCPM.POLY", "PValue.POLY",
               "FDR.POLY",    "KO1CPM.POLY",    "KO2CPM.POLY",    "WT1CPM.POLY",    "WT2CPM.POLY",    "GeneID","Length",
               "koRPKM.RNA", "wtRPKM.RNA", "koRPKM.POLY", "wtRPKM.POLY", "Biotype")


setwd("~/giagkas/RNAseq_POLYseq/")
write.table(pm, file = "RNA-Poly_2hpc.txt", sep = "\t", quote = F, row.names = F)

###
setwd("~/giagkas/RNAseq_POLYseq/plots/")
rpkm=m

rpkm$logKOrna= log2(rpkm$KOrna)
rpkm$logWTrna= log2(rpkm$WTrna)
rpkm$logKOpoly= log2(rpkm$KOpoly)
rpkm$logWTpoly= log2(rpkm$WTpoly)
m2=rpkm

fcru=subset(m2, (m2$WTrna>=1)&(m2$PValue.x<0.05)&(m2$logFC.x>0.584963))
fcrd=subset(m2, (m2$WTrna>=1)&(m2$PValue.x<0.05)&(m2$logFC.x<(-0.584963)))
fcpu=subset(m2, (m2$WTrna>=1)&(m2$WTpoly>=1)&(m2$PValue.y<0.05)&(m2$logFC.y>0.584963))
fcpd=subset(m2, (m2$WTrna>=1)&(m2$WTpoly>=1)&(m2$PValue.y<0.05)&(m2$logFC.y<(-0.584963)))


png(filename = "Abundance_Translation_2h_PC_wtThreshold.png", width =1200, height = 600)
par(mfrow=c(1,2),font.axis=2, cex.axis=1.5, cex.lab=1.5, lwd=2)
plot(m2$logWTrna,m2$logKOrna, pch=".", ylim = c(0,15),xlim = c(0,15), main= "Abundance 2hs", xlab = "log2RPKM WT",
     ylab = "log2RPKM KO")
points(fcru$logWTrna, fcru$logKOrna, pch="*", col="red")
points(fcrd$logWTrna, fcrd$logKOrna, pch="*", col="green")
legend("topleft", c(paste("FC>1.5","(n=",nrow(fcru),")", sep = " "),
                    paste("FC<-1.5","(n=",nrow(fcrd),")", sep = " ")), bty = "n", text.font=2, cex=1.5)
plot(m2$logWTpoly,m2$logKOpoly, pch=".", ylim = c(0,15), xlim = c(0,15),main= "Translation 2hs", xlab = "log2RPKM WT",
     ylab = "log2RPKM KO")
points(fcpu$logWTpoly, fcpu$logKOpoly, pch="*", col="red")
points(fcpd$logWTpoly, fcpd$logKOpoly, pch="*", col="green")
legend("topleft", c(paste("FC>1.5","(n=",nrow(fcpu),")", sep = " "),
                    paste("FC<-1.5","(n=",nrow(fcpd),")", sep = " ")), bty = "n", text.font=2, cex=1.5)
dev.off()



setwd("~/giagkas/RNAseq_POLYseq/plots/")
#m=subset(m, (m$WTrna>=1)&(m$WTpoly>=1))
m=subset(m, (m$WTrna>=1)&(m$WTpoly>=1)&(m$KOrna>=1)&(m$KOpoly>=1))
rpkm=m

rpkm$logKOrna= log2(rpkm$KOrna)
rpkm$logWTrna= log2(rpkm$WTrna)
rpkm$logKOpoly= log2(rpkm$KOpoly)
rpkm$logWTpoly= log2(rpkm$WTpoly)
m2=rpkm

fcru=subset(m2, (m2$PValue.x<0.05)&m2$logFC.x>0.584963)
fcrd=subset(m2, (m2$PValue.x<0.05)&m2$logFC.x<(-0.584963))
fcpu=subset(m2, (m2$PValue.y<0.05)&m2$logFC.y>0.584963)
fcpd=subset(m2, (m2$PValue.y<0.05)&m2$logFC.y<(-0.584963))

png(filename = "Abundance_Translation_2h_PC.png", width =1200, height = 600)
par(mfrow=c(1,2),font.axis=2, cex.axis=1.5, cex.lab=1.5, lwd=2)
plot(m2$logWTrna,m2$logKOrna, pch=".", ylim = c(0,15),xlim = c(0,15), main= "Abundance 2hs", xlab = "log2RPKM WT",
     ylab = "log2RPKM KO")
points(fcru$logWTrna, fcru$logKOrna, pch="*", col="red")
points(fcrd$logWTrna, fcrd$logKOrna, pch="*", col="green")
legend("topleft", c(paste("FC>1.5","(n=",nrow(fcru),")", sep = " "),
                    paste("FC<-1.5","(n=",nrow(fcrd),")", sep = " ")), bty = "n", text.font=2, cex=1.5)
plot(m2$logWTpoly,m2$logKOpoly, pch=".", ylim = c(0,15), xlim = c(0,15),main= "Translation 2hs", xlab = "log2RPKM WT",
     ylab = "log2RPKM KO")
points(fcpu$logWTpoly, fcpu$logKOpoly, pch="*", col="red")
points(fcpd$logWTpoly, fcpd$logKOpoly, pch="*", col="green")
legend("topleft", c(paste("FC>1.5","(n=",nrow(fcpu),")", sep = " "),
                    paste("FC<-1.5","(n=",nrow(fcpd),")", sep = " ")), bty = "n", text.font=2, cex=1.5)
dev.off()








##### M1-6hs  rpkm - Keep only RNAseq protein Coding RPKM>=1 
setwd("~/Documents/norm_Polys_to_RNAseq_DEGs/")
trl=read.table("Mus_musculus.NCBIM37.64-toMM9.gene-lenghts.txt")
colnames(trl)= c("GeneID", "GeneName", "Length")
trl=trl[!duplicated(trl["GeneName"]),]

setwd("~/giagkas/RNAseq_POLYseq/RNAseq/DEG/")
rn=read.table("edgeRgenes_featurecountExpr-6h.csv", header = T)
rn$GeneName=rownames(rn)
rn=subset(rn, select=c(9,1:8))
rownames(rn)=c()
colnames(rn)=c("GeneName", "logFC", "logCPM", "PValue", "FDR", "KO1", "KO2", "WT1", "WT2")
setwd("~/giagkas/RNAseq_POLYseq/Polysomes/DEG/")
po=read.table("edgeRgenes_featurecountExpr-6h.csv", header = T)
colnames(po)=c("GeneName", "logFC", "logCPM", "PValue", "FDR", "KO1", "KO2", "WT1", "WT2")
setwd("~/giagkas/RNAseq_POLYseq/")
m= merge(rn, po, by.x = "GeneName", by.y = "GeneName", all = T )

m=merge(m, trl, by= "GeneName", all.x = T)

#m[2:5]=NULL
#m[6:9]=NULL
#m[is.na(m)] <- 0

# RPKM=CPM/(transcript_length/1000)
rpkm=m
rpkm$KO1.x=m$KO1.x/(m$Length/1000)
rpkm$KO2.x=m$KO2.x/(m$Length/1000)
rpkm$WT1.x=m$WT1.x/(m$Length/1000)
rpkm$WT2.x=m$WT2.x/(m$Length/1000)
rpkm$KO1.y=m$KO1.y/(m$Length/1000)
rpkm$KO2.y=m$KO2.y/(m$Length/1000)
rpkm$WT1.y=m$WT1.y/(m$Length/1000)
rpkm$WT2.y=m$WT2.y/(m$Length/1000)

m=rpkm
m$KOrna=(m$KO1.x+m$KO2.x)/2
m$WTrna=(m$WT1.x+m$WT2.x)/2
m$KOpoly=(m$KO1.y+m$KO2.y)/2
m$WTpoly=(m$WT1.y+m$WT2.y)/2
am=m
#add gene biotype annotation
am=merge(am,z, by="GeneName", all.x=T)
colnames(am)=c("GeneName", "logFC.RNA",  "logCPM.RNA", "PValue.RNA", "FDR.RNA",    "KO1CPM.RNA",
               "KO2CPM.RNA",    "WT1CPM.RNA",    "WT2CPM.RNA",    "logFC.POLY",  "logCPM.POLY", "PValue.POLY",
               "FDR.POLY",    "KO1CPM.POLY",    "KO2CPM.POLY",    "WT1CPM.POLY",    "WT2CPM.POLY",    "GeneID","Length",
               "koRPKM.RNA", "wtRPKM.RNA", "koRPKM.POLY", "wtRPKM.POLY", "Biotype")

write.table(am, file = "RNA-Poly_6htotal.txt", sep = "\t", quote = F, row.names = F)


##protein coding
m=merge(m,z, by="GeneName", all.x=T)
m=subset(m, (m$Biotype=="protein_coding"))
pm=m
colnames(pm)=c("GeneName", "logFC.RNA",  "logCPM.RNA", "PValue.RNA", "FDR.RNA",    "KO1CPM.RNA",
               "KO2CPM.RNA",    "WT1CPM.RNA",    "WT2CPM.RNA",    "logFC.POLY",  "logCPM.POLY", "PValue.POLY",
               "FDR.POLY",    "KO1CPM.POLY",    "KO2CPM.POLY",    "WT1CPM.POLY",    "WT2CPM.POLY",    "GeneID","Length",
               "koRPKM.RNA", "wtRPKM.RNA", "koRPKM.POLY", "wtRPKM.POLY", "Biotype")



setwd("~/giagkas/RNAseq_POLYseq/")
write.table(pm, file = "RNA-Poly_6hpc.txt", sep = "\t", quote = F, row.names = F)


##

setwd("~/giagkas/RNAseq_POLYseq/plots/")
rpkm=m

rpkm$logKOrna= log2(rpkm$KOrna)
rpkm$logWTrna= log2(rpkm$WTrna)
rpkm$logKOpoly= log2(rpkm$KOpoly)
rpkm$logWTpoly= log2(rpkm$WTpoly)
m2=rpkm

fcru=subset(m2, (m2$WTrna>=1)&(m2$PValue.x<0.05)&(m2$logFC.x>0.584963))
fcrd=subset(m2, (m2$WTrna>=1)&(m2$PValue.x<0.05)&(m2$logFC.x<(-0.584963)))
fcpu=subset(m2, (m2$WTrna>=1)&(m2$WTpoly>=1)&(m2$PValue.y<0.05)&(m2$logFC.y>0.584963))
fcpd=subset(m2, (m2$WTrna>=1)&(m2$WTpoly>=1)&(m2$PValue.y<0.05)&(m2$logFC.y<(-0.584963)))

png(filename = "Abundance_Translation_6h_PC_wtThreshold.png", width =1200, height = 600)
par(mfrow=c(1,2),font.axis=2, cex.axis=1.5, cex.lab=1.5, lwd=2)
plot(m2$logWTrna,m2$logKOrna, pch=".", ylim = c(0,15),xlim = c(0,15), main= "Abundance 6hs", xlab = "log2RPKM WT",
     ylab = "log2RPKM KO")
points(fcru$logWTrna, fcru$logKOrna, pch="*", col="red")
points(fcrd$logWTrna, fcrd$logKOrna, pch="*", col="green")
legend("topleft", c(paste("FC>1.5","(n=",nrow(fcru),")", sep = " "),
                    paste("FC<-1.5","(n=",nrow(fcrd),")", sep = " ")), bty = "n", text.font=2, cex=1.5)
plot(m2$logWTpoly,m2$logKOpoly, pch=".", ylim = c(0,15), xlim = c(0,15),main= "Translation 6hs", xlab = "log2RPKM WT",
     ylab = "log2RPKM KO")
points(fcpu$logWTpoly, fcpu$logKOpoly, pch="*", col="red")
points(fcpd$logWTpoly, fcpd$logKOpoly, pch="*", col="green")
legend("topleft", c(paste("FC>1.5","(n=",nrow(fcpu),")", sep = " "),
                    paste("FC<-1.5","(n=",nrow(fcpd),")", sep = " ")), bty = "n", text.font=2, cex=1.5)

dev.off()


##

setwd("~/giagkas/RNAseq_POLYseq/plots/")
#m=subset(m, (m$WTrna>=1)&(m$WTpoly>=1))
m=subset(m, (m$WTrna>=1)&(m$WTpoly>=1)&(m$KOrna>=1)&(m$KOpoly>=1))
rpkm=m

rpkm$logKOrna= log2(rpkm$KOrna)
rpkm$logWTrna= log2(rpkm$WTrna)
rpkm$logKOpoly= log2(rpkm$KOpoly)
rpkm$logWTpoly= log2(rpkm$WTpoly)
m2=rpkm

fcru=subset(m2, (m2$PValue.x<0.05)&m2$logFC.x>0.584963)
fcrd=subset(m2, (m2$PValue.x<0.05)&m2$logFC.x<(-0.584963))
fcpu=subset(m2, (m2$PValue.y<0.05)&m2$logFC.y>0.584963)
fcpd=subset(m2, (m2$PValue.y<0.05)&m2$logFC.y<(-0.584963))

png(filename = "Abundance_Translation_6h_PC.png", width =1200, height = 600)
par(mfrow=c(1,2),font.axis=2, cex.axis=1.5, cex.lab=1.5, lwd=2)
plot(m2$logWTrna,m2$logKOrna, pch=".", ylim = c(0,15),xlim = c(0,15), main= "Abundance 6hs", xlab = "log2RPKM WT",
     ylab = "log2RPKM KO")
points(fcru$logWTrna, fcru$logKOrna, pch="*", col="red")
points(fcrd$logWTrna, fcrd$logKOrna, pch="*", col="green")
legend("topleft", c(paste("FC>1.5","(n=",nrow(fcru),")", sep = " "),
                    paste("FC<-1.5","(n=",nrow(fcrd),")", sep = " ")), bty = "n", text.font=2, cex=1.5)
plot(m2$logWTpoly,m2$logKOpoly, pch=".", ylim = c(0,15), xlim = c(0,15),main= "Translation 6hs", xlab = "log2RPKM WT",
     ylab = "log2RPKM KO")
points(fcpu$logWTpoly, fcpu$logKOpoly, pch="*", col="red")
points(fcpd$logWTpoly, fcpd$logKOpoly, pch="*", col="green")
legend("topleft", c(paste("FC>1.5","(n=",nrow(fcpu),")", sep = " "),
                    paste("FC<-1.5","(n=",nrow(fcpd),")", sep = " ")), bty = "n", text.font=2, cex=1.5)

dev.off()