segment.h 6.44 KB
Newer Older
Mikaël Salson's avatar
Mikaël Salson committed
1 2 3 4 5 6 7 8
#ifndef SEGMENT_H
#define SEGMENT_H

#include <string>
#include <iostream>
#include "fasta.h"
#include "dynprog.h"
#include "tools.h"
9
#include "germline.h"
Mikaël Salson's avatar
Mikaël Salson committed
10 11
#include "kmerstore.h"
#include "kmeraffect.h"
Mathieu Giraud's avatar
Mathieu Giraud committed
12 13
#include "affectanalyser.h"
#include "json.h"
Mikaël Salson's avatar
Mikaël Salson committed
14 15

#define EXTEND_D_ZONE 5
16 17 18

#define MIN_D_LENGTH 5          /* If a D-REGION is smaller than this threshold, it is not output */

Mikaël Salson's avatar
Mikaël Salson committed
19 20 21
#define RATIO_STRAND 5          /* The ratio between the affectations in one
                                   strand and the other, to safely attribute a
                                   segment to a given strand */
Mikaël Salson's avatar
Mikaël Salson committed
22

23
#define DETECT_THRESHOLD 5      /* If the number of both V and J affectations
24 25
                                   is above this threshold, then the sequence,
                                   if it not segmented, will be marked as AMBIGUOUS */
26

Mathieu Giraud's avatar
Mathieu Giraud committed
27 28
#define JSON_REMEMBER_BEST  4   /* The number of V/D/J predictions to keep  */

29 30
#define NO_LIMIT_VALUE  -1

31
#define THRESHOLD_NB_EXPECTED 1.0 /* Threshold of the accepted expected value for number of found k-mers */
32 33


34

Mikaël Salson's avatar
Mikaël Salson committed
35 36
using namespace std;

37 38 39 40
enum SEGMENTED { NOT_PROCESSED,
		 TOTAL_SEG_AND_WINDOW,
                 SEG_PLUS, SEG_MINUS,
                 UNSEG_TOO_SHORT, UNSEG_STRAND_NOT_CONSISTENT,
Mikaël Salson's avatar
Mikaël Salson committed
41
		 UNSEG_TOO_FEW_ZERO,  UNSEG_TOO_FEW_V, UNSEG_TOO_FEW_J, 
42
		 UNSEG_BAD_DELTA_MIN, UNSEG_BAD_DELTA_MAX, UNSEG_AMBIGUOUS, UNSEG_NOISY,
43
		 UNSEG_TOO_SHORT_FOR_WINDOW,
44

Mikaël Salson's avatar
Mikaël Salson committed
45
		 STATS_SIZE } ;
46 47 48 49
const char* const segmented_mesg[] = { "?",
                                       "SEG",
                                       "SEG_+", "SEG_-",
                                       "UNSEG too short", "UNSEG strand",
50
				       "UNSEG too few (0)", "UNSEG too few V", "UNSEG too few J",
51
				       "UNSEG < delta_min", "UNSEG > delta_max", "UNSEG ambiguous",
52
                                       "UNSEG noisy",
53
                                       "UNSEG too short w",
Mikaël Salson's avatar
Mikaël Salson committed
54
                                      } ;
Mikaël Salson's avatar
Mikaël Salson committed
55 56 57 58

class Segmenter {
protected:
  string sequence;
Mathieu Giraud's avatar
Mathieu Giraud committed
59 60
  int Vend, Jstart;
  int Dstart, Dend;
Marc Duez's avatar
Marc Duez committed
61
  int CDR3start, CDR3end;
Mikaël Salson's avatar
Mikaël Salson committed
62
  bool reversed, segmented, dSegmented;
63
  int because;
Mikaël Salson's avatar
Mikaël Salson committed
64 65 66 67 68 69

  string removeChevauchement();
  bool finishSegmentation();
  bool finishSegmentationD();

 public:
70
  Germline *segmented_germline;
71
  string label;
Mikaël Salson's avatar
Mikaël Salson committed
72
  string code;
Mikaël Salson's avatar
Mikaël Salson committed
73
  string code_short;
Mikaël Salson's avatar
Mikaël Salson committed
74
  string code_light;
75 76
  string info;        // .vdj.fa header, fixed fields
  string info_extra;  // .vdj.fa header, other information, at the end of the header
Mikaël Salson's avatar
Mikaël Salson committed
77 78
  int best_V, best_J ;
  int del_V, del_D_left, del_D_right, del_J ;
79
  string seg_V, seg_N, seg_J, system;
80

Mikaël Salson's avatar
Mikaël Salson committed
81
  int best_D;
82
  double evalue;
83 84
  double evalue_left;
  double evalue_right;
Mikaël Salson's avatar
Mikaël Salson committed
85 86 87
  string seg_N1, seg_D, seg_N2;
  Cost segment_cost;

Mathieu Giraud's avatar
Mathieu Giraud committed
88
  virtual ~Segmenter();
Mikaël Salson's avatar
Mikaël Salson committed
89

Mathieu Giraud's avatar
Mathieu Giraud committed
90
  /* Queries */
Mikaël Salson's avatar
Mikaël Salson committed
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 130 131 132 133 134 135 136 137 138 139 140 141



  Sequence getSequence() const ;

  /**
   * @param l: length around the junction
   * @return the string centered on the junction (ie. at position
   *         (getLeft() + getRight())/2).
   *         The string has length l unless the original string 
   *         is not long enough.
   *         The junction is revcomp-ed if the original string comes from reverse
   *         strand.
   */
  string getJunction(int l) const;

  /**
   * @return the left position (on forward strand) of the segmentation.
   */
  int getLeft() const;
  
  /**
   * @return the right position (on forward strand) of the segmentation
   */
  int getRight() const;
  
  /**
   * @return the left position (on forward strand) of the D segmentation.
   */
  int getLeftD() const;
  
  /**
   * @return the right position (on forward strand) of the D segmentation
   */
  int getRightD() const;

  /**
   * @return true iff the string comes from reverse strand
   */
  bool isReverse() const;

  /**
   * @return true iff the sequence has been successfully segmented
   */
  bool isSegmented() const;
  
  /**
   * @return true if a D gene was found in the N region
   */
  bool isDSegmented() const;

142 143 144 145 146 147 148 149 150
  /**
   * @return the status of the segmentation. Tells if the Sequence has been segmented
   *         of if it has not, what the reason is.
   * @assert getSegmentationStatus() == SEG_PLUS || getSegmentationStatus() == SEG_MINUS
   *         <==> isSegmented()
   */
  int getSegmentationStatus() const;

  string getInfoLine() const;
Mikaël Salson's avatar
Mikaël Salson committed
151

152 153 154 155 156
  /**
   * @post status == SEG_PLUS || status == SEG_MINUS <==> isSegmented()
   */
  void setSegmentationStatus(int status);

Mikaël Salson's avatar
Mikaël Salson committed
157 158 159 160 161 162 163 164 165 166 167
  friend ostream &operator<<(ostream &out, const Segmenter &s);
};



ostream &operator<<(ostream &out, const Segmenter &s);



class KmerSegmenter : public Segmenter
{
Mathieu Giraud's avatar
Mathieu Giraud committed
168
 private:
169
  KmerAffectAnalyser *kaa;
Mikaël Salson's avatar
Mikaël Salson committed
170 171 172 173
 protected:
  string affects;

 public:
174
  int score;
175 176
  int pvalue_left;
  int pvalue_right;
177
  KmerAffect before, after;
178 179

  KmerSegmenter();
Mikaël Salson's avatar
Mikaël Salson committed
180 181 182
  /**
   * Build a segmenter based on KmerSegmentation
   * @param seq: An object read from a FASTA/FASTQ file
183
   * @param germline: the germline
Mikaël Salson's avatar
Mikaël Salson committed
184
   */
185
  KmerSegmenter(Sequence seq, Germline *germline, double threshold = THRESHOLD_NB_EXPECTED, int multiplier=1);
Mathieu Giraud's avatar
Mathieu Giraud committed
186

187 188
  KmerSegmenter(const KmerSegmenter &seg);

Mathieu Giraud's avatar
Mathieu Giraud committed
189 190 191 192 193
  ~KmerSegmenter();

  /**
   * @return the KmerAffectAnalyser of the current sequence.
   */
194
  KmerAffectAnalyser *getKmerAffectAnalyser() const;
Mathieu Giraud's avatar
Mathieu Giraud committed
195

196
  string getInfoLineWithAffects() const;
197 198
  void toJsonList(JsonList *seg);

Mathieu Giraud's avatar
Mathieu Giraud committed
199
 private:
200 201 202
  void computeSegmentation(int strand, KmerAffect left, KmerAffect right,
                           int delta_min, int delta_max,
                           double threshold, int multiplier);
Mikaël Salson's avatar
Mikaël Salson committed
203 204
};

205 206 207

class KmerMultiSegmenter
{
208 209
 private:
  double threshold_nb_expected;
210 211 212
 public:
  /**
   * @param seq: An object read from a FASTA/FASTQ file
213 214 215 216
   * @param multigermline: the multigerm
   * @param threshold: threshold of randomly expected segmentation
   */
  KmerMultiSegmenter(Sequence seq, MultiGermline *multigermline, ostream *out_unsegmented,
217
                     double threshold = THRESHOLD_NB_EXPECTED, int nb_reads_for_evalue = 1);
218

219 220
  ~KmerMultiSegmenter();

221
  KmerSegmenter *the_kseg;
222
  MultiGermline *multi_germline;
223 224 225
};


Mikaël Salson's avatar
Mikaël Salson committed
226 227 228
class FineSegmenter : public Segmenter
{
 public:
Mikaël Salson's avatar
Mikaël Salson committed
229 230 231 232
   vector<pair<int, int> > score_V;
   vector<pair<int, int> > score_D;
   vector<pair<int, int> > score_J;
   
Mikaël Salson's avatar
Mikaël Salson committed
233 234 235
   /**
   * Build a fineSegmenter based on KmerSegmentation
   * @param seq: An object read from a FASTA/FASTQ file
236
   * @param germline: germline used
Mikaël Salson's avatar
Mikaël Salson committed
237
   */
238
  FineSegmenter(Sequence seq, Germline *germline, Cost segment_cost);
Mikaël Salson's avatar
Mikaël Salson committed
239 240 241
  
  /**
  * extend segmentation from VJ to VDJ
242
  * @param germline: germline used
Mikaël Salson's avatar
Mikaël Salson committed
243
  */
244
  void FineSegmentD(Germline *germline);
Marc Duez's avatar
Marc Duez committed
245
  void findCDR3();
Mathieu Giraud's avatar
Mathieu Giraud committed
246

247
  void toJsonList(JsonList *seg);
Mikaël Salson's avatar
Mikaël Salson committed
248
  
Mikaël Salson's avatar
Mikaël Salson committed
249 250 251 252 253
};



#endif