1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 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 128 129 130 131 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
|
#ifndef SEQLIB_BFC_H
#define SEQLIB_BFC_H
extern "C" {
#include <fml.h>
#include <fml/bfc.h>
}
#include "SeqLib/BamRecord.h"
#include "SeqLib/UnalignedSequence.h"
namespace SeqLib {
/** Class to perform error-correction using BFC algorithm
*
* BFC is designed and implemented by Heng Li (https://github.com/lh3/bfc).
* From Heng: It is a variant of the classical spectrum alignment algorithm introduced
* by Pevzner et al (2001). It uses an exhaustive search to find a k-mer path
* through a read that minimizeds a heuristic objective function jointly considering
* penalities on correction, quality and k-mer support.
*/
class BFC {
public:
/** Construct a new BFC engine */
BFC() {
m_idx = 0;
bfc_opt_init(&bfc_opt);
ch = NULL;
kmer = 0;
flt_uniq = 0;
n_seqs = 0;
m_seqs = NULL;
kcov = 0;
tot_k = 0;
sum_k = 0;
tot_len = 0;
m_seqs_size = 0;
}
~BFC() {
clear();
if (ch)
bfc_ch_destroy(ch);
}
/** Peform BFC error correction on the sequences stored in this object */
bool ErrorCorrect();
/** Train the error corrector using the reads stored in this object */
bool Train();
/** Add a sequence for either training or correction
* @param seq A sequence to be copied into this object (A, T, C, G)
*/
bool AddSequence(const char* seq, const char* qual, const char* name);
/** Set the k-mer size for training
* @note zero is auto
*/
void SetKmer(int k) { kmer = k; }
/** Correct a single new sequence not stored in object
* @param str Sequence of string to correct (ACTG)
* @param q Quality score of sequence to correct
* @value Returns true if corrected */
bool CorrectSequence(std::string& str, const std::string& q);
/** Clear the stored reads */
void clear();
/** Return the calculated kcov */
float GetKCov() const { return kcov; }
/** Return the calculated kcov */
int GetKMer() const { return kmer; }
/** Return the number of sequences controlled by this */
int NumSequences() const { return n_seqs; }
/** Return the next sequence stored in object
* @param s Empty string to be filled.
* @param q Empty string name to be filled.
* @value True if string was filled with sequence. False if no more sequences.
*/
bool GetSequence(std::string& s, std::string& q);
/** Reset the sequence iterator inside GetSequence
*/
void ResetGetSequence() { m_idx = 0; };
private:
size_t m_idx;
// the amount of memory allocated
size_t m_seqs_size;
void learn_correct();
bfc_opt_t bfc_opt;
// histogram of kmer occurences
uint64_t hist[256];
// diff histogram of kmers??
uint64_t hist_high[64];
uint64_t tot_len;
uint64_t sum_k; // total valid kmer count (kmers above min_count) ?
// total number of kmers?
uint64_t tot_k;
//
float kcov;
// reads to correct in place
fml_seq1_t * m_seqs;
// number of sequeces
size_t n_seqs;
// fermi lite options
fml_opt_t fml_opt;
// vector of names
std::vector<char*> m_names;
// assign names, qualities and seq to m_seqs
void allocate_sequences_from_reads(const BamRecordVector& brv);
// assign names, qualities and seq to m_seqs
void allocate_sequences_from_char(const std::vector<char*>& v);
void allocate_sequences_from_strings(const std::vector<std::string>& v);
// do the actual read correction
void correct_reads();
// 0 turns off filter uniq
int flt_uniq; // from fml_correct call
int l_pre;
// 0 is auto learn
int kmer;
// holds data after learning how to correct
bfc_ch_t *ch;
// holds data for actual error correction
ec_step_t es;
};
}
#endif
|