setwd("~/Documents/HuR_PARclip_withHeaders/mRNA-1kIntrons_v2/IL4/")
#no collapse to gene level. Transcript coordinates
a1=read.table(textConnection(gsub("_", "\t", readLines("sh-clusters-IL4.txt2-1kIntrons2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG.f.flt.anno1.bed"))))
head(a1)
colnames(a1)=c("TranscriptCoord","ClusterStart", "ClusterEnd", "ClusterID", "Sequence", "ReadCount",
               "ModeLocation","ConversionLocationCount","ConversionEventCount", "NonConversionEventCount","ModeScore",
               "AverageCnvPrc", "GroupConversionEvents", "SDev", "MaxConvPrc", "GeneID","CodingNonCoding", "ExonOrIntron")

#create custome ID
df1=as.data.frame(substr(a1$TranscriptCoord, 1, 18))
colnames(df1)="TranscriptID"
df=cbind(df1, a1)
a1 <- within(df,   ID<- paste(TranscriptID,ClusterID,ModeLocation,ConversionEventCount,  sep="_"))

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







#Genomic coordinates
a2=read.table(textConnection(gsub("_", "\t", readLines("sh-clusters-IL4.txt2-1kIntrons2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG.f.flt.anno1.bed.genomic2.bed"))))
colnames(a2)=c("Chr","ChrStart", "ChrEnd", "TranscriptID", "ClusterID", "Sequence", "ReadCount",
               "ModeLocation", "ConversionLocationCount", "ConversionEventCount", "NonConversionEventCount","ModeScore",
               "AverageCnvPrc", "GroupConversionEvents", "SDev", "MaxConvPrcScore", "Strand")
head(a2)
#create custome ID
a2 <- within(a2,   ID<- paste(TranscriptID,ClusterID,ModeLocation,ConversionEventCount,  sep="_"))
#find which are the duplicated values
n_occur <- data.frame(table(a2$ID))
n_occur[n_occur$Freq> 1,]




#3UTR
a3=read.table(textConnection(gsub("_", "\t", readLines("sh-clusters-IL4.txt2-1kIntrons2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG.f.flt.anno1.bed.genomic2.bed.anno2.exon.3utr.bed"))))
colnames(a3)=c("Chr","ChrStart", "ChrEnd", "TranscriptID", "ClusterID", "Sequence", "ReadCount",
               "ModeLocation", "ConversionLocationCount", "ConversionEventCount","NonConversionEventCount", "ModeScore",
               "AverageCnvPrc", "GroupConversionEvents", "SDev", "MaxConvPrcScore", "Strand","ExonOrIntron")

#create custome ID
a3 <- within(a3,   ID<- paste(TranscriptID,ClusterID,ModeLocation,ConversionEventCount,  sep="_"))
#find which are the duplicated values
n_occur <- data.frame(table(a3$ID))
n_occur[n_occur$Freq> 1,]

##5UTR
a4=read.table(textConnection(gsub("_", "\t", readLines("sh-clusters-IL4.txt2-1kIntrons2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG.f.flt.anno1.bed.genomic2.bed.anno2.exon.5utr.bed"))))
colnames(a4)=c("Chr","ChrStart", "ChrEnd", "TranscriptID", "ClusterID", "Sequence", "ReadCount",
               "ModeLocation", "ConversionLocationCount", "ConversionEventCount","NonConversionEventCount", "ModeScore",
               "AverageCnvPrc", "GroupConversionEvents", "SDev", "MaxConvPrcScore", "Strand","ExonOrIntron")

#create custome ID
a4 <- within(a4,   ID<- paste(TranscriptID,ClusterID,ModeLocation,ConversionEventCount,  sep="_"))
#find which are the duplicated values
n_occur <- data.frame(table(a4$ID))
n_occur[n_occur$Freq> 1,]


#CDS
a5=read.table(textConnection(gsub("_", "\t", readLines("sh-clusters-IL4.txt2-1kIntrons2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG.f.flt.anno1.bed.genomic2.bed.anno2.exon.cds.bed"))))
colnames(a5)=c("Chr","ChrStart", "ChrEnd", "TranscriptID", "ClusterID", "Sequence", "ReadCount",
               "ModeLocation", "ConversionLocationCount", "ConversionEventCount","NonConversionEventCount", "ModeScore",
               "AverageCnvPrc", "GroupConversionEvents", "SDev", "MaxConvPrcScore", "Strand","ExonOrIntron")
#create custome ID
a5 <- within(a5,   ID<- paste(TranscriptID,ClusterID,ModeLocation,ConversionEventCount,  sep="_"))
#find which are the duplicated values
n_occur <- data.frame(table(a5$ID))
n_occur[n_occur$Freq> 1,]




#ncRNA
a6=read.table(textConnection(gsub("_", "\t", readLines("sh-clusters-IL4.txt2-1kIntrons2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG.f.flt.anno1.bed.genomic2.bed.anno2.exon.ncRNA"))))
colnames(a6)=c("Chr","ChrStart", "ChrEnd", "TranscriptID", "ClusterID", "Sequence", "ReadCount",
               "ModeLocation", "ConversionLocationCount", "ConversionEventCount","NonConversionEventCount", "ModeScore",
               "AverageCnvPrc", "GroupConversionEvents", "SDev", "MaxConvPrcScore", "Strand","ExonOrIntron")

#create custome ID
a6 <- within(a6,   ID<- paste(TranscriptID,ClusterID,ModeLocation,ConversionEventCount,  sep="_"))
#find which are the duplicated values
n_occur <- data.frame(table(a6$ID))
n_occur[n_occur$Freq> 1,]




a2$X3UTR=a2$ID%in%a3$ID
a2$X5UTR=a2$ID%in%a4$ID
a2$CDS=a2$ID%in%a5$ID
a2$ncRNA=a2$ID%in%a6$ID
ann=read.table(textConnection(gsub(" ", "", readLines("~/Documents/HuR_PARclip_withHeaders/reference_gtf/Mus_musculus.NCBIM37.64-toMM9.transcript-gene_name.txt"))))
colnames(ann)=c("TranscriptID", "GeneID","GeneName","Biotype")
ma2=merge(ann, a2, by="TranscriptID", all.y = T)


#write.table(ma2, file="IL4_merged.txt", sep = "\t", quote = T, row.names = F)
test=read.table("IL4_merged.txt", header = T)


mf=a1
#test if the length of the sequence is in accordance with cluster start and end
tdf=list()
for(i in 1:dim(mf)[1]){
  x=nchar(as.character(mf[i, 5]))
  y=mf[i, 3]-mf[i, 2]
  test=ifelse((x-y)==1,1,2)
  tdf[i]=test
}
tdf=t(as.data.frame(tdf))
#if true, then seq length is in accordance with cluster start and end
summary(tdf[1:length(tdf),]==1)

mf=a2
#test if the length of the sequence is in accordance with genomic coordinates
cdf=list()
for(i in 1:dim(mf)[1]){
  x=nchar(as.character(mf[i, 6]))
  y=mf[i, 3]-mf[i, 2]
  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

#nimber of genes
n_g <- data.frame(table(ma2$GeneName))
n_g=n_g[n_g$Freq>= 1,]
















#??
a7=read.table("sh-clusters-IL4.txt2-1kIntrons2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG.f.flt.anno1.bed.genomic2.bed.anno2.intron")
a8=read.table("sh-clusters-IL4.txt2-1kIntrons2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG.f.flt.anno1.bed.genomic2.bed.exon-intron")
a9=read.table("sh-clusters-IL4.txt2-1kIntrons2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG.f.flt.anno1.bed.genomic2.bed.best.bed")
a10=read.table("sh-clusters-IL4.txt2-1kIntrons2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG.f.flt.anno1.bed.genomic2.bed.collapsed.best.per-gene-exons.bed")
a11=read.table("sh-clusters-IL4.txt2-1kIntrons2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG.f.flt.anno1.bed.genomic2.bed.anno2")