setwd("~/giagkas/parclip/1kIntron/IL4/")
a1=read.table("sh-clusters-IL4.txt2-1kIntrons2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG.f.flt.anno1.bed.genomic2.bed.per-gene.bed2",sep="\t")
a1$V13 <- gsub("_", "\t", a1$V4)
write.table(a1, file = "IL4-bed-per-gene.txt", quote = F, sep = "\t", row.names = F, col.names = F)
a1=read.delim("IL4-bed-per-gene.txt", header = F)
head(a1)

colnames(a1)=c("Chr","ChrStart", "ChrEnd","ID","MaxConvScore","Strand", "GeneName",  "Biotype", "2ryAnn_1",
               "2ryAnn_1_Occur","2ryAnn_2", "2ryAnn_2_Occur","TranscriptID", "ClusterID", "Sequence", "ReadCount",
               "ModeLocation","ConversionLocationCount","ConversionEventCount", "NonConversionEventCount","ModeScore",
               "AverageCnvPrc", "GroupConversionEvents", "SDev")


#find which are the duplicated values
n_occur <- data.frame(table(a1$ID))
n_occur[n_occur$Freq> 1,]

#3utr
a2=read.table("sh-clusters-IL4.txt2-1kIntrons2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG.f.flt.anno1.bed.genomic2.bed.3utr.best.bed",sep="\t")


#5utr
a3=read.table("sh-clusters-IL4.txt2-1kIntrons2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG.f.flt.anno1.bed.genomic2.bed.5utr.best.bed",sep="\t")

#cds
a4=read.table("sh-clusters-IL4.txt2-1kIntrons2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG.f.flt.anno1.bed.genomic2.bed.cds.best.bed",sep="\t")

#intron
a5=read.table("sh-clusters-IL4.txt2-1kIntrons2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG.f.flt.anno1.bed.genomic2.bed.intron.best.bed",sep="\t")


#ncRNA
a6=read.table("sh-clusters-IL4.txt2-1kIntrons2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG.f.flt.anno1.bed.genomic2.bed.ncRNA.best.bed",sep="\t")


#exon
a7=read.table("sh-clusters-IL4.txt2-1kIntrons2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG.f.flt.anno1.bed.genomic2.bed.exon.best.bed",sep="\t")

a1$X3UTR=a1$ID%in%a2$V4
summary(a1$X3UTR)
a1$X5UTR=a1$ID%in%a3$V4
summary(a1$X5UTR)
a1$CDS=a1$ID%in%a4$V4
summary(a1$CDS)
a1$intron=a1$ID%in%a5$V4
summary(a1$intron)
a1$ncRNA=a1$ID%in%a6$V4
summary(a1$ncRNA)
a1$exon=a1$ID%in%a7$V4
summary(a1$exon)

colnames(a1)
m=subset(a1, select = c(7,13,8,1,2,3,6,14,15,16,17,18,19,20,21,22,23,24,5,25,26,27,28,29,30,9,10,11,12,4))


#find which are the duplicated values
n_occur <- data.frame(table(m$ID))
n_occur=subset(n_occur, n_occur$Freq>1)
mt=subset(m, m$ID%in%n_occur$Var1)
which(mt$TranscriptID=="ENSMUST00000097424")
mt[74,]
mt[230,]
#write.table(mt, file = "IL4_duplicates.txt", col.names = T, row.names = F, quote = F, sep = "\t")


mf=m
#test if the length of the sequence is in accordance with genomic coordinates
colnames(mf)
cdf=list()
for(i in 1:dim(mf)[1]){
  x=nchar(as.character(mf[i, 9]))
  y=mf[i, 6]-mf[i, 5]
  test=ifelse((x-y)==1,1,2)
  cdf[i]=test
}
cdf=t(as.data.frame(cdf))
#if true, then seq length is in accordance with genomic coordinates
summary(cdf[1:length(cdf),]==1)
which(cdf[1:length(cdf),]!=1)
tr=as.integer(table(cdf[1:length(cdf),]==1)["TRUE"])
fs=as.integer(table(cdf[1:length(cdf),]==1)["FALSE"])
fs/tr

tst=3681
mf[tst,]
mf[tst,6]-mf[tst,5]
nchar(as.character(mf[tst, 9]))


###
t1=subset(m, (m$X3UTR=="FALSE")&(m$X5UTR=="FALSE")&(m$CDS=="FALSE")&(m$intron=="FALSE")&(m$ncRNA=="FALSE")&(m$exon=="FALSE"))

t2=subset(m, (m$X3UTR=="FALSE")&(m$X5UTR=="FALSE")&(m$CDS=="FALSE")&(m$intron=="FALSE")&(m$ncRNA=="TRUE")&(m$exon=="TRUE"))

t3=subset(m, (m$ncRNA=="TRUE"))
summary(t3$Biotype)
summary(t2$Biotype)

#3utr
t4=subset(m, (m$X3UTR=="TRUE"))
colnames(t4)
t4=subset(t4, select = c(1,2,3,4,5,6,30))
tm4=merge(t4, a2, by.x="ID", by.y="V4", all = T)

#5utr
t5=subset(m, (m$X5UTR=="TRUE"))
colnames(t5)
t5=subset(t5, select = c(1,2,3,4,5,6,30))
tm5=merge(t5, a3, by.x="ID", by.y="V4", all = T)


#cds
t6=subset(m, (m$CDS=="TRUE"))
colnames(t6)
t6=subset(t6, select = c(1,2,3,4,5,6,30))
tm6=merge(t6, a4, by.x="ID", by.y="V4", all = T)


#intron
t7=subset(m, (m$intron=="TRUE"))
colnames(t7)
t7=subset(t7, select = c(1,2,3,4,5,6,30))
tm7=merge(t7, a5, by.x="ID", by.y="V4", all = T)


#ncRNA
t8=subset(m, (m$ncRNA=="TRUE"))
colnames(t8)
t8=subset(t8, select = c(1,2,3,4,5,6,30))
tm8=merge(t8, a6, by.x="ID", by.y="V4", all = T)


#exon
t9=subset(m, (m$exon=="TRUE"))
colnames(t9)
t9=subset(t9, select = c(1,2,3,4,5,6,30))
tm9=merge(t9, a7, by.x="ID", by.y="V4", all = T)

t9=m[2503,]




###
setwd("~/Documents/norm_Polys_to_RNAseq_DEGs/RNAseq/")
rn=read.table("M2_IL4_RNAseq.txt", header = T)
colnames(rn)=c("GeneName", "logFC_RNA","logCPM_RNA","PValue_RNA", "FDR_RNA",
               "KO1CPM_RNA","KO2CPM_RNA","WT1CPM_RNA","WT2CPM_RNA", "GeneID",
               "GeneLength", "KORPKM_RNA", "WTRPKM_RNA")
rn$GeneID=NULL
rn$GeneLength=NULL
setwd("~/Documents/norm_Polys_to_RNAseq_DEGs/Polysomes/")
rp=read.table("M2_IL4_Polysomes.txt", header = T)
colnames(rp)=c("GeneName", "logFC_Poly","logCPM_Poly","PValue_Poly", "FDR_Poly",
               "KO1CPM_Poly","KO2CPM_Poly","WT1CPM_Poly","WT2CPM_Poly", "GeneID",
               "GeneLength", "KORPKM_Poly", "WTRPKM_Poly")
rp$GeneID=NULL
rp$GeneLength=NULL
setwd("~/giagkas/parclip/1kIntron/IL4/")


m1=merge(m, rn, by="GeneName", all = T)
m2=merge(m1, rp, by="GeneName", all = T)


m2$CLI=log2(1+m2$ConversionEventCount)
m2$nrmReadCount=m2$ReadCount/m2$WTRPKM_RNA
m2$nrmCnvEventCount=m2$ConversionEventCount/m2$WTRPKM_RNA
m2$nrmNonCnvEventCount=m2$NonConversionEventCount/m2$WTRPKM_RNA
m2$nrmCLI=log2(1+m2$nrmCnvEventCount)
colnames(m2)
fm=subset(m2, select = c(1:19,51:55, 20:50))
colnames(fm)

#fm$ModeLocation=NULL
setwd("~/giagkas/parclip/1kIntron/IL4/")
#write.table(fm, file="IL4_parclip_RNAseq_POLYseq.txt", quote = F, row.names = F, sep = "\t")