affectanalyser.h 11.5 KB
Newer Older
Mikaël Salson's avatar
Mikaël Salson committed
1 2 3 4 5 6 7
#ifndef AFFECT_ANALYSER_H
#define AFFECT_ANALYSER_H

#include "kmerstore.h"
#include <set>
#include <vector>
#include <cassert>
8
#include <map>
Mikaël Salson's avatar
Mikaël Salson committed
9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43

// Define two constant affectations: ambiguous and unknown.

/* Declaration of types */

typedef enum affect_options_e {
  AO_NONE, AO_NO_CONSECUTIVE, AO_NO_MULTIPLICITY
} affect_options_t;

/**
 * Class that records for every k-mer of a given sequence
 * in which sequences this k-mer was also seen.
 * It can either record one affectation per kmer (the only sequence where it
 * occurs or ambiguous case if there are several possibilities, or 
 * unknown otherwise), or all the possible affectations (ie. there is no ambiguous
 * case, all the possibilities for one k-mer are stored).
 *
 * Input: Index that constitutes the k-mer sequence repertoire
 * Input: Sequence whose k-mers must be affected
 */
template<class T>
class AffectAnalyser {
 public:

  /* Queries */

  /**
   * @return the total number of affectations
   */
  virtual int count() const = 0;

  /**
   * @param affect: An affectation
   * @return the number of times this affectation has been given in the read.
   */
44
  virtual int count(const T &affect) const = 0;
Mikaël Salson's avatar
Mikaël Salson committed
45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97

  /**
   * @param i: the position to consider
   * @pre i >= 0 && i < count()
   * @return the affectation of the k-mer at position i.
   */
  virtual const T&getAffectation(int i)  const = 0;

  /**
   * @param options: options can either be AO_NONE or AO_NO_CONSECUTIVE 
   * @return all the affectations contained in the read from left to right.
   *         if AO_NO_CONSECUTIVE is given: two consecutive elements in the vector
   *                                     will be different (we remove consecutive
   *                                     duplicates)
   */
  virtual vector<T> getAllAffectations(affect_options_t options) const = 0;

  /**
   * @return the distinct affectations
   */
  virtual set<T> getDistinctAffectations() const = 0;

  /**
   * @return the sequence we are analysing
   */
  virtual const string &getSequence() const = 0;

  /**
   * @param affect: an affectation
   * @return the first occurrence of this affectation in the read
   *         or string::npos if the affectation was not found
   * @post getAffectation(first(affect)) == affect 
   * ==>  getAffectation(1...first(affect)-1) != affect
   */
  virtual int first(const T &affect) const  = 0;

  /**
   * @param affect: an affectation
   * @return the last occurrence of this affectation in the read
   *         or string::npos if the affectation was not found
   * @post getAffectation(last(affect)) == affect 
   * ==> getAffectation(last(affect)+1 ... count() -1) != affect
   */
  virtual int last(const T &affect) const  = 0;

  /**
   * @return a string representation of the object
   */
  virtual string toString() const  = 0;
};

template <class T>
class KmerAffectAnalyser: public AffectAnalyser<T> {
98
 protected:
Mikaël Salson's avatar
Mikaël Salson committed
99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129
  IKmerStore<T> &kms;
  const string &seq;
  vector<T> affectations;
 public:
  /**
   * @param kms: the index storing the affectation for the k-mers
   *             (parameter is not copied)
   * @param seq: the sequence to analyse (parameter is not copied)
   */
  KmerAffectAnalyser(IKmerStore<T> &kms, const string &seq);
  ~KmerAffectAnalyser();

  int count() const;

  int count(const T &affect) const;

  const T&getAffectation(int i) const;

  vector<T> getAllAffectations(affect_options_t options) const;

  set<T> getDistinctAffectations() const;

  const string &getSequence() const;

  int first(const T &affect) const;

  int last(const T &affect) const ;

  string toString() const;
};

130 131 132 133 134 135 136 137
/**
 * Class that allows to count in constant time the number of affectations
 * before or after a given point.
 */
template <class T>
class CountKmerAffectAnalyser: public KmerAffectAnalyser<T> {
 private:
  map<T, int* >counts;
138
  int overlap;
139 140 141 142 143
 public:

  CountKmerAffectAnalyser(IKmerStore<T> &kms, const string &seq);
  ~CountKmerAffectAnalyser();

144 145
  int count() const;

146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164
  /**
   * @complexity constant time
   */
  int count(const T &affect) const;

  /**
   * Count the number of an affectation before (strictly) than a position
   * @complexity constant time
   */
  int countBefore(const T&affect, int pos) const;

  /**
   * Count the number of an affectation after (strictly) than a position)
   * @complexity constant time
   */
  int countAfter(const T&affect, int pos) const;

  /**
   * @return the first position pos in the sequence such that 
165 166
   *         countBefore(before, pos - overlap + 1 + overlap) 
             + countAfter(after, pos) is maximal
167
   *         and pos >= start, and the maximum is greater than min; 
168 169
   *         or -1 if such a position doesn't exist.
   *         Where overlap is getAllowedOverlap().
170 171
   * @complexity linear in getSequence().size() 
   */
172
  int firstMax(const T&before, const T&after, int start=0, int min=-1) const;
173 174 175

  /**
   * @return the last position pos in the sequence such that
176 177
   *         countBefore(before, pos - overlap + 1+ overlap) 
   *         + countAfter(after, pos) is maximal
178 179
   *         and pos <= end (if end == -1 considers end of sequence), and the 
   *         maximum is greater than min; or -1 if such a position doesn't exist.
180
   *         Where overlap is getAllowedOverlap().
181 182
   * @complexity linear in getSequence().size()
   */
183
  int lastMax(const T&before, const T&after, int end=-1, int min=-1) const;
184

185 186 187 188 189 190 191 192 193 194 195 196 197
  /**
   * @return the allowed overlap between two k-mers with distinct affectations
   * (default is 0)
   */
  int getAllowedOverlap();

  /**
   * Set the overlap allowed between two k-mers with two different affectations,
   * when looking for the maximum.
   * The overlap should not be greater than the span of the seed used.
   */
  void setAllowedOverlap(int overlap);

198 199 200 201 202 203 204 205
 private:
  /**
   * Build the counts map.
   */
  void buildCounts();

  /**
   * Search the maximum. Used by firstMax and lastMax.
206
   * @param iter: should be either 1 or -1
207 208
   */
  int searchMax(const T&before, const T &after,
209
                int start, int end, int iter, int min) const;
210
};
Mikaël Salson's avatar
Mikaël Salson committed
211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291

template <class T>
KmerAffectAnalyser<T>::KmerAffectAnalyser(IKmerStore<T> &kms, 
                                       const string &seq)
  :kms(kms), seq(seq) {
  assert(seq.length() >=  (size_t)kms.getS());
  affectations = kms.getResults(seq, true);
}

template <class T>
KmerAffectAnalyser<T>::~KmerAffectAnalyser(){}

template <class T>
int KmerAffectAnalyser<T>::count() const{
  return affectations.size();
}

template <class T>
int KmerAffectAnalyser<T>::count(const T &affect) const{
  int count = 0;
  for (typename vector<T>::const_iterator it = affectations.begin(); 
       it < affectations.end(); it++) {
    if (*it == affect)
      count++;
  }
  return count;
}

template <class T>
const T&KmerAffectAnalyser<T>::getAffectation(int i) const{
  assert(i >= 0 && i < count());
  return affectations[i];
}

template <class T>
vector<T> KmerAffectAnalyser<T>::getAllAffectations(affect_options_t options) const{
  if (options == AO_NONE)
    return affectations;
  vector<T> result;
  T previous = affectations[0];
  result.push_back(previous);
  for (size_t i = 1; i < affectations.size(); i++) {
    if (! (previous == affectations[i]))
      result.push_back(previous);
  }
  return result;
}

template <class T>
set<T> KmerAffectAnalyser<T>::getDistinctAffectations() const{
  set<T> result;
  for (size_t i = 0; i < affectations.size(); i++) {    
    result.insert(affectations[i]);
  }
  return result;
}

template <class T>
const string &KmerAffectAnalyser<T>::getSequence() const{
  return seq;
}

template <class T>
int KmerAffectAnalyser<T>::first(const T &affect) const{
  for (size_t i = 0; i < affectations.size(); i++) 
    if (affect == affectations[i])
      return i;
  return (int) string::npos;
}

template <class T>
int KmerAffectAnalyser<T>::last(const T &affect) const{
  for (size_t i = affectations.size(); i > 0;  i--) 
    if (affect == affectations[i-1])
      return i-1;
  return (int) string::npos;
}

template <class T>
string KmerAffectAnalyser<T>::toString() const{
  string kmer;
292 293 294 295 296 297
  for (size_t i = 0; i < affectations.size(); i++) {
    kmer += affectations[i].toString();
#ifdef DEBUG_KMERS
    kmer += ": "+spaced(seq.substr(i,kms.getS()), kms.getSeed())+"\n";
#endif
  }
Mikaël Salson's avatar
Mikaël Salson committed
298 299 300
  return kmer;
}

301 302 303 304 305
/* CountKmerAffectAnalyser */

template <class T>
CountKmerAffectAnalyser<T>::CountKmerAffectAnalyser(IKmerStore<T> &kms, const string &seq): KmerAffectAnalyser<T>(kms, seq) {
  buildCounts();
306
  overlap=0;
307 308 309 310 311 312 313 314 315 316 317 318 319
}

template <class T>
CountKmerAffectAnalyser<T>::~CountKmerAffectAnalyser() {
  set<T> affects = this->getDistinctAffectations();

  /* Initialize each key with a 0-integer array */
  for (typename set<T>::iterator it = affects.begin(); 
       it != affects.end(); it++) {
    delete [] counts[*it];
  }  
}

320 321 322 323 324
template <class T>
int CountKmerAffectAnalyser<T>::count() const {
  return KmerAffectAnalyser<T>::count();
}

325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350
template <class T>
int CountKmerAffectAnalyser<T>::count(const T &affect) const {
  if (counts.count(affect) == 0)
    return 0;

  return counts.find(affect)->second[KmerAffectAnalyser<T>::count() - 1];
}

template <class T>
int CountKmerAffectAnalyser<T>::countBefore(const T&affect, int pos) const {
  if (pos == 0 || counts.count(affect) == 0)
    return 0;
  return counts.find(affect)->second[pos-1];
}

template <class T>
int CountKmerAffectAnalyser<T>::countAfter(const T&affect, int pos) const {
  if (counts.count(affect) == 0)
    return 0;
  int length = KmerAffectAnalyser<T>::count();
  typename map<T, int*>::const_iterator it = counts.find(affect);
  return it->second[length-1] - it->second[pos];
}  

template <class T>
int CountKmerAffectAnalyser<T>::firstMax(const T&before, const T&after, 
351 352
                                         int start, int min) const {
  return searchMax(before, after, start, KmerAffectAnalyser<T>::count()-1,1, min);
353 354 355 356
}

template <class T>
int CountKmerAffectAnalyser<T>::lastMax(const T&before, const T&after, 
357
                                        int end, int min) const {
358 359
  if (end == -1)
    end = KmerAffectAnalyser<T>::count()-1;
360
  return searchMax(before, after, end, 0, -1, min);
361 362
}

363 364 365 366 367 368 369 370 371 372
template <class T>
int CountKmerAffectAnalyser<T>::getAllowedOverlap() {
  return overlap;
}

template <class T>
void CountKmerAffectAnalyser<T>::setAllowedOverlap(int overlap) {
  this->overlap = overlap;
}

373 374
template <class T>
int CountKmerAffectAnalyser<T>::searchMax(const T&before, const T& after,
375
                                          int start, int end, int iter, int min) const {
376
  int first_pos_max = -1;
377
  int max_value = min;
378
  int shift = KmerAffectAnalyser<T>::kms.getS() - overlap - 1;
379 380 381 382 383 384
  int shiftedStart = start, shiftedEnd = end;
  if (iter == 1)
    shiftedStart += shift;
  else
    shiftedEnd += shift;
  for (int i = shiftedStart; (i)*iter <= iter*shiftedEnd; i+=iter) {
385
    int value = countBefore(before, i - shift) + countAfter(after, i);
386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419
    if (value > max_value) {
      max_value = value;
      first_pos_max = i;
    }
  }
  return first_pos_max;
}

template <class T>
void CountKmerAffectAnalyser<T>::buildCounts() {
  int length = KmerAffectAnalyser<T>::count();
  set<T> affects = this->getDistinctAffectations();

  /* Initialize each key with a 0-integer array */
  for (typename set<T>::iterator it = affects.begin(); 
       it != affects.end(); it++) {
    counts[*it] = new int[length];
  }

  /* Fill in the counts arrays */
  for (int i = 0; i < length; i++) {
    T current = this->getAffectation(i);
    for (typename set<T>::iterator it = affects.begin(); 
         it != affects.end(); it++) {
      int value = (current == *it) ? 1 : 0;

      if (i == 0)
        counts[*it][i] = value;
      else
        counts[*it][i] = counts[*it][i-1] + value;
    }
  }
}

Mikaël Salson's avatar
Mikaël Salson committed
420
#endif