alg/seq_span.cpp

   1 #include <stdexcept>
   2 #include <set>
   3 #include <algorithm>
   4
   5 #include "seq_span.hpp"
   6 #include "mussa_exceptions.hpp"
   7
   8 SeqSpan::SeqSpan(const SeqSpan &o)
   9   : seq(o.seq),
  10     seq_start(o.seq_start),
  11     seq_count(o.seq_count),
  12     parent(o.parent),
  13     rc_seq(o.rc_seq)
  14 {
  15 }
  16
  17 SeqSpan::SeqSpan(const SeqSpan *p)
  18   : seq(p->seq),
  19     seq_start(p->seq_start),
  20     seq_count(p->seq_count),
  21     parent(p->parent),
  22     rc_seq(p->rc_seq)
  23 {
  24 }
  25
  26 SeqSpan::SeqSpan(const std::string &seq_,
  27                  AlphabetRef a,
  28                  strand_type strand_)
  29   : seq_start(0),
  30     seq_count(seq_.length()),
  31     parent()
  32 {
  33   switch (strand_) {
  34   case PlusStrand:
  35   case SingleStrand:
  36     // do nothing
  37     seq_strand = strand_;
  38     seq.reset(new SeqString(seq_, a));
  39   break;
  40   default:
  41     throw sequence_invalid_strand(
  42             "Can only initialize a Plus or Single strand from a string"
  43           );
  44     break;
  45   }
  46 }
  47
  48 SeqSpan::SeqSpan(const SeqSpanRef parent_,
  49                  size_type start_,
  50                  size_type count_,
  51                  strand_type strand_)
  52   : seq(parent_->seq),
  53     seq_start(parent_->seq_start + start_),
  54     parent(parent_)
  55 {
  56   if (count_ == npos)
  57     seq_count = parent_->seq_count;
  58   else
  59     seq_count = count_;
  60
  61   // blech this strand logic is too long
  62   switch(strand_) {
  63   case UnknownStrand:
  64   case PlusStrand:
  65   case MinusStrand:
  66   case BothStrand:
  67     seq_strand = strand_;
  68   break;
  69   case SameStrand:
  70     if (parent) {
  71       seq_strand = parent->strand();
  72     } else {
  73       throw sequence_invalid_strand("SameStrand is undefined with no parent");
  74     }
  75   break;
  76   case OppositeStrand:
  77     if (parent) {
  78       switch (parent->strand()) {
  79       case PlusStrand:
  80         seq_strand = MinusStrand;
  81       break;
  82       case MinusStrand:
  83         seq_strand = PlusStrand;
  84       break;
  85       case UnknownStrand:
  86       case BothStrand:
  87         seq_strand = parent->strand();
  88       break;
  89       case SingleStrand:
  90         throw sequence_invalid_strand(
  91                 "OppositeStrand is undefined with SingleStrands"
  92               );
  93         break;
  94       default:
  95         throw sequence_invalid_strand("Parent has an invalid strand type");
  96       }
  97     } else {
  98       throw sequence_invalid_strand("SameStrand is undefined with no parent");
  99     }
 100   break;
 101   case SingleStrand:
 102     if (parent) {
 103       if (parent->strand() == SingleStrand) {
 104         seq_strand = SingleStrand;
 105       } else {
 106         throw  sequence_invalid_strand("Can't change single strandedness");
 107       }
 108     } else {
 109       seq_strand = SingleStrand;
 110     }
 111   break;
 112   default:
 113     throw sequence_invalid_strand("unrecognized strand identifier");
 114     break;
 115   }
 116
 117   // Ack the complexity increases!
 118   // If our parent is on the minus strand, we need to adjust the start
 119   // and count to look like we're selecting from the right side of the
 120   // parent sequence (AKA the start of the minus strand)
 121   if (parent and parent->strand() == MinusStrand) {
 122     seq_start = parent->start() + parent->size() - (start_ + seq_count);
 123   }
 124 }
 125
 126 //////
 127 // Helpful operators
 128 SeqSpan &SeqSpan::operator=(const SeqSpan& s)
 129 {
 130   if (this != &s) {
 131     seq = s.seq;
 132     seq_start = s.seq_start;
 133     seq_count = s.seq_count;
 134     parent = s.parent;
 135     rc_seq = s.rc_seq;
 136   }
 137   return *this;
 138 }
 139
 140 std::ostream& operator<<(std::ostream& out, const SeqSpan& s)
 141 {
 142   out << s.sequence();
 143 }
 144
 145 /* Not implemented yet
 146 //! compare two spans
 147 //! \throws sequence_invalid_comparison
 148 friend bool operator<(const SeqSpan&, const SeqSpan&);
 149 bool operator<(const SeqSpan& a, const SeqSpan& b)
 150 {
 151   // are we subcomponents of the same sequence region?
 152   if (a.seq.get() == b.seq.get()) {
 153     if (a.seq_start < b.seq_start)
 154       return true;
 155     else
 156       return false;
 157   } else {
 158     throw mussa_error("can only compare two spans from the same sequence");
 159   }
 160 }
 161 */
 162
 163 bool operator==(const SeqSpan& a, const SeqSpan& b)
 164 {
 165   if (SeqSpan::isFamily(a, b)) {
 166     // can do fast comparison
 167     if (a.seq_start == b.seq_start and a.seq_count == b.seq_count) {
 168       // unknown strands match anything
 169       if (a.seq_strand == SeqSpan::UnknownStrand or
 170           b.seq_strand == SeqSpan::UnknownStrand) {
 171         return true;
 172       } else {
 173         // otherwise our stranded-ness needs to match
 174         return (a.seq_strand == b.seq_strand);
 175       }
 176     } else {
 177       return false;
 178     }
 179   }
 180   // if we're not part of the same seqspan heirarchy we're never equal
 181   // the string equality is handled off in Sequence
 182   return false;
 183 }
 184
 185 bool operator!=(const SeqSpan& a, const SeqSpan& b)
 186 {
 187   return not (a == b);
 188 }
 189
 190 SeqSpan::const_reference SeqSpan::operator[](SeqSpan::size_type i) const
 191 {
 192   return at(i);
 193 }
 194
 195 SeqSpan::const_reference SeqSpan::at(SeqSpan::size_type i) const
 196 {
 197   if (!seq) throw std::out_of_range("empty sequence");
 198   if (seq_strand == PlusStrand) {
 199     return seq->at(seq_start+i);
 200   } else {
 201     return seq->rc_map[seq->at((seq_start+seq_count-1)-i)];
 202   }
 203 }
 204
 205 const char *SeqSpan::data() const
 206 {
 207   if (seq) {
 208     return seq->c_str()+seq_start;
 209   } else
 210     return 0;
 211 }
 212
 213 SeqSpan::const_iterator SeqSpan::begin() const
 214 {
 215   if (seq and seq_count != 0) {
 216     if (seq_strand != MinusStrand) {
 217       return seq->begin()+seq_start;
 218     } else {
 219       if (not rc_seq) {
 220         initialize_rc_seq();
 221       }
 222       return rc_seq->begin();
 223     }
 224   }
 225   return SeqSpan::const_iterator(0);
 226 }
 227
 228 SeqSpan::const_iterator SeqSpan::end() const
 229 {
 230   if (seq and seq_count != 0) {
 231     if (seq_strand != MinusStrand) {
 232       return seq->begin() + seq_start + seq_count;
 233     } else {
 234       if (not rc_seq) {
 235         initialize_rc_seq();
 236       }
 237       return rc_seq->end();
 238     }
 239   }
 240   return SeqSpan::const_iterator(0);
 241 }
 242
 243 SeqSpan::const_reverse_iterator SeqSpan::rbegin() const
 244 {
 245   if (seq and seq_count != 0) {
 246     if (seq_strand != MinusStrand) {
 247       return seq->rbegin()+(seq->size()-(seq_start+seq_count));
 248     } else {
 249       if (not rc_seq) {
 250         initialize_rc_seq();
 251       }
 252       return rc_seq->rbegin();
 253     }
 254   }
 255   return SeqSpan::const_reverse_iterator();
 256 }
 257
 258 SeqSpan::const_reverse_iterator SeqSpan::rend() const
 259 {
 260   if (seq and seq_count != 0) {
 261     if (seq_strand != MinusStrand) {
 262       return rbegin() + seq_count;
 263     } else {
 264       if (not rc_seq) {
 265         initialize_rc_seq();
 266       }
 267       return rc_seq->rend();
 268     }
 269   }
 270   return SeqSpan::const_reverse_iterator();
 271 }
 272
 273 bool SeqSpan::empty() const
 274 {
 275   return (seq_count == 0) ? true : false;
 276 }
 277
 278 SeqSpan::size_type SeqSpan::find_first_not_of(
 279   const std::string& query,
 280   SeqSpan::size_type index) const
 281 {
 282   typedef std::set<std::string::value_type> sequence_set;
 283   sequence_set match_set;
 284
 285   for(std::string::const_iterator query_item = query.begin();
 286       query_item != query.end();
 287       ++query_item)
 288   {
 289     match_set.insert(*query_item);
 290   }
 291   for(const_iterator base = begin();
 292       base != end();
 293       ++base)
 294   {
 295     if(match_set.find(*base) == match_set.end()) {
 296       return base-begin();
 297     }
 298   }
 299   return SeqSpan::npos;
 300 }
 301
 302 void SeqSpan::setStart(SeqSpan::size_type v)
 303 {
 304   if (v > stop()) {
 305     // cry
 306     throw mussa_error("can't set Start > Stop");
 307   }
 308   seq_count += seq_start - v;
 309   seq_start = v;
 310 }
 311
 312 void SeqSpan::setStop(SeqSpan::size_type v)
 313 {
 314   if ( v < start() ) {
 315     // negative sized sequences are bad
 316     throw mussa_error("can't set Stop < Start");
 317   }
 318   seq_count = std::min<size_type>(v - seq_start, parentSize()-parentStart());
 319 }
 320
 321 SeqSpan::size_type SeqSpan::parentStart() const
 322 {
 323   if (!parent) {
 324     // no parent
 325     return start();
 326   } else {
 327     return start() - parent->start();
 328   }
 329 }
 330
 331 void SeqSpan::setParentStart(SeqSpan::size_type v)
 332 {
 333   setStart(parent->start() + v);
 334 }
 335
 336 SeqSpan::size_type SeqSpan::parentStop() const
 337 {
 338   if (!parent) {
 339     // no parent
 340     return stop();
 341   } else {
 342     return stop() - parent->start();
 343   }
 344 }
 345
 346 void SeqSpan::setParentStop(SeqSpan::size_type v)
 347 {
 348   setStop(parent->start() + v);
 349 }
 350
 351 SeqSpanRef SeqSpan::subseq(size_type start, size_type count, strand_type strand)
 352 {
 353   count = std::min<size_type>(count, seq_count - start);
 354
 355   SeqSpanRef new_span(new SeqSpan(this->shared_from_this(), start, count, strand));
 356   return new_span;
 357 }
 358
 359 std::string SeqSpan::sequence() const
 360 {
 361   if (seq) {
 362     if (seq_strand != MinusStrand)
 363       return seq->substr(seq_start, seq_count);
 364     else {
 365       initialize_rc_seq();
 366       return *rc_seq;
 367     }
 368   } else {
 369     return std::string();
 370   }
 371 }
 372
 373 bool SeqSpan::isFamily(const SeqSpan& a, const SeqSpan& b)
 374 {
 375   return a.seq.get() == b.seq.get();
 376 }
 377
 378 void SeqSpan::initialize_rc_seq() const
 379 {
 380   std::string new_rc_seq;
 381   std::string rc_map(seq->get_alphabet().get_complement_map());
 382
 383   new_rc_seq.reserve(size());
 384   size_type i = (seq_start + seq_count);
 385   while( i != seq_start ) {
 386     --i;
 387     new_rc_seq.append(1, rc_map[ seq->at(i) ]);
 388   }
 389
 390   SeqStringRef new_rc_ref(new SeqString(new_rc_seq, seq->get_alphabet_ref()));
 391   const_cast<SeqStringRef *>(&rc_seq)->swap(new_rc_ref);
 392 }