trapperDNP.cc

Go to the documentation of this file.

#include "MAl_addon.h"
#include "trapperDNP.h"
#include "dnp_types.h"
#include <cassert>

#include <cmath>

DNPclass::DNPclass( MAl_addon* mal, int minqual, int d_min, double p_max ) : 
  m(mal), num_candidates(0), qualmin(minqual), Dmin(d_min), p_tot_max(p_max)
  // {{{ 

{
  
  // Fill in the 'quality2error vector'

  for( int q = 0; q < 100; q++ ) {
    
    quality2error.push_back( std::pow(10.0, -(double)q/10.0)  );
    
  }
  base_trans['a'] = 'a';
  base_trans['A'] = 'a';

  base_trans['t'] = 't';
  base_trans['T'] = 't';
  
  base_trans['g'] = 'g';
  base_trans['G'] = 'g';

  base_trans['c'] = 'c';
  base_trans['C'] = 'c';
  
  base_trans['n'] = 'n';
  base_trans['N'] = 'n';

  base_trans['x'] = 'x';
  base_trans['X'] = 'x';
  
  base_trans['-'] = '-';

  base_trans['*'] = '*';

  cerr<<"Finding candidate positions...  ";
  find_DNP_candidate_positions();
  cerr<<"Done"<<endl;
  cerr<<"Evaluating candidates... ";
  evaluate_DNP_candidate_positions();
  cerr<<"Done"<<endl;
  cerr<<num_found<<" DNPs found at "<<num_pos<<" positions"<<endl;
}

// }}}
DNPclass::~DNPclass()
  // {{{ 

{

}

// }}}
//**************************************************************
void DNPclass::find_DNP_candidate_positions()
  // {{{ 

{
  
  
  int min = m->get_min_col();
  int max = m->get_max_col();
  
  
  
  for( int i = min; i < max; i++ ) {
    if ( column_is_candidate(i) ) {
      cand_vec.push_back(i);
    }
    int proc = 100 * i / max;
//     if ( proc % 5 == 0 )
//       cerr<<proc<<" % checked"<<endl;
  }
  assert( cand_vec.size() == cand_bases.size() );
  num_candidates += (int)cand_vec.size();
}

// }}}
int DNPclass::column_is_candidate(int col)
  // {{{ 

{
  /* 
     Definition of candidate base type:
     
     ¤ Must be present in at least Dmin reads
     ¤ Can be present in at most half the number of reads on column
     
     ¤ Extra: at present there can be only one candidate base type 
     on a column. Without this constraint, a more involved approach 
     has to be taken in order to identify candidate base types.
     
  */
 
  char cons = base_trans [ m->get_most_freq_base_on_col( col ) ];
  if ( cons != 'a' && cons != 'c' && 
       cons != 'g' && cons != 't'  ) return 0;
  


  char base_arr[] = {'a', 'c', 'g', 't'/*, '-'*/};
  int cand_num(0);
  int cand_base = 'u';
  int cand_found(0);
  
  for( int i = 0; i < 4; i++ ) {
    char base = base_trans [ base_arr[i] ];
    if ( base == cons /*|| base != m->get_entry( 0, col )*/ /*NOTA BENE*/) continue;
    
    if ( m->get_num_bases_on_col(col, base) >= Dmin ) {
      if ( m->get_num_bases_on_col(col, base) > cand_num ) {
        cand_num = m->get_num_bases_on_col(col, base);
        cand_base = base;
        cand_found = 1;
      }
    }
  }

  
  if ( cand_found ) {
    cand_bases.push_back( cand_base );
    cons_bases.push_back( cons );
    return 1;
  }
  return 0;
  
}

// }}}

void DNPclass::evaluate_DNP_candidate_positions()

  // {{{ 

{
  /*
    Tentative try: compare all candidate columns to each other
  
  */
  num_found = 0;
  num_pos = 0;
  
  assert( (int)confirmed_DNP_pos.size() == 0 );
  for( int i = 0; i < m->get_max_col(); i++ ) {
    confirmed_DNP_pos.push_back(0);
  }

  vector<char> first_base_col, second_base_col;
  vector<int> first_qual_col, second_qual_col;
  vector<int> rownum_vec;

  
  for( int i = 0; i < (int)cand_vec.size() - 1; i++ ) {
    int proc = 100 * i / cand_vec.size();
//     if ( proc % 5 == 0 )
//       cerr<<proc<<" % checked"<<endl;
    for( int j = i+1; j < (int)cand_vec.size(); j++ ) {
      
      pick_base_row_and_qual_columns(cand_vec[i], cand_vec[j], first_base_col, 
                                     second_base_col, first_qual_col, second_qual_col, 
                                     rownum_vec);
      
      if ( (int)first_base_col.size() < Dmin*2 ) {//Need sufficient number of bases on columns
        //cerr<<"Size too small! i: "<<i<<", j: "<<j<<endl;
        //cerr<<"cand_vec[i]: "<<cand_vec[i]<<", cand_vec[j]: "<<cand_vec[j]<<endl;
        break;
      }
      char cand_base_1 = cand_bases[i];
      char cand_base_2 = cand_bases[j];
      
      int num_cand_1(0), num_cand_2(0);
      double exp_err_first(0.0), exp_err_second(0.0);
      int num_corr(0);
      double expectation_val(0.0);
      look_for_corr(first_base_col, second_base_col, first_qual_col, 
                    second_qual_col, cand_base_1, cand_base_2, 
                    num_cand_1, num_cand_2, exp_err_first, exp_err_second,
                    num_corr, expectation_val);
      
//       double kvot;
//       if ( num_cand_1 > num_cand_2 ) {
//         kvot = double(num_cand_1)/double(num_cand_2);
//       }
//       else {
//         kvot = double(num_cand_2)/double(num_cand_1);
//       }
      
      if ( num_cand_1 < Dmin ) {
        //cerr<<"Too few candidates! i: "<<i<<", j: "<<j<<endl;
        break;
      }
      int DNPfound(0);
      double tot_prob=0;
      double col_prob=0;
      double p1(-1), p2(-1), corr_prob(-1);

      if ( num_corr >= Dmin ) { //NOTA BENE! Assumes Dmin = 3...
        //Computation of "the probability of observing any of the two columns by chance".
        
        p1 = poisson_probability(exp_err_first, num_cand_1);
        p2 = poisson_probability(exp_err_second, num_cand_1);
        col_prob = 1.0 - (1.0 - p1) * (1.0 - p2);
        
        assert( col_prob >= 0.0 && col_prob <= 1.0);//A bit shaky...
        
        corr_prob = poisson_probability(expectation_val, num_corr);
        tot_prob = col_prob * corr_prob;
        

        if (tot_prob <= p_tot_max || num_corr > Dmin )DNPfound = 1;
       
      }
      
      
      if ( DNPfound ) {

        if ( !confirmed_DNP_pos[ cand_vec[i] ] || 
             !confirmed_DNP_pos[ cand_vec[j] ]) {
          num_pos++;
        }
        

        assert( cand_vec[i] < (int)confirmed_DNP_pos.size() );
        assert( cand_vec[j] < (int)confirmed_DNP_pos.size() );

        confirmed_DNP_pos[ cand_vec[i] ] = 1;
        confirmed_DNP_pos[ cand_vec[j] ] = 1;
        set_DNPs( cand_vec[i], cand_vec[j], 
                  cand_bases[i], cand_bases[j], rownum_vec, 
                  i, j, cons_bases[i], cons_bases[j]);
      }

    }
  }
}

// }}}


void DNPclass::set_DNPs( int col1, int col2, char DNPbase1, char DNPbase2, 
                         vector<int> &rownum_vec, int ID1, int ID2, 
                         char cons1, char cons2)
// {{{ 

{

  for( int i = 0; i < (int)rownum_vec.size(); i++ ) {
    
    if ( m->get_entry( rownum_vec[i], col1 ) == DNPbase1 &&
         m->get_entry( rownum_vec[i], col2 ) == DNPbase2 &&
         m->get_qual(rownum_vec[i], col1) >= qualmin &&
         m->get_qual(rownum_vec[i], col2) >= qualmin) {
      //Set DNP
      if ( !m->is_DNP( rownum_vec[i] , col1 ) ) {
        num_found++;
      }
      if ( !m->is_DNP( rownum_vec[i] , col2 ) ) {
        num_found++;
      }
      

      m->set_DNP( rownum_vec[i] , col1, true, ID1, dnp_type(DNPbase1, cons1) );
      m->set_DNP( rownum_vec[i] , col2, true, ID2, dnp_type(DNPbase2, cons2));

      
    }
  }
}

// }}}

void DNPclass::look_for_corr( vector<char> &first_col, vector<char> &second_col, 
                              vector<int> &first_qual_col, vector<int> &second_qual_col, 
                              char cand_base_1, char cand_base_2, int &num_err_first, int &num_err_second, 
                              double &exp_err_first, double &exp_err_second, int &num_corr, 
                              double &expectation_val)
  // {{{ 

{
  //NOTA BENE! This function currently performs its analysis
  //with regard to the candidate base. That is, every base
  //present on the column is regarded as correct if it is of
  //another type than the candidate base.

  //The two results are stored in expectation_val and returned
  //in col_prob, respectively.

  num_err_first = num_err_second = 0;
  exp_err_first = exp_err_second = 0.0;
  expectation_val = 0;
  num_corr = 0;
  
  //Pre-computation of necessary params

  for( int i = 0; i < (int)first_col.size(); i++ ) {
    assert( first_qual_col[i] > -1 && first_qual_col[i] < 100 );
    assert( second_qual_col[i] > -1 && second_qual_col[i] < 100 );
    
    int no_1_cand(0), no_2_cand(0);

    exp_err_first += quality2error[ first_qual_col[i] ];
    exp_err_second += quality2error[ second_qual_col[i] ];
    
    
    if ( first_col[i] == cand_base_1 && first_qual_col[i] >= qualmin ) {
      num_err_first++;
      no_1_cand = 1;
    }
    if ( second_col[i] == cand_base_2 && second_qual_col[i] >= qualmin ) {
      num_err_second++;
      no_2_cand = 1;
    }
    if ( no_1_cand && no_2_cand )
      num_corr++;
  }
  

  //Computation of the expectation value for the number of coinciding deviations
  
  for( int i = 0; i < (int)first_col.size(); i++ ) {
    double err1 = quality2error[ first_qual_col[i] ];
    double err2 = quality2error[ second_qual_col[i] ];
    double new_sum = exp_err_first - err1;
    double Exp1 = num_err_first * err1 / ( num_err_first * err1 + new_sum * (1.0 - err1) );
    new_sum = exp_err_second - err2;
    double Exp2 = num_err_second * err2 / ( num_err_second * err2 + new_sum * (1.0 - err2) );
    
    expectation_val += Exp1 * Exp2;
  }

}

// }}}
double DNPclass::poisson_probability(double lambda, int obs)
  // {{{ 

{
  //   if ( lambda < 0.0 ) {
  //     cerr<<"lambda: "<<lambda<<endl;
  //   }
  assert( lambda >= 0.0);
  assert( obs >= 0 );
  double M = lambda;
  int N_in = obs;
  int N = N_in - 1;
  int J       (0);
  double S = exp(-M);
  double res(0.0);
  
  if(N < 0){
    return 1.0;
  }
  

  if(S==0){
    J = N;
    if(J > M){
      J = int(M);
    }

    double F(0.0);
    for( int i = 1; i <= J; i++ ) {
      F += log((double)i);
    }
    double L = -M-F+J*log(M);
    L = exp(L);
    F=1;

    double T(0.0);

    for( int i = 1; i <= J; i++ ) {
      F *= (J+1-i)/M;
      T += F;
    }
    
    if(J==N){
      T=(T+1)*L;
      res = 1.0 - T;

    }
    else{
      F=1;
      for( int i = 1; i <= N-J; i++ ) {
        F *= M/(i+J);
        T=T+F;
      }
      T=(T+1)*L;
      res = (1.0 - T);
    }

  }
  else{
    
    double T=S;
    for( int i = 0; i <= N-1; i++ ) {
      S = S*M/(i+1);
      T += S;
    }
    res = (1.0 - T);
  }

  if ( res < 0.0 ) {
    if ( -res < 0.0000000000001 )
      res = 0.0;
  }
  else if ( res > 1.0 ) {
    if ( res - 1.0 < 0.0000000000001 )
      res = 1.0;
  }
  assert( res <= 1.0 );
  assert( res >= 0.0 );

  return res;
}

// }}}
void DNPclass::pick_base_row_and_qual_columns(int firstDNPindex, int secondDNPindex, 
                                              vector<char> &firstDNPvec, vector<char> &secondDNPvec,
                                              vector<int> &first_qual_vec, vector<int> &second_qual_vec,
                                              vector<int> &rownum_vec) 
  // {{{ 

{
  firstDNPvec.erase( firstDNPvec.begin(), firstDNPvec.end() );
  secondDNPvec.erase( secondDNPvec.begin(), secondDNPvec.end() );
  first_qual_vec.erase( first_qual_vec.begin(), first_qual_vec.end() );
  second_qual_vec.erase( second_qual_vec.begin(), second_qual_vec.end() );
  rownum_vec.erase( rownum_vec.begin(), rownum_vec.end() );

  for( int j = 0; j < m->get_num_rows(); j++ ) {
    
    if ( firstDNPindex >= m->get_row_end(j) || firstDNPindex < m->get_row_start(j) )
      continue;
    if ( secondDNPindex >= m->get_row_end(j) || secondDNPindex < m->get_row_start(j) )
      continue;

    //Both bases come from the same read

    char baseOnPos1 = base_trans[ m->get_entry(j, firstDNPindex) ];      
    char baseOnPos2 = base_trans[ m->get_entry(j, secondDNPindex) ];

    //EK 030408
    if(baseOnPos1 == 'x' || baseOnPos2 == 'x')
      continue;
    
    assert(baseOnPos1 == 'a' || baseOnPos1 == 't' || baseOnPos1 == 'g' || 
           baseOnPos1 == 'c' || baseOnPos1 == 'n' || baseOnPos1 == '-' ||
           baseOnPos1 == 'x' || baseOnPos1 == '*' );
    assert(baseOnPos2 == 'a' || baseOnPos2 == 't' || baseOnPos2 == 'g' || 
           baseOnPos2 == 'c' || baseOnPos2 == 'n' || baseOnPos2 == '-' ||
           baseOnPos2 == 'x' || baseOnPos2 == '*');
    
    int qualOnPos1 = m->get_qual(j, firstDNPindex);      
    int qualOnPos2 = m->get_qual(j, secondDNPindex);
    
    firstDNPvec.push_back(baseOnPos1);
    secondDNPvec.push_back(baseOnPos2);

    first_qual_vec.push_back(qualOnPos1);
    second_qual_vec.push_back(qualOnPos2);

    assert( j < m->get_num_rows() );

    rownum_vec.push_back( j );
  }
}

// }}}

---