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 12.2 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
#include <iostream>
7
#include "bioreader.hpp"
Mikaël Salson's avatar
Mikaël Salson committed
8 9
#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"
15
#include "filter.h"
Mikaël Salson's avatar
Mikaël Salson committed
16

17 18
// #define DEBUG_EVALUE

Mikaël Salson's avatar
Mikaël Salson committed
19
#define EXTEND_D_ZONE 5
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
#define DETECT_THRESHOLD_STRAND 5   /* If the number of total affectations
                                       is above this threshold, then a sequence with no clearly attributed
27
                                       strand will be marked as STRAND_NOT_CONSISTENT */
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 BAD_EVALUE  1e10
32

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

36 37 38 39
#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). */

40 41 42
#define DEFAULT_WINDOW_SHIFT 5  /* Number of nucleotide to try shifting or
                                   reducing the window when it doesn't fit at
                                   its place */
43
#define MINIMAL_WINDOW_LENGTH 30 /* As we now dynamically adapt the window
44 45
                                    length we need to specify a minimum
                                    otherwise we could go as low as
46 47 48
                                    2*DEFAULT_WINDOW_SHIFT. Of course this can
                                    be overriden by the command line by
                                    providing a shorter -w*/
49

50
#define DEFAULT_MINIMIZE_WIDTH       12   /* Number of nucleotides on which the minimization integer is computed */
51

52 53 54 55 56
#define DEFAULT_MINIMIZE_ONE_MARGIN  20   /* Number of nucleotides on ends of the reads excluded from the minimization in SEG_METHOD_ONE.
                                             - A read smaller than (DEFAULT_MINIMIZE_ONE_MARGIN*2 + max(k, DEFAULT_MINIMIZE_WIDTH)), that is 52,
                                               will trigger UNSEG_TOO_SHORT_FOR_WINDOW.
                                             - Margins above the half of the window length may produce shortened or shifted windows
                                          */
57 58


Mikaël Salson's avatar
Mikaël Salson committed
59
using namespace std;
Marc Duez's avatar
Marc Duez committed
60
using json = nlohmann::json;
Mikaël Salson's avatar
Mikaël Salson committed
61

62 63 64
enum SEGMENTED { NOT_PROCESSED,
		 TOTAL_SEG_AND_WINDOW,
                 SEG_PLUS, SEG_MINUS,
65
                 SEG_CHANGED_WINDOW,
66
                 UNSEG_TOO_SHORT, UNSEG_STRAND_NOT_CONSISTENT,
67
		 UNSEG_TOO_FEW_ZERO,  UNSEG_ONLY_V, UNSEG_ONLY_J,
68
                 UNSEG_BAD_DELTA_MIN, UNSEG_AMBIGUOUS,
69
		 UNSEG_TOO_SHORT_FOR_WINDOW,
70

Mikaël Salson's avatar
Mikaël Salson committed
71
		 STATS_SIZE } ;
72 73 74

#define STATS_FIRST_UNSEG UNSEG_TOO_SHORT

75 76 77
const char* const segmented_mesg[] = { "?",
                                       "SEG",
                                       "SEG_+", "SEG_-",
78
                                       "SEG changed w",
79
                                       "UNSEG too short", "UNSEG strand",
80
				       "UNSEG too few V/J", "UNSEG only V/5'", "UNSEG only J/3'",
81
				       "UNSEG < delta_min", "UNSEG ambiguous",
82
                                       "UNSEG too short w",
Mikaël Salson's avatar
Mikaël Salson committed
83
                                      } ;
Mikaël Salson's avatar
Mikaël Salson committed
84

85 86


87 88 89 90
/**
 * An alignment box (AlignBox) gather all parameters for a recombined gene segment (V, D, J, other D...)
 **/

91 92 93
class AlignBox
{
 public:
94
  string key;
95
  string color;
96

97 98 99 100 101
  int del_left;
  int start;
  int end;
  int del_right;

102
  AlignBox(string key = "", string color="");
103
  string getSequence(string sequence);
104
  void addToJson(json &seg);
105

106 107 108 109 110
  /**
   * Returns 'V', 'D', 'J', or possibly '5', '4', '3', '?', depending on the ref_label and on the key
   */
  char getInitial();

111 112 113 114 115
  /**
   * Returns the length
   */
  int getLength();

116 117 118 119 120 121 122 123 124 125 126 127
  /**
   * Returns the position in the reference string corresponding to the position in the read
   * Preliminary implementation, only works for the start of V and J boxes
   */
  int posInRef(int i);

  /**
   * Format the reference string in a given range position, possibly adding the ANSI colors
   * where there is the alignment
   */
  string refToString(int from, int to);

128
  /* Identifier, label and sequence of the reference sequence (the best one) */
129 130 131 132
  int ref_nb;
  string ref_label;
  string ref;

133 134 135
  /* Marked position, for Cys104 and Phe118/Trp118 */
  int marked_pos;

136 137
  /* Identifiers and scores of other possible reference sequence */
  vector<pair<int, int> > score;
138 139
};

140
ostream &operator<<(ostream &out, const AlignBox &box);
141

142
/**
143
 * Check whether there is an overlap between two boxes,
144 145 146
 * 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
147
 * @param *box_left, *box_right            the two boxes
148 149
 * @param segment_cost:                    the cost used by the dynamic programing
 *
150 151
 * @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
152 153 154 155 156
 *
 * @return                                 the N segment
 */

string check_and_resolve_overlap(string seq, int seq_begin, int seq_end,
157 158
                                 AlignBox *box_left, AlignBox *box_right,
                                 Cost segment_cost);
159

Mikaël Salson's avatar
Mikaël Salson committed
160 161 162
class Segmenter {
protected:
  string sequence;
163
  string sequence_or_rc;
164

165
  // JUNCTIONstart/end and CDR3start/end are 1-based
166
  int JUNCTIONstart, JUNCTIONend;
167
  string JUNCTIONaa;
168
  bool JUNCTIONproductive;
169

Marc Duez's avatar
Marc Duez committed
170
  int CDR3start, CDR3end;
171 172
  string CDR3nuc;
  string CDR3aa;
173

Mikaël Salson's avatar
Mikaël Salson committed
174
  bool reversed, segmented, dSegmented;
175
  bool junctionChanged;
176
  int because;
Mikaël Salson's avatar
Mikaël Salson committed
177

178 179
  /**
   * Compares evalue_left, evalue_right and evalue against the provided threshold
180
   * @post some evalue is above the threshold ==> because is set to UNSEG_TOO_FEW_ZERO, UNSEG_ONLY_J or UNSEG_ONLY_V
181 182 183
   */
  void checkLeftRightEvaluesThreshold(double threshold, int strand);

Mikaël Salson's avatar
Mikaël Salson committed
184 185 186 187 188
  string removeChevauchement();
  bool finishSegmentation();
  bool finishSegmentationD();

 public:
189
  Germline *segmented_germline;
190
  string label;
Mikaël Salson's avatar
Mikaël Salson committed
191
  string code;
192 193
  string info;        // .vdj.fa header, fixed fields
  string info_extra;  // .vdj.fa header, other information, at the end of the header
194
  string seg_V, seg_N, seg_J, system;
195

196 197
  AlignBox *box_V, *box_D, *box_J;

198
  double evalue;
199 200
  double evalue_left;
  double evalue_right;
Mikaël Salson's avatar
Mikaël Salson committed
201 202 203
  string seg_N1, seg_D, seg_N2;
  Cost segment_cost;

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

Mathieu Giraud's avatar
Mathieu Giraud committed
206
  /* Queries */
Mikaël Salson's avatar
Mikaël Salson committed
207 208 209 210 211 212 213 214 215 216 217 218 219



  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.
220
   * @post If the size or position of the window had to be dynamically adapted to fit
221
   *       in the read, isJunctionChanged() will return true
Mikaël Salson's avatar
Mikaël Salson committed
222
   */
223
  string getJunction(int l, int shift=0);
Mikaël Salson's avatar
Mikaël Salson committed
224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259

  /**
   * @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;

260 261 262 263
  /**
   * @return true iff the junction has been shifted or shortened
   *              dynamically to fit into the sequence
   */
264
  bool isJunctionChanged() const;
265

266 267 268 269 270 271 272 273 274
  /**
   * @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
275

276 277 278 279 280
  /**
   * @post status == SEG_PLUS || status == SEG_MINUS <==> isSegmented()
   */
  void setSegmentationStatus(int status);

Mikaël Salson's avatar
Mikaël Salson committed
281 282 283 284 285 286 287 288 289 290 291
  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
292
 private:
293
  KmerAffectAnalyser *kaa;
Mikaël Salson's avatar
Mikaël Salson committed
294 295 296 297
 protected:
  string affects;

 public:
298
  int score;
299 300
  int pvalue_left;
  int pvalue_right;
301
  KmerAffect before, after;
302 303

  KmerSegmenter();
Mikaël Salson's avatar
Mikaël Salson committed
304 305 306
  /**
   * Build a segmenter based on KmerSegmentation
   * @param seq: An object read from a FASTA/FASTQ file
307
   * @param germline: the germline
Mikaël Salson's avatar
Mikaël Salson committed
308
   */
309
  KmerSegmenter(Sequence seq, Germline *germline, double threshold = THRESHOLD_NB_EXPECTED, double multiplier=1.0);
Mathieu Giraud's avatar
Mathieu Giraud committed
310

311 312
  KmerSegmenter(const KmerSegmenter &seg);

Mathieu Giraud's avatar
Mathieu Giraud committed
313 314 315 316 317
  ~KmerSegmenter();

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

320
  string getInfoLineWithAffects() const;
Marc Duez's avatar
Marc Duez committed
321
  json toJson();
322

Mathieu Giraud's avatar
Mathieu Giraud committed
323
 private:
324
  void computeSegmentation(int strand, KmerAffect left, KmerAffect right,
325
                           double threshold, double multiplier);
Mikaël Salson's avatar
Mikaël Salson committed
326 327
};

328 329 330

class KmerMultiSegmenter
{
331 332
 private:
  double threshold_nb_expected;
333 334 335
 public:
  /**
   * @param seq: An object read from a FASTA/FASTQ file
336 337 338 339
   * @param multigermline: the multigerm
   * @param threshold: threshold of randomly expected segmentation
   */
  KmerMultiSegmenter(Sequence seq, MultiGermline *multigermline, ostream *out_unsegmented,
340
                     double threshold = THRESHOLD_NB_EXPECTED, int nb_reads_for_evalue = 1);
341

342 343
  ~KmerMultiSegmenter();

344
  KmerSegmenter *the_kseg;
345
  MultiGermline *multi_germline;
346 347 348
};


Mikaël Salson's avatar
Mikaël Salson committed
349 350 351
class FineSegmenter : public Segmenter
{
 public:
Mikaël Salson's avatar
Mikaël Salson committed
352 353 354
   vector<pair<int, int> > score_V;
   vector<pair<int, int> > score_D;
   vector<pair<int, int> > score_J;
355 356

   vector <AlignBox*> boxes ;
357

Mikaël Salson's avatar
Mikaël Salson committed
358 359 360
   /**
   * Build a fineSegmenter based on KmerSegmentation
   * @param seq: An object read from a FASTA/FASTQ file
361
   * @param germline: germline used
362 363
   * @param kmer_threshold: The number of the best k-mers for filtering.
   * By default this parameter doesn't filter the germline.
Mikaël Salson's avatar
Mikaël Salson committed
364
   */
365
   FineSegmenter(Sequence seq, Germline *germline, Cost segment_cost,
366 367
                 double threshold = THRESHOLD_NB_EXPECTED, double multiplier=1.0,
                int kmer_threshold=NO_LIMIT_VALUE);
368 369 370

   ~FineSegmenter();

Mikaël Salson's avatar
Mikaël Salson committed
371 372
  /**
  * extend segmentation from VJ to VDJ
373
  * @param germline: germline used
Mikaël Salson's avatar
Mikaël Salson committed
374
  */
375
  void FineSegmentD(Germline *germline, bool several_D,
376
                    double threshold = THRESHOLD_NB_EXPECTED_D, double multiplier=1.0);
377 378

  bool FineSegmentD(Germline *germline,
379
                    AlignBox *box_Y, AlignBox *box_DD, AlignBox *box_Z,
380
                    int forbidden_id,
381
                    int extend_DD_on_Y, int extend_DD_on_Z,
382
                    double threshold = THRESHOLD_NB_EXPECTED_D, double multiplier=1.0);
383

Mathieu Giraud's avatar
Mathieu Giraud committed
384 385 386 387
  /**
   * 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
388
  void findCDR3();
Mathieu Giraud's avatar
Mathieu Giraud committed
389

Marc Duez's avatar
Marc Duez committed
390
  json toJson();
Mikaël Salson's avatar
Mikaël Salson committed
391
  
Mikaël Salson's avatar
Mikaël Salson committed
392 393 394
};


395
void align_against_collection(string &read, BioReader &rep, int forbidden_rep_id,
396
                              bool reverse_ref, bool reverse_both, bool local,
397
                              AlignBox *box, Cost segment_cost, bool banded_dp=true);
Mikaël Salson's avatar
Mikaël Salson committed
398 399

#endif