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
 118 //////
 119 // Helpful operators
 120 SeqSpan &SeqSpan::operator=(const SeqSpan& s)
 121 {
 122   if (this != &s) {
 123     seq = s.seq;
 124     seq_start = s.seq_start;
 125     seq_count = s.seq_count;
 126     parent = s.parent;
 127     rc_seq = s.rc_seq;
 128   }
 129   return *this;
 130 }
 131
 132 std::ostream& operator<<(std::ostream& out, const SeqSpan& s)
 133 {
 134   out << s.sequence();
 135 }
 136
 137 /* Not implemented yet
 138 //! compare two spans
 139 //! \throws sequence_invalid_comparison
 140 friend bool operator<(const SeqSpan&, const SeqSpan&);
 141 bool operator<(const SeqSpan& a, const SeqSpan& b)
 142 {
 143   // are we subcomponents of the same sequence region?
 144   if (a.seq.get() == b.seq.get()) {
 145     if (a.seq_start < b.seq_start)
 146       return true;
 147     else
 148       return false;
 149   } else {
 150     throw mussa_error("can only compare two spans from the same sequence");
 151   }
 152 }
 153 */
 154
 155 bool operator==(const SeqSpan& a, const SeqSpan& b)
 156 {
 157   if (SeqSpan::isFamily(a, b)) {
 158     // can do fast comparison
 159     if (a.seq_start == b.seq_start and a.seq_count == b.seq_count) {
 160       // unknown strands match anything
 161       if (a.seq_strand == SeqSpan::UnknownStrand or
 162           b.seq_strand == SeqSpan::UnknownStrand) {
 163         return true;
 164       } else {
 165         // otherwise our stranded-ness needs to match
 166         return (a.seq_strand == b.seq_strand);
 167       }
 168     } else {
 169       return false;
 170     }
 171   }
 172   // if we're not part of the same seqspan heirarchy we're never equal
 173   // the string equality is handled off in Sequence
 174   return false;
 175 }
 176
 177 bool operator!=(const SeqSpan& a, const SeqSpan& b)
 178 {
 179   return not (a == b);
 180 }
 181
 182 SeqSpan::const_reference SeqSpan::operator[](SeqSpan::size_type i) const
 183 {
 184   return at(i);
 185 }
 186
 187 SeqSpan::const_reference SeqSpan::at(SeqSpan::size_type i) const
 188 {
 189   if (!seq) throw std::out_of_range("empty sequence");
 190   if (seq_strand == PlusStrand) {
 191     return seq->at(seq_start+i);
 192   } else {
 193     return seq->rc_map[seq->at((seq_start+seq_count-1)-i)];
 194   }
 195 }
 196
 197 const char *SeqSpan::data() const
 198 {
 199   if (seq) {
 200     return seq->c_str()+seq_start;
 201   } else
 202     return 0;
 203 }
 204
 205 SeqSpan::const_iterator SeqSpan::begin() const
 206 {
 207   if (seq and seq_count != 0) {
 208     if (seq_strand != MinusStrand) {
 209       return seq->begin()+seq_start;
 210     } else {
 211       if (not rc_seq) {
 212         initialize_rc_seq();
 213       }
 214       return rc_seq->begin();
 215     }
 216   }
 217   return SeqSpan::const_iterator(0);
 218 }
 219
 220 SeqSpan::const_iterator SeqSpan::end() const
 221 {
 222   if (seq and seq_count != 0) {
 223     if (seq_strand != MinusStrand) {
 224       return seq->begin() + seq_start + seq_count;
 225     } else {
 226       if (not rc_seq) {
 227         initialize_rc_seq();
 228       }
 229       return rc_seq->end();
 230     }
 231   }
 232   return SeqSpan::const_iterator(0);
 233 }
 234
 235 SeqSpan::const_reverse_iterator SeqSpan::rbegin() const
 236 {
 237   if (seq and seq_count != 0) {
 238     if (seq_strand != MinusStrand) {
 239       return seq->rbegin()+(seq->size()-(seq_start+seq_count));
 240     } else {
 241       if (not rc_seq) {
 242         initialize_rc_seq();
 243       }
 244       return rc_seq->rbegin();
 245     }
 246   }
 247   return SeqSpan::const_reverse_iterator();
 248 }
 249
 250 SeqSpan::const_reverse_iterator SeqSpan::rend() const
 251 {
 252   if (seq and seq_count != 0) {
 253     if (seq_strand != MinusStrand) {
 254       return rbegin() + seq_count;
 255     } else {
 256       if (not rc_seq) {
 257         initialize_rc_seq();
 258       }
 259       return rc_seq->rend();
 260     }
 261   }
 262   return SeqSpan::const_reverse_iterator();
 263 }
 264
 265 bool SeqSpan::empty() const
 266 {
 267   return (seq_count == 0) ? true : false;
 268 }
 269
 270 SeqSpan::size_type SeqSpan::find_first_not_of(
 271   const std::string& query,
 272   SeqSpan::size_type index) const
 273 {
 274   typedef std::set<std::string::value_type> sequence_set;
 275   sequence_set match_set;
 276
 277   for(std::string::const_iterator query_item = query.begin();
 278       query_item != query.end();
 279       ++query_item)
 280   {
 281     match_set.insert(*query_item);
 282   }
 283   for(const_iterator base = begin();
 284       base != end();
 285       ++base)
 286   {
 287     if(match_set.find(*base) == match_set.end()) {
 288       return base-begin();
 289     }
 290   }
 291   return SeqSpan::npos;
 292 }
 293
 294 void SeqSpan::setStart(SeqSpan::size_type v)
 295 {
 296   if (v > stop()) {
 297     // cry
 298     throw mussa_error("can't set Start > Stop");
 299   }
 300   seq_count += seq_start - v;
 301   seq_start = v;
 302 }
 303
 304 void SeqSpan::setStop(SeqSpan::size_type v)
 305 {
 306   if ( v < start() ) {
 307     // negative sized sequences are bad
 308     throw mussa_error("can't set Stop < Start");
 309   }
 310   seq_count = std::min<size_type>(v - seq_start, parentSize()-parentStart());
 311 }
 312
 313 SeqSpan::size_type SeqSpan::parentStart() const
 314 {
 315   if (!parent) {
 316     // no parent
 317     return start();
 318   } else {
 319     return start() - parent->start();
 320   }
 321 }
 322
 323 void SeqSpan::setParentStart(SeqSpan::size_type v)
 324 {
 325   setStart(parent->start() + v);
 326 }
 327
 328 SeqSpan::size_type SeqSpan::parentStop() const
 329 {
 330   if (!parent) {
 331     // no parent
 332     return stop();
 333   } else {
 334     return stop() - parent->start();
 335   }
 336 }
 337
 338 void SeqSpan::setParentStop(SeqSpan::size_type v)
 339 {
 340   setStop(parent->start() + v);
 341 }
 342
 343 SeqSpanRef SeqSpan::subseq(size_type start, size_type count, strand_type strand)
 344 {
 345   count = std::min<size_type>(count, seq_count - start);
 346
 347   SeqSpanRef new_span(new SeqSpan(this->shared_from_this(), start, count, strand));
 348   return new_span;
 349 }
 350
 351 std::string SeqSpan::sequence() const
 352 {
 353   if (seq) {
 354     if (seq_strand != MinusStrand)
 355       return seq->substr(seq_start, seq_count);
 356     else {
 357       initialize_rc_seq();
 358       return *rc_seq;
 359     }
 360   } else {
 361     return std::string();
 362   }
 363 }
 364
 365 bool SeqSpan::isFamily(const SeqSpan& a, const SeqSpan& b)
 366 {
 367   return a.seq.get() == b.seq.get();
 368 }
 369
 370 void SeqSpan::initialize_rc_seq() const
 371 {
 372   std::string new_rc_seq;
 373   std::string rc_map(seq->get_alphabet().get_complement_map());
 374
 375   new_rc_seq.reserve(size());
 376   size_type i = (seq_start + seq_count);
 377   while( i != seq_start ) {
 378     --i;
 379     new_rc_seq.append(1, rc_map[ seq->at(i) ]);
 380   }
 381
 382   SeqStringRef new_rc_ref(new SeqString(new_rc_seq, seq->get_alphabet_ref()));
 383   const_cast<SeqStringRef *>(&rc_seq)->swap(new_rc_ref);
 384 }