affectanalyser.h 10.2 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 98 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

  /**
   * @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> {
 private:
  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 138 139 140 141 142
/**
 * 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;
 public:

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

143 144
  int count() const;

145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189
  /**
   * @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 
   *         countBefore(before, pos) + countAfter(after, pos) is maximal
   *         and pos >= start.
   * @complexity linear in getSequence().size() 
   */
  int firstMax(const T&before, const T&after, int start=0) const;

  /**
   * @return the last position pos in the sequence such that
   *         countBefore(before, pos) + countAfter(after, pos) is maximal
   *         and pos <= end (if end == -1 considers end of sequence)
   * @complexity linear in getSequence().size()
   */
  int lastMax(const T&before, const T&after, int end=-1) const;

 private:
  /**
   * Build the counts map.
   */
  void buildCounts();

  /**
   * Search the maximum. Used by firstMax and lastMax.
   */
  int searchMax(const T&before, const T &after,
                int start, int end, int iter) const;
};
Mikaël Salson's avatar
Mikaël Salson committed
190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 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

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;
271 272 273 274 275 276
  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
277 278 279
  return kmer;
}

280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297
/* CountKmerAffectAnalyser */

template <class T>
CountKmerAffectAnalyser<T>::CountKmerAffectAnalyser(IKmerStore<T> &kms, const string &seq): KmerAffectAnalyser<T>(kms, seq) {
  buildCounts();
}

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];
  }  
}

298 299 300 301 302
template <class T>
int CountKmerAffectAnalyser<T>::count() const {
  return KmerAffectAnalyser<T>::count();
}

303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 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 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381
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, 
                                      int start) const {
  return searchMax(before, after, start, KmerAffectAnalyser<T>::count()-1,1);
}

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

template <class T>
int CountKmerAffectAnalyser<T>::searchMax(const T&before, const T& after,
                                       int start, int end, int iter) const {
  int first_pos_max = -1;
  int max_value = -1;
  for (int i = start; i <= end; i+=iter) {
    int value = countBefore(before, i) + countAfter(after, i);
    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
382
#endif