segment.h 9.02 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 75 76 77 78 79 80 81 82 83 84 85 86 87 88
/**
 * Check whether there is an overlap between a Vend (*pos_seq_left) and a Jstart (*pos_seq_right),
 * 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
 * @param ref_left, ref_right:             the two reference sequences
 * @param *pos_seq_left, *pos_seq_right:   the initial end/start positions of the possibly overlapping sequences (Vend and Jstart)
 * @param *trim_ref_left, *trim_ref_right: reference to deletions
 * @param segment_cost:                    the cost used by the dynamic programing
 *
 * @post  trim_ref_left (at the end of ref_left) and trim_ref_right (at the beginning of ref_right)
 *        are set to the best number of nucleotides to trim in order to remove the overlap
 *        pos_seq_left and pos_seq_right are shifted by the good number of nucleotides
 *
 * @return                                 the N segment
 */

string check_and_resolve_overlap(string seq, int seq_begin, int seq_end,
89
                        string ref_left, string ref_right,
90 91
                        int *pos_seq_left, int *pos_seq_right,
                        int *trim_ref_left, int *trim_ref_right, Cost segment_cost);
92

Mikaël Salson's avatar
Mikaël Salson committed
93 94 95
class Segmenter {
protected:
  string sequence;
Mathieu Giraud's avatar
Mathieu Giraud committed
96 97
  int Vend, Jstart;
  int Dstart, Dend;
Marc Duez's avatar
Marc Duez committed
98
  int CDR3start, CDR3end;
Mikaël Salson's avatar
Mikaël Salson committed
99
  bool reversed, segmented, dSegmented;
100
  int because;
Mikaël Salson's avatar
Mikaël Salson committed
101

102 103 104 105 106 107
  /**
   * 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
108 109 110 111 112
  string removeChevauchement();
  bool finishSegmentation();
  bool finishSegmentationD();

 public:
113
  Germline *segmented_germline;
114
  string label;
Mikaël Salson's avatar
Mikaël Salson committed
115
  string code;
116 117
  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
118 119
  int best_V, best_J ;
  int del_V, del_D_left, del_D_right, del_J ;
120
  string seg_V, seg_N, seg_J, system;
121

Mikaël Salson's avatar
Mikaël Salson committed
122
  int best_D;
123
  double evalue;
124 125
  double evalue_left;
  double evalue_right;
Mikaël Salson's avatar
Mikaël Salson committed
126 127 128
  string seg_N1, seg_D, seg_N2;
  Cost segment_cost;

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

Mathieu Giraud's avatar
Mathieu Giraud committed
131
  /* Queries */
Mikaël Salson's avatar
Mikaël Salson committed
132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182



  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;

183 184 185 186 187 188 189 190 191
  /**
   * @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
192

193 194 195 196 197
  /**
   * @post status == SEG_PLUS || status == SEG_MINUS <==> isSegmented()
   */
  void setSegmentationStatus(int status);

Mikaël Salson's avatar
Mikaël Salson committed
198 199 200 201 202 203 204 205 206 207 208
  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
209
 private:
210
  KmerAffectAnalyser *kaa;
Mikaël Salson's avatar
Mikaël Salson committed
211 212 213 214
 protected:
  string affects;

 public:
215
  int score;
216 217
  int pvalue_left;
  int pvalue_right;
218
  KmerAffect before, after;
219 220

  KmerSegmenter();
Mikaël Salson's avatar
Mikaël Salson committed
221 222 223
  /**
   * Build a segmenter based on KmerSegmentation
   * @param seq: An object read from a FASTA/FASTQ file
224
   * @param germline: the germline
Mikaël Salson's avatar
Mikaël Salson committed
225
   */
226
  KmerSegmenter(Sequence seq, Germline *germline, double threshold = THRESHOLD_NB_EXPECTED, int multiplier=1);
Mathieu Giraud's avatar
Mathieu Giraud committed
227

228 229
  KmerSegmenter(const KmerSegmenter &seg);

Mathieu Giraud's avatar
Mathieu Giraud committed
230 231 232 233 234
  ~KmerSegmenter();

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

237
  string getInfoLineWithAffects() const;
Marc Duez's avatar
Marc Duez committed
238
  json toJson();
239

Mathieu Giraud's avatar
Mathieu Giraud committed
240
 private:
241 242
  void computeSegmentation(int strand, KmerAffect left, KmerAffect right,
                           double threshold, int multiplier);
Mikaël Salson's avatar
Mikaël Salson committed
243 244
};

245 246 247

class KmerMultiSegmenter
{
248 249
 private:
  double threshold_nb_expected;
250 251 252
 public:
  /**
   * @param seq: An object read from a FASTA/FASTQ file
253 254 255 256
   * @param multigermline: the multigerm
   * @param threshold: threshold of randomly expected segmentation
   */
  KmerMultiSegmenter(Sequence seq, MultiGermline *multigermline, ostream *out_unsegmented,
257
                     double threshold = THRESHOLD_NB_EXPECTED, int nb_reads_for_evalue = 1);
258

259 260
  ~KmerMultiSegmenter();

261
  KmerSegmenter *the_kseg;
262
  MultiGermline *multi_germline;
263 264 265
};


Mikaël Salson's avatar
Mikaël Salson committed
266 267 268
class FineSegmenter : public Segmenter
{
 public:
Mikaël Salson's avatar
Mikaël Salson committed
269 270 271 272
   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
273 274 275
   /**
   * Build a fineSegmenter based on KmerSegmentation
   * @param seq: An object read from a FASTA/FASTQ file
276
   * @param germline: germline used
Mikaël Salson's avatar
Mikaël Salson committed
277
   */
278 279
   FineSegmenter(Sequence seq, Germline *germline, Cost segment_cost,
                 double threshold = THRESHOLD_NB_EXPECTED, int multiplier=1);
Mikaël Salson's avatar
Mikaël Salson committed
280 281 282
  
  /**
  * extend segmentation from VJ to VDJ
283
  * @param germline: germline used
Mikaël Salson's avatar
Mikaël Salson committed
284
  */
285 286
  void FineSegmentD(Germline *germline,
                    double threshold = THRESHOLD_NB_EXPECTED_D, int multiplier=1);
Marc Duez's avatar
Marc Duez committed
287
  void findCDR3();
Mathieu Giraud's avatar
Mathieu Giraud committed
288

Marc Duez's avatar
Marc Duez committed
289
  json toJson();
Mikaël Salson's avatar
Mikaël Salson committed
290
  
Mikaël Salson's avatar
Mikaël Salson committed
291 292 293 294 295
};



#endif