Mise à jour terminée. Pour connaître les apports de la version 13.8.4 par rapport à notre ancienne version vous pouvez lire les "Release Notes" suivantes :
https://about.gitlab.com/releases/2021/02/11/security-release-gitlab-13-8-4-released/
https://about.gitlab.com/releases/2021/02/05/gitlab-13-8-3-released/

segment.h 6.4 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 26 27
                                   is above this threshold, then the sequence
                                   will be labeled as 'detected', and, if it
                                   not segmented, the remaining germlines will
                                   not be tested */
28

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

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

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

enum SEGMENTED { DONT_KNOW, SEG_PLUS, SEG_MINUS, UNSEG_TOO_SHORT, UNSEG_STRAND_NOT_CONSISTENT, 
		 UNSEG_TOO_FEW_ZERO,  UNSEG_TOO_FEW_V, UNSEG_TOO_FEW_J, 
37
		 UNSEG_BAD_DELTA_MIN, UNSEG_BAD_DELTA_MAX, UNSEG_AMBIGUOUS, UNSEG_NOISY,
Mikaël Salson's avatar
Mikaël Salson committed
38 39 40 41
		 TOTAL_SEG_AND_WINDOW, 
		 TOTAL_SEG_BUT_TOO_SHORT_FOR_THE_WINDOW,
		 STATS_SIZE } ;
const char* const segmented_mesg[] = { "?", "SEG_+", "SEG_-", "UNSEG too short", "UNSEG strand",  
42
				       "UNSEG too few (0)", "UNSEG too few V", "UNSEG too few J",
43
				       "UNSEG < delta_min", "UNSEG > delta_max", "UNSEG ambiguous",
44
                                       "UNSEG noisy",
Mikaël Salson's avatar
Mikaël Salson committed
45
                                       "= SEG, with window",
46
                                       "= SEG, no window",
Mikaël Salson's avatar
Mikaël Salson committed
47
                                      } ;
Mikaël Salson's avatar
Mikaël Salson committed
48 49 50 51 52

class Segmenter {
protected:
  string label;
  string sequence;
Mathieu Giraud's avatar
Mathieu Giraud committed
53 54
  int Vend, Jstart;
  int Dstart, Dend;
Marc Duez's avatar
Marc Duez committed
55
  int CDR3start, CDR3end;
Mikaël Salson's avatar
Mikaël Salson committed
56
  bool reversed, segmented, dSegmented;
57
  int because;
Mikaël Salson's avatar
Mikaël Salson committed
58 59 60 61 62 63

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

 public:
64
  Germline *segmented_germline;
Mikaël Salson's avatar
Mikaël Salson committed
65
  string code;
Mikaël Salson's avatar
Mikaël Salson committed
66
  string code_short;
Mikaël Salson's avatar
Mikaël Salson committed
67
  string code_light;
68 69
  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
70 71
  int best_V, best_J ;
  int del_V, del_D_left, del_D_right, del_J ;
72
  string seg_V, seg_N, seg_J, system;
73

Mikaël Salson's avatar
Mikaël Salson committed
74 75 76 77
  int best_D;
  string seg_N1, seg_D, seg_N2;
  Cost segment_cost;

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

Mathieu Giraud's avatar
Mathieu Giraud committed
80
  /* Queries */
Mikaël Salson's avatar
Mikaël Salson committed
81 82 83 84 85 86 87 88 89 90 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



  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;

132 133 134 135 136 137 138 139 140
  /**
   * @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
141

142 143 144 145 146
  /**
   * @post status == SEG_PLUS || status == SEG_MINUS <==> isSegmented()
   */
  void setSegmentationStatus(int status);

Mikaël Salson's avatar
Mikaël Salson committed
147 148 149 150 151 152 153 154 155 156 157
  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
158
 private:
159
  int detected;
160
  KmerAffectAnalyser *kaa;
Mikaël Salson's avatar
Mikaël Salson committed
161 162 163 164
 protected:
  string affects;

 public:
165
  bool isDetected() const;
166
  int score;
167 168

  KmerSegmenter();
Mikaël Salson's avatar
Mikaël Salson committed
169 170 171
  /**
   * Build a segmenter based on KmerSegmentation
   * @param seq: An object read from a FASTA/FASTQ file
172
   * @param germline: the germline
Mikaël Salson's avatar
Mikaël Salson committed
173
   */
174
  KmerSegmenter(Sequence seq, Germline *germline);
Mathieu Giraud's avatar
Mathieu Giraud committed
175

176 177
  KmerSegmenter(const KmerSegmenter &seg);

Mathieu Giraud's avatar
Mathieu Giraud committed
178 179
  ~KmerSegmenter();

180 181 182 183 184
  /**
   * @return probability that the score of the segmenter is 'at_least' or more
   */
  double getProbabilityAtLeastOrAbove(int at_least) const;

Mathieu Giraud's avatar
Mathieu Giraud committed
185 186 187
  /**
   * @return the KmerAffectAnalyser of the current sequence.
   */
188
  KmerAffectAnalyser *getKmerAffectAnalyser() const;
Mathieu Giraud's avatar
Mathieu Giraud committed
189

190 191
  void toJsonList(JsonList *seg);

Mathieu Giraud's avatar
Mathieu Giraud committed
192
 private:
193
  void computeSegmentation(int strand, KmerAffect left, KmerAffect right);
Mikaël Salson's avatar
Mikaël Salson committed
194 195
};

196 197 198

class KmerMultiSegmenter
{
199 200
 private:
  double threshold_nb_expected;
201 202 203
 public:
  /**
   * @param seq: An object read from a FASTA/FASTQ file
204 205 206 207 208 209 210
   * @param multigermline: the multigerm
   * @param threshold: threshold of randomly expected segmentation
   */
  KmerMultiSegmenter(Sequence seq, MultiGermline *multigermline, ostream *out_unsegmented,
                     double threshold = THRESHOLD_NB_EXPECTED);

  /**
211
   * @return expected number of Segmenter that would have yield the maximum score by chance
212
   */
213 214
  double getNbExpected() const;

215 216
  ~KmerMultiSegmenter();

217
  KmerSegmenter *the_kseg;
218
  MultiGermline *multi_germline;
219 220 221
};


Mikaël Salson's avatar
Mikaël Salson committed
222 223 224
class FineSegmenter : public Segmenter
{
 public:
Mikaël Salson's avatar
Mikaël Salson committed
225 226 227 228
   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
229 230 231
   /**
   * Build a fineSegmenter based on KmerSegmentation
   * @param seq: An object read from a FASTA/FASTQ file
232
   * @param germline: germline used
Mikaël Salson's avatar
Mikaël Salson committed
233
   */
234
  FineSegmenter(Sequence seq, Germline *germline, Cost segment_cost);
Mikaël Salson's avatar
Mikaël Salson committed
235 236 237
  
  /**
  * extend segmentation from VJ to VDJ
238
  * @param germline: germline used
Mikaël Salson's avatar
Mikaël Salson committed
239
  */
240
  void FineSegmentD(Germline *germline);
Marc Duez's avatar
Marc Duez committed
241
  void findCDR3();
Mathieu Giraud's avatar
Mathieu Giraud committed
242

243
  void toJsonList(JsonList *seg);
Mikaël Salson's avatar
Mikaël Salson committed
244
  
Mikaël Salson's avatar
Mikaël Salson committed
245 246 247 248 249
};



#endif