mal_readonly.cc

Go to the documentation of this file.

#include "mal_readonly.h"
#include "generaldata.h"
#include "readdata.h"
#include "featuredata.h"
#include <sstream>
#include <cctype>

using namespace std;


MAl_Readonly::MAl_Readonly( size_t bufsize, set<db_recno_t>& recnolist, TrapperDoc* pdoc ) :
  buff_size(bufsize + 1),// +1 because buffID == 0 is unused
  num_seq( recnolist.size() ),
  doc(pdoc),
  selectedReads(recnolist)
{
  //Init the bookkeeping arrays, the rest will be automatically filled as the reads are used

  //These guys are of the buffer size
  seqs.reserve( buff_size );
  quals.reserve( buff_size );
  DNPs.reserve( buff_size );//NOTA BENE, vector<bool> is special case...
  tags.reserve( buff_size );
  time_course.reserve( buff_size );

  names.reserve( buff_size );
  headers.reserve( buff_size );
  mates.reserve( buff_size );
  strands.reserve( buff_size );
  

  //New deal - NOT always kept in memory...
  seq_rows.reserve( buff_size );
  seq_begin_global.reserve( buff_size );
  seq_end_global.reserve( buff_size );//Unnecessary?? We have the sizes of seqs...
  seq_beginGood.reserve( buff_size );//NB, not global!
  seq_endGood.reserve( buff_size );//NB, not global!
  mate_lengths.reserve( buff_size );

  //Buffer stuff
  ID_to_dbID.reserve( num_seq );//Should be of actual data set size
  ID_to_buffID.reserve( num_seq );//Ditto
  buffID_to_dbID.reserve( buff_size );//buffer size
  buffID_to_ID.reserve( buff_size );//buffer size
  put_in_db.reserve( buff_size );//buffer size
  
  for( set<db_recno_t>::iterator it = recnolist.begin(); it != recnolist.end(); ++it ) {
    ID_to_dbID.push_back( *it );
    ID_to_buffID.push_back( 0 );
  }

  //Init all containers
  for( size_t i = 0; i < buff_size; i++ ) {

    buffID_to_dbID.push_back( 0 );
    buffID_to_ID.push_back( 0 );
    put_in_db.push_back(false);

    //NOTA BENE, ONLY START, END AND ROW so far...
    
    seq_rows.push_back( 0 );//Need better initval...
    seq_begin_global.push_back( 0 );//Ditto
    seq_end_global.push_back( 0 );//Ditto

    seqs.push_back( vector<base_t>() );
    quals.push_back( vector<qual_t>() );
    time_course.push_back( vector<double>() );
    DNPs.push_back( vector<dnp_struct>() );
    tags.push_back( multimap<size_t, tag_struct>() );
    names.push_back( string() );
    mates.push_back( string() );
    strands.push_back( string() );
    

    seq_beginGood.push_back( 0 );
    seq_endGood.push_back( 0 );
    mate_lengths.push_back( 0 );
    
  }
  cerr<<"buff_size: "<<buff_size<<endl;
}

MAl_Readonly::~MAl_Readonly()
{
  //Nothing to be done here...
}


//Public methods - common

size_t MAl_Readonly::get_num_seq()
{
  return num_seq;
}

std::string MAl_Readonly::get_name( size_t ID )
{
  ID = get_buffID( ID );
  
  assert( ID < names.size() );
  
  return names[ ID ];
  
}

std::string MAl_Readonly::get_header( size_t ID )
{
  ID = get_buffID( ID );
  
  assert( ID < headers.size() );
  
  return headers[ ID ];
  
}

std::string MAl_Readonly::get_seq( size_t ID )
{
  ID = get_buffID( ID );
  
  assert( ID < seqs.size() );
  
//   return seqs[ ID ];
  return string(seqs[ ID ].begin(), seqs[ ID ].end());
  
}

string MAl_Readonly::get_strand( size_t ID )
{
  ID = get_buffID( ID );
  
  assert( ID < strands.size() );
  
  return strands[ ID ];
  
}

size_t MAl_Readonly::get_len( size_t ID)
{
  ID = get_buffID( ID );
  
  assert( ID < seqs.size() );
  
  return seqs[ ID ].size();
  
}

void MAl_Readonly::select_read( size_t ID, bool status )
{
  
  assert( ID < ID_to_dbID.size() );
  
  db_recno_t recno = ID_to_dbID[ ID ];
  if ( recno == 0 ) return;

  set<db_recno_t>::iterator pos;
  if ( status == true )
    selectedReads.insert( recno );
  else if ( ( pos = selectedReads.find( recno )) != selectedReads.end() ) {
    selectedReads.erase( pos );
  }
  
  
}


//Public methods - separate interfaces

size_t MAl_Readonly::get_seq_row( size_t ID )
{
  ID = get_buffID( ID );

  assert( ID < seq_rows.size() );
  
  return seq_rows[ ID ];
  
}


size_t MAl_Readonly::get_seq_begin_global( size_t ID )
// int MAl_Readonly::get_seq_begin_global( size_t ID )
{
  ID = get_buffID( ID );

  assert( ID < seq_begin_global.size() );
  
  return seq_begin_global[ ID ];
                           
}


size_t MAl_Readonly::get_seq_end( size_t ID )
{
  
  return get_seq_end_global( ID ) - get_seq_begin_global( ID );
}

size_t MAl_Readonly::get_seq_end_global( size_t ID )
{
  ID = get_buffID( ID );

  assert( ID < seq_end_global.size() );
  
  return seq_end_global[ ID ];
  
}


size_t MAl_Readonly::get_beginGood( size_t ID )
{
  ID = get_buffID( ID );

  assert( ID < seq_beginGood.size() );
  
  return seq_beginGood[ ID ];
}

size_t MAl_Readonly::get_beginGood_global( size_t ID )
{
  return get_beginGood( ID ) + get_seq_begin_global( ID );
  
}


size_t MAl_Readonly::get_endGood( size_t ID )
{
  ID = get_buffID( ID );

  assert( ID < seq_endGood.size() );
  
  return seq_endGood[ ID ];
  
}

size_t MAl_Readonly::get_endGood_global( size_t ID )
{
  return get_endGood( ID ) + get_seq_begin_global( ID );
}

base_t MAl_Readonly::get_base( size_t ID, size_t index )
{
  
  ID = get_buffID( ID );
  
  assert( ID < seqs.size() );

  if (index >= seqs[ ID ].size() ) {
    cout << flush;
    cerr << "index (local) = " << index << ", seqs[ ID ].size() = " << seqs[ ID ].size() << endl;

  }
  assert( index < seqs[ ID ].size() );
  
  base_t tmp = seqs[ ID ][ index ];
  if (tmp == '*') {
    tmp = '-';
  }
  
  return tmp;
//  return seqs[ ID ][ index ];
}

base_t MAl_Readonly::get_base_global( size_t ID, size_t index )
{
  return get_base( ID, index - get_seq_begin_global( ID ) );
  
}

qual_t MAl_Readonly::get_qual( size_t ID, size_t index )
{
  ID = get_buffID( ID );
  
  assert( ID < quals.size() );
  assert( index < quals[ ID ].size() );
  
  return quals[ ID ][ index ];
  
}

qual_t MAl_Readonly::get_qual_global( size_t ID, size_t index )
{
  
  return get_qual( ID, index - get_seq_begin_global( ID ) );
  
}


bool MAl_Readonly::is_DNP(size_t ID, size_t index)
{
  ID = get_buffID( ID );

  if ( index >= DNPs[ ID ].size() ) {
    
  }

  assert( ID < DNPs.size() );
  assert( index < DNPs[ ID ].size() );
  
  return DNPs[ ID ][ index ].isDNP;
}

bool MAl_Readonly::is_DNP_global(size_t ID, size_t index)
{

  return is_DNP( ID, index - get_seq_begin_global( ID ) );
}

int MAl_Readonly::get_DNP_ID(size_t ID, size_t index)
{
  ID = get_buffID( ID );

  assert( ID < DNPs.size() );
  assert( index < DNPs[ ID ].size() );
  
  return DNPs[ ID ][ index ].ID;  
}

int MAl_Readonly::get_DNP_ID_global(size_t ID, size_t index)
{
  return get_DNP_ID( ID, index - get_seq_begin_global( ID ) ); 
}

int MAl_Readonly::get_DNP_type(size_t ID, size_t index)
{
  ID = get_buffID( ID );

  assert( ID < DNPs.size() );
  assert( index < DNPs[ ID ].size() );
  
  return DNPs[ ID ][ index ].type;
}

int MAl_Readonly::get_DNP_type_global(size_t ID, size_t index)
{
  return get_DNP_type( ID, index - get_seq_begin_global( ID ) );  
}


//Private methods

size_t MAl_Readonly::get_buffID(size_t ID)
{
  

  assert( ID < ID_to_buffID.size() );
  size_t buffID = ID_to_buffID[ ID ];

  if ( buffID == 0 ) {
    //Seq is not present in buffer and thus has to be read into buffer,
    //OR a new seq is being created
    buffID = next_buffID();

    if ( buffID >= put_in_db.size() ) {
      cerr<<"buffID: "<<buffID<<endl;
      cerr<<"put_in_db.size(): "<<put_in_db.size()<<endl;
    }
    
    assert( buffID < put_in_db.size() );
    
    if ( put_in_db[ buffID ] ) {
      
      flush_buffer( buffID );
    }
    
    //If dbID is 0, a new seq is being created and nothing has to be read from DB
    //This is taken care of in read_from_db()

    read_from_db( buffID, ID );

    
  }
  assert( buffID != 0 );
  return buffID;

}

size_t MAl_Readonly::next_buffID()
{
//   static size_t curr_buffID(static_cast<size_t>(-1));//Have to check this, but should work although size_t is unsigned
  static size_t curr_buffID(0);//ID == 0 is unused...
  ++curr_buffID;
  
  if ( curr_buffID == buff_size ) curr_buffID = 1;
  
  return curr_buffID;
  
}

void MAl_Readonly::flush_buffer( size_t buffID )
{
  //Not implemented in read only MAl
  
}

void MAl_Readonly::read_from_db(size_t buffID, size_t ID)
{

  db_recno_t recno = ID_to_dbID[ ID ];

  if ( recno != 0 ) {

    //NB: EXTREMELY important that DNP features and other features
    //that appear more than once on reads are read AFTER DnaStrData!!!

    read_seq_from_db( recno, buffID );
    read_feat_from_db( recno, buffID, "ReadMetaData");//Metadata
    read_feat_from_db( recno, buffID, "DnaStrData");//The seq
    read_feat_from_db( recno, buffID, "QualityData");//The qual
    read_feat_from_db( recno, buffID, "DnpData");//The dnps
    read_feat_from_db( recno, buffID, "TagData");//General tags
    
    
    /*
      read_feat_from_db( recno, buffID, "beginGood");//Maybe this should be added as data in the seq class...
      read_feat_from_db( recno, buffID, "endGood");//Ditto
      read_feat_from_db( recno, buffID, "name");//Ditto
      read_feat_from_db( recno, buffID, "header");//Ditto
      read_feat_from_db( recno, buffID, "strand");//Ditto
      
    */
    
  }
  

  ID_to_buffID[ ID ] = buffID;
  buffID_to_dbID[ buffID ] = recno;
  buffID_to_ID[ buffID ] = ID;
  put_in_db[ buffID ] = true;
}

void MAl_Readonly::read_seq_from_db( db_recno_t recno, size_t buffID )
{

  Database::PrimaryIterator<ReadData> read_it(doc, "ReadData");

  int ret_read = read_it.setFromRecno(recno);
  ReadData* r_test = (ret_read != DB_NOTFOUND) ? read_it.answer() : 0;
  assert( r_test );

  seq_begin_global[ buffID ] = r_test->startPos();
  seq_end_global[ buffID ] = r_test->endPos() + 1;//Index annoyance...
  seq_rows[ buffID ] = r_test->row();
//   names[ buffID ] = r_test->name();
//   mates[ buffID ] = r_test->mate();
//   strands[ buffID ] = r_test->strand();
//   seq_beginGood[ buffID ] = r_test->beginGood();
//   seq_endGood[ buffID ] = r_test->endGood() + 1;
//   mate_lengths[ buffID ] = r_test->mateLength();
//   time_course[ buffID ] = r_test->tc_vec.stlVector();
  
//   cerr<<"Reading seq: "<<endl;
//   cerr<<"recno: "<<recno<<endl;
//   cerr<<"r_test->startPos(): "<<r_test->startPos()<<endl;
//   cerr<<"r_test->endPos(): "<<r_test->endPos()<<endl;

}

void MAl_Readonly::read_feat_from_db( db_recno_t recno, size_t buffID, const string& data_type_name)
{
  Database::SecondaryIterator<FeatureData> feat_it("readRecno", doc, data_type_name);
  
  feat_it.key()->setReadRecno( recno );

  if ( feat_it.set() != 0 ) {

    return;
  }

  
  if ( DnaStrData* data = dynamic_cast<DnaStrData*>(feat_it.answer()) ) {
    seqs[buffID] = data->dnaVector.stlVector();
    //Have to init DNPs here to assert that DNP vec is of same size...
    DNPs[buffID] = vector<dnp_struct>(seqs[buffID].size(), dnp_struct());
  }
  else if ( QualityData* data = dynamic_cast<QualityData*>(feat_it.answer()) ) {
    quals[buffID] = data->qualityVector.stlVector();
  }
  else if ( DnpData* data = dynamic_cast<DnpData*>(feat_it.answer()) ) {
    
    do {

      DNPs[buffID][ data->startPos() ].isDNP = true;
      DNPs[buffID][ data->startPos() ].recno = data->getRecno();
      DNPs[buffID][ data->startPos() ].ID = data->get_dnpID();      
      DNPs[buffID][ data->startPos() ].type = data->get_dnp_type();      

    } while ( (feat_it.nextdup() == 0) && (data = dynamic_cast<DnpData*>(feat_it.answer())) );
  }
  else if ( TagData* data = dynamic_cast<TagData*>(feat_it.answer()) ) {
    do {
      ostringstream os;
      os << data->getScore() <<" "<< data->getInfo();
      tags[buffID].insert( pair<size_t, tag_struct>( data->endPos() + 1, tag_struct(data->startPos(), data->endPos() + 1, data_type_name, os.str(), data->getRecno()) ) );
      
    } while ( (feat_it.nextdup() == 0) && (data = dynamic_cast<TagData*>(feat_it.answer())) );
  }
  else if (ReadMetaData* data = dynamic_cast<ReadMetaData*>(feat_it.answer()) ) {
    
    names[ buffID ] = data->name();
    mates[ buffID ] = data->mate();
    strands[ buffID ] = data->strand();
    seq_beginGood[ buffID ] = data->beginGood();
    seq_endGood[ buffID ] = data->endGood() + 1;
    mate_lengths[ buffID ] = data->mateLength();
    time_course[ buffID ] = data->tc_vec.stlVector();
    
  }
  
}




base_t MAl_Readonly::comp_base(const char& base)
{
  char tmp = toupper(base);
  if ( tmp == 'A') {
    return base + 19;
  }
  else if ( tmp == 'T' ) {
    return base - 19;    
  }
  else if ( tmp == 'G' ) {
    return base - 4;    
  }
  else if ( tmp == 'C' ) {
    return base + 4;    
  }
  else {
    return base;
  }
}

double MAl_Readonly::get_max_expression(size_t ID)
{
  ID = get_buffID( ID );

  assert( ID < time_course.size() );
  
  double tc_max(0.000001);
  bool maxfound(false);
//   cerr<<names[ ID ]<<endl;
  
  for( size_t i = 0; i < time_course[ID].size(); i++ ) {
    if ( time_course[ID][i] > tc_max ) {
      tc_max = time_course[ID][i];
      maxfound = true;
    }
  }
  if ( maxfound ) {
//     cerr<<"yo, tc_max: "<<tc_max<<endl;
    return tc_max;
  }

  return 0.0;
}

size_t MAl_Readonly::get_num_expression_points(size_t ID)
{
  ID = get_buffID( ID );

  assert( ID < time_course.size() );

  return time_course[ID].size();
}

double MAl_Readonly::get_expression(size_t ID, size_t expression_index)
{
  ID = get_buffID( ID );

  assert( ID < time_course.size() );
  assert( expression_index < time_course[ID].size() );

  return time_course[ID][expression_index];
}

---