Attention une mise à jour du serveur va être effectuée le vendredi 16 avril entre 12h et 12h30. Cette mise à jour va générer une interruption du service de quelques minutes.

segment.h 5.61 KB
Newer Older
Mikaël Salson's avatar
Mikaël Salson committed
1 2 3 4 5 6 7 8 9 10
#ifndef SEGMENT_H
#define SEGMENT_H

#include <string>
#include <iostream>
#include "fasta.h"
#include "dynprog.h"
#include "tools.h"
#include "kmerstore.h"
#include "kmeraffect.h"
Mathieu Giraud's avatar
Mathieu Giraud committed
11 12
#include "affectanalyser.h"
#include "json.h"
Mikaël Salson's avatar
Mikaël Salson committed
13 14

#define EXTEND_D_ZONE 5
Mikaël Salson's avatar
Mikaël Salson committed
15 16 17
#define RATIO_STRAND 5          /* The ratio between the affectations in one
                                   strand and the other, to safely attribute a
                                   segment to a given strand */
Mikaël Salson's avatar
Mikaël Salson committed
18

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

Mikaël Salson's avatar
Mikaël Salson committed
21 22 23 24
using namespace std;

enum SEGMENTED { DONT_KNOW, SEG_PLUS, SEG_MINUS, UNSEG_TOO_SHORT, UNSEG_STRAND_NOT_CONSISTENT, 
		 UNSEG_TOO_FEW_ZERO,  UNSEG_TOO_FEW_V, UNSEG_TOO_FEW_J, 
Mikaël Salson's avatar
Mikaël Salson committed
25 26 27 28 29
		 UNSEG_BAD_DELTA_MIN, UNSEG_BAD_DELTA_MAX,
		 TOTAL_SEG_AND_WINDOW, 
		 TOTAL_SEG_BUT_TOO_SHORT_FOR_THE_WINDOW,
		 STATS_SIZE } ;
const char* const segmented_mesg[] = { "?", "SEG_+", "SEG_-", "UNSEG too short", "UNSEG strand",  
Mikaël Salson's avatar
Mikaël Salson committed
30
				       "UNSEG too few (zero)", "UNSEG too few V", "UNSEG too few J",
Mikaël Salson's avatar
Mikaël Salson committed
31 32 33 34
				       "UNSEG < delta_min", "UNSEG > delta_max",
                                       "= SEG, with window",
                                       "= SEG, but no window",
                                      } ;
Mikaël Salson's avatar
Mikaël Salson committed
35 36 37 38 39

class Segmenter {
protected:
  string label;
  string sequence;
Mathieu Giraud's avatar
Mathieu Giraud committed
40 41
  int Vend, Jstart;
  int Dstart, Dend;
Mikaël Salson's avatar
Mikaël Salson committed
42 43 44 45 46 47 48 49
  bool reversed, segmented, dSegmented;

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

 public:
  string code;
Mikaël Salson's avatar
Mikaël Salson committed
50
  string code_short;
Mikaël Salson's avatar
Mikaël Salson committed
51
  string code_light;
52 53
  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
54 55 56 57 58 59 60 61
  int best_V, best_J ;
  int del_V, del_D_left, del_D_right, del_J ;
  string seg_V, seg_N, seg_J;
  
  int best_D;
  string seg_N1, seg_D, seg_N2;
  Cost segment_cost;

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

Mathieu Giraud's avatar
Mathieu Giraud committed
64
  /* Queries */
Mikaël Salson's avatar
Mikaël Salson committed
65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 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



  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;


  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
128 129
 private:
  int because;                  
130
  KmerAffectAnalyser<KmerAffect> *kaa;
Mikaël Salson's avatar
Mikaël Salson committed
131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
 protected:
  string affects;

 public:
  /**
   * Build a segmenter based on KmerSegmentation
   * @param seq: An object read from a FASTA/FASTQ file
   * @param index: A Kmer index
   * @param delta_min: the minimal distance between the right bound and the left bound 
   *        so that the segmentation is accepted 
   *        (left bound: end of V, right bound : start of J)
   * @param delta_min: the maximal distance between the right bound and the left bound 
   *        so that the segmentation is accepted 
   *        (left bound: end of V, right bound : start of J)
   */
  KmerSegmenter(Sequence seq, IKmerStore<KmerAffect> *index, 
Mathieu Giraud's avatar
Mathieu Giraud committed
147 148 149 150 151 152 153
		int delta_min, int delta_max);

  ~KmerSegmenter();

  /**
   * @return the KmerAffectAnalyser of the current sequence.
   */
154
  KmerAffectAnalyser<KmerAffect> *getKmerAffectAnalyser() const;
Mathieu Giraud's avatar
Mathieu Giraud committed
155 156 157 158 159 160 161 162 163 164

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

 private:
165
  void computeSegmentation(int strand, int delta_min, int delta_max, int s);
Mikaël Salson's avatar
Mikaël Salson committed
166 167 168 169 170
};

class FineSegmenter : public Segmenter
{
 public:
Mikaël Salson's avatar
Mikaël Salson committed
171 172 173 174
   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
175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196
   /**
   * Build a fineSegmenter based on KmerSegmentation
   * @param seq: An object read from a FASTA/FASTQ file
   * @param rep_V: germline for V
   * @param rep_J: germline for J
   * @param delta_min: the minimal distance between the right bound and the left bound 
   *        so that the segmentation is accepted 
   *        (left bound: end of V, right bound : start of J)
   * @param delta_min: the maximal distance between the right bound and the left bound 
   *        so that the segmentation is accepted 
   *        (left bound: end of V, right bound : start of J)
   */
  FineSegmenter(Sequence seq, Fasta &rep_V, Fasta &rep_J,
		int delta_min, int delta_max, Cost segment_cost);
  
  /**
  * extend segmentation from VJ to VDJ
  * @param rep_J: germline for V
  * @param rep_J: germline for D
  * @param rep_J: germline for J
  */
  void FineSegmentD(Fasta &rep_V, Fasta &rep_D, Fasta &rep_J);
Mathieu Giraud's avatar
Mathieu Giraud committed
197 198

  JsonList toJsonList(Fasta &rep_V, Fasta &rep_D, Fasta &rep_J);
Mikaël Salson's avatar
Mikaël Salson committed
199
  
Mikaël Salson's avatar
Mikaël Salson committed
200 201 202 203 204
};



#endif