KmerMultiSegmenter: Number of expected hits randomly.

Should take into account the total number of sequences but we don't know it at first.
Maybe the value should be computed afterwards…

algo/core/segment.cpp

@@ -282,12 +282,13 @@ KmerSegmenter::~KmerSegmenter() {
     delete kaa;
 }
 
-KmerMultiSegmenter::KmerMultiSegmenter(Sequence seq, MultiGermline *multigermline, ostream *out_unsegmented)
+KmerMultiSegmenter::KmerMultiSegmenter(Sequence seq, MultiGermline *multigermline, ostream *out_unsegmented, double threshold)
 {
   int best_score_seg = 0 ; // Best score, segmented sequences
   int best_score_unseg = 0 ; // Best score, unsegmented sequences
   the_kseg = NULL;
   multi_germline = multigermline;
+  threshold_nb_expected = threshold;
   
   // Iterate over the germlines
   for (list<Germline*>::const_iterator it = multigermline->germlines.begin(); it != multigermline->germlines.end(); ++it)
@@ -345,7 +346,19 @@ KmerMultiSegmenter::KmerMultiSegmenter(Sequence seq, MultiGermline *multigermlin
       } else {
         delete kseg;
       }
-    } // end for (Germlines)  
+    } // end for (Germlines)
+  if (the_kseg->isSegmented() && getNbExpected() > threshold_nb_expected)
+    the_kseg->setSegmentationStatus(UNSEG_NOISY);
+}
+
+double KmerMultiSegmenter::getNbExpected() const {
+  int max = the_kseg->score;
+  double proba = 0;
+  int n = the_kseg->getSequence().sequence.size() - the_kseg->getKmerAffectAnalyser()->getIndex().getS();
+  for (int i = max; i<  n; i++) {
+    proba += nChoosek(n, i) * pow(the_kseg->getKmerAffectAnalyser()->getIndex().getIndexLoad(), i);
+  }
+  return multi_germline->germlines.size() * proba;
 }
 
 KmerMultiSegmenter::~KmerMultiSegmenter() {

algo/core/segment.h

@@ -28,6 +28,8 @@
 
 #define JSON_REMEMBER_BEST  4   /* The number of V/D/J predictions to keep  */
 
+#define THRESHOLD_NB_EXPECTED 1E-6 /* Threshold of the accepted expected value for number of found k-mers */
+
 using namespace std;
 
 enum SEGMENTED { DONT_KNOW, SEG_PLUS, SEG_MINUS, UNSEG_TOO_SHORT, UNSEG_STRAND_NOT_CONSISTENT, 
@@ -189,12 +191,23 @@ class KmerSegmenter : public Segmenter
 
 class KmerMultiSegmenter
 {
+ private:
+  double threshold_nb_expected;
  public:
   /**
    * @param seq: An object read from a FASTA/FASTQ file
-   * @param multigermline: the multigermline
+   * @param multigermline: the multigerm
+   * @param threshold: threshold of randomly expected segmentation
+   */
+  KmerMultiSegmenter(Sequence seq, MultiGermline *multigermline, ostream *out_unsegmented,
+                     double threshold = THRESHOLD_NB_EXPECTED);
+
+  /**
+   * @return true iff the best score is sufficiently different from the scores
+   * of the other segmenters.
    */
-  KmerMultiSegmenter(Sequence seq, MultiGermline *multigermline, ostream *out_unsegmented);
+  double getNbExpected() const;
+
   ~KmerMultiSegmenter();
 
   KmerSegmenter *the_kseg;