diff src/splicetrap.estimate.cpp @ 1:adc0f7765d85 draft

planemo upload

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/splicetrap.estimate.cpp	Thu Sep 07 15:06:58 2017 -0400
@@ -0,0 +1,854 @@
+//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;
+
+      
+}