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 10.6 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 17
// #define DEBUG_EVALUE

Mikaël Salson's avatar
Mikaël Salson committed
18
#define EXTEND_D_ZONE 5
19

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

24 25 26 27
#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 */

28 29 30 31 32
#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. */
33

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

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
#define THRESHOLD_NB_EXPECTED_D  .05 /* e-value threshold, D-REGION */
40

41 42 43 44
#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). */

45

46

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

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

Mikaël Salson's avatar
Mikaël Salson committed
58
		 STATS_SIZE } ;
59 60 61

#define STATS_FIRST_UNSEG UNSEG_TOO_SHORT

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

71 72


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

77 78 79
class AlignBox
{
 public:
80
  string key;
81
  string color;
82

83 84 85 86 87
  int del_left;
  int start;
  int end;
  int del_right;

88
  AlignBox(string key = "", string color="");
89
  string getSequence(string sequence);
90
  void addToJson(json &seg);
91

92 93 94 95 96 97 98 99 100 101 102 103
  /**
   * 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);

104
  /* Identifier, label and sequence of the reference sequence (the best one) */
105 106 107 108
  int ref_nb;
  string ref_label;
  string ref;

109 110 111
  /* Marked position, for Cys104 and Phe118/Trp118 */
  int marked_pos;

112 113
  /* Identifiers and scores of other possible reference sequence */
  vector<pair<int, int> > score;
114 115
};

116
ostream &operator<<(ostream &out, const AlignBox &box);
117

118
/**
119
 * Check whether there is an overlap between two boxes,
120 121 122
 * 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
123
 * @param *box_left, *box_right            the two boxes
124 125
 * @param segment_cost:                    the cost used by the dynamic programing
 *
126 127
 * @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
128 129 130 131 132
 *
 * @return                                 the N segment
 */

string check_and_resolve_overlap(string seq, int seq_begin, int seq_end,
133 134
                                 AlignBox *box_left, AlignBox *box_right,
                                 Cost segment_cost);
135

Mikaël Salson's avatar
Mikaël Salson committed
136 137 138
class Segmenter {
protected:
  string sequence;
139
  string sequence_or_rc;
140

141
  // JUNCTIONstart/end and CDR3start/end are 1-based
142
  int JUNCTIONstart, JUNCTIONend;
143
  string JUNCTIONaa;
144
  bool JUNCTIONproductive;
145

Marc Duez's avatar
Marc Duez committed
146
  int CDR3start, CDR3end;
147 148
  string CDR3nuc;
  string CDR3aa;
149

Mikaël Salson's avatar
Mikaël Salson committed
150
  bool reversed, segmented, dSegmented;
151
  int because;
Mikaël Salson's avatar
Mikaël Salson committed
152

153 154 155 156 157 158
  /**
   * 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
159 160 161 162 163
  string removeChevauchement();
  bool finishSegmentation();
  bool finishSegmentationD();

 public:
164
  Germline *segmented_germline;
165
  string label;
Mikaël Salson's avatar
Mikaël Salson committed
166
  string code;
167 168
  string info;        // .vdj.fa header, fixed fields
  string info_extra;  // .vdj.fa header, other information, at the end of the header
169
  string seg_V, seg_N, seg_J, system;
170

171 172
  AlignBox *box_V, *box_D, *box_J;

173
  double evalue;
174 175
  double evalue_left;
  double evalue_right;
Mikaël Salson's avatar
Mikaël Salson committed
176 177 178
  string seg_N1, seg_D, seg_N2;
  Cost segment_cost;

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

Mathieu Giraud's avatar
Mathieu Giraud committed
181
  /* Queries */
Mikaël Salson's avatar
Mikaël Salson committed
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 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232



  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;

233 234 235 236 237 238 239 240 241
  /**
   * @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
242

243 244 245 246 247
  /**
   * @post status == SEG_PLUS || status == SEG_MINUS <==> isSegmented()
   */
  void setSegmentationStatus(int status);

Mikaël Salson's avatar
Mikaël Salson committed
248 249 250 251 252 253 254 255 256 257 258
  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
259
 private:
260
  KmerAffectAnalyser *kaa;
Mikaël Salson's avatar
Mikaël Salson committed
261 262 263 264
 protected:
  string affects;

 public:
265
  int score;
266 267
  int pvalue_left;
  int pvalue_right;
268
  KmerAffect before, after;
269 270

  KmerSegmenter();
Mikaël Salson's avatar
Mikaël Salson committed
271 272 273
  /**
   * Build a segmenter based on KmerSegmentation
   * @param seq: An object read from a FASTA/FASTQ file
274
   * @param germline: the germline
Mikaël Salson's avatar
Mikaël Salson committed
275
   */
276
  KmerSegmenter(Sequence seq, Germline *germline, double threshold = THRESHOLD_NB_EXPECTED, double multiplier=1.0);
Mathieu Giraud's avatar
Mathieu Giraud committed
277

278 279
  KmerSegmenter(const KmerSegmenter &seg);

Mathieu Giraud's avatar
Mathieu Giraud committed
280 281 282 283 284
  ~KmerSegmenter();

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

287
  string getInfoLineWithAffects() const;
Marc Duez's avatar
Marc Duez committed
288
  json toJson();
289

Mathieu Giraud's avatar
Mathieu Giraud committed
290
 private:
291
  void computeSegmentation(int strand, KmerAffect left, KmerAffect right,
292
                           double threshold, double multiplier);
Mikaël Salson's avatar
Mikaël Salson committed
293 294
};

295 296 297

class KmerMultiSegmenter
{
298 299
 private:
  double threshold_nb_expected;
300 301 302
 public:
  /**
   * @param seq: An object read from a FASTA/FASTQ file
303 304 305 306
   * @param multigermline: the multigerm
   * @param threshold: threshold of randomly expected segmentation
   */
  KmerMultiSegmenter(Sequence seq, MultiGermline *multigermline, ostream *out_unsegmented,
307
                     double threshold = THRESHOLD_NB_EXPECTED, int nb_reads_for_evalue = 1);
308

309 310
  ~KmerMultiSegmenter();

311
  KmerSegmenter *the_kseg;
312
  MultiGermline *multi_germline;
313 314 315
};


Mikaël Salson's avatar
Mikaël Salson committed
316 317 318
class FineSegmenter : public Segmenter
{
 public:
Mikaël Salson's avatar
Mikaël Salson committed
319 320 321
   vector<pair<int, int> > score_V;
   vector<pair<int, int> > score_D;
   vector<pair<int, int> > score_J;
322 323

   vector <AlignBox*> boxes ;
Mikaël Salson's avatar
Mikaël Salson committed
324
   
Mikaël Salson's avatar
Mikaël Salson committed
325 326 327
   /**
   * Build a fineSegmenter based on KmerSegmentation
   * @param seq: An object read from a FASTA/FASTQ file
328
   * @param germline: germline used
Mikaël Salson's avatar
Mikaël Salson committed
329
   */
330
   FineSegmenter(Sequence seq, Germline *germline, Cost segment_cost,
331
                 double threshold = THRESHOLD_NB_EXPECTED, double multiplier=1.0);
332 333 334

   ~FineSegmenter();

Mikaël Salson's avatar
Mikaël Salson committed
335 336
  /**
  * extend segmentation from VJ to VDJ
337
  * @param germline: germline used
Mikaël Salson's avatar
Mikaël Salson committed
338
  */
339
  void FineSegmentD(Germline *germline, bool several_D,
340
                    double threshold = THRESHOLD_NB_EXPECTED_D, double multiplier=1.0);
341 342

  bool FineSegmentD(Germline *germline,
343
                    AlignBox *box_Y, AlignBox *box_DD, AlignBox *box_Z,
344
                    int forbidden_id,
345
                    int extend_DD_on_Y, int extend_DD_on_Z,
346
                    double threshold = THRESHOLD_NB_EXPECTED_D, double multiplier=1.0);
347

Mathieu Giraud's avatar
Mathieu Giraud committed
348 349 350 351
  /**
   * 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
352
  void findCDR3();
Mathieu Giraud's avatar
Mathieu Giraud committed
353

Marc Duez's avatar
Marc Duez committed
354
  json toJson();
Mikaël Salson's avatar
Mikaël Salson committed
355
  
Mikaël Salson's avatar
Mikaël Salson committed
356 357 358
};


359 360 361
void align_against_collection(string &read, Fasta &rep, int forbidden_rep_id,
                              bool reverse_ref, bool reverse_both, bool local,
                              AlignBox *box, Cost segment_cost);
Mikaël Salson's avatar
Mikaël Salson committed
362 363

#endif