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.56 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
#include "affectanalyser.h"
Marc Duez's avatar
Marc Duez committed
13
#include "../lib/json.hpp"
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 */

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

23 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
                                       stranf will be marked as STRAND_NOT_CONSISTEN */

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

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

33 34
#define NO_LIMIT_VALUE  -1

35
#define THRESHOLD_NB_EXPECTED 1.0 /* Threshold of the accepted expected value for number of found k-mers */
36 37


38

Mikaël Salson's avatar
Mikaël Salson committed
39
using namespace std;
Marc Duez's avatar
Marc Duez committed
40
using json = nlohmann::json;
Mikaël Salson's avatar
Mikaël Salson committed
41

42 43 44 45
enum SEGMENTED { NOT_PROCESSED,
		 TOTAL_SEG_AND_WINDOW,
                 SEG_PLUS, SEG_MINUS,
                 UNSEG_TOO_SHORT, UNSEG_STRAND_NOT_CONSISTENT,
Mikaël Salson's avatar
Mikaël Salson committed
46
		 UNSEG_TOO_FEW_ZERO,  UNSEG_TOO_FEW_V, UNSEG_TOO_FEW_J, 
47
                 UNSEG_BAD_DELTA_MIN, UNSEG_AMBIGUOUS,
48
		 UNSEG_TOO_SHORT_FOR_WINDOW,
49

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

class Segmenter {
protected:
  string sequence;
Mathieu Giraud's avatar
Mathieu Giraud committed
63 64
  int Vend, Jstart;
  int Dstart, Dend;
Marc Duez's avatar
Marc Duez committed
65
  int CDR3start, CDR3end;
Mikaël Salson's avatar
Mikaël Salson committed
66
  bool reversed, segmented, dSegmented;
67
  int because;
Mikaël Salson's avatar
Mikaël Salson committed
68 69 70 71 72 73

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

 public:
74
  Germline *segmented_germline;
75
  string label;
Mikaël Salson's avatar
Mikaël Salson committed
76
  string code;
Mikaël Salson's avatar
Mikaël Salson committed
77
  string code_short;
Mikaël Salson's avatar
Mikaël Salson committed
78
  string code_light;
79 80
  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
81 82
  int best_V, best_J ;
  int del_V, del_D_left, del_D_right, del_J ;
83
  string seg_V, seg_N, seg_J, system;
84

Mikaël Salson's avatar
Mikaël Salson committed
85
  int best_D;
86
  double evalue;
87 88
  double evalue_left;
  double evalue_right;
Mikaël Salson's avatar
Mikaël Salson committed
89 90 91
  string seg_N1, seg_D, seg_N2;
  Cost segment_cost;

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

Mathieu Giraud's avatar
Mathieu Giraud committed
94
  /* Queries */
Mikaël Salson's avatar
Mikaël Salson committed
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 132 133 134 135 136 137 138 139 140 141 142 143 144 145



  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;

146 147 148 149 150 151 152 153 154
  /**
   * @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
155

156 157 158 159 160
  /**
   * @post status == SEG_PLUS || status == SEG_MINUS <==> isSegmented()
   */
  void setSegmentationStatus(int status);

Mikaël Salson's avatar
Mikaël Salson committed
161 162 163 164 165 166 167 168 169 170 171
  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
172
 private:
173
  KmerAffectAnalyser *kaa;
Mikaël Salson's avatar
Mikaël Salson committed
174 175 176 177
 protected:
  string affects;

 public:
178
  int score;
179 180
  int pvalue_left;
  int pvalue_right;
181
  KmerAffect before, after;
182 183

  KmerSegmenter();
Mikaël Salson's avatar
Mikaël Salson committed
184 185 186
  /**
   * Build a segmenter based on KmerSegmentation
   * @param seq: An object read from a FASTA/FASTQ file
187
   * @param germline: the germline
Mikaël Salson's avatar
Mikaël Salson committed
188
   */
189
  KmerSegmenter(Sequence seq, Germline *germline, double threshold = THRESHOLD_NB_EXPECTED, int multiplier=1);
Mathieu Giraud's avatar
Mathieu Giraud committed
190

191 192
  KmerSegmenter(const KmerSegmenter &seg);

Mathieu Giraud's avatar
Mathieu Giraud committed
193 194 195 196 197
  ~KmerSegmenter();

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

200
  string getInfoLineWithAffects() const;
Marc Duez's avatar
Marc Duez committed
201
  json toJson();
202

Mathieu Giraud's avatar
Mathieu Giraud committed
203
 private:
204 205
  void computeSegmentation(int strand, KmerAffect left, KmerAffect right,
                           double threshold, int multiplier);
Mikaël Salson's avatar
Mikaël Salson committed
206 207
};

208 209 210

class KmerMultiSegmenter
{
211 212
 private:
  double threshold_nb_expected;
213 214 215
 public:
  /**
   * @param seq: An object read from a FASTA/FASTQ file
216 217 218 219
   * @param multigermline: the multigerm
   * @param threshold: threshold of randomly expected segmentation
   */
  KmerMultiSegmenter(Sequence seq, MultiGermline *multigermline, ostream *out_unsegmented,
220
                     double threshold = THRESHOLD_NB_EXPECTED, int nb_reads_for_evalue = 1);
221

222 223
  ~KmerMultiSegmenter();

224
  KmerSegmenter *the_kseg;
225
  MultiGermline *multi_germline;
226 227 228
};


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

Marc Duez's avatar
Marc Duez committed
250
  json toJson();
Mikaël Salson's avatar
Mikaël Salson committed
251
  
Mikaël Salson's avatar
Mikaël Salson committed
252 253 254 255 256
};



#endif