Commit ab7c5a7b authored by Mathieu Giraud's avatar Mathieu Giraud

Merge remote-tracking branch 'origin/master' into c++11-new

Conflicts:
	algo/tools/align.cpp
parents 06728f26 43b105e9
# 'Undefined name' for web2py globals 'db', 'auth', 'request'
pyflakes:
disable:
- F821
......@@ -13,6 +13,9 @@ endif
all:
make COVERAGE="$(COVERAGE_OPTION)" -C $(VIDJIL_ALGO_SRC)
static:
make all LDFLAGS="-static -static-libstdc++"
test:
make COVERAGE="$(COVERAGE)" unit
make functional
......@@ -22,6 +25,8 @@ test:
test_browser: unit_browser functional_browser
test_server: unit_server
test_tools:
make -C tools/tests
......@@ -83,6 +88,9 @@ functional_browser:
headless_browser:
make -C browser/test headless
unit_server:
make -C server/ unit
### Code coverage
coverage: unit_coverage should_coverage
......@@ -127,6 +135,7 @@ cleanall: clean
make -C data $^
make -C germline $^
make -C $(VIDJIL_ALGO_SRC) cleanall
make -C server cleanall
.PHONY: all test should clean cleanall distrib data germline unit_coverage should_coverage coverage
......@@ -136,7 +145,8 @@ RELEASE_SOURCE = $(wildcard $(VIDJIL_ALGO_SRC)/*.cpp) $(wildcard $(VIDJIL_ALGO_S
RELEASE_MAKE = ./Makefile $(VIDJIL_ALGO_SRC)/Makefile $(VIDJIL_ALGO_SRC)/core/Makefile $(VIDJIL_ALGO_SRC)/tests/Makefile $(VIDJIL_ALGO_SRC)/lib/Makefile germline/Makefile data/Makefile tools/tests/Makefile doc/Makefile
RELEASE_TESTS = doc/format-analysis.org data/get-sequences $(wildcard data/*.vidjil) $(wildcard data/*.analysis) $(wildcard data/*.fa) $(wildcard data/*.fq) $(VIDJIL_ALGO_SRC)/tests/should-vdj-to-tap.py $(wildcard $(VIDJIL_ALGO_SRC)/tests/should-vdj-tests/*.should-vdj.fa) $(VIDJIL_ALGO_SRC)/tests/should-to-tap.sh $(wildcard $(VIDJIL_ALGO_SRC)/tests/*.should_get) $(wildcard $(VIDJIL_ALGO_SRC)/tests/bugs/*.fa) $(wildcard $(VIDJIL_ALGO_SRC)/tests/bugs/*.should_get) $(VIDJIL_ALGO_SRC)/tests/format-json.sh $(wildcard doc/analysis-example*.vidjil) $(wildcard tools/tests/*.should_get) tools/tests/should-to-tap.sh tools/diff_json.sh
RELEASE_GERMLINES = germline/germline_id germline/get-saved-germline germline/get-germline germline/split-from-imgt.py
RELEASE_FILES = $(RELEASE_SOURCE) $(RELEASE_TESTS) $(RELEASE_MAKE) $(RELEASE_GERMLINES) doc/CHANGELOG doc/algo.org doc/LICENSE data/segmentation.fasta $(wildcard data/*.fa.gz) $(wildcard data/*.label)
RELEASE_HELP = doc/algo.org doc/locus.org doc/CHANGELOG doc/LICENSE
RELEASE_FILES = $(RELEASE_SOURCE) $(RELEASE_TESTS) $(RELEASE_MAKE) $(RELEASE_GERMLINES) $(RELEASE_HELP) data/segmentation.fasta $(wildcard data/*.fa.gz) $(wildcard data/*.label)
RELEASE_ARCHIVE = vidjil-$(RELEASE_TAG).tgz
CURRENT_DIR = vidjil
......
......@@ -28,14 +28,17 @@ EXEC=vidjil
MAINCORE=$(wildcard *.cpp)
LIBCORE=core/vidjil.a lib/lib.a
BINDIR=..
CGIDIR=../browser/cgi
BINDIR=../
CGIDIR=../browser/cgi/
VIDJIL=$(BINDIR)$(EXEC)
ALIGN_CGI=$(CGIDIR)align.cgi
CREATE_VERSION_GIT_H := $(shell test -x ./create-git-version-h.sh && ./create-git-version-h.sh)
.PHONY: all core lib clean forcedep
v: vidjil align.cgi
all: $(VIDJIL) $(ALIGN_CGI)
###
......@@ -66,22 +69,30 @@ cleanspaced:
###
align.cgi: cgi/align.o
$(ALIGN_CGI): cgi/align.o $(LIBCORE)
mkdir -p $(CGIDIR)
make -C core OPTIM="$(OPTIM)"
$(CC) -o $(CGIDIR)/align.cgi cgi/align.o $(LIBCORE) $(LDFLAGS) $(LDLIBS) $(CXXFLAGS)
$(CC) -o $@ $^ $(LDFLAGS) $(LDLIBS) $(CXXFLAGS)
###
all: $(EXEC)
debug:
make clean
make DEBUG="-ggdb"
$(EXEC): %: %.o
$(VIDJIL): $(BINDIR)%: %.o $(LIBCORE)
make -C core OPTIM="$(OPTIM)"
make -C lib OPTIM="$(OPTIM)"
$(CC) -o $(BINDIR)/$@ $^ $(LIBCORE) $(LDFLAGS) $(LDLIBS) $(CXXFLAGS)
$(CC) -o $@ $^ $(LDFLAGS) $(LDLIBS) $(CXXFLAGS)
###
# Subdirectories
###
core/%.a: FORCE
make -C core $(notdir $@)
lib/%.a: FORCE
make -C lib $(notdir $@)
clean:
make -C core $@
......@@ -103,6 +114,8 @@ forcedep:
make -C core forcedep
make -C lib forcedep
FORCE:
DEP=$(wildcard dep.mk)
ifeq (${DEP},)
......
......@@ -3,11 +3,13 @@
#include <ctype.h>
void Germline::init(string _code, char _shortcut,
int _delta_min, int _delta_max)
int _delta_min, int _delta_max,
int max_indexing)
{
code = _code ;
shortcut = _shortcut ;
index = 0 ;
this->max_indexing = max_indexing;
affect_5 = "V" ;
affect_4 = "" ;
......@@ -19,24 +21,23 @@ void Germline::init(string _code, char _shortcut,
delta_min = _delta_min ;
delta_max = _delta_max ;
stats_reads.setLabel(code);
stats_clones.setLabel("");
}
Germline::Germline(string _code, char _shortcut,
int _delta_min, int _delta_max)
int _delta_min, int _delta_max,
int max_indexing)
{
init(_code, _shortcut, _delta_min, _delta_max);
init(_code, _shortcut, _delta_min, _delta_max, max_indexing);
}
Germline::Germline(string _code, char _shortcut,
string f_rep_5, string f_rep_4, string f_rep_3,
int _delta_min, int _delta_max)
int _delta_min, int _delta_max,
int max_indexing)
{
init(_code, _shortcut, _delta_min, _delta_max);
init(_code, _shortcut, _delta_min, _delta_max, max_indexing);
f_reps_5.push_back(f_rep_5);
f_reps_4.push_back(f_rep_4);
......@@ -50,9 +51,10 @@ Germline::Germline(string _code, char _shortcut,
Germline::Germline(string _code, char _shortcut,
list <string> _f_reps_5, list <string> _f_reps_4, list <string> _f_reps_3,
int _delta_min, int _delta_max)
int _delta_min, int _delta_max,
int max_indexing)
{
init(_code, _shortcut, _delta_min, _delta_max);
init(_code, _shortcut, _delta_min, _delta_max, max_indexing);
f_reps_5 = _f_reps_5 ;
f_reps_4 = _f_reps_4 ;
......@@ -75,9 +77,10 @@ Germline::Germline(string _code, char _shortcut,
Germline::Germline(string _code, char _shortcut,
Fasta _rep_5, Fasta _rep_4, Fasta _rep_3,
int _delta_min, int _delta_max)
int _delta_min, int _delta_max,
int max_indexing)
{
init(_code, _shortcut, _delta_min, _delta_max);
init(_code, _shortcut, _delta_min, _delta_max, max_indexing);
rep_5 = _rep_5 ;
rep_4 = _rep_4 ;
......@@ -101,23 +104,22 @@ void Germline::set_index(IKmerStore<KmerAffect> *_index)
void Germline::update_index(IKmerStore<KmerAffect> *_index)
{
if (!_index) _index = index ;
_index->insert(rep_5, affect_5);
_index->insert(rep_5, affect_5, max_indexing);
if (affect_4.size())
_index->insert(rep_4, affect_4);
_index->insert(rep_3, affect_3);
_index->insert(rep_3, affect_3, -max_indexing);
}
void Germline::mark_as_ambiguous(Germline *other)
{
index->insert(other->rep_5, AFFECT_AMBIGUOUS_SYMBOL);
index->insert(other->rep_5, AFFECT_AMBIGUOUS_SYMBOL, max_indexing);
if (other->affect_4.size())
index->insert(other->rep_4, AFFECT_AMBIGUOUS_SYMBOL);
index->insert(other->rep_3, AFFECT_AMBIGUOUS_SYMBOL);
index->insert(other->rep_3, AFFECT_AMBIGUOUS_SYMBOL, -max_indexing);
}
......@@ -178,36 +180,37 @@ void MultiGermline::add_germline(Germline *germline, string seed)
germlines.push_back(germline);
}
void MultiGermline::build_default_set(string path)
void MultiGermline::build_default_set(string path, int max_indexing)
{
// Should parse 'data/germlines.data'
add_germline(new Germline("TRA", 'A', path + "/TRAV.fa", "", path + "/TRAJ.fa", -10, 20), SEED_S13);
add_germline(new Germline("TRB", 'B', path + "/TRBV.fa", path + "/TRBD.fa", path + "/TRBJ.fa", 0, 80), SEED_S12);
add_germline(new Germline("TRG", 'G', path + "/TRGV.fa", "", path + "/TRGJ.fa", -10, 30), SEED_S10);
add_germline(new Germline("TRD", 'D', path + "/TRDV.fa", path + "/TRDD.fa", path + "/TRDJ.fa", 0, 80), SEED_S10);
add_germline(new Germline("IGH", 'H', path + "/IGHV.fa", path + "/IGHD.fa", path + "/IGHJ.fa", 0, 80), SEED_S12);
add_germline(new Germline("IGK", 'K', path + "/IGKV.fa", "", path + "/IGKJ.fa", -10, 20), SEED_S10);
add_germline(new Germline("IGL", 'L', path + "/IGLV.fa", "", path + "/IGLJ.fa", -10, 20), SEED_S10);
add_germline(new Germline("TRA", 'A', path + "/TRAV.fa", "", path + "/TRAJ.fa", -10, 20, max_indexing), SEED_S13);
add_germline(new Germline("TRB", 'B', path + "/TRBV.fa", path + "/TRBD.fa", path + "/TRBJ.fa", 0, 80, max_indexing), SEED_S12);
add_germline(new Germline("TRG", 'G', path + "/TRGV.fa", "", path + "/TRGJ.fa", -10, 30, max_indexing), SEED_S10);
add_germline(new Germline("TRD", 'D', path + "/TRDV.fa", path + "/TRDD.fa", path + "/TRDJ.fa", 0, 80, max_indexing), SEED_S10);
add_germline(new Germline("IGH", 'H', path + "/IGHV.fa", path + "/IGHD.fa", path + "/IGHJ.fa", 0, 80, max_indexing), SEED_S12);
add_germline(new Germline("IGK", 'K', path + "/IGKV.fa", "", path + "/IGKJ.fa", -10, 20, max_indexing), SEED_S10);
add_germline(new Germline("IGL", 'L', path + "/IGLV.fa", "", path + "/IGLJ.fa", -10, 20, max_indexing), SEED_S10);
}
void MultiGermline::build_incomplete_set(string path)
void MultiGermline::build_incomplete_set(string path, int max_indexing)
{
// Should parse 'data/germlines.data'
// VdJa
add_germline(new Germline("VdJa", 'a', path + "/TRDV.fa", path + "/TRDD.fa", path + "/TRAJ.fa", -10, 80), SEED_S13);
// TRA+D
add_germline(new Germline("TRA+D", 'a', path + "/TRDV.fa", path + "/TRDD.fa", path + "/TRAJ.fa", -10, 80, max_indexing), SEED_S13);
add_germline(new Germline("TRA+D", 'a', path + "/TRDD_upstream.fa", "", path + "/TRAJ.fa", -10, 80, max_indexing), SEED_S13);
// DD2-DD3
add_germline(new Germline("TRD+", 'd', path + "/TRDD2-01.fa", "", path + "/TRDJ.fa", -10, 60), SEED_9);
add_germline(new Germline("TRD+", 'd', path + "/TRDV.fa", "", path + "/TRDD3-01.fa", -10, 50), SEED_9);
add_germline(new Germline("TRD+", 'd', path + "/TRDD2-01.fa", "", path + "/TRDD3-01.fa", -10, 50), SEED_9);
// DD-JD + DD2-DD3
add_germline(new Germline("TRD+", 'd', path + "/TRDD2_upstream.fa", "", path + "/TRDJ.fa", -10, 60, max_indexing), SEED_9);
add_germline(new Germline("TRD+", 'd', path + "/TRDV.fa", "", path + "/TRDD3_downstream.fa", -10, 50, max_indexing), SEED_9);
add_germline(new Germline("TRD+", 'd', path + "/TRDD2_upstream.fa", "", path + "/TRDD3_downstream.fa", -10, 50, max_indexing), SEED_9);
// DH-JH
add_germline(new Germline("IGH+", 'h', path + "/IGHD.fa", "", path + "/IGHJ.fa", -10, 20), SEED_S12);
add_germline(new Germline("IGH+", 'h', path + "/IGHD_upstream.fa", "", path + "/IGHJ.fa", -10, 20, max_indexing), SEED_S12);
// IGK: KDE, INTRON
add_germline(new Germline("IGK+", 'k', path + "/IGK-INTRON.fa", "", path + "/IGK-KDE.fa", -10, 80), SEED_S10);
add_germline(new Germline("IGK+", 'k', path + "/IGKV.fa", "", path + "/IGK-KDE.fa", -10, 80), SEED_S10);
add_germline(new Germline("IGK+", 'k', path + "/IGK-INTRON.fa", "", path + "/IGK-KDE.fa", -10, 80, max_indexing), SEED_S10);
add_germline(new Germline("IGK+", 'k', path + "/IGKV.fa", "", path + "/IGK-KDE.fa", -10, 80, max_indexing), SEED_S10);
}
/* if 'one_index_per_germline' was not set, this should be called once all germlines have been loaded */
......@@ -234,25 +237,10 @@ void MultiGermline::build_with_one_index(string seed, bool set_index)
insert_in_one_index(index, set_index);
}
void MultiGermline::out_stats(ostream &out)
{
out << " " ;
out << "reads av. len clones av. rds" ;
out << endl ;
for (list<Germline*>::const_iterator it = germlines.begin(); it != germlines.end(); ++it)
{
Germline *germline = *it ;
out << germline->stats_reads ;
out << germline->stats_clones ;
out << endl ;
}
}
/* Mark k-mers common to several germlines as ambiguous */
void MultiGermline::mark_cross_germlines_as_ambiguous()
{
string VdJa = "VdJa";
string VdJa = "TRA+D";
for (list<Germline*>::const_iterator it = germlines.begin(); it != germlines.end(); ++it)
{
......
......@@ -15,8 +15,11 @@ using namespace std;
class Germline {
private:
int max_indexing;
void init(string _code, char _shortcut,
int _delta_min, int _delta_max);
int _delta_min, int _delta_max,
int max_indexing);
public:
/*
......@@ -26,22 +29,27 @@ class Germline {
* @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)
* @param max_indexing: maximal length of the sequence to be indexed (0: all)
*/
Germline(string _code, char _shortcut,
list <string> f_rep_5, list <string> f_rep_4, list <string> f_rep_3,
int _delta_min, int _delta_max);
int _delta_min, int _delta_max,
int max_indexing=0);
Germline(string _code, char _shortcut,
string f_rep_5, string f_rep_4, string f_rep_3,
int _delta_min, int _delta_max);
int _delta_min, int _delta_max,
int max_indexing=0);
Germline(string _code, char _shortcut,
Fasta _rep_5, Fasta _rep_4, Fasta _rep_3,
int _delta_min, int _delta_max);
int _delta_min, int _delta_max,
int max_indexing=0);
Germline(string _code, char _shortcut,
int _delta_min, int _delta_max);
int _delta_min, int _delta_max,
int max_indexing=0);
~Germline();
......@@ -72,9 +80,6 @@ class Germline {
int delta_min;
int delta_max;
Stats stats_reads;
Stats stats_clones;
};
......@@ -96,8 +101,8 @@ class MultiGermline {
void insert(Germline *germline);
void add_germline(Germline *germline, string seed);
void build_default_set(string path);
void build_incomplete_set(string path);
void build_default_set(string path, int max_indexing);
void build_incomplete_set(string path, int max_indexing);
// Creates and update an unique index for all the germlines
// If 'set_index' is set, set this index as the index for all germlines
......@@ -105,8 +110,6 @@ class MultiGermline {
void build_with_one_index(string seed, bool set_index);
void mark_cross_germlines_as_ambiguous();
void out_stats(ostream &out);
};
ostream &operator<<(ostream &out, const MultiGermline &multigermline);
......
......@@ -56,23 +56,27 @@ public:
* @param label: label that must be associated to the given files
* @post All the sequences in the FASTA files have been indexed, and the label is stored in the list of labels
*/
void insert(Fasta& input, const string& label="");
void insert(Fasta& input, const string& label="", int keep_only = 0);
/**
* @param input: A list of FASTA files
* @param label: label that must be associated to the given files
* @post All the sequences in the FASTA files have been indexed, and the label is stored in the list of labels
*/
void insert(list<Fasta>& input, const string& label="");
void insert(list<Fasta>& input, const string& label="", int keep_only = 0);
/**
* @param input: A sequence to be cut in k-mers
* @param label: label that must be associated to the given files
* @param keep_only: if > 0 will keep at most the last keep_only nucleotides
* of the sequence. if < 0 will keep at most the first
* keep_only nucleotides of the sequence. if == 0,
* will keep all the sequence.
* @post All the k-mers in the sequence have been indexed.
*/
void insert(const seqtype &sequence,
const string &label,
bool ignore_extended_nucleotides=true);
bool ignore_extended_nucleotides=true, int keep_only = 0);
/**
* @param word: a k-mer
......@@ -105,6 +109,12 @@ public:
*/
string getSeed() const;
/**
* @param kmer: a kmer
* @return one label associated with the kmer
*/
string getLabel(T kmer) const;
/**
* @param seq: a sequence
* @param no_revcomp: force not to revcomp the sequence, even if
......@@ -187,28 +197,43 @@ IKmerStore<T>::~IKmerStore(){}
template<class T>
void IKmerStore<T>::insert(list<Fasta>& input,
const string &label){
const string &label,
int keep_only){
for(list<Fasta>::iterator it = input.begin() ; it != input.end() ; it++){
insert(*it, label);
insert(*it, label, keep_only);
}
}
template<class T>
void IKmerStore<T>::insert(Fasta& input,
const string &label){
const string &label,
int keep_only){
for (int r = 0; r < input.size(); r++) {
insert(input.sequence(r), label);
insert(input.sequence(r), label, true, keep_only);
}
labels.push_back(make_pair(T(label, 1), label)) ;
if (revcomp_indexed && ! T::hasRevcompSymetry()) {
labels.push_back(make_pair(T(label, -1), label)) ;
}
}
template<class T>
void IKmerStore<T>::insert(const seqtype &sequence,
const string &label,
bool ignore_extended_nucleotides){
for(size_t i = 0 ; i + s < sequence.length() + 1 ; i++) {
seqtype kmer = spaced(sequence.substr(i, s), seed);
bool ignore_extended_nucleotides,
int keep_only){
size_t start_indexing = 0;
size_t end_indexing = sequence.length();
if (keep_only > 0 && sequence.length() > (size_t)keep_only) {
start_indexing = sequence.length() - keep_only;
} else if (keep_only < 0 && sequence.length() > (size_t) -keep_only) {
end_indexing = -keep_only;
}
for(size_t i = start_indexing ; i + s < end_indexing + 1 ; i++) {
seqtype substr = sequence.substr(i, s);
seqtype kmer = spaced(substr, seed);
if (ignore_extended_nucleotides && has_extended_nucleotides(kmer))
continue;
......@@ -227,7 +252,7 @@ void IKmerStore<T>::insert(const seqtype &sequence,
}
this_kmer += T(label, strand);
if (revcomp_indexed && ! T::hasRevcompSymetry()) {
seqtype rc_kmer = revcomp(kmer);
seqtype rc_kmer = spaced(revcomp(substr), seed);
T &this_rc_kmer = this->get(rc_kmer);
if (this_rc_kmer.isNull())
nb_kmers_inserted++;
......@@ -278,6 +303,13 @@ string IKmerStore<T>::getSeed() const {
return seed;
}
template<class T>
string IKmerStore<T>::getLabel(T kmer) const {
for (typename list< pair<T, string> >::const_iterator it = labels.begin(); it != labels.end(); ++it)
if (it->first == kmer)
return it->second ;
return "" ;
}
// .getResults()
template<class T>
......
......@@ -5,13 +5,12 @@
#include "tools.h"
map <string, string> load_map(string map_file)
void load_into_map(map <string, string> &the_map, string map_file)
{
// Loads a simple file with key, values into a map
map <string, string> the_map ;
if (!map_file.size())
return the_map ;
return ;
cout << " <== " << map_file ;
......@@ -20,7 +19,6 @@ map <string, string> load_map(string map_file)
if (!f.is_open())
{
cout << " [failed] " << endl ;
return the_map ;
}
int nb_keys = 0 ;
......@@ -42,8 +40,6 @@ map <string, string> load_map(string map_file)
}
cout << ": " << nb_keys << " elements" << endl ;
return the_map ;
}
using namespace std;
#include <map>
#include <string>
#include <fstream>
......@@ -6,5 +6,5 @@ using namespace std;
#include <list>
#include "fasta.h"
map <string, string> load_map(string map_file);
void load_into_map(map <string, string> &the_map, string map_file);
#include "read_storage.h"
#include "tools.h"
size_t VirtualReadStorage::getMaxNbReadsStored() const{
return maxNbStored;
......@@ -12,24 +13,54 @@ void VirtualReadStorage::setMaxNbReadsStored(size_t nb_reads) {
//////////////////////////////////////////////////
BinReadStorage::BinReadStorage()
:nb_bins(0), max_score(0), nb_inserted(0), nb_stored(0), smallest_bin_not_empty(~0) {
bins = NULL;
}
:nb_bins(0), bins(NULL), score_bins(NULL), nb_scores(NULL), total_nb_scores(0), max_score(0),
nb_inserted(0), nb_stored(0), smallest_bin_not_empty(~0),label() {}
void BinReadStorage::init(size_t nb_bins, size_t max_score, const VirtualReadScore *vrs) {
void BinReadStorage::init(size_t nb_bins, size_t max_score, const VirtualReadScore *vrs, bool no_list) {
this->nb_bins = nb_bins;
this->max_score = max_score;
bins = new list<Sequence>[nb_bins+1];
if (no_list)
bins = NULL;
else
bins = new list<Sequence>[nb_bins+1];
score_bins = new double[nb_bins+1];
nb_scores = new size_t[nb_bins+1];
for (size_t i = 0; i <= nb_bins; i++) {
score_bins[i] = 0;
nb_scores[i] = 0;
}
scorer = vrs;
}
BinReadStorage::~BinReadStorage() {
if (bins)
delete [] bins;
if (score_bins) {
delete [] score_bins;
delete [] nb_scores;
}
}
void BinReadStorage::addScore(Sequence &s) {
addScore(scorer->getScore(s.sequence));
}
void BinReadStorage::addScore(float score) {
addScore(scoreToBin(score), score);
}
void BinReadStorage::addScore(size_t bin, float score) {
score_bins[bin] += score;
nb_scores[bin]++;
total_nb_scores++;
}
void BinReadStorage::add(Sequence &s) {
size_t bin = scoreToBin(scorer->getScore(s.sequence));
float score = scorer->getScore(s.sequence);
size_t bin = scoreToBin(score);
addScore(bin, score);
if (nb_stored < getMaxNbReadsStored()) {
bins[bin].push_back(s);
nb_stored++;
......@@ -48,10 +79,53 @@ void BinReadStorage::add(Sequence &s) {
nb_inserted++;
}
size_t BinReadStorage::getNbBins() const {
return nb_bins;
}
size_t BinReadStorage::getNbInserted() const {
return nb_inserted;
}
double BinReadStorage::getAverageScoreBySeq(Sequence &s) {
return getAverageScoreByScore(scorer->getScore(s.sequence));
}
double BinReadStorage::getAverageScoreByScore(float score) {
return getAverageScore(scoreToBin(score));
}
double BinReadStorage::getAverageScore(size_t bin) {
return getScore(bin) / getNbScores(bin);
}
double BinReadStorage::getInvertedAverageScore(size_t bin) {
return getNbScores(bin) / getScore(bin);
}
double BinReadStorage::getScoreBySeq(Sequence &s) {
return getScoreByScore(scorer->getScore(s.sequence));
}
double BinReadStorage::getScoreByScore(float score) {
return getScore(scoreToBin(score));
}
double BinReadStorage::getScore(size_t bin) {
if (bin > getNbBins()) {
double sum = 0;
for (size_t i = 0; i <= getNbBins(); i++)
sum += score_bins[i];
return sum;