windowExtractor.h 5.61 KB
Newer Older
Mathieu Giraud's avatar
Mathieu Giraud committed
1 2 3 4
#ifndef WINDOW_EXTRACTOR_H
#define WINDOW_EXTRACTOR_H

#include <iostream>
5 6
#include <string>
#include <map>
Mathieu Giraud's avatar
Mathieu Giraud committed
7
#include "segment.h"
8
#include "germline.h"
Mathieu Giraud's avatar
Mathieu Giraud committed
9 10 11
#include "kmerstore.h"
#include "kmeraffect.h"
#include "windows.h"
12
#include "read_storage.h"
13
#include "bioreader.hpp"
14 15 16

#define NB_BINS_CLONES 10
#define MAX_VALUE_BINS_CLONES 1000
Mathieu Giraud's avatar
Mathieu Giraud committed
17 18 19 20

using namespace std;

/**
21
 * This takes an OnlineBioReader reader as input and extract windows from the
Mathieu Giraud's avatar
Mathieu Giraud committed
22 23 24 25 26
 * sequences given in the input.
 */
class WindowExtractor {
 private:
  size_t nb_reads;
27 28
  map<string, BinReadStorage> stats_reads;
  map<string, BinReadStorage> stats_clones;
Mathieu Giraud's avatar
Mathieu Giraud committed
29 30 31

  ostream *out_segmented;
  ostream *out_unsegmented;
32
  ofstream **out_unsegmented_detail;
33
  bool unsegmented_detail_full;
34
  ostream *out_affects;
Mathieu Giraud's avatar
Mathieu Giraud committed
35

36
  Stats stats[STATS_SIZE];
37
  size_t max_reads_per_window;
38 39

  MultiGermline *multigermline;
Mathieu Giraud's avatar
Mathieu Giraud committed
40 41
 public:

42
  WindowExtractor(MultiGermline *multigermline);
Mathieu Giraud's avatar
Mathieu Giraud committed
43 44 45 46 47 48 49 50

  /**
   * Extract windows from the collection of input reads.
   * If (un)segmented sequences must be output, the functions 
   * set(Un)SegmentedOutput() must be called before.
   * @param reads: the collection of input reads
   * @param w: length of the window
   * @param windows_labels: Windows that must be kept and registered as such.
51
   * @param only_labeled_windows: remember only windows from windows_labels
52
   * @param nb_expected: maximal e-value of the segmentation
53
   * @param nb_reads_for_evalue: number of reads, used for e-value computation. Can be approximate or faked.
Mathieu Giraud's avatar
Mathieu Giraud committed
54 55 56 57 58
   * @return a pointer to a WindowsStorage that will contain all the windows.
   *         It is a pointer so that the WindowsStorage is not duplicated.
   * @post Statistics on segmentation will be provided through the getSegmentationStats() methods
   *       and getAverageSegmentationLength().
   */
59
  WindowsStorage *extract(OnlineBioReader *reads,
60
                          size_t w,
61
                          map<string, string> &windows_labels, bool only_labeled_windows=false,
62
                          bool keep_unsegmented_as_clone=false,
63
                          double nb_expected = THRESHOLD_NB_EXPECTED, int nb_reads_for_evalue = 1);
Mathieu Giraud's avatar
Mathieu Giraud committed
64 65 66 67 68 69 70 71

  /**
   * @return the average length of sequences whose segmentation has been classified as seg
   * @param seg: one of the segmentation available in the enum SEGMENTED 
   * @pre extract() must have been launched.
   */
  float getAverageSegmentationLength(SEGMENTED seg);

72 73 74 75 76
  /**
   * cf. WindowsStorage::getMaximalNbReadsPerWindow()
   */
  size_t getMaximalNbReadsPerWindow();

Mathieu Giraud's avatar
Mathieu Giraud committed
77 78 79 80 81 82 83 84 85 86 87 88
  /**
   * @return Total number of processed reads by the previous call to extract()
   */
  size_t getNbReads();

  /**
   * @return the number of sequences whose segmentation has been classified as seg
   * @param seg: one of the segmentation available in the enum SEGMENTED 
   * @pre extract() must have been launched.
   */
  size_t getNbSegmented(SEGMENTED seg);

89 90 91 92 93 94
   /**
   * @return the number of reads segmented from germline
   * @param germline_code: one of the germline code in multigermline
   * @pre extract() must have been launched.
   */
  size_t getNbReadsGermline(string germline_code);
95 96 97 98 99 100

  /**
   * cf. WindowsStorage::setMaximalNbReadsPerWindow()
   */
  void setMaximalNbReadsPerWindow(size_t max_reads);

Mathieu Giraud's avatar
Mathieu Giraud committed
101 102 103 104 105 106 107 108 109 110 111 112 113 114
  /**
   * Defines the output stream where the segmented sequences will be output.
   * Otherwise no output will be given.
   * @param out: The output stream
   */
  void setSegmentedOutput(ostream *out);

  /**
   * Defines the output stream where the unsegmented sequences will be output.
   * Otherwise no output will be given.
   * @param out: The output stream
   */
  void setUnsegmentedOutput(ostream *out);

115 116 117
  /**
   * Defines the output streams where the unsegmented sequences will be output, split by unsegmentation cause.
   * Otherwise no output will be given.
118 119
   * @param outs: The output streams
   * @param unsegmented_detail_full: Whether we should output UNSEG_TOO_FEW_ZERO reads
120
   */
121
  void setUnsegmentedDetailOutput(ofstream **outs, bool unsegmented_detail_full=false);
122

123 124 125 126 127 128 129
  /**
   * Defines the output stream where the detailed affects will be output.
   * Otherwise no output will be given.
   * @param out: The output stream
   */
  void setAffectsOutput(ostream *out);

130
  /**
131
   * Output the segmentation and germlines stats
132 133 134 135
   * @param out: The output stream
   */
  void out_stats(ostream &out);

136 137 138 139 140 141 142 143 144 145
  /**
   * Output segmentation stats
   */
  void out_stats_segmentation(ostream &out);

  /**
   * Output germlines statistics (read lengths per germline,
   * and number of reads per clones).
   */
  void out_stats_germlines(ostream &out);
146 147 148 149 150 151 152 153

  /**
   * Get best window length and shifts.
   * @param read_length: total length of the sequence
   * @param max_window_length: the maximal window length to start with
   * @param central_pos: start position of the second half of the window
   * @param shift: length of the shift
   * @return a pair (length, shift) giving the first window that fits at the central_pos
154
   *         (possibly shifted at most twice, as mentionned by shift) of the given length.
155 156 157 158 159 160 161
   *         shift == 0 iff length == max_window_length.
   *         In case of an odd window length, we assume that the second half of the
   *         window, starting at central_pos, is the longest.
   */
  static pair<int, int> get_best_length_shifts(size_t read_length, size_t max_window_length,
                                               int central_pos,
                                               int shift);
Mathieu Giraud's avatar
Mathieu Giraud committed
162 163 164 165 166
 private:
  /**
   * Initialize the statistics (put 0 everywhere).
   */
  void init_stats();
167 168 169 170 171

  /*
   * Fill the stats_clone member of the different Germlines
   */
  void fillStatsClones(WindowsStorage *storage);
Mathieu Giraud's avatar
Mathieu Giraud committed
172 173 174
};

#endif