testAffectAnalyser.cpp 26.7 KB
Newer Older
Mikaël Salson's avatar
Mikaël Salson committed
1 2 3 4 5
#include "tests.h"
#include <core/affectanalyser.h>
#include <core/kmerstore.h>
#include <core/kmeraffect.h>

6
template<template <class> class T>
Mikaël Salson's avatar
Mikaël Salson committed
7 8 9
void testAffectAnalyser1() {
  const int k = 4;
  const bool revcomp = false;
10 11

  T<KmerAffect> *index = createIndex<T<KmerAffect> >(k, revcomp);
Mikaël Salson's avatar
Mikaël Salson committed
12
  
13
  string sequence = "AAAACCCCCGGGGG";
14 15
  KmerAffectAnalyser kaa(*index, sequence);
  CountKmerAffectAnalyser ckaa(*index, sequence);
16
  TAP_TEST_EQUAL(ckaa.getAllowedOverlap(), 0, TEST_COUNT_AA_GET_OVERLAP, "");
17
  ckaa.setAllowedOverlap(k-1);
Mikaël Salson's avatar
Mikaël Salson committed
18

19 20 21
  set<KmerAffect> forbidden;
  forbidden.insert(KmerAffect::getAmbiguous());
  forbidden.insert(KmerAffect::getUnknown());
22

23 24
  TAP_TEST_EQUAL(ckaa.getAllowedOverlap(), k-1, TEST_COUNT_AA_GET_OVERLAP, "");
  TAP_TEST_EQUAL(ckaa.getSequence(), "AAAACCCCCGGGGG", TEST_AA_GET_SEQUENCE, "actual: " << ckaa.getSequence());
25

26 27
  bool full_length = (kaa.getAllAffectations(AO_NONE).size() == 14); // Hack for different lenghts
  int shift = 4/2 ; // getS(), See #3727
28 29 30 31
  TAP_TEST_EQUAL(kaa.minimize(KmerAffect::getAmbiguous(), 0, 4), 0+shift, TEST_AA_MINIMIZE, ""); // first k-mer, AAAA is ambiguous
  TAP_TEST_EQUAL(kaa.minimize(KmerAffect::getAmbiguous(), 1, 4), NO_MINIMIZING_POSITION, TEST_AA_MINIMIZE, ""); // too large margin (left side)
  TAP_TEST_EQUAL(kaa.minimize(KmerAffect("A", 1, k), 0, 3), NO_MINIMIZING_POSITION, TEST_AA_MINIMIZE, ""); // no non-ambiguous AAAA
  TAP_TEST_EQUAL(kaa.minimize(KmerAffect("C", 1, k), 3, 4), 4+shift, TEST_AA_MINIMIZE, ""); // margin = 3, does not affect C
32
  if (full_length)
33 34
  TAP_TEST_EQUAL(kaa.minimize(KmerAffect("C", 1, k), 5, 4), 5+shift, TEST_AA_MINIMIZE, ""); // margin = 5, second k-mer C exactly fits between both margins
  TAP_TEST_EQUAL(kaa.minimize(KmerAffect("G", 1, k), 5, 4), NO_MINIMIZING_POSITION, TEST_AA_MINIMIZE, ""); // too large margin (right side)
Mathieu Giraud's avatar
Mathieu Giraud committed
35

Mikaël Salson's avatar
Mikaël Salson committed
36
  for (int i = 2; i < nb_seq-1; i++) {
Mikaël Salson's avatar
Mikaël Salson committed
37 38
    // i starts at 2 because AAAA is not found: there is an ambiguity with
    // AAAA coming from AAAACAAAACAAAAC or AAAAAAAAAAAAAAA
39
    KmerAffect current_affect(seq[2*i+1], 1, k);
40 41
    TAP_TEST_EQUAL(kaa.count(current_affect), 0, TEST_AA_COUNT, "");
    TAP_TEST_EQUAL(ckaa.count(current_affect), 0, TEST_COUNT_AA_COUNT, ckaa.count(current_affect));
Mikaël Salson's avatar
Mikaël Salson committed
42 43 44 45
    TAP_TEST(kaa.first(current_affect) == (int)string::npos, TEST_AA_FIRST, "");
    TAP_TEST(kaa.last(current_affect) == (int)string::npos, TEST_AA_LAST, "");
  }
  for (int i = 0; i < 2; i++) {
46
    KmerAffect current_affect(seq[2*i+1], 1, k);
47 48
    TAP_TEST_EQUAL(kaa.count(current_affect), 2, TEST_AA_COUNT, kaa.count(current_affect));
    TAP_TEST_EQUAL(ckaa.count(current_affect), 2, TEST_COUNT_AA_COUNT, ckaa.count(current_affect));
Mikaël Salson's avatar
Mikaël Salson committed
49 50 51
    TAP_TEST(kaa.getAffectation(kaa.first(current_affect)) == current_affect, TEST_AA_GET_AFFECT, "");
    TAP_TEST(kaa.getAffectation(kaa.last(current_affect)) == current_affect, TEST_AA_GET_AFFECT, "");
  }
52 53 54
  TAP_TEST(kaa.count(KmerAffect(seq[2*(nb_seq-1)+1], 1, k)) == 1, TEST_AA_COUNT, "");
  TAP_TEST((kaa.first(KmerAffect(seq[2*(nb_seq-1)+1], 1, k))
           == kaa.last(KmerAffect(seq[2*(nb_seq-1)+1], 1, k)))
Mikaël Salson's avatar
Mikaël Salson committed
55
           == 1, TEST_AA_FIRST, "");
Mikaël Salson's avatar
Mikaël Salson committed
56

57 58
  TAP_TEST(ckaa.max(forbidden) == KmerAffect("C lots of", 1, k)
           || ckaa.max(forbidden) == KmerAffect("G lots of", 1, k),
59 60
           TEST_COUNT_AA_MAX, "max is " << ckaa.max(forbidden));

61
  TAP_TEST(ckaa.max() == KmerAffect::getUnknown(), 
62 63
           TEST_COUNT_AA_MAX, "max is " << ckaa.max());

64 65 66
  TAP_TEST(kaa.getAffectation(6  - k + 1).isUnknown(), TEST_AA_PREDICATES, "");
  TAP_TEST(kaa.getAffectation(11  - k + 1).isUnknown(), TEST_AA_PREDICATES, "");
  TAP_TEST(kaa.getAffectation(3  - k + 1).isAmbiguous(), TEST_AA_PREDICATES, "");
Mikaël Salson's avatar
Mikaël Salson committed
67 68 69
  
  TAP_TEST(kaa.getDistinctAffectations().size() == 5, TEST_AA_GET_DISTINCT_AFFECT, "");

70 71
  KmerAffect cAffect = KmerAffect(seq[1], 1, k);
  KmerAffect gAffect = KmerAffect(seq[3], 1, k);
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
  TAP_TEST_EQUAL(ckaa.countBefore(cAffect, 0), 0, TEST_COUNT_AA_COUNT_BEFORE, "");
  TAP_TEST_EQUAL(ckaa.countBefore(gAffect, 0), 0, TEST_COUNT_AA_COUNT_BEFORE, "");
  TAP_TEST_EQUAL(ckaa.countAfter(cAffect, 10), 0, TEST_COUNT_AA_COUNT_BEFORE, "");
  TAP_TEST_EQUAL(ckaa.countAfter(gAffect, 13  - k + 1), 0, TEST_COUNT_AA_COUNT_BEFORE, "");

  TAP_TEST_EQUAL(ckaa.countBefore(cAffect, 7  - k + 1), 0, TEST_COUNT_AA_COUNT_BEFORE, "");
  TAP_TEST_EQUAL(ckaa.countBefore(cAffect, 8  - k + 1), 1, TEST_COUNT_AA_COUNT_BEFORE, "");
  TAP_TEST_EQUAL(ckaa.countAfter(cAffect, 7  - k + 1), 1, TEST_COUNT_AA_COUNT_AFTER, "");
  TAP_TEST_EQUAL(ckaa.countAfter(cAffect, 8  - k + 1), 0, TEST_COUNT_AA_COUNT_AFTER, "");

  TAP_TEST_EQUAL(ckaa.countAfter(gAffect, 7  - k + 1), 2, TEST_COUNT_AA_COUNT_AFTER, "");
  TAP_TEST_EQUAL(ckaa.countAfter(gAffect, 8  - k + 1), 2, TEST_COUNT_AA_COUNT_AFTER, "");
  TAP_TEST_EQUAL(ckaa.countBefore(gAffect, 7  - k + 1), 0, TEST_COUNT_AA_COUNT_BEFORE, "");
  TAP_TEST_EQUAL(ckaa.countBefore(gAffect, 8  - k + 1 ), 0, TEST_COUNT_AA_COUNT_BEFORE, "");

  TAP_TEST_EQUAL(ckaa.countBefore(cAffect, 12  - k + 1), 2, TEST_COUNT_AA_COUNT_BEFORE, "");
  TAP_TEST_EQUAL(ckaa.countBefore(cAffect, 13  - k + 1), 2, TEST_COUNT_AA_COUNT_BEFORE, "");
  TAP_TEST_EQUAL(ckaa.countAfter(cAffect, 12  - k + 1), 0, TEST_COUNT_AA_COUNT_AFTER, "");
  TAP_TEST_EQUAL(ckaa.countAfter(cAffect, 13  - k + 1), 0, TEST_COUNT_AA_COUNT_AFTER, "");

  TAP_TEST_EQUAL(ckaa.countAfter(gAffect, 12  - k + 1), 1, TEST_COUNT_AA_COUNT_AFTER, "");
  TAP_TEST_EQUAL(ckaa.countAfter(gAffect, 13  - k + 1), 0, TEST_COUNT_AA_COUNT_AFTER, "");
  TAP_TEST_EQUAL(ckaa.countBefore(gAffect, 12  - k + 1), 0, TEST_COUNT_AA_COUNT_BEFORE, "");
  TAP_TEST_EQUAL(ckaa.countBefore(gAffect, 13  - k + 1), 1, TEST_COUNT_AA_COUNT_BEFORE, "");

  TAP_TEST_EQUAL(ckaa.firstMax(cAffect, gAffect), 9 - k + 1, TEST_COUNT_AA_FIRST_MAX, "");
  TAP_TEST_EQUAL(ckaa.lastMax(cAffect, gAffect), 11 -  k + 1, TEST_COUNT_AA_LAST_MAX, ckaa.lastMax(cAffect, gAffect));
Mikaël Salson's avatar
Mikaël Salson committed
99 100

  // Test affectation with two affects that are not in the sequence
101 102
  KmerAffect aAffect = KmerAffect(seq[5], 1, k);
  KmerAffect tAffect = KmerAffect(seq[7], 1, k);
103 104
  TAP_TEST_EQUAL(ckaa.firstMax(aAffect, tAffect), -1, TEST_COUNT_AA_FIRST_MAX, "");
  TAP_TEST_EQUAL(ckaa.lastMax(aAffect, tAffect), - 1, 
Mikaël Salson's avatar
Mikaël Salson committed
105
           TEST_COUNT_AA_LAST_MAX, "");
106 107
  TAP_TEST_EQUAL(ckaa.countAfter(tAffect, 7  - k + 1), 0, TEST_COUNT_AA_COUNT_AFTER, "");
  TAP_TEST_EQUAL(ckaa.countBefore(tAffect, 7  - k + 1), 0, TEST_COUNT_AA_COUNT_BEFORE, "");
Mikaël Salson's avatar
Mikaël Salson committed
108

109 110
  // Test affectation with one affect not in the sequence

111 112
  TAP_TEST_EQUAL(ckaa.firstMax(cAffect, tAffect), -1, TEST_COUNT_AA_FIRST_MAX, "");
  TAP_TEST_EQUAL(ckaa.lastMax(cAffect, tAffect), -1, 
113 114
           TEST_COUNT_AA_LAST_MAX, "");

115 116
  TAP_TEST_EQUAL(ckaa.firstMax(aAffect, gAffect), -1, TEST_COUNT_AA_FIRST_MAX, "");
  TAP_TEST_EQUAL(ckaa.lastMax(aAffect, gAffect), -1, TEST_COUNT_AA_LAST_MAX, "");
Mikaël Salson's avatar
Mikaël Salson committed
117 118 119
  delete index;
}

Mikaël Salson's avatar
Mikaël Salson committed
120
// Test with revcomp
121
template<template <class> class T>
Mikaël Salson's avatar
Mikaël Salson committed
122 123 124
void testAffectAnalyser2() {
  const int k = 5;
  const bool revcomp = true;
125
  T<KmerAffect> *index = createIndex<T<KmerAffect> >(k, revcomp);
126
  string sequence = "TTTTTGGGGG";
127 128
  KmerAffectAnalyser kaa(*index, sequence);
  CountKmerAffectAnalyser ckaa(*index, sequence);
129
  ckaa.setAllowedOverlap(k-1);
130

131 132 133
  set<KmerAffect> forbidden;
  forbidden.insert(KmerAffect::getAmbiguous());
  forbidden.insert(KmerAffect::getUnknown());
Mikaël Salson's avatar
Mikaël Salson committed
134
  
135 136
  TAP_TEST_EQUAL(kaa.getSequence(), "TTTTTGGGGG", TEST_AA_GET_SEQUENCE, "actual: ");
  TAP_TEST_EQUAL(ckaa.getSequence(), "TTTTTGGGGG", TEST_AA_GET_SEQUENCE, "actual: " << ckaa.getSequence());
137

138 139 140 141 142
  TAP_TEST(kaa.getAffectation(1) == KmerAffect(seq[2*(nb_seq-1)+1], -1, k), TEST_AA_GET_AFFECT, "");
  TAP_TEST(kaa.count(kaa.getAffectation(1)) == 1, TEST_AA_GET_AFFECT, "");
  TAP_TEST(ckaa.count(kaa.getAffectation(1)) == 1, TEST_COUNT_AA_COUNT, "");
  TAP_TEST(kaa.getAffectation(0) == kaa.getAffectation(10 - k), TEST_AA_GET_AFFECT, "");
  TAP_TEST(kaa.getAffectation(0).isAmbiguous(), TEST_AA_PREDICATES, "");
Mikaël Salson's avatar
Mikaël Salson committed
143

144
  for (int i = 6; i < 14 - k; i++)
145
    TAP_TEST(kaa.getAffectation(i  - k + 1).isUnknown(), TEST_AA_PREDICATES, "");
Mikaël Salson's avatar
Mikaël Salson committed
146 147 148

  TAP_TEST(kaa.getDistinctAffectations().size() == 3, TEST_AA_GET_DISTINCT_AFFECT, "");

149
  TAP_TEST(ckaa.max(forbidden) == KmerAffect(seq[2*(nb_seq-1)+1], -1, k),
150 151
           TEST_COUNT_AA_MAX, "max is " << ckaa.max(forbidden));

152
  TAP_TEST(ckaa.max() == KmerAffect::getUnknown(), 
153 154
           TEST_COUNT_AA_MAX, "max is " << ckaa.max());

155

156
  for (int i = 4; i < 14 - k; i++)
157
    TAP_TEST(kaa.getAffectation(i  - k + 1) == kaa.getAllAffectations(AO_NONE)[i  - k + 1], TEST_AA_GET_ALL_AO_NONE, "");
158 159

  TAP_TEST(kaa.getAffectation(0) == kaa.getAllAffectations(AO_NO_CONSECUTIVE)[0], TEST_AA_GET_ALL_AO_NO_CONSECUTIVE, "");
160
  if ((size_t)kaa.count() == sequence.length() - k + 1) {
161
    TAP_TEST(kaa.getAllAffectations(AO_NO_CONSECUTIVE).size() == 4, TEST_AA_GET_ALL_AO_NO_CONSECUTIVE, "size = " << kaa.getAllAffectations(AO_NO_CONSECUTIVE).size());
162
  } else if ((size_t) kaa.count() == sequence.length()) {
163 164
    TAP_TEST(kaa.getAllAffectations(AO_NO_CONSECUTIVE).size() == 5, TEST_AA_GET_ALL_AO_NO_CONSECUTIVE, "size = " << kaa.getAllAffectations(AO_NO_CONSECUTIVE).size());
  }
165 166 167
  TAP_TEST(kaa.getAffectation(1) == kaa.getAllAffectations(AO_NO_CONSECUTIVE)[1], TEST_AA_GET_ALL_AO_NO_CONSECUTIVE, "actual: " << kaa.getAllAffectations(AO_NO_CONSECUTIVE)[1] << ", expected: " << kaa.getAffectation(1));
  TAP_TEST(kaa.getAffectation(2) == kaa.getAllAffectations(AO_NO_CONSECUTIVE)[2], TEST_AA_GET_ALL_AO_NO_CONSECUTIVE, kaa.getAllAffectations(AO_NO_CONSECUTIVE)[2] << ", expected: " << kaa.getAffectation(2));
  TAP_TEST(kaa.getAllAffectations(AO_NO_CONSECUTIVE)[3] == kaa.getAffectation(10-k), TEST_AA_GET_ALL_AO_NO_CONSECUTIVE, kaa.getAllAffectations(AO_NO_CONSECUTIVE)[3] << ", expected: " << kaa.getAffectation(10-k));
168

Mikaël Salson's avatar
Mikaël Salson committed
169 170 171
  delete index;
}

172 173 174 175 176 177 178 179 180 181 182

template<template <class> class T>
void testAffectAnalyserMaxes() {
  const int k = 4;
  const bool revcomp = false;

  T<KmerAffect> *index = createIndex<T<KmerAffect> >(k, revcomp);
  
  string sequence = "ACCCCAGGGGGA";
  CountKmerAffectAnalyser ckaa(*index, sequence);

183 184
  KmerAffect cAffect = KmerAffect("C", 1, k);
  KmerAffect gAffect = KmerAffect("G", 1, k);
185 186 187 188 189 190 191 192 193

  set<KmerAffect> forbidden;
  forbidden.insert(KmerAffect::getAmbiguous());
  forbidden.insert(KmerAffect::getUnknown());

  TAP_TEST(ckaa.max(forbidden) == gAffect, TEST_COUNT_AA_MAX, "max is " << ckaa.max(forbidden));
  TAP_TEST(ckaa.max12(forbidden).first == gAffect, TEST_COUNT_AA_MAX12, "max1 is " << ckaa.max12(forbidden).first);
  TAP_TEST(ckaa.max12(forbidden).second == cAffect, TEST_COUNT_AA_MAX12, "max2 is " << ckaa.max12(forbidden).second);

194 195 196
  TAP_TEST(ckaa.sortLeftRight(ckaa.max12(forbidden)).first == cAffect, TEST_AA_SORT_LEFT_RIGHT, "bad max12, left");
  TAP_TEST(ckaa.sortLeftRight(ckaa.max12(forbidden)).second == gAffect, TEST_AA_SORT_LEFT_RIGHT, "bad max12, right");

197 198 199 200 201 202 203 204 205 206 207 208 209
  forbidden.insert(gAffect);
  TAP_TEST(ckaa.max(forbidden) == cAffect, TEST_COUNT_AA_MAX, "max is " << ckaa.max(forbidden));
  TAP_TEST(ckaa.max12(forbidden).first == cAffect, TEST_COUNT_AA_MAX12, "max1 is " << ckaa.max12(forbidden).first);
  TAP_TEST(ckaa.max12(forbidden).second == KmerAffect::getUnknown(), TEST_COUNT_AA_MAX12, "max2 is " << ckaa.max12(forbidden).second);

  forbidden.insert(cAffect);
  TAP_TEST(ckaa.max(forbidden) == KmerAffect::getUnknown(), TEST_COUNT_AA_MAX, "max is " << ckaa.max(forbidden));
  TAP_TEST(ckaa.max12(forbidden).first == KmerAffect::getUnknown(), TEST_COUNT_AA_MAX12, "max1 is " << ckaa.max12(forbidden).first);
  TAP_TEST(ckaa.max12(forbidden).second == KmerAffect::getUnknown(), TEST_COUNT_AA_MAX12, "max2 is " << ckaa.max12(forbidden).second);

  delete index;
}

210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230
template<template <class> class T>
void testAffectAnalyserSpaced() {
  const string seed = "##-##";
  const bool revcomp = false;

  T<KmerAffect> *index = createIndex<T<KmerAffect> >(seed, revcomp);

  string sequence = "AAAACCCCCGGGGG";
  // AA-AC, AA-CC, AA-CC, AC-CC, CC-CC, CC-CG, CC-GG, CG-GG, GG-GG
  KmerAffectAnalyser kaa(*index, sequence);
  CountKmerAffectAnalyser ckaa(*index, sequence);

  KmerAffect affect_AAAC(seq[11], 1, seed.size());
  KmerAffect affect_AAAC_bad(seq[11], 1, seed_weight(seed));
  KmerAffect affect_CCCC(seq[1], 1, seed.size());
  KmerAffect affect_CCCC_bad(seq[1], 1, seed_weight(seed));
  KmerAffect affect_GGGG(seq[3], 1, seed.size());
  KmerAffect affect_GGGG_bad(seq[3], 1, seed_weight(seed));

  TAP_TEST(kaa.toString().find("+V _ _ _+C _ _ _ _+G") != string::npos,
           TEST_AA_TO_STRING, "");
231 232 233
  TAP_TEST_EQUAL(kaa.count(affect_AAAC), 1, TEST_AA_COUNT, "Expected 1, got " << kaa.count(affect_AAAC) << " in " << __PRETTY_FUNCTION__);
  TAP_TEST_EQUAL(kaa.count(affect_CCCC), 1, TEST_AA_COUNT, "Expected 1, got " << kaa.count(affect_CCCC) << " in " << __PRETTY_FUNCTION__);
  TAP_TEST_EQUAL(kaa.count(affect_GGGG), 1, TEST_AA_COUNT, "Expected 1, got " << kaa.count(affect_GGGG) << " in " << __PRETTY_FUNCTION__);
234

235 236 237
  TAP_TEST_EQUAL(kaa.count(affect_AAAC_bad), 0, TEST_AA_COUNT, "Expected 0, got " << kaa.count(affect_AAAC_bad) << " in " << __PRETTY_FUNCTION__);
  TAP_TEST_EQUAL(kaa.count(affect_CCCC_bad), 0, TEST_AA_COUNT, "Expected 0, got " << kaa.count(affect_CCCC_bad) << " in " << __PRETTY_FUNCTION__);
  TAP_TEST_EQUAL(kaa.count(affect_GGGG_bad), 0, TEST_AA_COUNT, "Expected 0, got " << kaa.count(affect_GGGG_bad) << " in " << __PRETTY_FUNCTION__);
238
  delete index;
239 240
}

241 242 243 244 245 246 247 248 249
template<template <class> class T>
void testGetMaximum() {
  const int k = 4;
  const bool revcomp = true;
  T<KmerAffect> *index = createIndex<T<KmerAffect> >(k, revcomp);

  KmerAffect a[] = {AFFECT_J_BWD, AFFECT_J_BWD, AFFECT_UNKNOWN, AFFECT_UNKNOWN, AFFECT_UNKNOWN, AFFECT_UNKNOWN, AFFECT_V_BWD, AFFECT_V_BWD, AFFECT_V_BWD, AFFECT_UNKNOWN, AFFECT_UNKNOWN, AFFECT_UNKNOWN, AFFECT_J_BWD, AFFECT_J_BWD};
  //  0 1 2 3 4 5 6 7 8 9  11  13
  // J-J- _ _ _ _V-V-V- _ _ _J-J-
250
  //   ^^^^^^^^^^
251 252
  vector<KmerAffect> affectations(a, a+sizeof(a)/sizeof(KmerAffect));

253 254 255
  string seq = "AAAAAAAAAAAAAA";
  KmerAffectAnalyser kaa(*index, seq, affectations);
  TAP_TEST_EQUAL(kaa.getSequence(), seq, TEST_AA_GET_SEQUENCE, "");
256

257
  affect_infos results = kaa.getMaximum(AFFECT_J_BWD, AFFECT_V_BWD, 2., 0);
258 259 260 261
  TAP_TEST(! results.max_found, TEST_AA_GET_MAXIMUM_MAX_FOUND, 
           "(" << results.first_pos_max
           << ", " << results.last_pos_max << ")");

Mathieu Giraud's avatar
Mathieu Giraud committed
262
  results = kaa.getMaximum(AFFECT_J_BWD, AFFECT_V_BWD, 0.9, 0);
263 264 265 266 267 268
  TAP_TEST(results.max_found , 
           TEST_AA_GET_MAXIMUM_MAX_FOUND, "(" << results.first_pos_max
           << ", " << results.last_pos_max << ")");
  TAP_TEST(results.first_pos_max == 5 && results.last_pos_max == 5,
           TEST_AA_GET_MAXIMUM_POSITIONS,"(" << results.first_pos_max
           << ", " << results.last_pos_max << ")");
269
  TAP_TEST_EQUAL(results.max_value, 2, TEST_AA_GET_MAXIMUM_VALUE,
270 271
           "max = " << results.max_value);

Mathieu Giraud's avatar
Mathieu Giraud committed
272
  results = kaa.getMaximum(AFFECT_J_BWD, AFFECT_V_BWD, 0.9, k);
273 274 275 276 277 278
  TAP_TEST(results.max_found, 
           TEST_AA_GET_MAXIMUM_MAX_FOUND, "(" << results.first_pos_max
           << ", " << results.last_pos_max << ")");
  TAP_TEST(results.first_pos_max == 1 && results.last_pos_max == 5,
           TEST_AA_GET_MAXIMUM_POSITIONS,"(" << results.first_pos_max
           << ", " << results.last_pos_max << ")");
279
  TAP_TEST_EQUAL(results.max_value, 2, TEST_AA_GET_MAXIMUM_VALUE, "");
280

Mathieu Giraud's avatar
Mathieu Giraud committed
281
  affect_infos results2 = kaa.getMaximum(AFFECT_J_BWD, AFFECT_V_BWD, 0.9, k+5);
282 283
  TAP_TEST(results == results2, TEST_AA_GET_MAXIMUM_VALUE, "");

284 285 286 287 288
  KmerAffect a2[] = {AFFECT_V, AFFECT_V, AFFECT_V, AFFECT_V, AFFECT_V, 
                     AFFECT_V, AFFECT_V, AFFECT_V, AFFECT_V, AFFECT_V,
                     AFFECT_J, AFFECT_J, AFFECT_J,
                     AFFECT_V, AFFECT_V, AFFECT_V};
  //  0 1 2 3 4 5 6 7 8 9  11  13  15
289 290
  //  V V V V V V V V V V J J J V V V
  //                    ^^^^^^^^^^^^
291
  vector<KmerAffect> affectations2(a2, a2+sizeof(a2)/sizeof(KmerAffect));
292 293
  seq = "AAAAAAAAAAAAAAAA";
  KmerAffectAnalyser kaa2(*index, seq, affectations2);
294
  results = kaa2.getMaximum(AFFECT_V, AFFECT_J, 2., 0);
295 296 297 298 299 300 301 302
  TAP_TEST(! results.max_found, 
           TEST_AA_GET_MAXIMUM_MAX_FOUND, "(" << results.first_pos_max
           << ", " << results.last_pos_max << ")");

  results = kaa2.getMaximum(AFFECT_V, AFFECT_J, 1., k);
  TAP_TEST(! results.max_found, 
           TEST_AA_GET_MAXIMUM_MAX_FOUND, "(" << results.first_pos_max
           << ", " << results.last_pos_max << ")");
303
  TAP_TEST_EQUAL(results.max_value, 10, TEST_AA_GET_MAXIMUM_VALUE,
304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320
           "max = " << results.max_value);
  TAP_TEST(results.first_pos_max == 9
           && results.last_pos_max == 15, TEST_AA_GET_MAXIMUM_POSITIONS,
           "max positions: [" << results.first_pos_max << ", " 
           << results.last_pos_max << "]");
  TAP_TEST(results.nb_before_right == 0 && results.nb_after_right == 0, 
           TEST_AA_GET_MAXIMUM_COUNTS, 
           "before right: " << results.nb_before_right
           << ", after right: " << results.nb_after_right);

  results = kaa2.getMaximum(AFFECT_V_BWD, AFFECT_J_BWD);
  // No result

  TAP_TEST(! results.max_found,
           TEST_AA_GET_MAXIMUM_MAX_FOUND, 
           "(" << results.first_pos_max << ", " 
           << results.last_pos_max << ")");
321
  TAP_TEST_EQUAL(results.max_value, 0, TEST_AA_GET_MAXIMUM_VALUE,
322 323 324 325 326 327 328 329 330
           "max = " << results.max_value);

  
  KmerAffect a3[] = {AFFECT_V, AFFECT_V, AFFECT_V, AFFECT_V, AFFECT_V, 
                     AFFECT_V, AFFECT_V, AFFECT_V, AFFECT_V, AFFECT_V,
                     AFFECT_UNKNOWN, AFFECT_UNKNOWN, AFFECT_UNKNOWN,
                     AFFECT_UNKNOWN, AFFECT_UNKNOWN, AFFECT_UNKNOWN,
                     AFFECT_J, AFFECT_J, AFFECT_V, AFFECT_J};
  //  0 1 2 3 4 5 6 7 8 9  11  13  15  17  19
331
  //  V V V V V V V V V V _ _ _ _ _ _ J J V J
332 333
  //0 0 0 0 0 1 2 3 4 5 6 7 8 9101010 9 8 7 6
  //                           ^^^^^^
334 335
  //  V V V V V V                     J J   J    // nb_before_left (6), nb_after_right (3)
  //                                      V      // nb_before_right (1)
336
  vector<KmerAffect> affectations3(a3, a3+sizeof(a3)/sizeof(KmerAffect));
337 338
  seq = "AAAAAAAAAAAAAAAAAAAAA";
  KmerAffectAnalyser kaa3(*index, seq, affectations3);
339
  // span = 4, maxOverlap = 0
340
  results = kaa3.getMaximum(AFFECT_V, AFFECT_J, 2., 0);
341 342 343

  TAP_TEST(results.max_found, TEST_AA_GET_MAXIMUM_MAX_FOUND,
           "max_found = " << results.max_found);
344
  TAP_TEST(results.max_value == 10, TEST_AA_GET_MAXIMUM_VALUE,
345 346 347 348 349 350 351 352 353 354 355 356
           "max = " << results.max_value);
  TAP_TEST(results.first_pos_max == 13 && results.last_pos_max == 15,
           TEST_AA_GET_MAXIMUM_POSITIONS, 
           "first = " << results.first_pos_max 
           << ", last = " << results.last_pos_max);
  TAP_TEST(results.nb_before_left == 10 && results.nb_before_right == 1
           && results.nb_after_left == 0 && results.nb_after_right == 3,
           TEST_AA_GET_MAXIMUM_COUNTS, 
           "before:: left: " << results.nb_before_left <<", right: " 
           << results.nb_before_right << "\nafter:: left: " 
           << results.nb_after_left << ", right: " 
           << results.nb_after_right);
357

358 359 360 361 362

  KmerAffect a4[] = {AFFECT_V, AFFECT_V, AFFECT_V, AFFECT_V, AFFECT_V,
                     AFFECT_V, AFFECT_V, AFFECT_V, AFFECT_V, AFFECT_V,
                     AFFECT_J, AFFECT_J, AFFECT_J, AFFECT_J, AFFECT_J,
                     AFFECT_J, AFFECT_J, AFFECT_J, AFFECT_J, AFFECT_J};
363
  //  0 1 2 3 4 5 6 7 8 9  11  13  15  17  19    // i
364
  //  V V V V V V V V V V J J J J J J J J J J
365
  //0 0 0 0 0 1 2 3 4 5 6 6 6 6 6 5 4 3 2 1 0    // currentValue, after iteration i
366
  //                    ^^^^^^^^^
367
  //  V V V V V V                 J J J J J J    // nb_before_left (6), nb_after_right (6)
368
  vector<KmerAffect> affectations4(a4, a4+sizeof(a4)/sizeof(KmerAffect));
369
  KmerAffectAnalyser kaa4(*index, seq, affectations4);
370
  // span = 4, maxOverlap = 0
371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393
  results = kaa4.getMaximum(AFFECT_V, AFFECT_J, 2., 0);

  TAP_TEST(results.max_found, TEST_AA_GET_MAXIMUM_MAX_FOUND,
           "max_found = " << results.max_found);
  TAP_TEST(results.max_value == 6, TEST_AA_GET_MAXIMUM_VALUE,
           "max = " << results.max_value);
  TAP_TEST(results.first_pos_max == 9 && results.last_pos_max == 13,
           TEST_AA_GET_MAXIMUM_POSITIONS,
           "first = " << results.first_pos_max
           << ", last = " << results.last_pos_max);
  TAP_TEST(results.nb_before_left == 6 && results.nb_before_right == 0
           && results.nb_after_left == 0 && results.nb_after_right == 6,
           TEST_AA_GET_MAXIMUM_COUNTS,
           "before:: left: " << results.nb_before_left <<", right: "
           << results.nb_before_right << "\nafter:: left: "
           << results.nb_after_left << ", right: "
           << results.nb_after_right);


  KmerAffect a5[] = {AFFECT_V, AFFECT_V, AFFECT_V, AFFECT_V, AFFECT_J,
                     AFFECT_J, AFFECT_J, AFFECT_V, AFFECT_V, AFFECT_V,
                     AFFECT_V, AFFECT_V, AFFECT_J, AFFECT_J, AFFECT_J,
                     AFFECT_J, AFFECT_J, AFFECT_J, AFFECT_J, AFFECT_J};
394
  //  0 1 2 3 4 5 6 7 8 9  11  13  15  17  19   // i
395
  //  V V V V J J J V V V V V J J J J J J J J
396
  //0 0 0 0 0 0 0 0 1 1 1 1 2 2 2 2 2 1 0-1-2   // currentValue, after iteration i
397
  //                        ^^^^^^^^^
398 399 400
  //  V V V V       V                 J J J J   // nb_before_left (5), nb_after_right (4)
  //          J J J                             // nb_after_left (3)

401
  vector<KmerAffect> affectations5(a5, a5+sizeof(a5)/sizeof(KmerAffect));
402
  KmerAffectAnalyser kaa5(*index, seq, affectations5);
403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420
  results = kaa5.getMaximum(AFFECT_V, AFFECT_J, 2., 0);

  TAP_TEST(! results.max_found, TEST_AA_GET_MAXIMUM_MAX_FOUND,
           "max_found = " << results.max_found);
  TAP_TEST(results.max_value == 2, TEST_AA_GET_MAXIMUM_VALUE,
           "max = " << results.max_value);
  TAP_TEST(results.first_pos_max == 11 && results.last_pos_max == 15,
           TEST_AA_GET_MAXIMUM_POSITIONS,
           "first = " << results.first_pos_max
           << ", last = " << results.last_pos_max);
  TAP_TEST(results.nb_before_left == 5 && results.nb_before_right == 0
           && results.nb_after_left == 3 && results.nb_after_right == 4,
           TEST_AA_GET_MAXIMUM_COUNTS,
           "before:: left: " << results.nb_before_left <<", right: "
           << results.nb_before_right << "\nafter:: left: "
           << results.nb_after_left << ", right: "
           << results.nb_after_right);

421 422 423 424 425 426 427 428

  KmerAffect a6[] = {AFFECT_J_BWD, AFFECT_J_BWD, AFFECT_J_BWD,
                     AFFECT_J_BWD, AFFECT_J_BWD, AFFECT_J_BWD, AFFECT_J_BWD, AFFECT_J_BWD,
                     AFFECT_V_BWD, AFFECT_V_BWD, AFFECT_V_BWD,
                     AFFECT_V_BWD, AFFECT_V_BWD, AFFECT_J_BWD,
                     AFFECT_J_BWD, AFFECT_J_BWD, AFFECT_V_BWD, AFFECT_V_BWD, AFFECT_V_BWD, AFFECT_V_BWD};
  //  0 1 2 3 4 5 6 7 8 9  11  13  15  17  19
  // J-J-J-J-J-J-J-J-V-V-V-V-V-J-J-J-V-V-V-V-
429
  //               ^^^^^^^^^^
430
  vector<KmerAffect> affectations6(a6, a6+sizeof(a6)/sizeof(KmerAffect));
431
  KmerAffectAnalyser kaa6(*index, seq, affectations6);
432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449
  results = kaa6.getMaximum(AFFECT_J_BWD, AFFECT_V_BWD, 2., 0);

  TAP_TEST(! results.max_found, TEST_AA_GET_MAXIMUM_MAX_FOUND,
           "max_found = " << results.max_found);
  TAP_TEST(results.max_value == 4, TEST_AA_GET_MAXIMUM_VALUE,
           "max = " << results.max_value);
  TAP_TEST(results.first_pos_max == 7 && results.last_pos_max == 11,
           TEST_AA_GET_MAXIMUM_POSITIONS,
           "first = " << results.first_pos_max
           << ", last = " << results.last_pos_max);
  TAP_TEST(results.nb_before_left == 4 && results.nb_before_right == 3
           && results.nb_after_left == 0 && results.nb_after_right == 5,
           TEST_AA_GET_MAXIMUM_COUNTS,
           "before:: left: " << results.nb_before_left <<", right: "
           << results.nb_before_right << "\nafter:: left: "
           << results.nb_after_left << ", right: "
           << results.nb_after_right);

450
  delete index;
451 452
}

453 454 455
/**
 * A sequence and its revcomp are not affected in the same way.
 */
456
template<template <class> class T>
457
void testBugAffectAnalyser() {
458 459
  BioReader seqV("../../germline/homo-sapiens/IGHV.fa", 2);
  BioReader seqJ("../../germline/homo-sapiens/IGHJ.fa", 2);
460
  BioReader data("data/bug-revcomp.fa", 1, " ");
461

462 463
  int k = 9;
  T<KmerAffect> index(k, true);
464 465
  index.insert(seqV, "V");
  index.insert(seqJ, "J");
466
  index.finish_building();
467

468
  TAP_TEST_EQUAL(data.size(), 2, TEST_FASTA_SIZE, 
469 470 471 472 473 474 475
           "Should have 2 sequences (one seq and its revcomp), " 
           << data.size() << " instead");

  TAP_TEST(data.read(0).sequence.size() == data.read(1).sequence.size(),
           TEST_FASTA_SEQUENCE, 
           "Sequences should be of same length: sequence and its revcomp");

476 477
  KmerAffectAnalyser fwdAffect(index, data.read(0).sequence);
  KmerAffectAnalyser bwdAffect(index, data.read(1).sequence);
478

479
  int total = data.read(0).sequence.size() -k + 1;
480

481 482 483
  TAP_TEST(fwdAffect.getSequence() == data.read(0).sequence, TEST_AA_GET_SEQUENCE, "actual: " << fwdAffect.getSequence() << ", expected: " << data.read(0).sequence);
  TAP_TEST(bwdAffect.getSequence() == data.read(1).sequence, TEST_AA_GET_SEQUENCE, "actual: " << bwdAffect.getSequence() << ", expected: " << data.read(1).sequence);

484 485 486 487 488 489 490
  TAP_TEST(fwdAffect.count() == bwdAffect.count(),
           TEST_AA_COUNT,
           "Both sequences should have the same amount of affectations. "
           << fwdAffect.count() << " for the fwd, and " << bwdAffect.count()
           << " for the bwd instead");

  for (int i = 0; i < total; i++) {
491 492
    const KmerAffect fwd = fwdAffect.getAffectation(i);
    const KmerAffect bwd = bwdAffect.getAffectation(total - 1 - i);
493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512

    if (fwd.isAmbiguous() || bwd.isAmbiguous()
        || fwd.isUnknown() || bwd.isUnknown()) {
      TAP_TEST(fwd == bwd, TEST_AA_PREDICATES, 
               "If ambiguous or unknown, both affectations should be the same (fwd="
               << fwd << ", bwd=" << bwd << ", i= " << i << ") "
               << __PRETTY_FUNCTION__);
    } else {
      TAP_TEST(fwd.getLabel() == bwd.getLabel(), TEST_AA_REVCOMP_LABEL,
               "Label should be the same, instead: fwd=" << fwd.getLabel()
               << ", bwd=" << bwd.getLabel() << ", i=" <<i
               << " " << __PRETTY_FUNCTION__);
      TAP_TEST(-1*fwd.getStrand() == bwd.getStrand(), TEST_AA_REVCOMP_STRAND,
               "Strands should be the opposite, instead: fwd=" << fwd.getStrand()
               << ", bwd=" << bwd.getStrand() << ", i=" << i << " "
               << __PRETTY_FUNCTION__);
    }
  }
}

Mikaël Salson's avatar
Mikaël Salson committed
513
void testAffectAnalyser() {
514 515 516
  testAffectAnalyser1<MapKmerStore>();
  testAffectAnalyser2<MapKmerStore>();
  testBugAffectAnalyser<MapKmerStore>();
517
  testGetMaximum<MapKmerStore>();
518
  testAffectAnalyserSpaced<MapKmerStore>();
519

520 521
  testAffectAnalyser1<ArrayKmerStore>();
  testAffectAnalyser2<ArrayKmerStore>();
522
  testAffectAnalyserMaxes<ArrayKmerStore>();
523
  testBugAffectAnalyser<ArrayKmerStore>();
524
  testGetMaximum<ArrayKmerStore>();
525 526
  testAffectAnalyserSpaced<ArrayKmerStore>();

527 528 529 530 531

  testAffectAnalyser1<PointerACAutomaton>();
  testAffectAnalyser2<PointerACAutomaton>();
  testAffectAnalyserMaxes<PointerACAutomaton>();
  testBugAffectAnalyser<PointerACAutomaton>();
532
  testGetMaximum<PointerACAutomaton>();
533
  testAffectAnalyserSpaced<PointerACAutomaton>();
Mikaël Salson's avatar
Mikaël Salson committed
534
}