windowExtractor.h 5.92 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
#include "read_score.h"
15
#include "output.h"
16 17 18

#define NB_BINS_CLONES 10
#define MAX_VALUE_BINS_CLONES 1000
Mathieu Giraud's avatar
Mathieu Giraud committed
19 20 21 22

using namespace std;

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

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

38
  Stats stats[STATS_SIZE];
39
  size_t max_reads_per_window;
40 41

  MultiGermline *multigermline;
Mathieu Giraud's avatar
Mathieu Giraud committed
42 43
 public:

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

  /**
   * 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.
53
   * @param only_labeled_windows: remember only windows from windows_labels
54
   * @param nb_expected: maximal e-value of the segmentation
55
   * @param nb_reads_for_evalue: number of reads, used for e-value computation. Can be approximate or faked.
56
   * @param scorer: how reads are scored (only the best ones are keeped for large clones)
57
   * @param output: global output, used here for warnings
Mathieu Giraud's avatar
Mathieu Giraud committed
58 59 60 61 62
   * @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().
   */
63
  WindowsStorage *extract(OnlineBioReader *reads,
64
                          size_t w,
65
                          map<string, string> &windows_labels, bool only_labeled_windows=false,
66
                          bool keep_unsegmented_as_clone=false,
67
                          double nb_expected = THRESHOLD_NB_EXPECTED, int nb_reads_for_evalue = 1,
68 69
                          VirtualReadScore *scorer = &DEFAULT_READ_SCORE,
                          SampleOutput *output = NULL);
Mathieu Giraud's avatar
Mathieu Giraud committed
70 71 72 73 74 75 76 77

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

78 79 80 81 82
  /**
   * cf. WindowsStorage::getMaximalNbReadsPerWindow()
   */
  size_t getMaximalNbReadsPerWindow();

Mathieu Giraud's avatar
Mathieu Giraud committed
83 84 85 86 87 88 89 90 91 92 93 94
  /**
   * @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);

95 96 97 98 99 100
   /**
   * @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);
101 102 103 104 105 106

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

Mathieu Giraud's avatar
Mathieu Giraud committed
107 108 109 110 111 112 113 114 115 116 117 118 119 120
  /**
   * 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);

121 122 123
  /**
   * Defines the output streams where the unsegmented sequences will be output, split by unsegmentation cause.
   * Otherwise no output will be given.
124 125
   * @param outs: The output streams
   * @param unsegmented_detail_full: Whether we should output UNSEG_TOO_FEW_ZERO reads
126
   */
127
  void setUnsegmentedDetailOutput(ofstream **outs, bool unsegmented_detail_full=false);
128

129 130 131 132 133 134 135
  /**
   * 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);

136
  /**
137
   * Output the segmentation and germlines stats
138 139 140 141
   * @param out: The output stream
   */
  void out_stats(ostream &out);

142 143 144 145 146 147 148 149 150 151
  /**
   * 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);
152 153 154 155 156 157 158 159

  /**
   * 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
160
   *         (possibly shifted at most twice, as mentionned by shift) of the given length.
161 162 163 164 165 166 167
   *         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
168 169 170 171 172
 private:
  /**
   * Initialize the statistics (put 0 everywhere).
   */
  void init_stats();
173 174 175 176 177

  /*
   * Fill the stats_clone member of the different Germlines
   */
  void fillStatsClones(WindowsStorage *storage);
Mathieu Giraud's avatar
Mathieu Giraud committed
178 179 180
};

#endif