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 9.87 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 27 28 29 30
#define DETECT_THRESHOLD 5      /* This threshold only affects the classification of unsegmented sequences.
                                   If the numbers of both V and J affectations are above this threshold,
                                   then the unsegmented sequence will be marked as AMBIGUOUS.
                                   Otherwise, if the number of either V or J affectations if above this threshold,
                                   then the unsegmented sequence will be marked as either ONLY V or ONLY J. */
31

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

34
#define BAD_EVALUE  1e10
35

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

39 40 41 42
#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). */

43

44

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

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

Mikaël Salson's avatar
Mikaël Salson committed
56
		 STATS_SIZE } ;
57 58 59

#define STATS_FIRST_UNSEG UNSEG_TOO_SHORT

60 61 62 63
const char* const segmented_mesg[] = { "?",
                                       "SEG",
                                       "SEG_+", "SEG_-",
                                       "UNSEG too short", "UNSEG strand",
64
				       "UNSEG too few V/J", "UNSEG only V/5'", "UNSEG only J/3'",
65
				       "UNSEG < delta_min", "UNSEG ambiguous",
66
                                       "UNSEG too short w",
Mikaël Salson's avatar
Mikaël Salson committed
67
                                      } ;
Mikaël Salson's avatar
Mikaël Salson committed
68

69 70


71 72 73 74
/**
 * An alignment box (AlignBox) gather all parameters for a recombined gene segment (V, D, J, other D...)
 **/

75 76 77
class AlignBox
{
 public:
78 79
  string key;

80 81 82 83 84
  int del_left;
  int start;
  int end;
  int del_right;

85
  AlignBox(string key = "");
86
  string getSequence(string sequence);
87
  void addToJson(json &seg);
88 89

  /* Identifier, label and sequence of the reference sequence (the best one) */
90 91 92 93
  int ref_nb;
  string ref_label;
  string ref;

94 95 96
  /* Marked position, for Cys104 and Phe118/Trp118 */
  int marked_pos;

97 98
  /* Identifiers and scores of other possible reference sequence */
  vector<pair<int, int> > score;
99 100
};

101
ostream &operator<<(ostream &out, const AlignBox &box);
102

103
/**
104
 * Check whether there is an overlap between two boxes,
105 106 107
 * If this is the case, fix this overlap (finding the best split point), and update segmentation accordingly
 * @param seq:                             the read
 * @param seq_begin, seq_end:              the positions to consider on 'seq' for the two sequences that may overlap
108
 * @param *box_left, *box_right            the two boxes
109 110
 * @param segment_cost:                    the cost used by the dynamic programing
 *
111 112
 * @post  box_left->del_left and box_right->del_right are set to the best number of nucleotides to trim in order to remove the overlap.
 *        box_left->end and box_right->start are shifted by the good number of nucleotides
113 114 115 116 117
 *
 * @return                                 the N segment
 */

string check_and_resolve_overlap(string seq, int seq_begin, int seq_end,
118 119
                                 AlignBox *box_left, AlignBox *box_right,
                                 Cost segment_cost);
120

Mikaël Salson's avatar
Mikaël Salson committed
121 122 123
class Segmenter {
protected:
  string sequence;
124
  string sequence_or_rc;
125 126

  int JUNCTIONstart, JUNCTIONend;
127
  string JUNCTIONaa;
128
  bool JUNCTIONproductive;
129

Marc Duez's avatar
Marc Duez committed
130
  int CDR3start, CDR3end;
131 132
  string CDR3nuc;
  string CDR3aa;
133

Mikaël Salson's avatar
Mikaël Salson committed
134
  bool reversed, segmented, dSegmented;
135
  int because;
Mikaël Salson's avatar
Mikaël Salson committed
136

137 138 139 140 141 142
  /**
   * 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
143 144 145 146 147
  string removeChevauchement();
  bool finishSegmentation();
  bool finishSegmentationD();

 public:
148
  Germline *segmented_germline;
149
  string label;
Mikaël Salson's avatar
Mikaël Salson committed
150
  string code;
151 152
  string info;        // .vdj.fa header, fixed fields
  string info_extra;  // .vdj.fa header, other information, at the end of the header
153
  string seg_V, seg_N, seg_J, system;
154

155 156
  AlignBox *box_V, *box_D, *box_J;

157
  double evalue;
158 159
  double evalue_left;
  double evalue_right;
Mikaël Salson's avatar
Mikaël Salson committed
160 161 162
  string seg_N1, seg_D, seg_N2;
  Cost segment_cost;

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

Mathieu Giraud's avatar
Mathieu Giraud committed
165
  /* Queries */
Mikaël Salson's avatar
Mikaël Salson committed
166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216



  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;

217 218 219 220 221 222 223 224 225
  /**
   * @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
226

227 228 229 230 231
  /**
   * @post status == SEG_PLUS || status == SEG_MINUS <==> isSegmented()
   */
  void setSegmentationStatus(int status);

Mikaël Salson's avatar
Mikaël Salson committed
232 233 234 235 236 237 238 239 240 241 242
  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
243
 private:
244
  KmerAffectAnalyser *kaa;
Mikaël Salson's avatar
Mikaël Salson committed
245 246 247 248
 protected:
  string affects;

 public:
249
  int score;
250 251
  int pvalue_left;
  int pvalue_right;
252
  KmerAffect before, after;
253 254

  KmerSegmenter();
Mikaël Salson's avatar
Mikaël Salson committed
255 256 257
  /**
   * Build a segmenter based on KmerSegmentation
   * @param seq: An object read from a FASTA/FASTQ file
258
   * @param germline: the germline
Mikaël Salson's avatar
Mikaël Salson committed
259
   */
260
  KmerSegmenter(Sequence seq, Germline *germline, double threshold = THRESHOLD_NB_EXPECTED, int multiplier=1);
Mathieu Giraud's avatar
Mathieu Giraud committed
261

262 263
  KmerSegmenter(const KmerSegmenter &seg);

Mathieu Giraud's avatar
Mathieu Giraud committed
264 265 266 267 268
  ~KmerSegmenter();

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

271
  string getInfoLineWithAffects() const;
Marc Duez's avatar
Marc Duez committed
272
  json toJson();
273

Mathieu Giraud's avatar
Mathieu Giraud committed
274
 private:
275 276
  void computeSegmentation(int strand, KmerAffect left, KmerAffect right,
                           double threshold, int multiplier);
Mikaël Salson's avatar
Mikaël Salson committed
277 278
};

279 280 281

class KmerMultiSegmenter
{
282 283
 private:
  double threshold_nb_expected;
284 285 286
 public:
  /**
   * @param seq: An object read from a FASTA/FASTQ file
287 288 289 290
   * @param multigermline: the multigerm
   * @param threshold: threshold of randomly expected segmentation
   */
  KmerMultiSegmenter(Sequence seq, MultiGermline *multigermline, ostream *out_unsegmented,
291
                     double threshold = THRESHOLD_NB_EXPECTED, int nb_reads_for_evalue = 1);
292

293 294
  ~KmerMultiSegmenter();

295
  KmerSegmenter *the_kseg;
296
  MultiGermline *multi_germline;
297 298 299
};


Mikaël Salson's avatar
Mikaël Salson committed
300 301 302
class FineSegmenter : public Segmenter
{
 public:
Mikaël Salson's avatar
Mikaël Salson committed
303 304 305
   vector<pair<int, int> > score_V;
   vector<pair<int, int> > score_D;
   vector<pair<int, int> > score_J;
306 307

   vector <AlignBox*> boxes ;
Mikaël Salson's avatar
Mikaël Salson committed
308
   
Mikaël Salson's avatar
Mikaël Salson committed
309 310 311
   /**
   * Build a fineSegmenter based on KmerSegmentation
   * @param seq: An object read from a FASTA/FASTQ file
312
   * @param germline: germline used
Mikaël Salson's avatar
Mikaël Salson committed
313
   */
314 315
   FineSegmenter(Sequence seq, Germline *germline, Cost segment_cost,
                 double threshold = THRESHOLD_NB_EXPECTED, int multiplier=1);
316 317 318

   ~FineSegmenter();

Mikaël Salson's avatar
Mikaël Salson committed
319 320
  /**
  * extend segmentation from VJ to VDJ
321
  * @param germline: germline used
Mikaël Salson's avatar
Mikaël Salson committed
322
  */
323
  void FineSegmentD(Germline *germline, bool several_D,
324
                    double threshold = THRESHOLD_NB_EXPECTED_D, int multiplier=1);
325 326

  bool FineSegmentD(Germline *germline,
327
                    AlignBox *box_Y, AlignBox *box_DD, AlignBox *box_Z,
328
                    int forbidden_id,
329
                    int extend_DD_on_Y, int extend_DD_on_Z,
330 331
                    double threshold = THRESHOLD_NB_EXPECTED_D, int multiplier=1);

Mathieu Giraud's avatar
Mathieu Giraud committed
332 333 334 335
  /**
   * find JUNCTION/CDR3, by using marked Cys104 and Phe118/Trp118 positions
   * in the germline V and J genes and the backtrack of the DP matrix
   */
Marc Duez's avatar
Marc Duez committed
336
  void findCDR3();
Mathieu Giraud's avatar
Mathieu Giraud committed
337

Marc Duez's avatar
Marc Duez committed
338
  json toJson();
Mikaël Salson's avatar
Mikaël Salson committed
339
  
Mikaël Salson's avatar
Mikaël Salson committed
340 341 342 343 344
};



#endif