trash.cc

Go to the documentation of this file.

/* FILE: trash_revcomp.cc
   Copyright Erik Arner and Martti Tammi, University of Uppsala,
   Dept. of Genetics and Pathology, 2000.

   This program uses the files produced by PHRED. The reads that don't
   have acceptable quality, will be deleted from the PHD and fasta files.
   An additional trash.mtas file is created and the deleted sequences
   are stored there.

   This program is a module in the TRAP program package.

   MOdified in 2004 to fit with the Trapper system.
*/

#include <iostream>
#include <cmath>
#include <cassert>
#include "trash.h"


using namespace std;

void TrashAlgo::start()
  // {{{ 
{



  TRASH_WINDOW_LEN = 35;
  TRASH_WINDOW_STEP = 3;
  MIN_NUM_ADJ_WIN_GOOD = 4;
  MIN_LEN_GOOD = 50;
  GOOD_PERCENTAGE = 0.89;
  MAX_VARIANCE_IN_QUALWIN = 0.0035;
  

  //Note: this assumes that quality values are 99 at the most!
  for( int i = 0; i < 100; i++ ) {
    logQual_to_quality.push_back( 1 - pow(10.0,double(-i)/10.0) );
    decode_qualval_square.push_back( logQual_to_quality[i] * logQual_to_quality[i] );
  }


   
   
  for( size_t i = 0; i < getMAl()->get_num_seq(); i++ ) {
    
    trash(i);
    
    //Temp test
//     for( size_t j = 0; j < getMAl()->get_seq_end( i ); j++ ) {
//       if ( j < getMAl()->get_beginGood( i ) ||  j >= getMAl()->get_endGood( i ) )
//         getMAl()->set_base( i, j, 'x' );
//     }
    
  }

    
}

// }}}
void TrashAlgo::trash(size_t seqNumber)
  // {{{ 

{
  /* This function computes the quality regions of the
     sequences, and marks those sequences with too poor
     quality to be deleted from further investigation.

  */


  int num_adj_win_good(0);
  int begin_good(-1), end_good(-1);
  int bg_found(0), eg_found(0);
  int  reset_good(1), curr_win_good(0);
    
  //Left end
  bool break_for(false);
  int start = 0;
  int stop = getMAl()->get_seq_end(seqNumber) - TRASH_WINDOW_LEN + 1;
 
  for( int j = start; j < stop; j += TRASH_WINDOW_STEP ) {
      
    if ( reset_good ) {
      begin_good = j;
    }
    reset_good = 0;
    
    //Adjust window size for gaps
    int num_bases_in_win(0);
    int num_gaps_in_win(0);
    while ( num_bases_in_win < TRASH_WINDOW_LEN ) {
      if ( j + TRASH_WINDOW_LEN + num_gaps_in_win > getMAl()->get_seq_end(seqNumber) ) {
        
        break_for = true;
        break;
      }
//       cerr<<getMAl()->get_base( seqNumber, j + num_bases_in_win + num_gaps_in_win );
      if ( getMAl()->get_base( seqNumber, j + num_bases_in_win + num_gaps_in_win ) == '-' ) {

        num_gaps_in_win++;
      }
      else {
        num_bases_in_win++;
      }
    }
//     cerr<<endl;
    if ( break_for ) break;

    j += check_qual_in_win(j, j + TRASH_WINDOW_LEN + num_gaps_in_win, seqNumber, curr_win_good);

    /* Window status is now computed and it may be
       possible to set begin_good */

    if ( curr_win_good ) {//good or hq
      num_adj_win_good++;
      
      if ( num_adj_win_good >= MIN_NUM_ADJ_WIN_GOOD ) {
        //Found begin_good
        bg_found = 1;
      }
    }
    else {
      if ( !bg_found ) {
        reset_good = 1;
        num_adj_win_good = 0;
      }
    }
      
    if ( bg_found ) {
      break;
    }
  }//End left

  /* Right end -- same thing, but backwards.
     Unnecessary if begin good not found. */

  if ( bg_found ) {
    reset_good = 1;
    curr_win_good = 0;
    num_adj_win_good = 0;

    break_for = false;
    
    start = getMAl()->get_seq_end(seqNumber) - 1;
    stop = TRASH_WINDOW_LEN - 2;
      
    for( int j = start; j > stop; j -= TRASH_WINDOW_STEP ) {
        
      if ( reset_good ) {
        end_good = j + 1;
      }

      reset_good = 0;

      //Adjust window size for gaps
      int num_bases_in_win(0);
      int num_gaps_in_win(0);
      while ( num_bases_in_win < TRASH_WINDOW_LEN ) {
//         cerr<<getMAl()->get_base( seqNumber, j - num_gaps_in_win - num_bases_in_win );
        if ( j - TRASH_WINDOW_LEN + 1 - num_gaps_in_win < 0 ) {
          
          break_for = true;
          break;
        }
        if ( getMAl()->get_base( seqNumber, j - num_gaps_in_win - num_bases_in_win ) == '-' ) {
          num_gaps_in_win++;
        }
        else {
          num_bases_in_win++;
        }
      }
//       cerr<<endl;
      if ( break_for ) break;
      
      
      j -= check_qual_in_win(j - TRASH_WINDOW_LEN + 1 - num_gaps_in_win, j + 1, seqNumber, curr_win_good);
        
      /* Window status is now computed and it may be
         possible to set end_good */
        
      if ( curr_win_good ) {
        num_adj_win_good++;
          
        if ( num_adj_win_good >= MIN_NUM_ADJ_WIN_GOOD ) {
          //Found begin_good
          eg_found = 1;
        }
      }
      else {
        if ( !eg_found ) {
          reset_good = 1;
          num_adj_win_good = 0;
        }
      }
        
      if ( eg_found ) {
        break;
      }
    }//End right forj
  }//End right
    
  //Set the parameters for the read
    
  int len = end_good - begin_good;
  if ( ( bg_found && eg_found ) && ( len >= MIN_LEN_GOOD ) ) {

    getMAl()->set_beginGood(seqNumber, begin_good);
    getMAl()->set_endGood(seqNumber, end_good);

  }
  else {

    kill_read(seqNumber);
  }

  //}//End fori loop over all seqs


}

// }}}
int TrashAlgo::check_qual_in_win(int start, int stop, int ID, int &curr_win_good)
  // {{{ 

{

//   cerr<<"stop: "<<stop<<endl;
//   cerr<<"start: "<<start<<endl;
  assert( stop <= getMAl()->get_seq_end( ID) );

  double sum_qual(0.0), sum_qual_square(0.0), average_qual(0.0);
  int minLogQual(0);
  int num_gap(0);

  //Check quality in window
  for( int k = start; k < stop; k++ ) {
//     cerr<<getMAl()->get_base( ID, k );
//     if ( stop >= getMAl()->get_seq_end( ID) ) {
//       assert( stop == getMAl()->get_seq_end( ID) );
//       //Should not happen...
// //       assert(false);
//       break;
//     }
    
    if ( getMAl()->get_base( ID, k ) == '-' ) {

      num_gap++;
      continue;
    }
    
    
    sum_qual += logQual_to_quality[ getMAl()->get_qual(ID, k) ];
    sum_qual_square += decode_qualval_square[ getMAl()->get_qual(ID, k) ];
    int tempQual = getMAl()->get_qual(ID, k);
    if(tempQual >= minLogQual)
      minLogQual = tempQual;
  }
//   cerr<<endl;
  average_qual = sum_qual / (double)TRASH_WINDOW_LEN;

  if ( average_qual >= GOOD_PERCENTAGE ) 
    {
      
      //Window may qualify as good and high quality
      double Q = sum_qual_square - sum_qual*sum_qual/(double)TRASH_WINDOW_LEN;
      assert(TRASH_WINDOW_LEN > 1);
      double variance = Q / (double)(TRASH_WINDOW_LEN - 1);
      if ( variance <= MAX_VARIANCE_IN_QUALWIN || minLogQual >= 27) {
        //Window qualifies as good
        curr_win_good = 1;
      }
      else {
        //Variance too high for good and hq
        curr_win_good = 0;
      }
      
    }
  else 
    {
      //Window does not qualify as good or high quality
      curr_win_good = 0;
    }

  if ( stop - start != TRASH_WINDOW_LEN + num_gap ) {
    cerr<<"stop - start: "<<stop - start<<endl;
    cerr<<"TRASH_WINDOW_LEN + num_gap: "<<TRASH_WINDOW_LEN + num_gap<<endl;
    cerr<<"stop: "<<stop<<endl;
    cerr<<"start: "<<start<<endl;
    cerr<<"TRASH_WINDOW_LEN: "<<TRASH_WINDOW_LEN<<endl;
    cerr<<"num_gap: "<<num_gap<<endl;
  }
//   cerr<<"num_gap: "<<num_gap<<endl;
  
  assert( stop - start == TRASH_WINDOW_LEN + num_gap );
  return num_gap;
}

// }}}
void TrashAlgo::kill_read(int ID)
  // {{{ 

{

  getMAl()->set_beginGood( ID, getMAl()->get_len( ID )/2 );
  getMAl()->set_endGood( ID, getMAl()->get_len( ID )/2 );
}

// }}}

---