c++/EEDB/Chrom.cpp
/**** get_subsequence
  Description  : uses ChromChunk objects and the sequence in the database to 
                 return the actual sequence in this region. Since the
                 Chrom is assigned to a specific Assembly one only needs to specify
                 the start/end and an optional strand to fetch the sequence.
  Parameter[1] : chrom_start 
                 the chromosome start of the region to fetch
  Parameter[2] : chrom_end  
                 the chromosome end of the region to fetch
  Parameter[3] : strand <optional> as "-" or "+"
                 the strand of the sequence. if "-" then it will return the sequence 
                 on the reverse strand by reverse complementing the sequence
  Returntype   : Bio::Seq instance or undef if a data error happens
  Errors       : if the region is not valid or if data is not present it will return undef
  Exceptions   : none
***/
get_subsequence
Description  : uses ChromChunk objects and the sequence in the database to
return the actual sequence in this region. Since the
Chrom is assigned to a specific Assembly one only needs to specify
the start/end and an optional strand to fetch the sequence.
Parameter[1] : chrom_start
the chromosome start of the region to fetch
Parameter[2] : chrom_end
the chromosome end of the region to fetch
Parameter[3] : strand  as "-" or "+"
the strand of the sequence.
get_subsequence Description : uses ChromChunk objects and the sequence in the database to return the actual sequence in this region . Since the Chrom is assigned to a specific Assembly one only needs to specify the start / end and an optional strand to fetch the sequence . Parameter [ 1 ] : chrom_start the chromosome start of the region to fetch Parameter [ 2 ] : chrom_end the chromosome end of the region to fetch Parameter [ 3 ] : strand as " - " or " + " the strand of the sequence .
string  EEDB::Chrom::get_subsequence(long int chrom_start, long int chrom_end, string strand) {
  string sequence;
  if(assembly() == NULL) { return ""; }
  if(!(assembly()->sequence_loaded())) { return ""; }
  if(chrom_start < 0) { chrom_start=1; }                 
  if(chrom_end < 0)   { chrom_end = _chrom_length -1; }  
  if(chrom_start > chrom_end) { 
    long int t1 = chrom_start;
    chrom_start = chrom_end;
    chrom_end = t1;;
  }
  string asm_name =  assembly()->assembly_name();
  vector<DBObject*> chunks;
  if(_zdxstream!=NULL) {
    _zdxstream->stream_by_named_region(asm_name, _chrom_name, chrom_start, chrom_end);
    while(MQDB::DBObject *obj = _zdxstream->next_in_stream()) {
      if(obj->classname() == EEDB::ChromChunk::class_name) {
        chunks.push_back(obj);
      } else {
        obj->release();
      }
    }
  }
  else if(_database!=NULL) {
    chunks = EEDB::ChromChunk::fetch_all_by_named_region(_database, asm_name, _chrom_name, chrom_start, chrom_end);
  }
  if(chunks.empty()) { return ""; }
  for(unsigned int i=0; i<chunks.size(); i++) {
    EEDB::ChromChunk *chunk = (EEDB::ChromChunk*)chunks[i];
    if((chrom_start >= chunk->chrom_start()) and (chrom_end <= chunk->chrom_end())) {
      sequence = chunk->get_subsequence(chrom_start, chrom_end);
      if(strand == "-") {
        std::reverse(sequence.begin(),sequence.end());
        string::iterator it1;
        string rosette1 ="acgtrymkswhbvdnxACGTRYMKSWHBVDNX";
        string rosette2 ="tgcayrkmswdvbhnxTGCAYRKMSWDVBHNX";
        for(it1=sequence.begin(); it1!=sequence.end(); it1++) {
          size_t found = rosette1.find(*it1);
          if (found!=std::string::npos) {
            *it1 = rosette2[found];
          }
        }
      }
      return sequence;
    }
  }
  long int ts = chrom_start; 
  for(unsigned int i=0; i<chunks.size(); i++) {
    EEDB::ChromChunk *chunk = (EEDB::ChromChunk*)chunks[i];
    long int te = chunk->chrom_end();
    if(chrom_end < te) { te = chrom_end; }  
    sequence += chunk->get_subsequence(ts, te);
    ts = te+1;
    if(ts > chrom_end) { break; } 
  }
  if(strand == "-") {
  }
  return sequence;
}
string EEDB :: Chrom :: get_subsequence ( long int chrom_start , long int chrom_end , string strand ) { string sequence ; if ( assembly ( ) == NULL ) { return " " ; } if ( ! ( assembly ( ) -> sequence_loaded ( ) ) ) { return " " ; } if ( chrom_start < 0 ) { chrom_start = 1 ; } if ( chrom_end < 0 ) { chrom_end = _chrom_length - 1 ; } if ( chrom_start > chrom_end ) { long int t1 = chrom_start ; chrom_start = chrom_end ; chrom_end = t1 ; ; } string asm_name = assembly ( ) -> assembly_name ( ) ; vector < DBObject * > chunks ; if ( _zdxstream != NULL ) { _zdxstream -> stream_by_named_region ( asm_name , _chrom_name , chrom_start , chrom_end ) ; while ( MQDB :: DBObject * obj = _zdxstream -> next_in_stream ( ) ) { if ( obj -> classname ( ) == EEDB :: ChromChunk :: class_name ) { chunks . push_back ( obj ) ; } else { obj -> release ( ) ; } } } else if ( _database != NULL ) { chunks = EEDB :: ChromChunk :: fetch_all_by_named_region ( _database , asm_name , _chrom_name , chrom_start , chrom_end ) ; } if ( chunks . empty ( ) ) { return " " ; } for ( unsigned int i = 0 ; i < chunks . size ( ) ; i ++ ) { EEDB :: ChromChunk * chunk = ( EEDB :: ChromChunk * ) chunks [ i ] ; if ( ( chrom_start >= chunk -> chrom_start ( ) ) and ( chrom_end <= chunk -> chrom_end ( ) ) ) { sequence = chunk -> get_subsequence ( chrom_start , chrom_end ) ; if ( strand == " " ) { std :: reverse ( sequence . begin ( ) , sequence . end ( ) ) ; string :: iterator it1 ; string rosette1 = " " ; string rosette2 = " " ; for ( it1 = sequence . begin ( ) ; it1 != sequence . end ( ) ; it1 ++ ) { size_t found = rosette1 . find ( * it1 ) ; if ( found != std :: string :: npos ) { * it1 = rosette2 [ found ] ; } } } return sequence ; } } long int ts = chrom_start ; for ( unsigned int i = 0 ; i < chunks . size ( ) ; i ++ ) { EEDB :: ChromChunk * chunk = ( EEDB :: ChromChunk * ) chunks [ i ] ; long int te = chunk -> chrom_end ( ) ; if ( chrom_end < te ) { te = chrom_end ; } sequence += chunk -> get_subsequence ( ts , te ) ; ts = te + 1 ; if ( ts > chrom_end ) { break ; } } if ( strand == " " ) { } return sequence ; }
get_subsequence
Description  : uses ChromChunk objects and the sequence in the database to
return the actual sequence in this region.
get_subsequence Description : uses ChromChunk objects and the sequence in the database to return the actual sequence in this region .
//must provide valid coordinates //fprintf(stderr, "get_subsequence %s %ld %ld %s\n", xml().c_str(), chrom_start, chrom_end, strand.c_str()); //default to start of chromosome //default to end of chromosome //fetch_all_named_region returns chunks in sorted order by chrom_start //first see if region fits entirely inside one chunk //coord compare is much faster than string concats //do rev-complement here //// more than one chunk required to cover entire region so concat //// keep track of how much has already been concatonated by $ts and loop //// then do reverse complement at end //keep running track of position //end is in middle of this chunk //we are done //do rev-complement here