library(GenomicRanges)
library(GenomicAlignments)
library(GenomicFeatures)

source("/data/results/tools/r/GenomicFeatures/R/my-coordinate-mapping-methods.R")
setwd("/data/results/reference/mmu")

txdb <- makeTxDbFromGFF(file ="Mus_musculus.NCBIM37.64.gtf", format= "gtf", dataSource="Mus_musculus.NCBIM37.64", organism= "Mus musculus")
read_csv<-read.csv(file="/data/images/proton/DKlab/mr/parclip/paralyzer/PARalyzer_v1_1b/res/recommended_settings_gt5tc_gt0.25TtoC2/sh-clusters-IFN.txt2.csv.bed.gt5tc.gt0.25TtoC2.plus.noIGG.anno.csv", sep="", header = FALSE)

m_cols<-data.frame(GeneName=read_csv[,1], CDS_start=read_csv[,2],CDS_end=read_csv[,3])
df<-data.frame(seqnames=read_csv[,4], start=read_csv[,5], end=read_csv[,6], mcols=m_cols )
df<-data.frame(seqnames=read_csv[,4], start=read_csv[,5], end=read_csv[,6])
gr_ranges<-makeGRangesFromDataFrame(df)

mapTrRegion2Genome=function(ts,te,trid){
 ex=exons(txdb, filter=list(tx_name = trid ));
 r=ranges(ex);
 s=start(r);
 e=end(r);rl=NULL;
 st=as.logical(strand(ex)[1]=="-");
 st2=as.character(strand(ex)[1]);
 sn=as.character(seqnames(ex)[1]);
 if (sn=="MT"){sn="M";}
 sn=paste("chr",sn,sep="");
 ne=length(ex);
 sok=0;
 if (st) {
  for (j in (ne):1) {
   start = s[j]
   end = e[j];
   width = end - start + 1;
   if (sok==1){
    rl=append(rl,end);
   }
   if (ts <= width) {
     if (sok==0){
       mts=end-ts;
       mtsa=start;
       rl=c(mts+1);
       sok=1;sei=j;
     }
   } else{
     ts = ts- width;
   }
   if (te <= width)  {
     rl=append(rl,end-te+1);	
     break;
   } else {
     te=te-width;
     if (sok==1){
      rl=append(rl,start);
     }
   }
  }   
  return(c(st2,sn,rl));
} else {
   for (j in 1:ne) {
   start = s[j]
   end = e[j];
   width = end - start + 1;
   if (sok==1){
    rl=append(rl,start);
   }
#   print (c(j,ts,te,start,end,sok))
   if (ts <= width) {
      if (sok==0){
        mts=start+ts;
        rl=c(mts-1);
	mtsa=end;
	sok=1;sei=j;
      }
   } else{
    ts = ts - width;
   }
   if (te <= width) {
     rl=append(rl,start+te-1);	   
     break;
   } else {     
     te=te-width;
     if (sok==1){
      rl=append(rl,end);
     }
   }
  }   
  return(c(st2,sn,rl));
 }
}

mapAllTrRegions2Genome=function(gr_ranges){
 all_s=NULL;
 all_e=NULL;
 all_st=NULL;
 all_id=NULL;
 all_sn=NULL;
 all_i=NULL;
 ids=as.vector(gr_ranges@seqnames)
 for (i in 1:length(ids)){
    m=mapTrRegion2Genome(start(gr_ranges)[i],end(gr_ranges)[i],ids[i]);
    for (r in seq(3,length(m),2)){
      all_id=append(all_id,ids[i]);
      all_st=append(all_st,m[1]);
      all_i=append(all_i,i);
      all_sn=append(all_sn,m[2]);
      all_s=append(all_s,m[r]);
      all_e=append(all_e,m[r+1]);
    }
 }
 return(cbind(all_sn,all_s,all_e,all_st,all_id,all_i))
}


mm=mapAllTrRegions2Genome(gr_ranges)
write.table(mm,"IFN-mapped2.csv", sep="\t", row.names=FALSE)