segment.h 7.13 KB
Newer Older
Mikaël Salson's avatar
Mikaël Salson committed
1 2 3 4
#ifndef SEGMENT_H
#define SEGMENT_H

#include <string>
5
#include <fstream>
Mikaël Salson's avatar
Mikaël Salson committed
6 7 8 9
#include <iostream>
#include "fasta.h"
#include "dynprog.h"
#include "tools.h"
10
#include "germline.h"
Mikaël Salson's avatar
Mikaël Salson committed
11 12
#include "kmerstore.h"
#include "kmeraffect.h"
Mathieu Giraud's avatar
Mathieu Giraud committed
13
#include "affectanalyser.h"
Marc Duez's avatar
Marc Duez committed
14
#include "../lib/json.hpp"
Mikaël Salson's avatar
Mikaël Salson committed
15 16

#define EXTEND_D_ZONE 5
17 18

#define MIN_D_LENGTH 5          /* If a D-REGION is smaller than this threshold, it is not output */
19
#define MIN_MATCHES 10          /* If a V/J-REGION does not give an alignment score with at least this number of matches, the FineSegmenter does not segment the sequence */
20

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

25 26 27 28
#define DETECT_THRESHOLD_STRAND 5   /* If the number of total affectations
                                       is above this threshold, then a sequence with no clearly attributed
                                       stranf will be marked as STRAND_NOT_CONSISTEN */

29
#define DETECT_THRESHOLD 5      /* If the number of both V and J affectations
30 31
                                   is above this threshold, then the sequence,
                                   if it not segmented, will be marked as AMBIGUOUS */
32

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

35
#define NO_LIMIT_VALUE  -1
36
#define BAD_EVALUE  1e10
37

38
#define THRESHOLD_NB_EXPECTED 1.0 /* Threshold of the accepted expected value for number of found k-mers */
39 40


41

Mikaël Salson's avatar
Mikaël Salson committed
42
using namespace std;
Marc Duez's avatar
Marc Duez committed
43
using json = nlohmann::json;
Mikaël Salson's avatar
Mikaël Salson committed
44

45 46 47 48
enum SEGMENTED { NOT_PROCESSED,
		 TOTAL_SEG_AND_WINDOW,
                 SEG_PLUS, SEG_MINUS,
                 UNSEG_TOO_SHORT, UNSEG_STRAND_NOT_CONSISTENT,
49
		 UNSEG_TOO_FEW_ZERO,  UNSEG_TOO_FEW_J, UNSEG_TOO_FEW_V,
50
                 UNSEG_BAD_DELTA_MIN, UNSEG_AMBIGUOUS,
51
		 UNSEG_TOO_SHORT_FOR_WINDOW,
52

Mikaël Salson's avatar
Mikaël Salson committed
53
		 STATS_SIZE } ;
54 55 56 57
const char* const segmented_mesg[] = { "?",
                                       "SEG",
                                       "SEG_+", "SEG_-",
                                       "UNSEG too short", "UNSEG strand",
58
				       "UNSEG too few V/J", "UNSEG only V", "UNSEG only J",
59
				       "UNSEG < delta_min", "UNSEG ambiguous",
60
                                       "UNSEG too short w",
Mikaël Salson's avatar
Mikaël Salson committed
61
                                      } ;
Mikaël Salson's avatar
Mikaël Salson committed
62 63 64 65

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

72 73 74 75 76 77
  /**
   * Compares evalue_left, evalue_right and evalue against the provided threshold
   * @post some evalue is above the threshold ==> because is set to UNSEG_TOO_FEW_ZERO, UNSEG_TOO_FEW_V or UNSEG_TOO_FEW_J
   */
  void checkLeftRightEvaluesThreshold(double threshold, int strand);

Mikaël Salson's avatar
Mikaël Salson committed
78 79 80 81 82
  string removeChevauchement();
  bool finishSegmentation();
  bool finishSegmentationD();

 public:
83
  Germline *segmented_germline;
84
  string label;
Mikaël Salson's avatar
Mikaël Salson committed
85
  string code;
Mikaël Salson's avatar
Mikaël Salson committed
86
  string code_short;
Mikaël Salson's avatar
Mikaël Salson committed
87
  string code_light;
88 89
  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
90 91
  int best_V, best_J ;
  int del_V, del_D_left, del_D_right, del_J ;
92
  string seg_V, seg_N, seg_J, system;
93

Mikaël Salson's avatar
Mikaël Salson committed
94
  int best_D;
95
  double evalue;
96 97
  double evalue_left;
  double evalue_right;
Mikaël Salson's avatar
Mikaël Salson committed
98 99 100
  string seg_N1, seg_D, seg_N2;
  Cost segment_cost;

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

Mathieu Giraud's avatar
Mathieu Giraud committed
103
  /* Queries */
Mikaël Salson's avatar
Mikaël Salson committed
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 142 143 144 145 146 147 148 149 150 151 152 153 154



  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;

155 156 157 158 159 160 161 162 163
  /**
   * @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
164

165 166 167 168 169
  /**
   * @post status == SEG_PLUS || status == SEG_MINUS <==> isSegmented()
   */
  void setSegmentationStatus(int status);

Mikaël Salson's avatar
Mikaël Salson committed
170 171 172 173 174 175 176 177 178 179 180
  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
181
 private:
182
  KmerAffectAnalyser *kaa;
Mikaël Salson's avatar
Mikaël Salson committed
183 184 185 186
 protected:
  string affects;

 public:
187
  int score;
188 189
  int pvalue_left;
  int pvalue_right;
190
  KmerAffect before, after;
191 192

  KmerSegmenter();
Mikaël Salson's avatar
Mikaël Salson committed
193 194 195
  /**
   * Build a segmenter based on KmerSegmentation
   * @param seq: An object read from a FASTA/FASTQ file
196
   * @param germline: the germline
Mikaël Salson's avatar
Mikaël Salson committed
197
   */
198
  KmerSegmenter(Sequence seq, Germline *germline, double threshold = THRESHOLD_NB_EXPECTED, int multiplier=1);
Mathieu Giraud's avatar
Mathieu Giraud committed
199

200 201
  KmerSegmenter(const KmerSegmenter &seg);

Mathieu Giraud's avatar
Mathieu Giraud committed
202 203 204 205 206
  ~KmerSegmenter();

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

209
  string getInfoLineWithAffects() const;
Marc Duez's avatar
Marc Duez committed
210
  json toJson();
211

Mathieu Giraud's avatar
Mathieu Giraud committed
212
 private:
213 214
  void computeSegmentation(int strand, KmerAffect left, KmerAffect right,
                           double threshold, int multiplier);
Mikaël Salson's avatar
Mikaël Salson committed
215 216
};

217 218 219

class KmerMultiSegmenter
{
220 221
 private:
  double threshold_nb_expected;
222 223 224
 public:
  /**
   * @param seq: An object read from a FASTA/FASTQ file
225 226 227 228
   * @param multigermline: the multigerm
   * @param threshold: threshold of randomly expected segmentation
   */
  KmerMultiSegmenter(Sequence seq, MultiGermline *multigermline, ostream *out_unsegmented,
229
                     double threshold = THRESHOLD_NB_EXPECTED, int nb_reads_for_evalue = 1);
230

231 232
  ~KmerMultiSegmenter();

233
  KmerSegmenter *the_kseg;
234
  MultiGermline *multi_germline;
235 236 237
};


Mikaël Salson's avatar
Mikaël Salson committed
238 239 240
class FineSegmenter : public Segmenter
{
 public:
Mikaël Salson's avatar
Mikaël Salson committed
241 242 243 244
   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
245 246 247
   /**
   * Build a fineSegmenter based on KmerSegmentation
   * @param seq: An object read from a FASTA/FASTQ file
248
   * @param germline: germline used
Mikaël Salson's avatar
Mikaël Salson committed
249
   */
250 251
   FineSegmenter(Sequence seq, Germline *germline, Cost segment_cost,
                 double threshold = THRESHOLD_NB_EXPECTED, int multiplier=1);
Mikaël Salson's avatar
Mikaël Salson committed
252 253 254
  
  /**
  * extend segmentation from VJ to VDJ
255
  * @param germline: germline used
Mikaël Salson's avatar
Mikaël Salson committed
256
  */
257
  void FineSegmentD(Germline *germline);
Marc Duez's avatar
Marc Duez committed
258
  void findCDR3();
Mathieu Giraud's avatar
Mathieu Giraud committed
259

Marc Duez's avatar
Marc Duez committed
260
  json toJson();
Mikaël Salson's avatar
Mikaël Salson committed
261
  
Mikaël Salson's avatar
Mikaël Salson committed
262 263 264 265 266
};



#endif