segment.h 11.4 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"
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
#define DETECT_THRESHOLD_STRAND 5   /* If the number of total affectations
                                       is above this threshold, then a sequence with no clearly attributed
26
                                       strand will be marked as STRAND_NOT_CONSISTENT */
27

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

30
#define BAD_EVALUE  1e10
31

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

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

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

49 50 51
#define DEFAULT_MINIMIZE_WIDTH       12

#define DEFAULT_MINIMIZE_ONE_MARGIN  20
52

53

Mikaël Salson's avatar
Mikaël Salson committed
54
using namespace std;
Marc Duez's avatar
Marc Duez committed
55
using json = nlohmann::json;
Mikaël Salson's avatar
Mikaël Salson committed
56

57 58 59
enum SEGMENTED { NOT_PROCESSED,
		 TOTAL_SEG_AND_WINDOW,
                 SEG_PLUS, SEG_MINUS,
60
                 SEG_SHORTER_WINDOW,
61
                 UNSEG_TOO_SHORT, UNSEG_STRAND_NOT_CONSISTENT,
62
		 UNSEG_TOO_FEW_ZERO,  UNSEG_ONLY_V, UNSEG_ONLY_J,
63
                 UNSEG_BAD_DELTA_MIN, UNSEG_AMBIGUOUS,
64
		 UNSEG_TOO_SHORT_FOR_WINDOW,
65

Mikaël Salson's avatar
Mikaël Salson committed
66
		 STATS_SIZE } ;
67 68 69

#define STATS_FIRST_UNSEG UNSEG_TOO_SHORT

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

80 81


82 83 84 85
/**
 * An alignment box (AlignBox) gather all parameters for a recombined gene segment (V, D, J, other D...)
 **/

86 87 88
class AlignBox
{
 public:
89
  string key;
90
  string color;
91

92 93 94 95 96
  int del_left;
  int start;
  int end;
  int del_right;

97
  AlignBox(string key = "", string color="");
98
  string getSequence(string sequence);
99
  void addToJson(json &seg);
100

101 102 103 104 105
  /**
   * Returns 'V', 'D', 'J', or possibly '5', '4', '3', '?', depending on the ref_label and on the key
   */
  char getInitial();

106 107 108 109 110
  /**
   * Returns the length
   */
  int getLength();

111 112 113 114 115 116 117 118 119 120 121 122
  /**
   * 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);

123
  /* Identifier, label and sequence of the reference sequence (the best one) */
124 125 126 127
  int ref_nb;
  string ref_label;
  string ref;

128 129 130
  /* Marked position, for Cys104 and Phe118/Trp118 */
  int marked_pos;

131 132
  /* Identifiers and scores of other possible reference sequence */
  vector<pair<int, int> > score;
133 134
};

135
ostream &operator<<(ostream &out, const AlignBox &box);
136

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

string check_and_resolve_overlap(string seq, int seq_begin, int seq_end,
152 153
                                 AlignBox *box_left, AlignBox *box_right,
                                 Cost segment_cost);
154

Mikaël Salson's avatar
Mikaël Salson committed
155 156 157
class Segmenter {
protected:
  string sequence;
158
  string sequence_or_rc;
159

160
  // JUNCTIONstart/end and CDR3start/end are 1-based
161
  int JUNCTIONstart, JUNCTIONend;
162
  string JUNCTIONaa;
163
  bool JUNCTIONproductive;
164

Marc Duez's avatar
Marc Duez committed
165
  int CDR3start, CDR3end;
166 167
  string CDR3nuc;
  string CDR3aa;
168

Mikaël Salson's avatar
Mikaël Salson committed
169
  bool reversed, segmented, dSegmented;
170
  bool shiftedJunction;
171
  int because;
Mikaël Salson's avatar
Mikaël Salson committed
172

173 174
  /**
   * Compares evalue_left, evalue_right and evalue against the provided threshold
175
   * @post some evalue is above the threshold ==> because is set to UNSEG_TOO_FEW_ZERO, UNSEG_ONLY_J or UNSEG_ONLY_V
176 177 178
   */
  void checkLeftRightEvaluesThreshold(double threshold, int strand);

Mikaël Salson's avatar
Mikaël Salson committed
179 180 181 182 183
  string removeChevauchement();
  bool finishSegmentation();
  bool finishSegmentationD();

 public:
184
  Germline *segmented_germline;
185
  string label;
Mikaël Salson's avatar
Mikaël Salson committed
186
  string code;
187 188
  string info;        // .vdj.fa header, fixed fields
  string info_extra;  // .vdj.fa header, other information, at the end of the header
189
  string seg_V, seg_N, seg_J, system;
190

191 192
  AlignBox *box_V, *box_D, *box_J;

193
  double evalue;
194 195
  double evalue_left;
  double evalue_right;
Mikaël Salson's avatar
Mikaël Salson committed
196 197 198
  string seg_N1, seg_D, seg_N2;
  Cost segment_cost;

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

Mathieu Giraud's avatar
Mathieu Giraud committed
201
  /* Queries */
Mikaël Salson's avatar
Mikaël Salson committed
202 203 204 205 206 207 208 209 210 211 212 213 214



  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.
215 216
   * @post If the size or position of the window had to be dynamically adapted to fit
   *       in the read, isJunctionShifted() will return true
Mikaël Salson's avatar
Mikaël Salson committed
217
   */
218
  string getJunction(int l, int shift=0);
Mikaël Salson's avatar
Mikaël Salson committed
219 220 221 222 223 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

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

255 256 257 258 259 260
  /**
   * @return true iff the junction has been shifted or shortened
   *              dynamically to fit into the sequence
   */
  bool isJunctionShifted() const;

261 262 263 264 265 266 267 268 269
  /**
   * @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
270

271 272 273 274 275
  /**
   * @post status == SEG_PLUS || status == SEG_MINUS <==> isSegmented()
   */
  void setSegmentationStatus(int status);

Mikaël Salson's avatar
Mikaël Salson committed
276 277 278 279 280 281 282 283 284 285 286
  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
287
 private:
288
  KmerAffectAnalyser *kaa;
Mikaël Salson's avatar
Mikaël Salson committed
289 290 291 292
 protected:
  string affects;

 public:
293
  int score;
294 295
  int pvalue_left;
  int pvalue_right;
296
  KmerAffect before, after;
297 298

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

306 307
  KmerSegmenter(const KmerSegmenter &seg);

Mathieu Giraud's avatar
Mathieu Giraud committed
308 309 310 311 312
  ~KmerSegmenter();

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

315
  string getInfoLineWithAffects() const;
Marc Duez's avatar
Marc Duez committed
316
  json toJson();
317

Mathieu Giraud's avatar
Mathieu Giraud committed
318
 private:
319
  void computeSegmentation(int strand, KmerAffect left, KmerAffect right,
320
                           double threshold, double multiplier);
Mikaël Salson's avatar
Mikaël Salson committed
321 322
};

323 324 325

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

337 338
  ~KmerMultiSegmenter();

339
  KmerSegmenter *the_kseg;
340
  MultiGermline *multi_germline;
341 342 343
};


Mikaël Salson's avatar
Mikaël Salson committed
344 345 346
class FineSegmenter : public Segmenter
{
 public:
Mikaël Salson's avatar
Mikaël Salson committed
347 348 349
   vector<pair<int, int> > score_V;
   vector<pair<int, int> > score_D;
   vector<pair<int, int> > score_J;
350 351

   vector <AlignBox*> boxes ;
Mikaël Salson's avatar
Mikaël Salson committed
352
   
Mikaël Salson's avatar
Mikaël Salson committed
353 354 355
   /**
   * Build a fineSegmenter based on KmerSegmentation
   * @param seq: An object read from a FASTA/FASTQ file
356
   * @param germline: germline used
Mikaël Salson's avatar
Mikaël Salson committed
357
   */
358
   FineSegmenter(Sequence seq, Germline *germline, Cost segment_cost,
359
                 double threshold = THRESHOLD_NB_EXPECTED, double multiplier=1.0);
360 361 362

   ~FineSegmenter();

Mikaël Salson's avatar
Mikaël Salson committed
363 364
  /**
  * extend segmentation from VJ to VDJ
365
  * @param germline: germline used
Mikaël Salson's avatar
Mikaël Salson committed
366
  */
367
  void FineSegmentD(Germline *germline, bool several_D,
368
                    double threshold = THRESHOLD_NB_EXPECTED_D, double multiplier=1.0);
369 370

  bool FineSegmentD(Germline *germline,
371
                    AlignBox *box_Y, AlignBox *box_DD, AlignBox *box_Z,
372
                    int forbidden_id,
373
                    int extend_DD_on_Y, int extend_DD_on_Z,
374
                    double threshold = THRESHOLD_NB_EXPECTED_D, double multiplier=1.0);
375

Mathieu Giraud's avatar
Mathieu Giraud committed
376 377 378 379
  /**
   * 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
380
  void findCDR3();
Mathieu Giraud's avatar
Mathieu Giraud committed
381

Marc Duez's avatar
Marc Duez committed
382
  json toJson();
Mikaël Salson's avatar
Mikaël Salson committed
383
  
Mikaël Salson's avatar
Mikaël Salson committed
384 385 386
};


387
void align_against_collection(string &read, BioReader &rep, int forbidden_rep_id,
388 389
                              bool reverse_ref, bool reverse_both, bool local,
                              AlignBox *box, Cost segment_cost);
Mikaël Salson's avatar
Mikaël Salson committed
390 391

#endif