Attention une mise à jour du serveur va être effectuée le vendredi 16 avril entre 12h et 12h30. Cette mise à jour va générer une interruption du service de quelques minutes.

segment.h 8.39 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 RATIO_STRAND 2          /* The ratio between the affectations in one
Mikaël Salson's avatar
Mikaël Salson committed
19 20
                                   strand and the other, to safely attribute a
                                   segment to a given strand */
Mikaël Salson's avatar
Mikaël Salson committed
21

22 23 24 25
#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 */

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

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

32
#define BAD_EVALUE  1e10
33

34
#define THRESHOLD_NB_EXPECTED 1.0 /* Threshold of the accepted expected value for number of found k-mers */
35
#define THRESHOLD_NB_EXPECTED_D  .05 /* e-value threshold, D-REGION */
36

37 38 39 40
#define BOTTOM_TRIANGLE_SHIFT  20   /* Should equal to (max allowed 'k-band') + (max allowed number of V/J deletions) - (min size to recognize facing J/V)
                                       As we need ~10 bp to recognize the facing V/J, this value should be large enough to handle V/J deletions until ~30 bp,
                                       (and even larger V/J deletions if there is a large facing J or V in the read). */

41

42

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

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

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

64 65 66 67 68 69 70 71



 /**
  * Find the best split point when there may be a overlap between a left and a right region
  * @param overlap: the length of the overlap
  * @param seq_left, seq_right: the two sequences that overlap
  * @param ref_left, ref_right: the two reference sequences
72
  * @param pos_seq_left, pos_seq_right: the initial end/start position of two overlaping sequences
73
  * @param segment_cost: dp cost
74 75 76 77
  *
  * @post  trim_ref_left (at the end of ref_left) and trim_ref_right (at the beginning of ref_right)
  *        are set to the best number of nucleotides to trim in order to remove the overlap
  *        pos_seq_left and pos_seq_right are shifted by the good number of nucleotides (relatively to seq_{left,right})
78 79 80 81
  */

void best_overlap_split(int overlap, string seq_left, string seq_right,
                        string ref_left, string ref_right,
82 83
                        int *pos_seq_left, int *pos_seq_right,
                        int *trim_ref_left, int *trim_ref_right, Cost segment_cost);
84

Mikaël Salson's avatar
Mikaël Salson committed
85 86 87
class Segmenter {
protected:
  string sequence;
Mathieu Giraud's avatar
Mathieu Giraud committed
88 89
  int Vend, Jstart;
  int Dstart, Dend;
Marc Duez's avatar
Marc Duez committed
90
  int CDR3start, CDR3end;
Mikaël Salson's avatar
Mikaël Salson committed
91
  bool reversed, segmented, dSegmented;
92
  int because;
Mikaël Salson's avatar
Mikaël Salson committed
93

94 95 96 97 98 99
  /**
   * 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
100 101 102 103 104
  string removeChevauchement();
  bool finishSegmentation();
  bool finishSegmentationD();

 public:
105
  Germline *segmented_germline;
106
  string label;
Mikaël Salson's avatar
Mikaël Salson committed
107
  string code;
Mikaël Salson's avatar
Mikaël Salson committed
108
  string code_short;
Mikaël Salson's avatar
Mikaël Salson committed
109
  string code_light;
110 111
  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
112 113
  int best_V, best_J ;
  int del_V, del_D_left, del_D_right, del_J ;
114
  string seg_V, seg_N, seg_J, system;
115

Mikaël Salson's avatar
Mikaël Salson committed
116
  int best_D;
117
  double evalue;
118 119
  double evalue_left;
  double evalue_right;
Mikaël Salson's avatar
Mikaël Salson committed
120 121 122
  string seg_N1, seg_D, seg_N2;
  Cost segment_cost;

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

Mathieu Giraud's avatar
Mathieu Giraud committed
125
  /* Queries */
Mikaël Salson's avatar
Mikaël Salson committed
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 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176



  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;

177 178 179 180 181 182 183 184 185
  /**
   * @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
186

187 188 189 190 191
  /**
   * @post status == SEG_PLUS || status == SEG_MINUS <==> isSegmented()
   */
  void setSegmentationStatus(int status);

Mikaël Salson's avatar
Mikaël Salson committed
192 193 194 195 196 197 198 199 200 201 202
  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
203
 private:
204
  KmerAffectAnalyser *kaa;
Mikaël Salson's avatar
Mikaël Salson committed
205 206 207 208
 protected:
  string affects;

 public:
209
  int score;
210 211
  int pvalue_left;
  int pvalue_right;
212
  KmerAffect before, after;
213 214

  KmerSegmenter();
Mikaël Salson's avatar
Mikaël Salson committed
215 216 217
  /**
   * Build a segmenter based on KmerSegmentation
   * @param seq: An object read from a FASTA/FASTQ file
218
   * @param germline: the germline
Mikaël Salson's avatar
Mikaël Salson committed
219
   */
220
  KmerSegmenter(Sequence seq, Germline *germline, double threshold = THRESHOLD_NB_EXPECTED, int multiplier=1);
Mathieu Giraud's avatar
Mathieu Giraud committed
221

222 223
  KmerSegmenter(const KmerSegmenter &seg);

Mathieu Giraud's avatar
Mathieu Giraud committed
224 225 226 227 228
  ~KmerSegmenter();

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

231
  string getInfoLineWithAffects() const;
Marc Duez's avatar
Marc Duez committed
232
  json toJson();
233

Mathieu Giraud's avatar
Mathieu Giraud committed
234
 private:
235 236
  void computeSegmentation(int strand, KmerAffect left, KmerAffect right,
                           double threshold, int multiplier);
Mikaël Salson's avatar
Mikaël Salson committed
237 238
};

239 240 241

class KmerMultiSegmenter
{
242 243
 private:
  double threshold_nb_expected;
244 245 246
 public:
  /**
   * @param seq: An object read from a FASTA/FASTQ file
247 248 249 250
   * @param multigermline: the multigerm
   * @param threshold: threshold of randomly expected segmentation
   */
  KmerMultiSegmenter(Sequence seq, MultiGermline *multigermline, ostream *out_unsegmented,
251
                     double threshold = THRESHOLD_NB_EXPECTED, int nb_reads_for_evalue = 1);
252

253 254
  ~KmerMultiSegmenter();

255
  KmerSegmenter *the_kseg;
256
  MultiGermline *multi_germline;
257 258 259
};


Mikaël Salson's avatar
Mikaël Salson committed
260 261 262
class FineSegmenter : public Segmenter
{
 public:
Mikaël Salson's avatar
Mikaël Salson committed
263 264 265 266
   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
267 268 269
   /**
   * Build a fineSegmenter based on KmerSegmentation
   * @param seq: An object read from a FASTA/FASTQ file
270
   * @param germline: germline used
Mikaël Salson's avatar
Mikaël Salson committed
271
   */
272 273
   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
274 275 276
  
  /**
  * extend segmentation from VJ to VDJ
277
  * @param germline: germline used
Mikaël Salson's avatar
Mikaël Salson committed
278
  */
279 280
  void FineSegmentD(Germline *germline,
                    double threshold = THRESHOLD_NB_EXPECTED_D, int multiplier=1);
Marc Duez's avatar
Marc Duez committed
281
  void findCDR3();
Mathieu Giraud's avatar
Mathieu Giraud committed
282

Marc Duez's avatar
Marc Duez committed
283
  json toJson();
Mikaël Salson's avatar
Mikaël Salson committed
284
  
Mikaël Salson's avatar
Mikaël Salson committed
285 286 287 288 289
};



#endif