###compute the expression of transrcripts
###20141227

args <- commandArgs(TRUE)
library(rjags)
library(coda)

data_count= read.table(args)
data_count=data_count[,-c(9,11)]
colnames(data_count)=c("chromosome","source","feature","start","end","score","strand","frame","gene_id","transcript_id","counts","non.zero.base.number","length","fraction")


data_count[,"gene_id"]=as.character(as.matrix(data_count[,"gene_id"]))
data_count[,"transcript_id"]=as.character(as.matrix(data_count[,"transcript_id"]))

data_count[,"gene_id"]=as.character(strsplit(data_count[,"gene_id"],";"))
data_count[,"transcript_id"]=as.character(strsplit(data_count[,"transcript_id"],";"))

gene_list=unique(data_count[,"gene_id"])

##################################estimate expression ####################################################################

multi_result=data.frame('Mean'=0,'SD'=0,'Naive SE'=0,'Time_series SE'=0,id="a",'confi_2.5'=0,'confi_25'=0,'confi_50'=0,'confi_75'=0,'confi_97.5'=0)
single_result=data.frame('Mean'=0,'SD'=0,'Naive SE'=0,'Time_series SE'=0,id="a",'confi_2.5'=0,'confi_25'=0,'confi_50'=0,'confi_75'=0,'confi_97.5'=0)

for(i in 1:length(gene_list)){
   gene_level=data_count[which(data_count[,"gene_id"]==gene_list[i]),]
   start=gene_level[,"start"]
   exon=unique(start)
   exon_length=rep(0,length(exon))
   for(j in 1:length(exon)){
       exon_length[j]=gene_level[which(gene_level[,"start"]==exon[j])[1],"length"]
   }
   trans=unique(gene_level[,"transcript_id"])
   A=matrix(0,length(trans),length(exon))##row corresponding to transcript, col is exon
   row.names(A)=trans
   exon_counts=rep(0,length(exon))
   for(k in 1:length(trans)){
       tmp=gene_level[which(gene_level[,"transcript_id"]==trans[k]),]
       for(m in 1:dim(tmp)[1]){
            index=which(exon==tmp[m,"start"])
            A[k,index]=1
            exon_counts[index]=tmp[m,"counts"]
       }   
   }
   length_transcript=as.numeric(as.character(A%*%exon_length))
   zeros=rep(0,length(exon))
   
   if(length(trans)>1){
        data_list=list('E'=length(exon),'T'=length(trans),'Sample'=exon_counts,'zeros'=zeros,'l'=exon_length,'A'=t(A),'r'=0,'mu'=1500,'alpha'=rep(1,length(trans)))
         jags <- jags.model('GPD_multi_iso.bug',data = data_list,n.chains = 1,n.adapt = 10000) 
         update(jags, 10000)
         parameter=c('N','p')
         gene.sample=coda.samples(jags,parameter,n.iter=5000,thin=1)
          statis=summary(gene.sample)
          id=c(gene_list[i],trans)
          Empi=as.data.frame(statis$statistics)
          Quan=as.data.frame(statis$quantiles)
          gene.result=cbind(cbind(Empi,id),Quan)
          colnames(gene.result)=c("Mean","SD","Naive.SE","Time_series.SE","id",'confi_2.5','confi_25','confi_50','confi_75','confi_97.5')
          multi_result=rbind(multi_result,gene.result)
          
   }
   else{
      
         data_list=list('E'=length(exon),'Sample'=exon_counts,'zeros'=zeros,'l'=exon_length,'A'=t(A),'r'=0,'mu'=1500) 
         jags <- jags.model('GPD_single_iso.bug',data = data_list,n.chains = 1,n.adapt = 10000) 
         update(jags, 10000)
         parameter=c('f')
         gene.sample=coda.samples(jags,parameter,n.iter=5000,thin=1)
         statis=summary(gene.sample)
         id=gene_list[i]
         Empi=t(as.data.frame(statis$statistics))
         Quan=t(as.data.frame(statis$quantiles))
         gene.result=cbind(cbind(Empi,id),Quan)
         colnames(gene.result)=c("Mean","SD","Naive.SE","Time_series.SE","id",'confi_2.5','confi_25','confi_50','confi_75','confi_97.5')
         single_result=rbind(single_result,gene.result)
         
   }
   
}
multi_result=multi_result[-1,]
write.table(multi_result,file=paste(args,"_GIExp_multi_result.txt",sep=""),quote=FALSE,row.names=FALSE,col.names=TRUE)
single_result=single_result[-1,]
write.table(single_result,file=paste(args,"_GIExp_single_result.txt",sep=""),quote=FALSE,row.names=FALSE,col.names=TRUE)

















