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