annotate src/splicetrap.estimate.cpp @ 5:2ebca9da5e42 draft default tip

planemo upload

//Author: Jie Wu@CSHL
//TO replace the original Perl code

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <string.h>
#include <unistd.h>
#include <sstream>
#include <map>
#include <vector>
#include <math.h>
#include <time.h>
#include <limits>
#include <iomanip>

using namespace std;

int MAX_LINE_LEN = 1024*1024;

void printusage()
{
    cout<< "\tUsage:"<<endl;
    cout<<"\t-s\tread_size"<<endl;
    cout<<"\t-b\tIRM file"<<endl;
    cout<<"\t-f\tFZM file"<<endl;
    cout<<"\t-1\tmapping result 1"<<endl;
    cout<<"\t-2\tmapping result 2"<<endl;
    cout<<"\t-d\tTXdb database file"<<endl;
    cout<<"\t-o\toutput prefix"<<endl;
    
}

struct ReadPairInfo
{
	bool isPaired;
	vector<int> IsoFlag; //flag as 1 if mapped to the isoform
	vector<int> Size; //fragment size in the isoform
	vector<int> Pos;
};

class TXdb_entry{
	public:
	 string chrid;
	 //int start;
	 //int end;
	 string id;
	 bool iscomplete;
	 int start;
	 int end;
	 int exonstarts[3];
	 char strand;
	 int exonsize[3];
	 int len1;
	 int len2;
	 vector< ReadPairInfo > pairs;
	 int snum1; //single-end reads
	 int snum2;
	 int snum12;
	 float ir, bir;
	 int enum1; //exon body reads
	 int enum2;
	 int enum3;
	 
	 int jnum12;
	 int jnum23;
	 int jnum13;

	 int totalreadnum;
	 void estimate(vector<double> &FZM, map<string, vector<double> > &IRM, int &read_size)
	 {
		 float Maxe=0;
		 float Max = -numeric_limits<float>::max();
		 float BMaxe=0;
		 float BMax=-numeric_limits<float>::max();
		 for(float e1=0.001; e1<1;e1=e1+0.001)
		 {
			 float efix1 = e1*len1/
			               (e1*exonsize[1]+len2);
			 float efix2 = 1-efix1;
			 float LL = snum12*log( efix1/(len1-read_size+1)
			          + efix2/(len2-read_size+1) )
			          + snum1*log(efix1/(len1-read_size+1)) 
				  + snum2*log(efix2/(len2-read_size+1));
				  //cout<<LL<<endl;
				  //cout<<efix1<<"\t"  <<pairs.size()<<"\t"  <<len2<<"\t"	<<endl;
		         //int num1=0, num2=0,num12=0;
		         for(int i=0;i<pairs.size();i++)
			  {
                            if(pairs[i].IsoFlag[0] == 1)
			    {
				    if(pairs[i].IsoFlag[1] == 0)
				    {
					   // num1++;
					    float tmp=FZM[pairs[i].Size[0]]*efix1/(len1-pairs[i].Size[0]+1);
					    if (tmp == 0)
					    {
					    	 LL = LL - 308;
					    }
					    else
					    {
						 LL = LL + log(tmp);
					    }
				    }
				    else
				    {
					   // num12++;
					    float tmp=(FZM[pairs[i].Size[0]]*efix1/(len1-pairs[i].Size[0]+1)) 
					          + (FZM[pairs[i].Size[1]]*efix2/(len2-pairs[i].Size[1]+1));
				  	    if(tmp ==0)
					    {
					    LL = LL -308;
					    }
					    else
					    {
						    LL = LL + log(tmp);
					    }
				    }
			    }
			    else
			    {
				    //num2++;
				    float tmp = FZM[pairs[i].Size[1]]*efix2/(len2-pairs[i].Size[1]+1);
				    if(tmp ==0)
				    {
				    LL= LL - 308; 
				    }
				    else
				    {
					    LL = LL +log(tmp);
				    }
			    }
			  }
			  //cout<<num1<<"\t"<<num2<<"\t"<<num12<<"\n";
			  //cout<<LL<<endl;
			  
			  if(IRM["CA"].size() > 0)
			  {
				  float BLL;
				  string eventtype =id.substr(0,2);
				  if(!eventtype.compare("ME"))
				  {
					  eventtype = "CA";
				  }
				  BLL = LL + log(IRM[eventtype][int(e1/0.001)]);
				  if(BLL > BMax)
				  { 
					  BMax = BLL;
					  BMaxe = e1;
				  }
			  }
				  //cout<<LL<<endl;
			  if(LL > Max)
			  {
				  Max= LL;
				  Maxe = e1;

			  }
				  
		 }
		 ir= Maxe;
		 bir= BMaxe;
	 }
		
	 TXdb_entry(char *line)
	 {
		 size_t found;
		 string linestr(line);
		 found = linestr.find("[L]");
		 if(found != string::npos)
		 {
			 iscomplete=1;
		
		 	char *token;

		 	token = strtok(line,"\t");
		 	if(token==NULL)
		 	{
				 cout<<"error1"<<endl;
				 exit(0);
			 }
			 else
			 {
				 chrid = token;
			 }
			 token = strtok(NULL,"\t");
			 if(token==NULL)
			 {
				 cout<<"error2"<<endl;
				 exit(0);
			 }
			 else
			 {
				 sscanf(token,"%d",&start);
			 }
			 token = strtok(NULL,"\t");
        	         if(token==NULL)
                	 {
                        	 cout<<"error3"<<endl;
	                         exit(0);
        	         }
                	 else
	                 {
        	                 sscanf(token,"%d",&end);
	                 }
			 token = strtok(NULL,"\t");
			 if(token==NULL)
			 {
				 cout<<"error4"<<endl;exit(0);
			 }
			 else
			 {
				 id = token;
				 id = id.substr(0,id.size()-3);
			 }
			 token = strtok(NULL,"\t");
			 token = strtok(NULL,"\t");
			 if(token==NULL)
			 {
				  cout<<"error5"<<endl;exit(0);
			 }
			 else
			 {
				 sscanf(token,"%c",&strand);
			 }
			 for(int i=0;i<5;i++)
			 {
				 token=strtok(NULL,"\t");
			 }
			 if(token==NULL){
				 cout<<"error6"<<endl;exit(0);
			 }
			 else
			 {
				 string str_tmp(token);
				 size_t found_tmp;
				 int *pexonsize = exonsize;
				 while((found_tmp = str_tmp.find(",")) != string::npos){
					 sscanf((str_tmp.substr(0,found_tmp)).c_str(),"%d",pexonsize);
					 pexonsize++;
					 str_tmp.erase(0,found_tmp+1);
				 }
				 sscanf(str_tmp.c_str(),"%d",pexonsize); 
				 
			 }
			 token=strtok(NULL,"\t");
			 
			if(token==NULL){
				 cout<<"error7"<<endl;exit(0);
			 }
			 else
			 {
				 string str_tmp(token);
				 size_t found_tmp;
				 int *pexonstarts = exonstarts;
				 while((found_tmp = str_tmp.find(",")) != string::npos){
					 sscanf((str_tmp.substr(0,found_tmp)).c_str(),"%d",pexonstarts);
					 pexonstarts++;
					 str_tmp.erase(0,found_tmp+1);
				 }
				 sscanf(str_tmp.c_str(),"%d",pexonstarts); 
				 
			 }
		 len1 = exonsize[0]+exonsize[1]+exonsize[2];
			 len2 = exonsize[0]+exonsize[2];
			 snum12=0;
			 snum1 =0;
			 snum2 =0;
			 enum1=0;
			 enum2=0;
			 enum3=0;
			 jnum12=0;
			 jnum13=0;
			 jnum23=0;
			 totalreadnum=0;			
		 }
		else
		{
			iscomplete = 0;
		}
	 }
};

class ReadInfo
{
	public:
	bool isUseful;
	char *pline;
	void scan(map<string, vector<int> > &map_pos )
	//eventname, pos in L/S iso 
	{
		char *token;
		token = strtok(pline, "\t");
		token = strtok(NULL, "\t");
		token = strtok(NULL, "\t");
		token = strtok(NULL, "\t");
		string mapinfo(token);
		size_t  beg, end;
		end = 0;
		while( (beg= mapinfo.find_first_of("/", end))!=string::npos)
		{
			beg ++;
			end = mapinfo.find_first_of("[",beg);
			string eventid = mapinfo.substr( beg, end-beg ) ;
			vector<int> pos(2,-1);
			if(!map_pos.count(eventid))
			{
				map_pos.insert(pair<string, vector<int> >
					(eventid, pos) );
			}
			end ++;
			string LorS = mapinfo.substr( end,1 );
			beg = mapinfo.find_first_of(":",end );
			beg ++;
			end = mapinfo.find_first_of(",",beg );
			string pos_str(mapinfo.substr(beg, end-2-beg));
			
			if(LorS.compare("L") == 0)
			{
				sscanf(pos_str.c_str(),"%d",&(map_pos[eventid][0]) );
			}
			else
			{
				sscanf(pos_str.c_str(),"%d",&(map_pos[eventid][1]) );
			}
			
			//cout<<eventid<<"\t"<<map_pos[eventid][0]<<"\t"<<map_pos[eventid][1]<<endl;
		}
		
	}
	ReadInfo(char *line)
	{
		string line_str(line);
		string line_str_last2=line_str.substr(line_str.size()-2,2);
		if(line_str_last2.compare("NM") == 0 )
		{
			isUseful = 0;
		}
		else if(line_str_last2.compare("MT") == 0 )
		{
			isUseful = 0;
		}
		else
		{
			isUseful = 1;
			pline = line;
		}
	}
};



class ReadPair
{
	public:
	map<string, ReadPairInfo> mappinginfo;
	ReadPair(map<string, vector<int> > &read1_mappos, map<string, vector<int> > &read2_mappos, int &read_size)
	{
		map<string, vector<int> >::iterator it1 = read1_mappos.begin();
		for(;it1!=read1_mappos.end();it1++)
		{
			ReadPairInfo ReadPairInfo_tmp;
			ReadPairInfo_tmp.IsoFlag.resize(2);
			ReadPairInfo_tmp.Size.resize(2);
			//ReadPairInfo_tmp.Pos.resize(4);
			ReadPairInfo_tmp.Pos.assign(4,-1);
			string eventid = it1->first;
			if(read2_mappos.count(eventid) > 0)
			{
				ReadPairInfo_tmp.isPaired = 0;
				for(int i=0;i<2;i++)
				{
					//cout<<read1_mappos[eventid][i]<<endl;
					if(read2_mappos[eventid][i] == -1 || read1_mappos[eventid][i] == -1)
					{
						ReadPairInfo_tmp.IsoFlag[i] = 0;
						ReadPairInfo_tmp.Size[i]    = 0;
						//discard reads like this
					}
					else
					{

						ReadPairInfo_tmp.isPaired = 1;
						ReadPairInfo_tmp.IsoFlag[i] = 1;
						ReadPairInfo_tmp.Size[i]    = abs(read2_mappos[eventid][i]
								-read1_mappos[eventid][i])
								+read_size;
						ReadPairInfo_tmp.Pos[2*i]    = read1_mappos[eventid][i];
						ReadPairInfo_tmp.Pos[2*i+1]    = read2_mappos[eventid][i];
					}
				}
				mappinginfo.insert (pair<string, ReadPairInfo>(eventid, ReadPairInfo_tmp));
			}
			else
			{
				ReadPairInfo_tmp.isPaired = 0;
				for(int i=0;i<2;i++)
				{
					ReadPairInfo_tmp.Pos[2*i]    = read1_mappos[eventid][i];
					ReadPairInfo_tmp.Size[i] = 0;
					if(read1_mappos[eventid][i] == -1)
					{
						ReadPairInfo_tmp.IsoFlag[i] = 0;
						
					}
					else
					{
						ReadPairInfo_tmp.IsoFlag[i] = 1;
					}
				}
				mappinginfo.insert (pair<string, ReadPairInfo>(eventid, ReadPairInfo_tmp));
			}
		}
		map<string, vector<int> >::iterator it2 = read2_mappos.begin();
		for(;it2!=read2_mappos.end();it2++)
		{
			string eventid = it2->first;
			if(mappinginfo.count(eventid)>0)
			{
				continue;
			}
			ReadPairInfo ReadPairInfo_tmp;
			ReadPairInfo_tmp.IsoFlag.resize(2);
			ReadPairInfo_tmp.Size.resize(2);
			//ReadPairInfo_tmp.Pos.resize(4);
			ReadPairInfo_tmp.Pos.assign(4,-1);
			if(read1_mappos.count(eventid) > 0)
			{
				ReadPairInfo_tmp.isPaired = 0;
				for(int i=0;i<2;i++)
				{
					if(read2_mappos[eventid][i] == -1 || read1_mappos[eventid][i] == -1)
					{
						ReadPairInfo_tmp.IsoFlag[i] = 0;
						ReadPairInfo_tmp.Size[i]    = 0;
					}
					else
					{
						ReadPairInfo_tmp.isPaired = 1;
						ReadPairInfo_tmp.IsoFlag[i] = 1;
						ReadPairInfo_tmp.Size[i]    = abs(read2_mappos[eventid][i]
								-read1_mappos[eventid][i]
								+read_size);
						ReadPairInfo_tmp.Pos[2*i] = read1_mappos[eventid][i];
						ReadPairInfo_tmp.Pos[2*i+1]=read2_mappos[eventid][i];
					}
				}
				mappinginfo.insert (pair<string, ReadPairInfo>(eventid, ReadPairInfo_tmp));
			}
			else
			{
				ReadPairInfo_tmp.isPaired = 0;
				for(int i=0;i<2;i++)
				{
					ReadPairInfo_tmp.Pos[2*i+1]=read2_mappos[eventid][i];
					ReadPairInfo_tmp.Size[i] = 0;
					if(read2_mappos[eventid][i] == -1)
					{
						ReadPairInfo_tmp.IsoFlag[i] = 0;
					}
					else
					{
						ReadPairInfo_tmp.IsoFlag[i] = 1;
					}
				}
				mappinginfo.insert (pair<string, ReadPairInfo>(eventid, ReadPairInfo_tmp));
			}

		}
	}
};

int main (int argc, char **argv) {
    int next_option;
    int digit_optind = 0;
    int read_size = 36;
    char *dbfile_name = 0,
    	 *IRM_file_name = 0,
	 *FZM_file_name = 0,
	 *output_file_prefix = 0,
	 *read1_nomt_file_name = 0,
	 *read2_nomt_file_name = 0;

    while ( (next_option = getopt(argc, argv, "s:d:b:f:o:1:2:")) != -1) {
        int this_option_optind = optind ? optind : 1;
        switch (next_option) {
	case 's':	sscanf (optarg,"%d",&read_size);	break;
        case 'd':	dbfile_name = optarg;	break;
        case 'b':	IRM_file_name = optarg;	break;
        case 'f':	FZM_file_name  = optarg;	break;
	case 'o':	output_file_prefix = optarg;	break;
        case '1':	read1_nomt_file_name = optarg;	break;
        case '2':	read2_nomt_file_name = optarg;	break;
	default:	break;
	}
      }
    if (!read_size||
	    !dbfile_name||
	    !FZM_file_name||
	    !output_file_prefix||
	    !read1_nomt_file_name||
	    !read2_nomt_file_name)
      {
	  printusage();
	  return 0;
      }
      clock_t clock_start, clock_finish;
      clock_start=clock();
      string ratio_file_name = output_file_prefix;
      ratio_file_name.append(".ratio");
      string log_file_name = output_file_prefix;
      log_file_name.append(".log");
      string num_file_name = output_file_prefix;
      num_file_name.append(".nums");
      
      ofstream log_file(log_file_name.c_str());
      cout<<"Log file name:\t"<<log_file_name<<endl;
      cout<<"Parameters:"<<endl;
      cout<<"Read size:\t"<<read_size<<endl;
      log_file<<"#read_size:\t"<<read_size<<endl;
      cout<<"Database:\t"<<dbfile_name<<endl;
      log_file<<"#dbfile_name:\t"<<dbfile_name<<endl;
      cout<<"Fragment size file name:\t"<<FZM_file_name<<endl;
      log_file<<"#FZM_file_name:\t"<<FZM_file_name<<endl;
      cout<<"Output file prefix:\t"<<output_file_prefix<<endl;
      cout<<"Read end1 file:\t"<<read1_nomt_file_name<<endl;
      cout<<"Read end2 file:\t"<<read2_nomt_file_name<<endl;
      log_file<<"#read1_nomt_file_name:\t"<<read1_nomt_file_name<<endl;
      log_file<<"#read2_nomt_file_name:\t"<<read2_nomt_file_name<<endl;
      if(IRM_file_name)
      {
	      cout<<"Inclusion ratio distn. file:\t"<<IRM_file_name<<endl;
	      log_file<<"#IRM_file_name:\t"<<IRM_file_name<<endl;
	      
      }
      else
      {
	      cout<<"Inclusion ratio distn. file:\tNA"<<endl;
	      log_file<<"#IRM_file_name:\tNA"<<endl;
      }

      //start reading files, dbfile first
      cout<<"Reading database file..."<<endl;
      map<string, TXdb_entry> TXdb_entries;
      ifstream dbfile(dbfile_name);
      char dbfile_line[MAX_LINE_LEN];
      while(dbfile.getline(dbfile_line,MAX_LINE_LEN)){
	      TXdb_entry entry_tmp(dbfile_line);
	      if(entry_tmp.iscomplete)
	      {
		TXdb_entries.insert(pair<string, TXdb_entry>
				(entry_tmp.id,entry_tmp)
				);
	      	//cout<<entry_tmp.exonsize[1]<<endl;
	      }
	      
      }
      dbfile.close();
      
      //start reading bayes prior if there is
      map<string, vector< double > >  IRM;
      IRM.insert(pair<string, vector<double> > ("CA", vector<double>()));
      IRM.insert(pair<string, vector<double> > ("CS", vector<double>()));
      IRM.insert(pair<string, vector<double> > ("IR", vector<double>()));
      IRM.insert(pair<string, vector<double> > ("AA", vector<double>()));
      IRM.insert(pair<string, vector<double> > ("AD", vector<double>()));
      //IRM.insert(pair<string, vector<double> > ("AI", vector<double>()));
      if(IRM_file_name)
      {
	cout<<"Reading inclusion ratio distn. file..."<<endl;
	ifstream IRM_file(IRM_file_name);
	char IRM_line[MAX_LINE_LEN];
	IRM_file.getline(IRM_line, MAX_LINE_LEN);
	while(IRM_file.getline(IRM_line, MAX_LINE_LEN))
	{
		char *token;
		token = strtok(IRM_line, "\t");
		if (token != NULL) {
			double hist_tmp;
			sscanf(token,"%lf",&hist_tmp);
			IRM["CA"].push_back(hist_tmp);
		}
		token = strtok(NULL,"\t");
		if (token != NULL) {
			double hist_tmp;
			sscanf(token,"%lf",&hist_tmp);
			IRM["AD"].push_back(hist_tmp);
		}
		token = strtok(NULL,"\t");
		if (token != NULL) {
			double hist_tmp;
			sscanf(token,"%lf",&hist_tmp);
			IRM["AA"].push_back(hist_tmp);
		}
/*		token = strtok(NULL,"\t");
		if (token != NULL) {
			double hist_tmp;
			sscanf(token,"%lf",&hist_tmp);
			IRM["AI"].push_back(hist_tmp);
		}*/
		token = strtok(NULL,"\t");
		if (token != NULL) {
			double hist_tmp;
			sscanf(token,"%lf",&hist_tmp);
			IRM["IR"].push_back(hist_tmp);
		}
		token = strtok(NULL,"\t");
		if (token != NULL) {
			double hist_tmp;
			sscanf(token,"%lf",&hist_tmp);
			IRM["CS"].push_back(hist_tmp);
		}

	}
	IRM_file.close();
      }
      
      //read FZM
      cout<<"Reading fragment size distn. file..."<<endl;
      ifstream FZM_file(FZM_file_name);
      vector<double> FZM;
      char FZM_file_line[MAX_LINE_LEN];
      FZM_file.getline(FZM_file_line,MAX_LINE_LEN);
      while(FZM_file.getline (FZM_file_line,MAX_LINE_LEN) )
      {
	     double hist_tmp; 
	     sscanf(FZM_file_line, "%lf", &hist_tmp) ;
	     FZM.push_back(hist_tmp);
      }
      FZM_file.close();
      //start to read reads
      cout<<"Reading reads from both ends...."<<endl;
      ifstream read1_nomt_file(read1_nomt_file_name);
      ifstream read2_nomt_file(read2_nomt_file_name);
      char read1_line[MAX_LINE_LEN];
      char read2_line[MAX_LINE_LEN];
      int total_single_read_num=0;
      int total_pair_read_num=0;
      int total_read_num=0;
      while(read1_nomt_file.getline(read1_line, MAX_LINE_LEN))
      {
	      total_read_num++;
	      ReadInfo read1(read1_line);
	      map<string, vector<int> > read1_mappos;
	      if(read1.isUseful)
	      {
	        read1.scan(read1_mappos);
	      }
	      read2_nomt_file.getline(read2_line, MAX_LINE_LEN);
	      ReadInfo read2(read2_line);
	      map<string, vector<int> > read2_mappos;
	      if(read2.isUseful)
	      {
		read2.scan(read2_mappos);
	      }
	      //pairing
	      ReadPair read_pair_ins (read1_mappos, read2_mappos, read_size);
	      //feed them into TXdb_entries
	      map<string , ReadPairInfo>::iterator it = read_pair_ins.mappinginfo.begin();
	      for(; it!=read_pair_ins.mappinginfo.end();it++)
	      {
		      string eventid = it->first;
		      map<string, TXdb_entry>::iterator txdb_it = TXdb_entries.find(eventid);
		      if(txdb_it == TXdb_entries.end() )
		      {
			      continue;
		      }
		      //get numbers
		      for(int i=0; i<4;i++)
		      {
			      if(read_pair_ins.mappinginfo[eventid].Pos[i]!=-1)
			      {
				      int pos_tmp1 = read_pair_ins.mappinginfo[eventid].Pos[i];
				      int pos_tmp2;
				      if(i<2)
				      { 
					      if (txdb_it->second.strand == '-')
						      pos_tmp1 = txdb_it->second.len1 -pos_tmp1 +1- read_size;
					      pos_tmp2 = pos_tmp1 + read_size;
					      if (pos_tmp2 <= txdb_it->second.exonsize[0])
					      {   txdb_it->second.enum1++;continue;}
					      if (pos_tmp1 <= txdb_it->second.exonsize[0] )
					      {
						     if(read_pair_ins.mappinginfo[eventid].IsoFlag[1] !=1)
					             {   
						      txdb_it->second.jnum12++;
					             }
					             else
					             {
						      txdb_it->second.enum1++;
					             }
						     continue;
					      }
					      if (pos_tmp2 <= txdb_it->second.exonsize[0]+txdb_it->second.exonsize[1])
                                              {   txdb_it->second.enum2++;continue;}
					      if (pos_tmp1 <= txdb_it->second.exonsize[0]+txdb_it->second.exonsize[1])
					      {
						     if(read_pair_ins.mappinginfo[eventid].IsoFlag[1] !=1)
					             {   
						      txdb_it->second.jnum23++;
					             }
					             else
					             {
						      txdb_it->second.enum2++;
					             }
						     continue;
					        
					      }
					      txdb_it->second.enum3++;
				      }
				      else
				      {
					      if (read_pair_ins.mappinginfo[eventid].IsoFlag[0] == 1 
						      && read_pair_ins.mappinginfo[eventid].IsoFlag[1] == 1)
						     continue; 
					      if (txdb_it->second.strand == '-')
						      pos_tmp1 = txdb_it->second.len2 -pos_tmp1 +1- read_size;
					      pos_tmp2 = pos_tmp1 + read_size;
					      if (pos_tmp2 <= txdb_it->second.exonsize[0])
					      {   txdb_it->second.enum1++;continue;}
					      if (pos_tmp1 <= txdb_it->second.exonsize[0] )
					      {
						     if(read_pair_ins.mappinginfo[eventid].IsoFlag[0] !=1)
					             {   
						      txdb_it->second.jnum13++;
					             }
					             else
					             {
						      txdb_it->second.enum1++;
					             }
						     continue;
					      }
					      txdb_it->second.enum3++;
						      
					      
				      }
			      }
		      }

		      if(read_pair_ins.mappinginfo[eventid].isPaired)
		      {
			      total_pair_read_num += 2;
			      txdb_it->second.totalreadnum += 2;
			      (txdb_it->second).pairs.push_back(read_pair_ins.mappinginfo[eventid]);
		      }
		      else
		      {
			      total_single_read_num += 2;
			      txdb_it->second.totalreadnum += 1;
			      if(read_pair_ins.mappinginfo[eventid].IsoFlag[0] == 1)
			      {
				      if(read_pair_ins.mappinginfo[eventid].IsoFlag[1] == 0)
				      {
			                (txdb_it->second).snum1++;
				      }
				      else
				      {
					(txdb_it->second).snum12++;
				      }
			      }
			      else
			      {
				      (txdb_it->second).snum2++;
			      }
		      }
	      }

	      
      }

      read1_nomt_file.close();
      read2_nomt_file.close();
      cout<<"Total read number:\t"<<total_read_num*2<<endl;
      cout<<"Singled reads used:\t"<<total_single_read_num<<endl;
      cout<<"Paired reads used:\t"<<total_pair_read_num<<endl;

      log_file<<"#total_read_num:\t"<<total_read_num*2<<endl;
      log_file<<"#total_single_read_num:\t"<<total_single_read_num<<endl;
      log_file<<"#total_pair_read_num:\t"<<total_pair_read_num<<endl;
      
      //Calculate ratios and output
      cout<<"Now doing estimations and output to "<<ratio_file_name<<"\t"<<num_file_name<<endl;
      ofstream ratio_file(ratio_file_name.c_str());
      ofstream num_file(num_file_name.c_str());
      map<string, TXdb_entry>::iterator txdb_it = TXdb_entries.begin();
      for(;txdb_it!=TXdb_entries.end();txdb_it++)
      {
	      string eventid = txdb_it->first;
	      //cout<< eventid <<endl; 
	      txdb_it->second.estimate(FZM,IRM,read_size);
/*	      ratio_file<< eventid<<"\t"
	      		<<setprecision(4)<<fixed
	                <<txdb_it->second.ir<<"\t"
			<<txdb_it->second.bir<<"\t"
			<<"NA\tNA\tNA\t"
			<<txdb_it->second.jnum12<<"\t"
			<<txdb_it->second.jnum23<<"\t"
			<<txdb_it->second.jnum13<<"\t"
			<<(1.0*read_size*( txdb_it->second.enum1+txdb_it->second.jnum12 +txdb_it->second.jnum13 )/txdb_it->second.exonsize[0] )<<"\t"
			<<(1.0*read_size*( txdb_it->second.enum2+txdb_it->second.jnum12 +txdb_it->second.jnum23 )/txdb_it->second.exonsize[1] )<<"\t"
			<<(1.0*read_size*( txdb_it->second.enum3+txdb_it->second.jnum23 +txdb_it->second.jnum13 )/txdb_it->second.exonsize[2] )<<"\t"
			<<txdb_it->second.exonsize[0]<<"\t"
			<<txdb_it->second.exonsize[1]<<"\t"
			<<txdb_it->second.exonsize[2]<<"\t"
		<<(1.0*read_size*txdb_it->second.totalreadnum/txdb_it->second.len1)<<endl;
		*/
	      ratio_file<<eventid<<"\t"
	      		<<setprecision(4)<<fixed
			<<txdb_it->second.ir<<"\t"  //1
			<<txdb_it->second.bir<<"\t" //2
			<<txdb_it->second.chrid<<"\t" //3
			<<txdb_it->second.start<<"," //
			<<txdb_it->second.start+txdb_it->second.exonsize[0]<<"\t" //4
			<<txdb_it->second.start+txdb_it->second.exonstarts[1]<<","
			<<txdb_it->second.start+txdb_it->second.exonstarts[1]+txdb_it->second.exonsize[1]<<"\t"//5
			<<txdb_it->second.start+txdb_it->second.exonstarts[2]<<","
			<<txdb_it->second.start+txdb_it->second.exonstarts[2]+txdb_it->second.exonsize[2]<<"\t"//6
			<<txdb_it->second.strand<<"\t"//7
			<<txdb_it->second.jnum12<<"\t"//8
			<<txdb_it->second.jnum23<<"\t"//9
			<<txdb_it->second.jnum13<<"\t"//10
			<<(1.0*read_size*( txdb_it->second.enum1+txdb_it->second.jnum12 +txdb_it->second.jnum13 )/txdb_it->second.exonsize[0] )<<"\t"//11
			<<(1.0*read_size*( txdb_it->second.enum2+txdb_it->second.jnum12 +txdb_it->second.jnum23 )/txdb_it->second.exonsize[1] )<<"\t"//12
			<<(1.0*read_size*( txdb_it->second.enum3+txdb_it->second.jnum23 +txdb_it->second.jnum13 )/txdb_it->second.exonsize[2] )<<"\t"//13`
			<<(1.0*read_size*txdb_it->second.totalreadnum/txdb_it->second.len1)<<"\t"//14
			<<txdb_it->second.exonsize[0]<<"\t"//15
			<<txdb_it->second.exonsize[1]<<"\t"//16
			<<txdb_it->second.exonsize[2]<<endl;//17

	      num_file<< eventid<<"\t"
			<<txdb_it->second.chrid<<"\t"
			<<txdb_it->second.start<<","
			<<txdb_it->second.start+txdb_it->second.exonsize[0]<<"\t"
			<<txdb_it->second.start+txdb_it->second.exonstarts[1]<<","
			<<txdb_it->second.start+txdb_it->second.exonstarts[1]+txdb_it->second.exonsize[1]<<"\t"
			<<txdb_it->second.start+txdb_it->second.exonstarts[2]<<","
			<<txdb_it->second.start+txdb_it->second.exonstarts[2]+txdb_it->second.exonsize[2]<<"\t"
			<<txdb_it->second.strand<<"\t"
                <<txdb_it->second.enum1<<"\t"
			<<txdb_it->second.enum2<<"\t"
			<<txdb_it->second.enum3<<"\t"
			<<txdb_it->second.jnum12<<"\t"
			<<txdb_it->second.jnum23<<"\t"
			<<txdb_it->second.jnum13<<endl;
      }
      log_file.close();
      ratio_file.close();
      num_file.close();
      clock_finish = clock();
      cout<<"Done! time used:"<<(double)(clock_finish-clock_start)/CLOCKS_PER_SEC<<" seconds"<<endl;

      
}