[med-svn] [Git][med-team/gatb-core][upstream] New upstream version 1.4.1+git20190813.a73b6dd+dfsg
Andreas Tille
gitlab at salsa.debian.org
Tue Aug 20 12:10:18 BST 2019
Andreas Tille pushed to branch upstream at Debian Med / gatb-core
Commits:
8febe06b by Andreas Tille at 2019-08-17T20:09:04Z
New upstream version 1.4.1+git20190813.a73b6dd+dfsg
- - - - -
28 changed files:
- gatb-core/.gitignore
- gatb-core/CMakeLists.txt
- gatb-core/examples/Makefile
- gatb-core/src/gatb/bank/impl/BankFasta.cpp
- gatb-core/src/gatb/bcalm2/bglue_algo.cpp
- gatb-core/src/gatb/bcalm2/ograph.cpp
- gatb-core/src/gatb/debruijn/impl/BranchingAlgorithm.cpp
- gatb-core/src/gatb/debruijn/impl/GraphUnitigs.cpp
- gatb-core/src/gatb/debruijn/impl/IterativeExtensions.cpp
- gatb-core/src/gatb/debruijn/impl/IterativeExtensions.hpp
- gatb-core/src/gatb/debruijn/impl/LinkTigs.cpp
- gatb-core/src/gatb/debruijn/impl/LinkTigs.hpp
- gatb-core/src/gatb/debruijn/impl/Simplifications.cpp
- gatb-core/src/gatb/debruijn/impl/Simplifications.hpp
- gatb-core/src/gatb/debruijn/impl/Terminator.cpp
- gatb-core/src/gatb/kmer/impl/ConfigurationAlgorithm.cpp
- gatb-core/src/gatb/kmer/impl/DebloomAlgorithm.cpp
- gatb-core/src/gatb/kmer/impl/MPHFAlgorithm.cpp
- gatb-core/src/gatb/kmer/impl/Model.hpp
- gatb-core/src/gatb/kmer/impl/PartiInfo.cpp
- gatb-core/src/gatb/kmer/impl/PartiInfo.hpp
- gatb-core/src/gatb/kmer/impl/RepartitionAlgorithm.cpp
- gatb-core/src/gatb/kmer/impl/SortingCountAlgorithm.cpp
- gatb-core/src/gatb/template/TemplateSpecialization10.cpp.in
- gatb-core/src/gatb/tools/collections/impl/MapMPHF.hpp
- gatb-core/src/gatb/tools/misc/impl/Histogram.hpp
- gatb-core/src/gatb/tools/misc/impl/Pool.hpp
- gatb-core/test/jenkins/test-suite-fedora20-gcc-4.8.sh
Changes:
=====================================
gatb-core/.gitignore
=====================================
@@ -4,4 +4,8 @@ CMakeLists.txt.user
/.cproject
*.leon
.DS_Store
+thirdparty/hdf5/src/H5Edefin.h
+thirdparty/hdf5/src/H5Einit.h
+thirdparty/hdf5/src/H5Epubgen.h
+thirdparty/hdf5/src/H5Eterm.h
=====================================
gatb-core/CMakeLists.txt
=====================================
@@ -1,7 +1,7 @@
project(gatb-core)
# We set the required version
-cmake_minimum_required (VERSION 3.1.0)
+cmake_minimum_required (VERSION 3.10.0)
################################################################################
# The version number.
@@ -211,7 +211,7 @@ set (gatb-core-flags ${LIBRARY_COMPILE_DEFINITIONS})
set (gatb-core-includes ${PROJECT_BINARY_DIR}/include ${PROJECT_BINARY_DIR}/include/${CMAKE_BUILD_TYPE} ${PROJECT_SOURCE_DIR}/src ${PROJECT_SOURCE_DIR}/thirdparty ${gatb-core-extra-libraries-inc})
# We define the libraries used for linking binary based on gatb core
-set (gatb-core-libraries gatbcore-static dl pthread z hdf5 ${gatb-core-extra-libraries})
+set (gatb-core-libraries gatbcore-static dl pthread z hdf5-static ${gatb-core-extra-libraries})
# We define the directory where to find cmake helpers
set (gatb-core-cmake ${CMAKE_MODULE_PATH})
@@ -257,7 +257,7 @@ ADD_SUBDIRECTORY(thirdparty)
# DEPENDENCIES
################################################################################
# we must be sure that hdf5 is built and installed before building gatb-core
-ADD_DEPENDENCIES (gatbcore-static hdf5 hdf5_postbuild)
+ADD_DEPENDENCIES (gatbcore-static hdf5-static hdf5_postbuild)
################################################################################
# DOCUMENTATION GENERATION
=====================================
gatb-core/examples/Makefile
=====================================
@@ -8,7 +8,7 @@ GATB=$(shell pwd)/../
GATB_LIB=$(GATB)/build/lib
CXXFLAGS = -std=c++0x -O3 -I$(GATB)/src -I$(GATB)/build/include -I$(GATB)/thirdparty
-LDFLAGS= -L$(GATB_LIB) -lgatbcore -lpthread -lz -lhdf5 -std=c++0x -ldl -static
+LDFLAGS= -L$(GATB_LIB) -lgatbcore -lpthread -lhdf5 -lz -std=c++0x -ldl -static
SRCS = $(wildcard *.cpp)
=====================================
gatb-core/src/gatb/bank/impl/BankFasta.cpp
=====================================
@@ -58,7 +58,7 @@ typedef struct
{
gzFile stream;
unsigned char *buffer;
- int buffer_start, buffer_end;
+ uint64_t buffer_start, buffer_end;
bool eof;
char last_char;
@@ -79,7 +79,7 @@ struct variable_string_t
variable_string_t () : length(0), max(0), string(0) {}
~variable_string_t () { if (string!=0) { FREE(string); } }
- int length, max;
+ uint64_t length, max;
char *string;
};
@@ -433,7 +433,7 @@ inline signed int buffered_gets (
if (bf->buffer_start >= bf->buffer_end && bf->eof) return -1;
while (1)
{
- int i;
+ uint64_t i;
if (bf->buffer_start >= bf->buffer_end) if (!rebuffer (bf)) break;
if (allow_spaces)
{
@@ -452,6 +452,7 @@ inline signed int buffered_gets (
if (is_power_of_2(s->max)) { s->max ++; }
nearest_power_of_2(s->max);
s->string = (char*) REALLOC (s->string, s->max);
+ //std::cout << "realloc of size " << s->max << " res: " << (uint64_t)(s->string) << std::endl;
}
memcpy (s->string + s->length, bf->buffer + bf->buffer_start, i - bf->buffer_start);
s->length += i - bf->buffer_start;
@@ -765,7 +766,7 @@ void BankFasta::Iterator::estimate (u_int64_t& number, u_int64_t& totalSize, u_i
{
// linear extrapolation
number = (number * _ref.getSize()) / actualPosition;
- totalSize = (totalSize * _ref.getSize()) / actualPosition;
+ totalSize = totalSize * ( (float)_ref.getSize() / (float)actualPosition) ;
}
}
=====================================
gatb-core/src/gatb/bcalm2/bglue_algo.cpp
=====================================
@@ -60,10 +60,18 @@ using namespace gatb::core::tools::dp::impl;
using namespace gatb::core::tools::collections;
using namespace gatb::core::tools::collections::impl;
-
-
using namespace std;
+// let's be clear here:
+// UF hashes will be stored in 32 bits for efficiency (as I don't want to have a 64-bits UF for memory reasons, also, would require to modify unionFind.hpp)
+typedef uint32_t uf_hashes_t;
+// but there can be more than 2^{32}Â sequences in the glue file
+typedef uint64_t seq_idx_t;
+// so, potentially, more than 2^{32} UF hashes (but not necessarily, consider that some sequences don't need to be glued)
+// what will happen is that more one UF class won't be linked to a single unitig, but multiple unitigs
+// let's hope that there won't be saturation (only 1 UF class with all unitigs)
+// if this happens, then "Top 10 glue partitions by size:" will show only one entry and BCALM will blow up in memory
+// a fix would be to use a 64 bits UF (to be coded later)
namespace gatb { namespace core { namespace debruijn { namespace impl {
@@ -189,18 +197,17 @@ static string skip_first_abundance(const string& list)
return res;
}
-
template<int SPAN>
struct markedSeq
{
// there used to be "string seq; string abundance" but i noticed that i did not need that info for determining the chain of glues. not much space saved though (like 10-20%). I suppose the biggest memory-hog is the ks/ke unordered_map
- uint64_t index;
+ seq_idx_t index;
bool rc;
bool lmark, rmark;
typedef typename Kmer<SPAN>::Type Type;
Type ks, ke; // [start,end] kmers of seq, in canonical form (redundant information with seq, but helpful)
- markedSeq(uint64_t index, bool lmark, bool rmark, const Type &ks, const Type &ke) : index(index), rc(false), lmark(lmark), rmark(rmark), ks(ks), ke(ke) {};
+ markedSeq(seq_idx_t index, bool lmark, bool rmark, const Type &ks, const Type &ke) : index(index), rc(false), lmark(lmark), rmark(rmark), ks(ks), ke(ke) {};
void revcomp()
{
@@ -213,22 +220,22 @@ struct markedSeq
// hack to refer to a sequences in msInPart as reverse complemented
-static uint32_t is_rev_index(uint32_t index)
+static seq_idx_t is_rev_index(seq_idx_t index)
{
- return ((index >> 31 & 1) == 1);
+ return (((index >> (uint64_t)63) & (uint64_t)1) == 1);
}
-static uint32_t rev_index(uint32_t index)
+static seq_idx_t rev_index(seq_idx_t index)
{
if (is_rev_index(index))
{ std::cout << "Error: glue sequence index too large " << index << std::endl; exit(1);}
- return index | (1<<31);
+ return index | ((uint64_t)1<<(uint64_t)63);
}
-static uint32_t no_rev_index(uint32_t index)
+static seq_idx_t no_rev_index(seq_idx_t index)
{
- return index & ((1LL<<31) - 1LL);
+ return index & (((uint64_t)1<<(uint64_t)63) - (uint64_t)1);
}
//typedef lazy::memory::buffer_allocator<markedSeq> custom_allocator_t;
@@ -238,25 +245,27 @@ static uint32_t no_rev_index(uint32_t index)
* output: res, a list of lists of sequences that will be glued together
*/
template<int SPAN>
-static void determine_order_sequences(vector<vector<uint32_t>> &res, const vector<markedSeq<SPAN>> &markedSequences, int kmerSize, bool debug=false)
+static void determine_order_sequences(vector<vector<seq_idx_t>> &res, vector<bool> &res_is_circular, const vector<markedSeq<SPAN>> &markedSequences, int kmerSize, bool debug=false)
{
typedef typename Kmer<SPAN>::Type Type;
- unordered_map<Type, set<uint32_t> > kmerIndex;
- set<uint32_t> usedSeq;
- unsigned int nb_chained = 0;
+ unordered_map<Type, set<seq_idx_t> > kmerIndex;
+ set<seq_idx_t> usedSeq;
+ uint64_t nb_chained = 0;
// index kmers to their seq
// kmerIndex associates a kmer extremity to its index in markedSequences
- for (uint32_t i = 0; i < markedSequences.size(); i++)
+ for (seq_idx_t i = 0; i < markedSequences.size(); i++)
{
kmerIndex[markedSequences[i].ks].insert(i);
kmerIndex[markedSequences[i].ke].insert(i);
}
- auto glue_from_extremity = [&](markedSeq<SPAN> current, uint32_t chain_index, uint32_t markedSequence_index, bool expect_circular=false)
+ auto glue_from_extremity = [&](markedSeq<SPAN> current, seq_idx_t chain_index, seq_idx_t markedSequence_index, bool expect_circular=false)
{
- uint32_t first_index = markedSequence_index;
- vector<uint32_t> chain;
+#ifndef NDEBUG
+ seq_idx_t first_index = markedSequence_index;
+#endif
+ vector<seq_idx_t> chain;
chain.push_back(chain_index);
bool rmark = current.rmark;
@@ -268,14 +277,14 @@ static void determine_order_sequences(vector<vector<uint32_t>> &res, const vecto
std::cout << "current ke " << current.ke << " index " << no_rev_index(chain_index) << " markings: " << current.lmark << current.rmark <<std::endl;
// this sequence has a rmark, so necessarily there is another sequence to glue it with. find it here.
- set<uint32_t> candidateSuccessors = kmerIndex[current.ke];
+ set<seq_idx_t> candidateSuccessors = kmerIndex[current.ke];
assert(candidateSuccessors.find(markedSequence_index) != candidateSuccessors.end()); // remove the current seq from our indexing data structure
candidateSuccessors.erase(markedSequence_index);
assert(candidateSuccessors.size() == 1); // normally there is exactly one sequence to glue with
- uint32_t successor_index = *candidateSuccessors.begin(); // pop()
+ seq_idx_t successor_index = *candidateSuccessors.begin(); // pop()
assert(successor_index != markedSequence_index);
markedSeq<SPAN> successor = markedSequences[successor_index];
@@ -315,6 +324,11 @@ static void determine_order_sequences(vector<vector<uint32_t>> &res, const vecto
if (usedSeq.find(markedSequence_index) != usedSeq.end())
{
assert(markedSequence_index == first_index);
+ if (debug)
+ std::cout << "breaking at circular unitig" << std::endl;
+ // here it would be tempting to remove the last element of the chain to remove the last kmer
+ // but this strategy doesn't work as that element might be longer than a kmer
+ // so the preferred strategy is to cut the last nucleotide of a circular unitig
break;
}
}
@@ -329,11 +343,12 @@ static void determine_order_sequences(vector<vector<uint32_t>> &res, const vecto
}
res.push_back(chain);
+ res_is_circular.push_back(expect_circular);
nb_chained += chain.size();
};
- // iterated markedSequences, and picks extremities of a chain
+ // iterate markedSequences, and picks extremities of a chain
for (unsigned int i = 0; i < markedSequences.size(); i++)
{
markedSeq<SPAN> current = markedSequences[i];
@@ -352,7 +367,7 @@ static void determine_order_sequences(vector<vector<uint32_t>> &res, const vecto
}
/* normalize sequence so that lmark is false */
- uint32_t chain_index = markedSequences[i].index;
+ seq_idx_t chain_index = markedSequences[i].index;
if (current.lmark)
{
current.revcomp();
@@ -365,13 +380,13 @@ static void determine_order_sequences(vector<vector<uint32_t>> &res, const vecto
}
- // special case: a circular unitig, to be glued with multiple sequences all containing doubled kmers at extremities
- // my fix plan: we pick an extremity at random, and chop the last nucleotide and mark it to not be glued. also find the corresponding kmer in other extremity, and mark it as not to be glued
+ // handle the special cases undetected in the previous loop:
+ // they are circular unitigs, to be glued with other sequences all containing doubled kmers at extremities
while (nb_chained < markedSequences.size())
{
//std::cout << "nb chained " << nb_chained << " markedseq size" << markedSequences.size() << std::endl;
- vector<int> remaining_indices;
- for (uint32_t i = 0; i < markedSequences.size(); i++)
+ vector<seq_idx_t> remaining_indices;
+ for (seq_idx_t i = 0; i < markedSequences.size(); i++)
{
if (usedSeq.find(i) == usedSeq.end())
remaining_indices.push_back(i);
@@ -380,23 +395,17 @@ static void determine_order_sequences(vector<vector<uint32_t>> &res, const vecto
assert(remaining_indices.size()>0);
markedSeq<SPAN> current = markedSequences[remaining_indices[0]];
- uint32_t chain_index = markedSequences[remaining_indices[0]].index;
+ seq_idx_t chain_index = markedSequences[remaining_indices[0]].index;
glue_from_extremity(current, chain_index, remaining_indices[0], true);
}
-
- if (nb_chained < markedSequences.size())
- {
- std::cout << " Note: " << markedSequences.size() - nb_chained << " sequence chunks not returned in output unitigs (likely small circular contigs)" << std::endl;
- }
- // assert(sequences.size() == nb_chained); // make sure we've scheduled to glue all sequences in this partition
}
/* straightforward glueing of a chain
* sequences should be ordered and in the right orientation
- * so, it' just a matter of chopping of the first kmer
+ * so, it's just a matter of chopping of the first kmer of elements i>1 of each chain
*/
-static void glue_sequences(vector<uint32_t> &chain, std::vector<std::string> &sequences, std::vector<std::string> &abundances, int kmerSize, string &res_seq, string &res_abundances)
+static void glue_sequences(vector<seq_idx_t> &chain, bool is_circular, std::vector<std::string> &sequences, std::vector<std::string> &abundances, int kmerSize, string &res_seq, string &res_abundances)
{
bool debug=false;
@@ -406,7 +415,7 @@ static void glue_sequences(vector<uint32_t> &chain, std::vector<std::string> &se
if (debug) std::cout << "glueing new chain: ";
for (auto it = chain.begin(); it != chain.end(); it++)
{
- uint32_t idx = *it;
+ seq_idx_t idx = *it;
string seq = sequences[no_rev_index(idx)];
string abs = abundances[no_rev_index(idx)];
@@ -434,7 +443,24 @@ static void glue_sequences(vector<uint32_t> &chain, std::vector<std::string> &se
previous_kmer = seq.substr(seq.size() - k);
assert(previous_kmer.size() == k);
}
- if (debug) std::cout << std::endl;
+ if (is_circular)
+ {
+ if (debug) std::cout << "chopping off last nucleotide" << std::endl;
+ if (debug) std::cout << res_seq << std::endl;
+ if (debug) std::cout << res_abundances << std::endl;
+ // trick: do it with the first kmer instead
+ res_seq = rc(res_seq);
+ res_abundances = reverse_abundances(res_abundances);
+
+ res_seq = res_seq.substr(1);
+ res_abundances = skip_first_abundance(res_abundances);
+
+ res_seq = rc(res_seq);
+ res_abundances = reverse_abundances(res_abundances);
+ if (debug) std::cout << res_seq << std::endl;
+ if (debug) std::cout << res_abundances << std::endl;
+ }
+ if (debug) std::cout << std::endl;
}
@@ -479,8 +505,6 @@ class hasher_t
}
};
-typedef uint64_t partition_t;
-
/* computes and uniquifies the hashes of marked kmers at extremities of all to-be-glued sequences */
template <int SPAN>
void prepare_uf(std::string prefix, IBank *in, const int nb_threads, int& kmerSize, int pass, int nb_passes, uint64_t &nb_elts, uint64_t estimated_nb_glue_sequences)
@@ -489,7 +513,7 @@ void prepare_uf(std::string prefix, IBank *in, const int nb_threads, int& kmerSi
std::atomic<unsigned long> nb_marked_extremities, nb_unmarked_extremities;
nb_marked_extremities = 0; nb_unmarked_extremities = 0;
- std::vector<std::vector<partition_t >> uf_hashes_vectors(nb_threads);
+ std::vector<std::vector<uf_hashes_t>> uf_hashes_vectors(nb_threads);
// relatively accurate number of sequences to be inserted
for (int i = 0; i < nb_threads; i++)
@@ -497,7 +521,7 @@ void prepare_uf(std::string prefix, IBank *in, const int nb_threads, int& kmerSi
/* class (formerly a simple lambda function) to process a kmer and decide which bucket(s) it should go to */
/* needed to make it a class because i want it to remember its thread index */
- class UniquifyKeys
+ class RepartHashes
{
typedef typename Kmer<SPAN>::ModelCanonical ModelCanon;
@@ -506,13 +530,13 @@ void prepare_uf(std::string prefix, IBank *in, const int nb_threads, int& kmerSi
ModelCanon modelCanon;
Hasher_T<ModelCanon> hasher;
std::atomic<unsigned long> &nb_marked_extremities, &nb_unmarked_extremities;
- std::vector<std::vector<partition_t >> &uf_hashes_vectors;
+ std::vector<std::vector<uf_hashes_t>> &uf_hashes_vectors;
int _currentThreadIndex;
public:
- UniquifyKeys(int k, int pass, int nb_passes, int nb_threads,
+ RepartHashes(int k, int pass, int nb_passes, int nb_threads,
std::atomic<unsigned long> &nb_marked_extremities, std::atomic<unsigned long> & nb_unmarked_extremities,
- std::vector<std::vector<partition_t >> &uf_hashes_vectors
+ std::vector<std::vector<uf_hashes_t>> &uf_hashes_vectors
) : k(k), pass(pass), nb_passes(nb_passes), nb_threads(nb_threads), modelCanon(k), hasher(modelCanon),
nb_marked_extremities(nb_marked_extremities), nb_unmarked_extremities(nb_unmarked_extremities),
uf_hashes_vectors(uf_hashes_vectors), _currentThreadIndex(-1)
@@ -529,12 +553,13 @@ void prepare_uf(std::string prefix, IBank *in, const int nb_threads, int& kmerSi
if (lmark)
{
const string kmerBegin = seq.substr(0, k );
- // UF of canonical kmers in ModelCanon form, then hashed
+ // canonical kmers in ModelCanon form, then hashed
const typename ModelCanon::Kmer kmmerBegin = modelCanon.codeSeed(kmerBegin.c_str(), Data::ASCII);
const uint64_t h1 = hasher(kmmerBegin);
if (h1 % (uint64_t)nb_passes == (uint64_t)pass)
{
- uf_hashes_vectors[thread].push_back(h1);
+ uf_hashes_vectors[thread].push_back(h1); // this is where the 64 bits "hasher() result" to the 32 bits "uf_hash_t" conversion is taking place
+ // but hopefully it's okay
nb_marked_extremities++;
}
}
@@ -579,10 +604,10 @@ void prepare_uf(std::string prefix, IBank *in, const int nb_threads, int& kmerSi
Dispatcher dispatcher (nb_threads);
- UniquifyKeys uniquifyKeys(kmerSize, pass, nb_passes, nb_threads,
+ RepartHashes repartHashes(kmerSize, pass, nb_passes, nb_threads,
nb_marked_extremities, nb_unmarked_extremities,
uf_hashes_vectors);
- dispatcher.iterate (in->iterator(), uniquifyKeys);
+ dispatcher.iterate (in->iterator(), repartHashes);
logging( std::to_string(nb_marked_extremities.load()) + " marked kmers, " + std::to_string(nb_unmarked_extremities.load()) + " unmarked kmers");
@@ -593,7 +618,7 @@ void prepare_uf(std::string prefix, IBank *in, const int nb_threads, int& kmerSi
prepareUF(it->item());*/
- logging("created vector of hashes, size approx " + std::to_string( sizeof(partition_t)*nb_marked_extremities.load()/1024/1024) + " MB)");
+ logging("created vector of hashes, size approx " + std::to_string( sizeof(uf_hashes_t)*nb_marked_extremities.load()/1024/1024) + " MB)");
ThreadPool uf_sort_pool(nb_threads); // ThreadPool
// ctpl::thread_pool uf_merge_pool(nb_threads);
@@ -602,7 +627,7 @@ void prepare_uf(std::string prefix, IBank *in, const int nb_threads, int& kmerSi
{
auto sortuniq = [&uf_hashes_vectors, i] (int thread_id)
{
- std::vector<partition_t> &vec = uf_hashes_vectors[i];
+ std::vector<uf_hashes_t> &vec = uf_hashes_vectors[i];
sort( vec.begin(), vec.end() );
vec.erase( unique( vec.begin(), vec.end() ), vec.end() );
};
@@ -613,13 +638,16 @@ void prepare_uf(std::string prefix, IBank *in, const int nb_threads, int& kmerSi
uf_sort_pool.join(); // ThreadPool
-
// a single-threaded merge and write to file before they're loaded again in bglue
- BagFile<uint64_t> * bagf = new BagFile<uint64_t>( prefix+".glue.hashes."+ to_string(pass)); LOCAL(bagf);
- Bag<uint64_t> * currentbag = new BagCache<uint64_t> ( bagf, 10000 ); LOCAL(currentbag);// really? we have to through these hoops to do a simple binary file in gatb? gotta change this.
- uint64_t nb_elts_pass = 0;
+ BagFile<uf_hashes_t> * bagf = new BagFile<uf_hashes_t>( prefix+".glue.hashes."+ to_string(pass)); LOCAL(bagf);
+ Bag<uf_hashes_t> * currentbag = new BagCache<uf_hashes_t> ( bagf, 10000 ); LOCAL(currentbag);// really? we have to through these hoops to do a simple binary file in gatb? gotta change this.
+
+ uint64_t nb_elts_pass = 0; // it's a counter
- priority_queue<std::tuple<uint64_t,int>, std::vector<std::tuple<uint64_t,int>>, std::greater<std::tuple<uint64_t,int>> > pq; // http://stackoverflow.com/questions/2439283/how-can-i-create-min-stl-priority-queue
+ // tuple is uf_hash_t, thread_id(=int)
+ priority_queue<std::tuple<uf_hashes_t,int>, std::vector<std::tuple<uf_hashes_t,int>>, std::greater<std::tuple<uf_hashes_t,int>> > pq; // http://stackoverflow.com/questions/2439283/how-can-i-create-min-stl-priority-queue
+
+ // auxiliary info for threads
vector<uint64_t> hash_vector_idx(nb_threads);
vector<uint64_t> hash_vector_size(nb_threads);
@@ -632,11 +660,11 @@ void prepare_uf(std::string prefix, IBank *in, const int nb_threads, int& kmerSi
pq.emplace(make_tuple(uf_hashes_vectors[i][hash_vector_idx[i]++], i));
}
- uint64_t prev = 0;
+ uint64_t prev = UINT64_MAX; // should be uf_hash_t but let's have it larger to mark a special "prev" initial value
while (pq.size() > 0)
{
- std::tuple<uint64_t, int> elt = pq.top(); pq.pop();
- uint64_t cur = get<0>(elt);
+ std::tuple<uf_hashes_t, int> elt = pq.top(); pq.pop();
+ uf_hashes_t cur = get<0>(elt);
//std::cout << "got " << cur << " queue " << get<1>(elt) << std::endl;
if (cur != prev)
{
@@ -658,7 +686,7 @@ void prepare_uf(std::string prefix, IBank *in, const int nb_threads, int& kmerSi
free_memory_vector(uf_hashes_vectors);
- logging("pass " + to_string(pass+1) + "/" + to_string(nb_passes) + ", " + std::to_string(nb_elts_pass) + " unique hashes written to disk, size " + to_string(nb_elts_pass* sizeof(partition_t) / 1024/1024) + " MB");
+ logging("pass " + to_string(pass+1) + "/" + to_string(nb_passes) + ", " + std::to_string(nb_elts_pass) + " unique hashes written to disk, size " + to_string(nb_elts_pass* sizeof(uf_hashes_t) / 1024/1024) + " MB");
nb_elts += nb_elts_pass;
}
@@ -739,11 +767,11 @@ void bglue(Storage *storage,
prepare_uf<SPAN>(prefix, in, nb_threads, kmerSize, pass, nb_prepare_passes, nb_elts, nb_glue_sequences);
// load uf hashes from disk
- std::vector<partition_t> uf_hashes;
+ std::vector<uf_hashes_t> uf_hashes;
uf_hashes.reserve(nb_elts);
for (int pass = 0; pass < nb_prepare_passes; pass++)
{
- IteratorFile<uint64_t> file(prefix+".glue.hashes." + to_string(pass));
+ IteratorFile<uf_hashes_t> file(prefix+".glue.hashes." + to_string(pass));
for (file.first(); !file.isDone(); file.next())
uf_hashes.push_back(file.item());
}
@@ -751,12 +779,15 @@ void bglue(Storage *storage,
uf_hashes.push_back(0);
unsigned long nb_uf_keys = uf_hashes.size();
- logging("loaded all unique UF elements (" + std::to_string(nb_uf_keys) + ") into a single vector of size " + to_string(nb_uf_keys* sizeof(partition_t) / 1024/1024) + " MB");
+ logging("loaded all unique UF elements (" + std::to_string(nb_uf_keys) + ") into a single vector of size " + to_string(nb_uf_keys* sizeof(uf_hashes_t) / 1024/1024) + " MB");
int gamma = 3; // make it even faster.
- // TODO: why not iterate the UF hashes from disk instead of loading them in memory?
- boomphf::mphf<partition_t , /*TODO we don't need hasher_t here now that we're not hashing kmers, but I forgot to change*/ hasher_t< partition_t> > uf_mphf(nb_uf_keys, uf_hashes, nb_threads, gamma, verbose);
+ // side question:
+ // why not iterate the UF hashes from disk instead of loading them in memory (to construt that mphf)?
+ // I believe it's because the UF will anyway be loaded in memory after and will also take as much space as those hashes
+
+ boomphf::mphf<uf_hashes_t, /*TODO we don't need hasher_t here now that we're not hashing kmers, but I forgot to change*/ hasher_t< uf_hashes_t> > uf_mphf(nb_uf_keys, uf_hashes, nb_threads, gamma, verbose);
free_memory_vector(uf_hashes);
@@ -766,6 +797,12 @@ void bglue(Storage *storage,
logging("UF MPHF constructed (" + std::to_string(uf_mphf_memory/8/1024/1024) + " MB)" );
}
+ if (nb_uf_keys > (uint64_t)UINT32_MAX)
+ {
+ std::cout << "cannot create a union-find data structure, too many elements. This should in fact not even happen. Please contact a BCALM developer" << std::endl;
+ exit(1);
+ }
+
// create a UF data structure
// this one stores nb_uf_keys * uint64_t (actually, atomic's). so it's bigger than uf_hashes
unionFind ufkmers(nb_uf_keys);
@@ -779,7 +816,7 @@ void bglue(Storage *storage,
// those were toy one, here is the real one:
// instead of UF of kmers, we do a union find of hashes of kmers. less memory. will have collisions, but that's okay i think. let's see.
- // actually, in the current implementation, partition_t is not used, but values are indeed hardcoded in 32 bits (the UF implementation uses a 64 bits hash table for internal stuff)
+ // actually, in the current implementation, values are indeed hardcoded in 32 bits (the UF implementation uses a 64 bits hash table but don't get fooled, it's 32 bits of values and 32 bits of internal stuff)
// We loop over sequences.
@@ -875,11 +912,12 @@ void bglue(Storage *storage,
return;
/* now we're mirroring the UF to a vector of uint32_t's, it will take less space, and strictly same information
- * this is to get rid of the rank (one uint32) per element in the current UF implementation
- * we're using the disk to save space of populating one vector from the other in memory. */
+ * this is to get rid of the rank (one uint32) per element in the current UF implementation.
+ * To do this, we're using the disk to save space of populating one vector from the other in memory.
+ * (saves having to allocate both vectors at the same time) */
- BagFile<uint64_t> *ufkmers_bagf = new BagFile<uint64_t>(prefix+".glue.uf"); LOCAL(ufkmers_bagf);
- BagCache<uint64_t> *ufkmers_bag = new BagCache<uint64_t>( ufkmers_bagf, 10000 ); LOCAL(ufkmers_bag);
+ BagFile<uf_hashes_t> *ufkmers_bagf = new BagFile<uf_hashes_t>(prefix+".glue.uf"); LOCAL(ufkmers_bagf);
+ BagCache<uf_hashes_t> *ufkmers_bag = new BagCache<uf_hashes_t>( ufkmers_bagf, 10000 ); LOCAL(ufkmers_bag);
for (unsigned long i = 0; i < nb_uf_keys; i++)
//ufkmers_vector[i] = ufkmers.find(i); // just in-memory without the disk
@@ -892,15 +930,15 @@ void bglue(Storage *storage,
ufkmers_bag->flush();
- std::vector<uint32_t > ufkmers_vector(nb_uf_keys);
- IteratorFile<uint64_t> ufkmers_file(prefix+".glue.uf");
+ std::vector<uf_hashes_t> ufkmers_vector(nb_uf_keys);
+ IteratorFile<uf_hashes_t> ufkmers_file(prefix+".glue.uf");
unsigned long i = 0;
for (ufkmers_file.first(); !ufkmers_file.isDone(); ufkmers_file.next())
ufkmers_vector[i++] = ufkmers_file.item();
System::file().remove (prefix+".glue.uf");
- logging("loaded 32-bit UF (" + to_string(nb_uf_keys*sizeof(uint32_t)/1024/1024) + " MB)");
+ logging("loaded 32-bit UF (" + to_string(nb_uf_keys*sizeof(uf_hashes_t)/1024/1024) + " MB)");
// setup output file
string output_prefix = prefix;
@@ -987,7 +1025,7 @@ void bglue(Storage *storage,
{
const string abundances = comment.substr(3);
float mean_abundance = get_mean_abundance(abundances);
- uint32_t sum_abundances = get_sum_abundance(abundances);
+ uint64_t sum_abundances = get_sum_abundance(abundances);
// km is not a standard GFA field so i'm putting it in lower case as per the spec
output(seq, out, "LN:i:" + to_string(seq.size()) + " KC:i:" + to_string(sum_abundances) + " km:f:" + to_string_with_precision(mean_abundance));
@@ -1061,7 +1099,7 @@ void bglue(Storage *storage,
outLock.unlock();
unordered_map<int, vector< markedSeq<SPAN> >> msInPart;
- uint64_t seq_index = 0;
+ seq_idx_t seq_index = 0;
for (it.first(); !it.isDone(); it.next()) // BankFasta
{
@@ -1097,14 +1135,16 @@ void bglue(Storage *storage,
seq_index++;
}
- // now iterates all sequences in a partition to determine the order in which they're going to be glues (avoid intermediate gluing)
- vector<vector<uint32_t>> ordered_sequences_idxs ;
+ vector<vector<seq_idx_t>> seqs_to_glue;
+ vector<bool> seqs_to_glue_is_circular;
+
+ // now iterates all sequences in a partition to determine the order in which they're going to be glued
for (auto it = msInPart.begin(); it != msInPart.end(); it++)
{
bool debug = false; //debug = it->first == 38145; // debug specific partition
//std::cout << "1.processing partition " << it->first << std::endl;
- determine_order_sequences<SPAN>(ordered_sequences_idxs, it->second, kmerSize, debug); // return indices of markedSeq's inside it->second
- //std::cout << "2.processing partition " << it->first << " nb ordered sequences: " << ordered_sequences_idxs.size() << std::endl;
+ determine_order_sequences<SPAN>(seqs_to_glue, seqs_to_glue_is_circular, it->second, kmerSize, debug); // return indices of markedSeq's inside it->second
+ //std::cout << "2.processing partition " << it->first << " nb of final sequences: " << seqs_to_glue.size() << std::endl;
free_memory_vector(it->second);
}
@@ -1125,10 +1165,12 @@ void bglue(Storage *storage,
abundances.push_back(abundance_str);
}
- for (auto itO = ordered_sequences_idxs.begin(); itO != ordered_sequences_idxs.end(); itO++)
+ uint64_t nb_seqs_to_glue = seqs_to_glue.size();
+ assert(seqs_to_glue_is_circular.size() == nb_seqs_to_glue);
+ for (uint64_t i = 0; i < nb_seqs_to_glue; i++)
{
string seq, abs;
- glue_sequences(*itO, sequences, abundances, kmerSize, seq, abs); // takes as input the indices of ordered sequences, and the markedSeq's themselves
+ glue_sequences(seqs_to_glue[i], seqs_to_glue_is_circular[i], sequences, abundances, kmerSize, seq, abs); // takes as input the indices of ordered sequences, whether that sequence is circular, and the markedSeq's themselves along with their abundances
float mean_abundance = get_mean_abundance(abs);
uint32_t sum_abundances = get_sum_abundance(abs);
@@ -1137,7 +1179,8 @@ void bglue(Storage *storage,
}
}
- free_memory_vector(ordered_sequences_idxs);
+ free_memory_vector(seqs_to_glue);
+ free_memory_vector(seqs_to_glue_is_circular);
partitionBank.finalize(); // BankFasta
=====================================
gatb-core/src/gatb/bcalm2/ograph.cpp
=====================================
@@ -233,7 +233,7 @@ void graph3<span>::debruijn(){
minusone.setVal(-1);
left.push_back({0,minusone, SEQ_LEFT}); // dummy kmer so that we dont need to check bounds.. clever..
right.push_back({0,minusone, SEQ_LEFT});
- uint debug_index = 0;
+ //uint debug_index = 0;
for (uint32_t i = 0; i< indiceUnitigs; i++)
{
@@ -255,7 +255,7 @@ void graph3<span>::debruijn(){
//~ if (debug_index > 0) if (kL.indice == debug_index || kR.indice == debug_index ) std::cout << " identical, kl / kR " << kL.indice << " " << kR.indice << " unitigs " << unitigs[kL.indice] << " " << unitigs[kR.indice] << " positions " << kL.position << " " << kR.position << std::endl;
//~ if(isNumber (unitigs[kL.indice][0])){
//~ }
- if(not kL.indice==kR.indice){
+ if(kL.indice!=kR.indice){
update_connected(kL);
update_connected(kR);
}
@@ -265,7 +265,7 @@ void graph3<span>::debruijn(){
++iL;++iR;
if(left[iL].kmmer==kL.kmmer){
go=false;
- if(not left[iL].indice==right[iR].indice){
+ if(left[iL].indice!=right[iR].indice){
update_connected(left[iL]);
}
//~ while(left[++iL].kmmer<=kR.kmmer ){}
@@ -277,7 +277,7 @@ void graph3<span>::debruijn(){
}
if(right[iR].kmmer==kL.kmmer){
go=false;
- if(not left[iL].indice==right[iR].indice){
+ if(left[iL].indice!=right[iR].indice){
update_connected(right[iR]);
}
//~ while(right[++iR].kmmer<=kL.kmmer ){}
=====================================
gatb-core/src/gatb/debruijn/impl/BranchingAlgorithm.cpp
=====================================
@@ -151,8 +151,8 @@ struct FunctorNodes
void operator() (Node& node)
{
// We get branching nodes neighbors for the current node.
- GraphVector<Node> successors = graph->template successors (node);
- GraphVector<Node> predecessors = graph->template predecessors (node);
+ GraphVector<Node> successors = graph->successors (node);
+ GraphVector<Node> predecessors = graph->predecessors (node);
if ( ! (successors.size()==1 && predecessors.size()==1) )
{
=====================================
gatb-core/src/gatb/debruijn/impl/GraphUnitigs.cpp
=====================================
@@ -888,7 +888,9 @@ GraphUnitigsTemplate<span>::GraphUnitigsTemplate (tools::misc::IProperties* para
/** We create a storage instance. */
/* (this is actually loading, not creating, the storage at "uri") */
BaseGraph::_storageMode = load_from_hdf5 ? STORAGE_HDF5 : STORAGE_FILE;
- BaseGraph::setStorage (StorageFactory(BaseGraph::_storageMode).create (input, false, false));
+
+ bool append = true; // in principle we wouldn't need to modify the .h5 file when building unitigs, but at some very rare locations in the code, we do (like nb_unitigs)
+ BaseGraph::setStorage (StorageFactory(BaseGraph::_storageMode).create (input, false, false, false, append));
/** We get some properties. */
BaseGraph::_state = (typename GraphUnitigsTemplate<span>::StateMask) atol (BaseGraph::getGroup().getProperty ("state").c_str());
@@ -1742,6 +1744,10 @@ simplePathLongest_avance(const NodeGU& node, Direction dir, int& seqLength, int&
GraphVector<EdgeGU> in_neighbors_vec = this->neighborsEdge (neighbors[0].to, reverse(dir));
int in_neighbors = in_neighbors_vec.size();
+
+ if (in_neighbors == 0)
+ std:: cout << "simplePathLongest_avance stopped at " << toString(cur_node) << " because of surprising lack of in-neighbor given that we had to come from somewhere" << std::endl;
+
assert(in_neighbors >= 1);
/** We check we have no in-branching. */
if (in_neighbors > 1)
=====================================
gatb-core/src/gatb/debruijn/impl/IterativeExtensions.cpp
=====================================
@@ -80,7 +80,7 @@ struct NodeDepth
template<size_t span, typename Node, typename Edge, typename Graph>
bool IterativeExtensions<span, Node, Edge, Graph>::compare_and_mark_last_k_minus_one_mer (const string& node, set<kmer_type>& kmers_set)
{
- KmerModel leftKmer = modelMinusOne.codeSeed (node.c_str(), Data::ASCII, node.size() - modelMinusOne.getKmerSize());
+ KmerModelDirect leftKmer = modelMinusOne.codeSeed (node.c_str(), Data::ASCII, node.size() - modelMinusOne.getKmerSize());
kmer_type kmer = leftKmer.value();
if (kmers_set.find(kmer) != kmers_set.end())
@@ -115,7 +115,7 @@ IterativeExtensions<span, Node, Edge, Graph>::IterativeExtensions (
max_depth(max_depth), max_nodes(max_nodes)
{
model = Model (graph.getKmerSize() );
- modelMinusOne = Model (graph.getKmerSize()-1);
+ modelMinusOne = ModelDirect (graph.getKmerSize()-1);
}
/*********************************************************************
@@ -168,7 +168,7 @@ void IterativeExtensions<span, Node, Edge, Graph>::construct_linear_seqs (
#ifndef DONTMARK
set<kmer_type> already_extended_from;
- // compare_and_mark_last_k_minus_one_mer(L, already_extended_from); // mark first kmer to never extend from it again, // L will be marked at first iteration below
+ compare_and_mark_last_k_minus_one_mer(L, already_extended_from); // mark first kmer to never extend from it again
#endif
/** We will need a Path object and a Sequence object during the extension. */
@@ -271,12 +271,14 @@ void IterativeExtensions<span, Node, Edge, Graph>::construct_linear_seqs (
}
#ifndef DONTMARK
- // make sure this is the only time we see this (k-1)-overlap
- bool already_seen = compare_and_mark_last_k_minus_one_mer (graph.toString(ksd.node), already_extended_from);
- if (already_seen)
- {
- INFO (("... XXX not extending node %s becaues last k-1-mer was already seen\n", seq.toString().c_str()));
- continue;
+ if(len_right>0){ //==0 particular case of the first contig, contig already marked (before the loop), do not mark again, this will stop extension
+ // make sure this is the only time we see this (k-1)-overlap (end of the just built contig)
+ bool already_seen = compare_and_mark_last_k_minus_one_mer (graph.toString(endNode), already_extended_from);
+ if (already_seen)
+ {
+ INFO (("... XXX not extending node %s becaues last k-1-mer was already seen\n", seq.toString().c_str()));
+ continue;
+ }
}
#endif
=====================================
gatb-core/src/gatb/debruijn/impl/IterativeExtensions.hpp
=====================================
@@ -46,6 +46,8 @@ public:
/** ShortcutS. */
typedef typename kmer::impl::Kmer<span>::ModelCanonical Model;
typedef typename kmer::impl::Kmer<span>::ModelCanonical::Kmer KmerModel;
+ typedef typename kmer::impl::Kmer<span>::ModelDirect ModelDirect;
+ typedef typename kmer::impl::Kmer<span>::ModelDirect::Kmer KmerModelDirect;
typedef typename kmer::impl::Kmer<span>::Type kmer_type;
/** Constructor.
@@ -109,7 +111,7 @@ private:
int max_nodes;
Model model;
- Model modelMinusOne;
+ ModelDirect modelMinusOne;
/** Fill a Sequence instance from the results of the current graph traversal.
* \param[in] node : starting node of the path
=====================================
gatb-core/src/gatb/debruijn/impl/LinkTigs.cpp
=====================================
@@ -30,7 +30,7 @@ using namespace gatb::core::system::impl;
namespace gatb { namespace core { namespace debruijn { namespace impl {
static constexpr int nb_passes = 8;
- static void write_final_output(const string& unitigs_filename, bool verbose, BankFasta* out, uint64_t &nb_unitigs);
+ static void write_final_output(const string& unitigs_filename, bool verbose, BankFasta* out, uint64_t &nb_unitigs, bool renumber_unitigs);
static bool get_link_from_file(std::ifstream& input, std::string &link, uint64_t &unitig_id);
/* this procedure finds the overlaps between unitigs, using a hash table of all extremity (k-1)-mers
@@ -42,30 +42,68 @@ namespace gatb { namespace core { namespace debruijn { namespace impl {
*
* so the memory usage is just that of the hash tables that record kmers, not of the links
*
- * Assumption: FASTA header of unitigs starts with a unique number (unitig ID).
- * Normally bcalm outputs consecutive unitig ID's but LinkTigs could also work with non-consecutive, non-sorted IDs
+ * Two modes of operation:
+ *
+ * renumber_unitigs == true: FASTA header can be anything. Useful for any program that has removed some unitigs, e.g. merci.
+ * LinkTigs will take the header and split it into space-separated fields, remove the first field and keep the remaining ones.
+ * The first field will be replaced by numbered IDs in consecutive order, between 0 and |nb_unitigs|-1.
+ *
+ * renumber_unitigs == false: then FASTA headers of unitigs _needs_ to start with a unique number (unitig ID).
+ * Normally bcalm outputs consecutive unitig ID's but LinkTigs can also work with non-consecutive, non-sorted IDs
*/
template<size_t span>
-void link_tigs(string unitigs_filename, int kmerSize, int nb_threads, uint64_t &nb_unitigs, bool verbose)
+void link_tigs(string unitigs_filename, int kmerSize, int nb_threads, uint64_t &nb_unitigs, bool verbose, bool renumber_unitigs)
{
bcalm_logging = verbose;
- BankFasta* out = new BankFasta(unitigs_filename+".indexed");
- if (kmerSize < 4) { std::cout << "error, recent optimizations (specifically link_unitigs) don't support k<5 for now" << std::endl; exit(1); }
+ BankFasta* out = new BankFasta(unitigs_filename+".linked");
+ if (kmerSize < 4) { std::cout << "error, link_unitigs doesn't support k<5, sorry. Contact a developer if you really need k<4 support (alternatively: construct that tiny dBG using Python :)" << std::endl; exit(1); }
logging("Finding links between unitigs");
for (int pass = 0; pass < nb_passes; pass++)
- link_unitigs_pass<span>(unitigs_filename, verbose, pass, kmerSize);
+ link_unitigs_pass<span>(unitigs_filename, verbose, pass, kmerSize, renumber_unitigs);
- write_final_output(unitigs_filename, verbose, out, nb_unitigs);
+ write_final_output(unitigs_filename, verbose, out, nb_unitigs, renumber_unitigs);
delete out;
system::impl::System::file().remove (unitigs_filename);
- system::impl::System::file().rename (unitigs_filename+".indexed", unitigs_filename);
+ system::impl::System::file().rename (unitigs_filename+".linked", unitigs_filename);
logging("Done finding links between unitigs");
}
+/* strip L:'s from a comment line*/
+static string remove_previous_links(string &header)
+{
+ bool debug = false;
+ if (debug) std::cout << "parsing unitig links for " << header << std::endl;
+ string res = "";
+ std::stringstream stream(header);
+ while(1) {
+ string tok;
+ stream >> tok;
+ if(!stream)
+ break;
+
+ string field = tok.substr(0,2);
+
+ if (field != "L:")
+ {
+ res += tok + " ";
+ }
+ }
+ //res = res.substr(0,res.size()-2);
+ if (debug) std::cout << "returning " << res<< std::endl;
+ return res;
+}
+
+/* keep all the comment line except the first field*/
+static string strip_first_field(string &header)
+{
+ return header.substr(header.find(' ')+1);
+}
+
+
/*
* takes all the prefix.links.* files, sorted by unitigs.
* do a n-way merge to gather the links for each unitig in unitig order
@@ -73,7 +111,7 @@ void link_tigs(string unitigs_filename, int kmerSize, int nb_threads, uint64_t &
*/
-static void write_final_output(const string& unitigs_filename, bool verbose, BankFasta* out, uint64_t &nb_unitigs)
+static void write_final_output(const string& unitigs_filename, bool verbose, BankFasta* out, uint64_t &nb_unitigs, bool renumber_unitigs)
{
logging("gathering links from disk");
std::ifstream* inputLinks[nb_passes];
@@ -88,7 +126,8 @@ static void write_final_output(const string& unitigs_filename, bool verbose, Ban
string cur_links, seq, comment;
seq = itSeq->toString();
comment = itSeq->getComment();
-
+ comment = remove_previous_links(comment);
+
for (int pass = 0; pass < nb_passes; pass++)
{
string link; uint64_t unitig;
@@ -102,6 +141,10 @@ static void write_final_output(const string& unitigs_filename, bool verbose, Ban
uint64_t last_unitig = 0;
nb_unitigs = 0; // passed variable
+ bool first_one = true;
+ if (renumber_unitigs)
+ comment = to_string(nb_unitigs) + " " + strip_first_field(comment);
+
// nb_passes-way merge sort
while ((!finished.all()) || pq.size() > 0)
@@ -110,6 +153,12 @@ static void write_final_output(const string& unitigs_filename, bool verbose, Ban
int pass = get<1>(cur);
uint64_t unitig = get<0>(cur);
+ if (first_one) // handles the case where first unitig isn't labeled 0
+ {
+ last_unitig = unitig;
+ first_one = false;
+ }
+
if (unitig != last_unitig)
{
Sequence s (Data::ASCII);
@@ -123,6 +172,9 @@ static void write_final_output(const string& unitigs_filename, bool verbose, Ban
itSeq.next();
seq = itSeq->toString();
comment = itSeq->getComment();
+ comment = remove_previous_links(comment);
+ if (renumber_unitigs)
+ comment = to_string(nb_unitigs) + " " + strip_first_field(comment);
}
cur_links += get<2>(cur);
@@ -213,12 +265,13 @@ static void record_links(uint64_t utig_id, int pass, const string &link, std::of
template<size_t span>
-void link_unitigs_pass(const string unitigs_filename, bool verbose, const int pass, const int kmerSize)
+void link_unitigs_pass(const string unitigs_filename, bool verbose, const int pass, const int kmerSize, const bool renumber_unitigs)
{
typedef typename kmer::impl::Kmer<span>::ModelCanonical Model;
typedef typename kmer::impl::Kmer<span>::Type Type;
bool debug = false;
+ uint64_t utig_counter = 0;
BankFasta inputBank (unitigs_filename);
BankFasta::Iterator itSeq (inputBank);
@@ -239,6 +292,8 @@ void link_unitigs_pass(const string unitigs_filename, bool verbose, const int pa
const string& comment = itSeq->getComment();
unsigned long utig_id = std::stoul(comment.substr(0, comment.find(' ')));
+
+ if (renumber_unitigs) utig_id = utig_counter;
if (is_in_pass(seq, pass, UNITIG_BEGIN, kmerSize))
{
@@ -257,6 +312,8 @@ void link_unitigs_pass(const string unitigs_filename, bool verbose, const int pa
utigs_links_map[kmerEnd.value()].push_back(eEnd.pack());
// there is no UNITIG_BOTH here because we're taking (k-1)-mers.
}
+
+ utig_counter++;
}
std::ofstream links_file(unitigs_filename+".links." +to_string(pass));
@@ -264,6 +321,9 @@ void link_unitigs_pass(const string unitigs_filename, bool verbose, const int pa
uint64_t nb_hashed_entries = 0;
for (auto v : utigs_links_map)
nb_hashed_entries += v.second.size();
+
+ utig_counter = 0;
+
logging("step 2 (" + to_string(utigs_links_map.size()) + "kmers/" + to_string(nb_hashed_entries) + "extremities)");
for (itSeq.first(); !itSeq.isDone(); itSeq.next())
@@ -271,6 +331,8 @@ void link_unitigs_pass(const string unitigs_filename, bool verbose, const int pa
const string& seq = itSeq->toString();
const string& comment = itSeq->getComment();
unsigned long utig_id = std::stoul(comment.substr(0, comment.find(' ')));
+
+ if (renumber_unitigs) utig_id = utig_counter;
if (debug) std::cout << "unitig " << std::to_string(utig_id) << " : " << seq << std::endl;
@@ -378,6 +440,8 @@ void link_unitigs_pass(const string unitigs_filename, bool verbose, const int pa
}
record_links(utig_id, pass, out_links, links_file);
}
+
+ utig_counter++;
}
}
=====================================
gatb-core/src/gatb/debruijn/impl/LinkTigs.hpp
=====================================
@@ -30,10 +30,10 @@ namespace gatb { namespace core { namespace debruijn { namespace impl {
template<size_t SPAN>
- void link_tigs( std::string prefix, int kmerSize, int nb_threads, uint64_t &nb_unitigs, bool verbose);
+ void link_tigs( std::string prefix, int kmerSize, int nb_threads, uint64_t &nb_unitigs, bool verbose, bool renumber_unitigs = false);
template<size_t span>
- void link_unitigs_pass(const std::string unitigs_filename, bool verbose, const int pass, const int kmerSize);
+ void link_unitigs_pass(const std::string unitigs_filename, bool verbose, const int pass, const int kmerSize, const bool renumber_unitigs);
}}}}
=====================================
gatb-core/src/gatb/debruijn/impl/Simplifications.cpp
=====================================
@@ -1108,7 +1108,8 @@ void Simplifications<GraphType,Node,Edge>::heuristic_most_covered_path_unitigs(
current_node = neighbors[0].to;
if (_graph.degree(current_node, reverse(dir)) <= 1)
{
- std::cout << "Weird, there was supposed to be an in-neighbor. Maybe there's a loop. Remove this print if it never happens" << std::endl;
+ // actually this happens for kmers that are linked to themselves (e.g. k=5, AAATT, is connected to its own revcomp AATTT, so it's fine to end the simple path here on those cases)
+ //std::cout << "Weird, there was supposed to be an in-neighbor. Maybe there's a loop. Remove this print if it never happens" << "in degree " << _graph.degree(current_node, reverse(dir)) << " node" << _graph.toString(current_node) << std::endl;
return;
}
@@ -1280,7 +1281,7 @@ unsigned long Simplifications<GraphType,Node,Edge>::removeBulges()
unsigned int k = _graph.getKmerSize();
unsigned int maxBulgeLength = std::max((unsigned int)((double)k * _bulgeLen_kMult), (unsigned int)(k + _bulgeLen_kAdd)); // SPAdes, exactly
unsigned int backtrackingLimit = k+_bulgeAltPath_kAdd;//maxBulgeLength; // arbitrary, but if too high it will take much time; // with unitigs, no reason that it has to depend on k, but for some reason, setting it to just "k" doesnt remove nearly as many bulges as k=20. todo investigate that someday.
- unsigned int altPathCovMult = _bulgeAltPath_covMult;
+ double altPathCovMult = _bulgeAltPath_covMult;
// stats
//
=====================================
gatb-core/src/gatb/debruijn/impl/Simplifications.hpp
=====================================
@@ -66,7 +66,7 @@ public:
double _bulgeLen_kMult;
unsigned int _bulgeLen_kAdd;
unsigned int _bulgeAltPath_kAdd;
- unsigned int _bulgeAltPath_covMult;
+ double _bulgeAltPath_covMult;
double _ecLen_kMult;
double _ecRCTCcutoff;
=====================================
gatb-core/src/gatb/debruijn/impl/Terminator.cpp
=====================================
@@ -46,7 +46,7 @@ BranchingTerminatorTemplate<Node,Edge,Graph>::BranchingTerminatorTemplate (const
: TerminatorTemplate<Node,Edge,Graph> (graph)
{
/** We loop over the branching nodes. */
- GraphIterator<BranchingNode_t<Node> > itBranching = this->_graph.Graph::template iteratorBranching();
+ GraphIterator<BranchingNode_t<Node> > itBranching = this->_graph.iteratorBranching();
for (itBranching.first(); !itBranching.isDone(); itBranching.next())
{
/** We add the current branching node into the map. */
@@ -196,7 +196,7 @@ void BranchingTerminatorTemplate<Node,Edge,Graph>::mark (Node& node)
}
/** We loop the neighbors edges of the current node. */
- GraphVector<Edge> neighbors = this->_graph.template neighborsEdge(node.kmer);
+ GraphVector<Edge> neighbors = this->_graph.neighborsEdge(node.kmer);
/** We loop the branching neighbors. */
for (size_t i=0; i<neighbors.size(); i++)
=====================================
gatb-core/src/gatb/kmer/impl/ConfigurationAlgorithm.cpp
=====================================
@@ -252,10 +252,17 @@ void ConfigurationAlgorithm<span>::execute ()
/** We get some information about the bank. */
_bank->estimate (_config._estimateSeqNb, _config._estimateSeqTotalSize, _config._estimateSeqMaxSize);
+
+ //printf("_estimateSeqNb %llu _estimateSeqTotalSize %llu _estimateSeqMaxSize %llu \n", _config._estimateSeqNb, _config._estimateSeqTotalSize, _config._estimateSeqMaxSize);
/** We get the number of sub banks. */
_config._nb_banks = _bank->getCompositionNb();
+ if(_config._nb_banks == 1 )
+ {
+ _config._solidityKind = KMER_SOLIDITY_SUM;
+ }
+
/** We memorize the number of abundance min values set by the user.
* Note that it can be lower than the number of banks. */
_config._abundanceUserNb = _config._abundance.size();
@@ -311,11 +318,14 @@ void ConfigurationAlgorithm<span>::execute ()
_config._volume = _config._kmersNb * sizeof(Type) / MBYTE; // in MBytes
+ //printf("estimated usedSeqLen %llu _estimateSeqNb %llu _kmersNb %llu _volume %llu \n", usedSeqLen, _config._estimateSeqNb, kmersNb, _config._volume);
+
if (_config._volume == 0) { _config._volume = 1; } // tiny files fix
u_int64_t volume_minim = _config._volume * 0.5 *1.2 ; //0.5 for using kxmers 1.2 if bad repartition of minimizers ( todo sampling to assert ram usage)
if (volume_minim == 0) { volume_minim = 1; } // tiny files fix
+ // volume_minim is used a bit later
/** We get max(75%, 100% - X GB) */
if (_config._max_disk_space == 0) { _config._max_disk_space = std::max ((75*_config._available_space)/100, _config._available_space-available_space_min); }
@@ -341,7 +351,7 @@ void ConfigurationAlgorithm<span>::execute ()
//_nb_passes = 1; //do not constrain nb passes on disk space anymore (anyway with minim, not very big)
//increase it only if ram issue
- //printf("_volume %lli volume_minim %lli _max_disk_space %lli _nb_passes init %i \n", _volume,volume_minim,_max_disk_space,_nb_passes);
+ //printf("_volume %lli volume_minim %lli _max_disk_space %lli _nb_passes init %i \n", _config._volume,volume_minim, _config._max_disk_space, _config._nb_passes);
size_t max_open_files = System::file().getMaxFilesNumber() / 2;
=====================================
gatb-core/src/gatb/kmer/impl/DebloomAlgorithm.cpp
=====================================
@@ -337,6 +337,9 @@ void DebloomAlgorithm<span>::execute_aux (
/** We may have reach the maximum number of items in the partition. */
if (partition.size() >= partition.getMaxNbItems())
{
+ //first erase the destination file because BagFile creator does not erase previous file anymore
+ system::impl::System::file().remove (outputUri);
+
/** We exclude the partition content from the critical false positive file. */
end_debloom_partition (
partition,
=====================================
gatb-core/src/gatb/kmer/impl/MPHFAlgorithm.cpp
=====================================
@@ -249,17 +249,22 @@ void MPHFAlgorithm<span,Abundance_t,NodeState_t>::populate ()
/** We get the abundance of the current kmer. */
int abundance = itKmers->item().abundance;
- if (abundance > max_abundance_discrete)
- {
- _nb_abundances_above_precision++;
- //std::cout << "found abundance larger than discrete: " << abundance << std::endl;
- abundance = max_abundance_discrete;
- }
+
+ int idx ;
+ if (abundance >= max_abundance_discrete)
+ {
+ _nb_abundances_above_precision++;
+ //std::cout << "found abundance larger than discrete: " << abundance << std::endl;
+ idx = _abundanceDiscretization.size() -2 ;
+ }
+ else
+ {
+ //get first cell strictly greater than abundance
+ std::vector<int>::iterator up = std::upper_bound(_abundanceDiscretization.begin(), _abundanceDiscretization.end(), abundance);
+ up--; // get previous cell
+ idx = up- _abundanceDiscretization.begin() ;
+ }
- //get first cell strictly greater than abundance
- std::vector<int>::iterator up = std::upper_bound(_abundanceDiscretization.begin(), _abundanceDiscretization.end(), abundance);
- up--; // get previous cell
- int idx = up- _abundanceDiscretization.begin() ;
/** We set the abundance of the current kmer. */
_abundanceMap->at (h) = idx;
=====================================
gatb-core/src/gatb/kmer/impl/Model.hpp
=====================================
@@ -1434,7 +1434,7 @@ struct Kmer
//reste du newbyte avec le superk
- int skid =1;
+ uint skid =1;
while(true)
{
=====================================
gatb-core/src/gatb/kmer/impl/PartiInfo.cpp
=====================================
@@ -24,7 +24,9 @@
using namespace std;
#define DEBUG(a) //printf a
-// wanted to debug separately, DEBUG has more output than debug2
+// wanted to debug things separately:
+// now there are many different and mutually exclusive debug messages
+// DEBUG has more output than debug2
// to see top20 repartitions of bins, uncomment IFDEBUG and DEBUG2
#define IFDEBUG(a) //a
#define DEBUG2(a) //printf a
@@ -229,12 +231,11 @@ void Repartitor::load (tools::storage::impl::Group& group)
tools::storage::impl::Storage::istream is (group, "minimRepart");
is.read ((char*)&_nbpart, sizeof(_nbpart));
- is.read ((char*)&_mm, sizeof(_mm));
is.read ((char*)&_nb_minims, sizeof(_nb_minims));
is.read ((char*)&_nbPass, sizeof(_nbPass));
- DEBUG (("[Repartitor::load] : _nbpart=%d _mm=%d _nb_minims=%d _nbPass=%d \n",
- _nbpart, _mm, _nb_minims, _nbPass
+ DEBUG (("[Repartitor::load] : _nbpart=%d _nb_minims=%d _nbPass=%d \n",
+ _nbpart, _nb_minims, _nbPass
));
/** We allocate a table whose size is the number of possible minimizers. */
@@ -269,15 +270,14 @@ void Repartitor::load (tools::storage::impl::Group& group)
*********************************************************************/
void Repartitor::save (tools::storage::impl::Group& group)
{
- DEBUG (("[Repartitor::save] : _nbpart=%d _mm=%d _nb_minims=%d _nbPass=%d \n",
- _nbpart, _mm, _nb_minims, _nbPass
+ DEBUG (("[Repartitor::save] : _nbpart=%d _nb_minims=%d _nbPass=%d \n",
+ _nbpart, _nb_minims, _nbPass
));
bool hasMinimizerFrequencies = _freq_order != NULL;
tools::storage::impl::Storage::ostream os (group, "minimRepart");
os.write ((const char*)&_nbpart, sizeof(_nbpart));
- os.write ((const char*)&_mm, sizeof(_mm));
os.write ((const char*)&_nb_minims, sizeof(_nb_minims));
os.write ((const char*)&_nbPass, sizeof(_nbPass));
os.write ((const char*)_repart_table.data(), sizeof(Value) * _nb_minims);
@@ -304,9 +304,8 @@ void Repartitor::save (tools::storage::impl::Group& group)
*********************************************************************/
void Repartitor::printInfo ()
{
- size_t nbMinimizers = 1 << (_mm*2);
- printf("Repartitor : nbMinimizers=%ld\n", nbMinimizers);
- for (u_int64_t ii=0; ii<nbMinimizers; ii++ ) { printf(" table[%ld] = %d \n",ii,_repart_table[ii]); }
+ printf("Repartitor : nbMinimizers=%ld\n", _nb_minims);
+ for (u_int64_t ii=0; ii<_nb_minims; ii++ ) { printf(" table[%ld] = %d \n",ii,_repart_table[ii]); }
}
/********************************************************************************/
=====================================
gatb-core/src/gatb/kmer/impl/PartiInfo.hpp
=====================================
@@ -24,6 +24,7 @@
#include <gatb/kmer/impl/Model.hpp>
#include <gatb/tools/storage/impl/Storage.hpp>
+#include <mutex>
#include <queue>
/********************************************************************************/
@@ -42,16 +43,43 @@ namespace impl {
template <size_t xmer>
class PartiInfo
{
+ struct mmer_bin_record {
+ uint64_t nb_superks, nb_kmers, nb_kxmers;
+ mmer_bin_record& operator+=(const mmer_bin_record& other) {
+ nb_superks += other.nb_superks;
+ nb_kmers += other.nb_kmers;
+ nb_kxmers += other.nb_kxmers;
+ return *this;
+ }
+ };
+
+ struct parti_record {
+ uint64_t nb_kmers, nb_kxmers;
+ uint64_t nbk_per_radix[xmer * 256];
+
+ uint64_t& getNbKmer(int radix, int xx)
+ { return nbk_per_radix[xx * 256 + radix]; }
+
+ parti_record& operator+=(const parti_record& other) {
+ nb_kmers += other.nb_kmers;
+ nb_kxmers += other.nb_kxmers;
+ for(int i = 0 ; i < xmer * 256 ; i++) {
+ nbk_per_radix[i] += other.nbk_per_radix[i];
+ }
+ return *this;
+ }
+ };
+
public:
inline void incKmer (int numpart, u_int64_t val=1)
{
- _nb_kmers_per_parti[numpart]+=val;
+ _parti_records[numpart].nb_kmers+=val;
}
inline void incKxmer(int numpart, u_int64_t val=1)
{
- _nb_kxmers_per_parti[numpart]+=val; //now used to store number of kxmers per parti
+ _parti_records[numpart].nb_kxmers+=val; //now used to store number of kxmers per parti
}
//superksize in number of kmer inside this superk
@@ -60,73 +88,64 @@ public:
{
_nb_superk_total+=val;
_nb_kmer_total += val*superksize;
- _superk_per_mmer_bin[numbin]+= val ;
- _kmer_per_mmer_bin[numbin]+= val *superksize;
+ _mmer_bin_records[numbin].nb_superks+= val ;
+ _mmer_bin_records[numbin].nb_kmers+= val *superksize;
}
//kxmer count (regardless of x size), used for ram
inline void incKxmer_per_minimBin(int numbin, u_int64_t val=1)
{
- _kxmer_per_mmer_bin[numbin]+= val ;
+ _mmer_bin_records[numbin].nb_kxmers+= val ;
}
//numaprt, radix, size of kx mer
inline void incKmer_and_rad(int numpart, int radix,int x, u_int64_t val=1) //appele ds vrai loop
{
- _nb_kxmers_per_parti[numpart] += val; // now used to store number of kxmers per parti
- _nb_kmers_per_parti [numpart] += (val * (x+1)); // number of 'real' kmers
- _nbk_per_radix_per_part[x][radix][numpart]+=val; // contains number of kx mer per part per radix per x
+ parti_record& record = _parti_records[numpart];
+ record.nb_kxmers += val; // now used to store number of kxmers per parti
+ record.nb_kmers += (val * (x+1)); // number of 'real' kmers
+ record.getNbKmer(radix, x)+=val; // contains number of kx mer per part per radix per x
}
- /** */
- PartiInfo& operator+= (const PartiInfo& other)
- {
- //add other parti info , synced
-
- for (int np=0; np<_nbpart; np++)
- {
- for (size_t xx=0; xx<xmer; xx++)
- {
- for (int rad=0; rad<256; rad++)
- {
- __sync_fetch_and_add( & (_nbk_per_radix_per_part[xx][rad][np]), other.getNbKmer(np,rad,xx) );
- }
- }
-
- __sync_fetch_and_add (_nb_kmers_per_parti + np, other.getNbKmer (np));
- __sync_fetch_and_add (_nb_kxmers_per_parti + np, other.getNbSuperKmer (np));
+ /** Adds count to *this from other */
+ PartiInfo &add(const PartiInfo &other) {
+ for (int np = 0; np < _nbpart; np++) {
+ _parti_records[np] += other._parti_records[np];
}
- for (u_int64_t ii=0; ii< _num_mm_bins; ii++)
- {
- __sync_fetch_and_add (_superk_per_mmer_bin + ii, other.getNbSuperKmer_per_minim (ii));
- __sync_fetch_and_add (_kmer_per_mmer_bin + ii, other.getNbKmer_per_minim (ii));
- __sync_fetch_and_add (_kxmer_per_mmer_bin + ii, other.getNbKxmer_per_minim (ii));
+ for (int ii = 0; ii < _num_mm_bins; ii++) {
+ _mmer_bin_records[ii] += other._mmer_bin_records[ii];
}
- __sync_fetch_and_add(&_nb_superk_total , other._nb_superk_total);
- __sync_fetch_and_add(&_nb_kmer_total , other._nb_kmer_total);
-
-
+ _nb_superk_total += other._nb_superk_total;
+ _nb_kmer_total += other._nb_kmer_total;
+
return *this;
}
+ /** Adds count to *this from other while preventing concurrent access to *this */
+ PartiInfo &add_sync(const PartiInfo &other) {
+ // Locking is way faster than doing individuals "lock add"s for all the fields
+ std::lock_guard<std::mutex> guard{_mut};
+ return add(other);
+ }
+
/** */
inline u_int64_t getNbKmer(int numpart) const
{
- return _nb_kmers_per_parti[numpart];
+ return _parti_records[numpart].nb_kmers;
}
/** Get nbk in bin radix of parti numpart */
inline u_int64_t getNbKmer(int numpart, int radix, int xx) const
{
- return _nbk_per_radix_per_part[xx][radix][numpart];
+ return _parti_records[numpart].getNbKmer(radix, xx);
}
/** */
inline u_int64_t getNbSuperKmer(int numpart) const //now used for number of kxmers (indistinctive of xsize)
{
- return _nb_kxmers_per_parti[numpart];
+ return _parti_records[numpart].nb_kxmers;
}
inline u_int64_t getNbSuperKmerTotal() const
@@ -141,37 +160,26 @@ public:
/** */
inline u_int64_t getNbSuperKmer_per_minim(int numbin) const
{
- return _superk_per_mmer_bin[numbin];
+ return _mmer_bin_records[numbin].nb_superks;
}
/** */
inline u_int64_t getNbKmer_per_minim(int numbin) const
{
- return _kmer_per_mmer_bin[numbin];
+ return _mmer_bin_records[numbin].nb_kmers;
}
/** */
inline u_int64_t getNbKxmer_per_minim(int numbin) const
{
- return _kxmer_per_mmer_bin[numbin];
+ return _mmer_bin_records[numbin].nb_kxmers;
}
/** */
void clear()
{
- memset (_nb_kmers_per_parti, 0, _nbpart * sizeof(u_int64_t));
- memset (_nb_kxmers_per_parti, 0, _nbpart * sizeof(u_int64_t));
- memset (_superk_per_mmer_bin, 0, _num_mm_bins * sizeof(u_int64_t));
- memset (_kmer_per_mmer_bin, 0, _num_mm_bins * sizeof(u_int64_t));
- memset (_kxmer_per_mmer_bin, 0, _num_mm_bins * sizeof(u_int64_t));
-
- for (size_t xx=0; xx<xmer; xx++)
- {
- for(int ii=0; ii<256; ii++)
- {
- memset(_nbk_per_radix_per_part[xx][ii], 0, _nbpart*sizeof(u_int64_t));
- }
- }
+ memset (_parti_records, 0, _nbpart * sizeof(parti_record));
+ memset (_mmer_bin_records, 0, _num_mm_bins * sizeof(mmer_bin_record));
}
/** */
@@ -209,12 +217,8 @@ public:
for (int np=0; np<_num_mm_bins; np++)
{
- typedef typename Kmer<31>::Type Typem; //should be kmer size
- Typem cur;
- cur.setVal(np);
-
if(this->getNbSuperKmer_per_minim(np)!=0 || this->getNbKmer_per_minim(np) !=0 )
- printf("Bin[%5i (%s) ]= %lli %lli\n",np,cur.toString(_mm).c_str(), this->getNbSuperKmer_per_minim(np),this->getNbKmer_per_minim(np) );
+ printf("Bin[%5i]= %lli %lli\n",np, this->getNbSuperKmer_per_minim(np),this->getNbKmer_per_minim(np) );
sumk += this->getNbKmer_per_minim(np);
sumsuperk += this->getNbSuperKmer_per_minim(np);
@@ -225,49 +229,27 @@ public:
}
/** Constructor. */
- PartiInfo(int nbpart, int minimsize) : _nbpart(nbpart), _mm(minimsize)
+ PartiInfo(int nbpart, int minimsize) : _nbpart(nbpart)
{
_nb_superk_total =0;
_nb_kmer_total =0;
- _nb_kmers_per_parti = (u_int64_t*) CALLOC (nbpart, sizeof(u_int64_t));
- _nb_kxmers_per_parti = (u_int64_t*) CALLOC (nbpart, sizeof(u_int64_t));
- _num_mm_bins = 1 << (2*_mm);
- _superk_per_mmer_bin = (u_int64_t*) CALLOC (_num_mm_bins, sizeof(u_int64_t));
- _kmer_per_mmer_bin = (u_int64_t*) CALLOC (_num_mm_bins, sizeof(u_int64_t));
- _kxmer_per_mmer_bin = (u_int64_t*) CALLOC (_num_mm_bins, sizeof(u_int64_t));
-
- for(size_t xx=0; xx<xmer; xx++)
- {
- for(int ii=0; ii<256; ii++)
- {
- _nbk_per_radix_per_part[xx][ii] = (u_int64_t*) CALLOC (nbpart, sizeof(u_int64_t));
- }
- }
+ _parti_records = (parti_record*) CALLOC (_nbpart, sizeof(parti_record));
+ _num_mm_bins = 1 << (2*minimsize);
+ _mmer_bin_records = (mmer_bin_record*) CALLOC (_num_mm_bins, sizeof(mmer_bin_record));
}
+
/** Constructor (copy). Needed for fillparti class. */
PartiInfo(const PartiInfo& cr)
//the copy contr realloc its own arrays, zero init
{
_num_mm_bins = cr._num_mm_bins;
_nbpart = cr._nbpart;
- _mm = cr._mm;
_nb_superk_total = cr._nb_superk_total;
_nb_kmer_total = cr._nb_kmer_total;
- _nb_kmers_per_parti = (u_int64_t*) CALLOC (_nbpart, sizeof(u_int64_t));
- _nb_kxmers_per_parti = (u_int64_t*) CALLOC (_nbpart, sizeof(u_int64_t));
- _superk_per_mmer_bin = (u_int64_t*) CALLOC (_num_mm_bins, sizeof(u_int64_t));
- _kmer_per_mmer_bin = (u_int64_t*) CALLOC (_num_mm_bins, sizeof(u_int64_t));
- _kxmer_per_mmer_bin = (u_int64_t*) CALLOC (_num_mm_bins, sizeof(u_int64_t));
-
- for(size_t xx=0; xx<xmer; xx++)
- {
- for(int ii=0; ii<256; ii++)
- {
- _nbk_per_radix_per_part[xx][ii] = (u_int64_t *) CALLOC(_nbpart,sizeof(u_int64_t));
- }
- }
+ _parti_records = (parti_record*) CALLOC (_nbpart, sizeof(parti_record));
+ _mmer_bin_records = (mmer_bin_record*) CALLOC (_num_mm_bins, sizeof(mmer_bin_record));
// printf("PartiInfo copy constr %p _nb_kmers_per_parti %p \n",this,_nb_kmers_per_parti);
@@ -276,32 +258,18 @@ public:
/** Destructor. */
~PartiInfo()
{
- FREE (_nb_kmers_per_parti);
- FREE (_nb_kxmers_per_parti);
- FREE (_superk_per_mmer_bin);
- FREE (_kmer_per_mmer_bin);
- FREE (_kxmer_per_mmer_bin);
-
- for(size_t xx=0; xx<xmer; xx++) { for(int ii=0; ii<256; ii++) { FREE(_nbk_per_radix_per_part[xx][ii]); } }
-
- // printf("print info destroyed %p _nb_kmers_per_parti %p \n",this,_nb_kmers_per_parti);
+ FREE (_parti_records);
+ FREE (_mmer_bin_records);
}
private:
-
- u_int64_t* _nb_kmers_per_parti;
- u_int64_t* _nb_kxmers_per_parti; //now used to store number of kxmers per parti
- u_int64_t* _superk_per_mmer_bin;
- u_int64_t _nb_superk_total;
- u_int64_t _nb_kmer_total;
- u_int64_t* _kmer_per_mmer_bin;
- u_int64_t* _kxmer_per_mmer_bin;
-
- u_int64_t* _nbk_per_radix_per_part[xmer][256];//number of kxmer per parti per rad
+ std::mutex _mut = {}; // Locking for merging with main PartiInfo
+ u_int64_t _nbpart;
u_int64_t _num_mm_bins;
-
- int _nbpart;
- int _mm;
+ u_int64_t _nb_superk_total;
+ u_int64_t _nb_kmer_total;
+ parti_record* _parti_records;
+ mmer_bin_record* _mmer_bin_records;
};
/********************************************************************************/
@@ -331,13 +299,13 @@ public:
/** Constructor
* \param[in] nbpart : hash value will be in range [0..nbpart-1]
* \param[in] minimsize : size of the minimizers. */
- Repartitor (int nbpart=0, int minimsize=0, int nbPass=1) : _nbpart(nbpart), _mm(minimsize), _nb_minims(1 << (_mm*2)), _nbPass(nbPass), _freq_order(0)
+ Repartitor (int nbpart=0, int minimsize=0, int nbPass=1) : _nbpart(nbpart), _nb_minims(1 << (minimsize*2)), _nbPass(nbPass), _freq_order(0)
{
if (nbpart <= 0) { system::Exception("Repartitor: nbpart (%d) should be > 0", nbpart); }
}
/** Constructor */
- Repartitor (tools::storage::impl::Group& group) : _nbpart(0), _mm(0), _nb_minims(0), _nbPass(0), _freq_order(0) { this->load (group); }
+ Repartitor (tools::storage::impl::Group& group) : _nbpart(0), _nb_minims(0), _nbPass(0), _freq_order(0) { this->load (group); }
/** Destructor */
~Repartitor () { if (_freq_order) { delete[] _freq_order; } }
@@ -411,7 +379,6 @@ private:
}
u_int16_t _nbpart;
- u_int16_t _mm;
u_int64_t _nb_minims;
u_int16_t _nbPass;
std::vector<Value> _repart_table ;
=====================================
gatb-core/src/gatb/kmer/impl/RepartitionAlgorithm.cpp
=====================================
@@ -182,7 +182,7 @@ public:
size_t superKmerLen = superKmer.size();
/** We increase superkmer counter the current minimizer. */
- _local_pInfo.incSuperKmer_per_minimBin (superKmer.minimizer, superKmerLen);
+ _pInfo.incSuperKmer_per_minimBin (superKmer.minimizer, superKmerLen);
/** We loop over the kmer of the superkmer (except the first one).
* We update the pInfo each time we find a kxmer in the superkmer. */
@@ -193,7 +193,7 @@ public:
if (superKmer[ii].which() != prev_which || kx_size >= _kx) // kxmer_size = 1 //cost should diminish with larger kxmer
{
/** We increase the number of kxmer found for the current minimizer. */
- _local_pInfo.incKxmer_per_minimBin (superKmer.minimizer);
+ _pInfo.incKxmer_per_minimBin (superKmer.minimizer);
kx_size = 0;
}
else
@@ -205,7 +205,7 @@ public:
}
/** We add the pending kxmer to the bin. */
- _local_pInfo.incKxmer_per_minimBin (superKmer.minimizer);
+ _pInfo.incKxmer_per_minimBin (superKmer.minimizer);
// see if we need to stop
_nbSuperKmersSeenSoFar ++;
@@ -221,8 +221,6 @@ public:
SampleRepart (
Model& model,
Configuration& config,
- size_t nbPasses,
- size_t currentPass,
size_t nbPartitions,
IteratorListener* progress,
bool * cancelIterator,
@@ -230,24 +228,17 @@ public:
BankStats& bankStats,
PartiInfo<5>& pInfo
)
- : Sequence2SuperKmer<span> (model, nbPasses, currentPass, nbPartitions, progress, bankStats)
- ,_kx(4), _extern_pInfo(pInfo), _local_pInfo(config._nb_partitions, model.getMmersModel().getKmerSize()),
+ : Sequence2SuperKmer<span> (model, 1, 0, nbPartitions, progress, bankStats)
+ ,_kx(4), _pInfo(pInfo),
_cancelIterator(cancelIterator), _nbSeqsToSee(nbSeqsToSee), _nbSuperKmersSeenSoFar(0)
{
}
- /** Destructor. */
- ~SampleRepart ()
- {
- //add to global parti_info
- _extern_pInfo += _local_pInfo;
- }
-
private:
size_t _kx;
- PartiInfo<5>& _extern_pInfo;
- PartiInfo<5> _local_pInfo;
+ PartiInfo<5>& _pInfo;
+
bool* _cancelIterator;
size_t _nbSeqsToSee;
size_t _nbSuperKmersSeenSoFar;
@@ -332,6 +323,9 @@ void RepartitorAlgorithm<span>::computeFrequencies (Repartitor& repartitor)
u_int64_t nbseq_sample = std::min ( u_int64_t (estimateSeqNb * 0.05) ,u_int64_t( 50000000ULL) ) ;
// TODO would be better to just stop estimating minimizer frequency when it becomes stable. not after a fixed number of reads
+ if (nbseq_sample == 0)
+ nbseq_sample = 1;
+
u_int64_t rg = ((u_int64_t)1 << (2*_config._minim_size));
//cout << "\nAllocating " << ((rg*sizeof(uint32_t))/1024) << " KB for " << _minim_size <<"-mers frequency counting (" << rg << " elements total)" << endl;
uint32_t *m_mer_counts = new uint32_t[rg];
@@ -426,52 +420,14 @@ void RepartitorAlgorithm<span>::computeRepartition (Repartitor& repartitor)
Iterator<Sequence>* it = _bank->iterator(); LOCAL (it);
std::vector<Iterator<Sequence>*> itBanks = it->getComposition();
- /*
- bool dummyBoolean = false;
- SampleRepart<span> sampleRepart (
- model,
- _config,
- 1, // we don't care about the actual number of passes, we just use 1
- 0, // we don't care about the actual number of passes, the current one is 0
- _config._nb_partitions,
- NULL,
- &dummyBoolean, // will be set to true when iteration needs to be stopped
- nbseq_sample, // how many sequences we need to see
- bstatsDummy,
- sample_info
- );*/
-
- //SerialDispatcher serialDispatcher;
- //BankStats bstatsDummy;
-
for(size_t i=0; i<_config._nb_banks; i++){
- //cout << i << endl;
- //u_int64_t nbseq_sample_ = min(nbseq_sample, u_int64_t (500000));
- //TruncateIterator<Sequence>* it = new TruncateIterator<Sequence> (*(itBanks[i]), nbseq_sample_);
- //LOCAL(it);
- //serialDispatcher.iterate (it, sampleRepart);
- //Iterator<Sequence>* bankit = itBanks[i];
- //LOCAL(bankit);
CancellableIterator<Sequence>* cancellable_it = new CancellableIterator<Sequence> (*itBanks[i]);
LOCAL(cancellable_it);
- //cout << nbseq_sample << endl;
- /** We create a sequence iterator and give it a progress message */
- //Iterator<Sequence>* it_all_reads = createIterator<Sequence> (
- // cancellable_it,
- // _bank->getNbItems(),
- // Stringify::format (progressFormat0, bankShortName.c_str()).c_str()
- // );
- //LOCAL (it_all_reads);
-
-
-
/** We compute a distribution of Superkmers from a part of the bank. */
serialDispatcher.iterate (cancellable_it, SampleRepart<span> (
model,
_config,
- 1, // we don't care about the actual number of passes, we just use 1
- 0, // we don't care about the actual number of passes, the current one is 0
_config._nb_partitions,
NULL,
&(cancellable_it->_cancel), // will be set to true when iteration needs to be stopped
@@ -493,11 +449,9 @@ void RepartitorAlgorithm<span>::computeRepartition (Repartitor& repartitor)
CancellableIterator<Sequence>* cancellable_it = new CancellableIterator<Sequence> (*(it));
LOCAL(cancellable_it);
-
// how many seqs we need to see
u_int64_t nbseq_sample = std::max ( u_int64_t (_config._estimateSeqNb * 0.05) ,u_int64_t( 1000000ULL) ) ;
- //cout << nbseq_sample << endl;
/** We create a sequence iterator and give it a progress message */
Iterator<Sequence>* it_all_reads = createIterator<Sequence> (
cancellable_it,
@@ -506,7 +460,6 @@ void RepartitorAlgorithm<span>::computeRepartition (Repartitor& repartitor)
);
LOCAL (it_all_reads);
-
BankStats bstatsDummy;
/** We compute a distribution of Superkmers from a part of the bank. */
@@ -514,8 +467,6 @@ void RepartitorAlgorithm<span>::computeRepartition (Repartitor& repartitor)
serialDispatcher.iterate (it_all_reads, SampleRepart<span> (
model,
_config,
- 1, // we don't care about the actual number of passes, we just use 1
- 0, // we don't care about the actual number of passes, the current one is 0
_config._nb_partitions,
NULL,
&(cancellable_it->_cancel), // will be set to true when iteration needs to be stopped
=====================================
gatb-core/src/gatb/kmer/impl/SortingCountAlgorithm.cpp
=====================================
@@ -867,12 +867,12 @@ public:
/** Destructor. */
- ~FillPartitions ()
+ virtual ~FillPartitions ()
{
//printf("destruc fillparti _superkmerFiles %p \n",_superkmerFiles);
//add to global parti_info
- _extern_pInfo += _local_pInfo;
+ _extern_pInfo.add_sync(_local_pInfo);
}
private:
@@ -1007,12 +1007,12 @@ public:
/** Destructor. */
- ~FillPartitions ()
+ virtual ~FillPartitions ()
{
//printf("destruc fillparti _superkmerFiles %p \n",_superkmerFiles);
//add to global parti_info
- _extern_pInfo += _local_pInfo;
+ _extern_pInfo.add_sync(_local_pInfo);
}
private:
@@ -1143,12 +1143,12 @@ private:
/** Destructor. */
- ~FillPartitions ()
+ virtual ~FillPartitions ()
{
//printf("destruc fillparti _superkmerFiles %p \n",_superkmerFiles);
//add to global parti_info
- _extern_pInfo += _local_pInfo;
+ _extern_pInfo.add_sync(_local_pInfo);
}
private:
@@ -1183,6 +1183,7 @@ void SortingCountAlgorithm<span>::fillPartitions (size_t pass, Iterator<Sequence
DEBUG (("SortingCountAlgorithm<span>::fillPartitions _kmerSize=%d _minim_size=%d \n", _config._kmerSize, _config._minim_size));
+ _nbKmersPerPartitionPerBank.clear();
if(_config._solidityKind != KMER_SOLIDITY_SUM)
{
/** We delete the previous partitions storage. */
@@ -1226,72 +1227,79 @@ void SortingCountAlgorithm<span>::fillPartitions (size_t pass, Iterator<Sequence
/** We have to reinit the progress instance since it may have been used by SampleRepart before. */
_progress->init();
- /** We may have several input banks instead of a single one. */
- std::vector<Iterator<Sequence>*> itBanks = itSeq->getComposition();
-
- /** We first reset the vector holding the kmers number for each partition and for each bank.
- * It can be seen as the following matrix:
- *
- * part0 part1 part2 ... partJ
- * bank0 xxx xxx xxx xxx
- * bank1 xxx xxx xxx xxx
- * ...
- * bankI xxx xxx xxx xxx
- *
- * Here xxx is the number of items found for the bank I in the partition J
- */
- _nbKmersPerPartitionPerBank.clear();
-
- /** We launch the iteration of the sequences iterator with the created functors. */
- for (size_t i=0; i<itBanks.size(); i++)
- {
- size_t groupSize = 1000;
+ if (_config._solidityKind == KMER_SOLIDITY_SUM) {
+ /** We launch the iteration of the sequences iterator with the created
+ * functors. */
+
+ size_t groupSize = 1000;
bool deleteSynchro = true;
-
- /** We fill the partitions. Each thread will read synchronously and will call FillPartitions
- * in a synchronous way (in order to have global BanksStats correctly computed). */
-
- if(_config._solidityKind == KMER_SOLIDITY_SUM)
- {
- getDispatcher()->iterate (itBanks[i], FillPartitions<span,true> (
- model, _config._nb_passes, pass, _config._nb_partitions, _config._nb_cached_items_per_core_per_part, _progress, _bankStats, _tmpPartitions, *_repartitor, pInfo,_superKstorage
- ), groupSize, deleteSynchro);
- }
- else
- {
- getDispatcher()->iterate (itBanks[i], FillPartitions<span,false> (
- model, _config._nb_passes, pass, _config._nb_partitions, _config._nb_cached_items_per_core_per_part, _progress, _bankStats, _tmpPartitions, *_repartitor, pInfo,_superKstorage
- ), groupSize, deleteSynchro);
- }
-
-
- /** We flush the partitions in order to be sure to have the exact number of items per partition. */
- if(_config._solidityKind != KMER_SOLIDITY_SUM)
- {
- _tmpPartitions->flush();
-
- /** We get a snapshot of items number in each partition. */
- vector<size_t> nbItems;
- for (size_t p=0; p<_config._nb_partitions; p++)
- {
- nbItems.push_back ((*_tmpPartitions)[p].getNbItems()); //todo for multi count
- }
-
- /** We add the current number of kmers in each partition for the reached ith bank. */
- _nbKmersPerPartitionPerBank.push_back (nbItems);
- }
-
- //GR: close the input bank here with call to finalize
- itBanks[i]->finalize();
- }
-
- if(_config._solidityKind == KMER_SOLIDITY_SUM)
- {
+
+ /** We fill the partitions. Each thread will read synchronously and will
+ * call FillPartitions in a synchronous way (in order to have global
+ * BanksStats correctly computed). */
+ getDispatcher()->iterate(
+ itSeq,
+ FillPartitions<span, true>(
+ model, _config._nb_passes, pass, _config._nb_partitions,
+ _config._nb_cached_items_per_core_per_part, _progress, _bankStats,
+ _tmpPartitions, *_repartitor, pInfo, _superKstorage),
+ groupSize, deleteSynchro);
+
+ // GR: close the input bank here with call to finalize
+ itSeq->finalize();
+
_superKstorage->flushFiles();
_superKstorage->closeFiles();
+ } else {
+ /** We may have several input banks instead of a single one. */
+ std::vector<Iterator<Sequence>*> itBanks = itSeq->getComposition();
+
+ /** We first reset the vector holding the kmers number for each partition and for each bank.
+ * It can be seen as the following matrix:
+ *
+ * part0 part1 part2 ... partJ
+ * bank0 xxx xxx xxx xxx
+ * bank1 xxx xxx xxx xxx
+ * ...
+ * bankI xxx xxx xxx xxx
+ *
+ * Here xxx is the number of items found for the bank I in the partition J
+ */
+
+ /** We launch the iteration of the sequences iterator with the created functors. */
+ for (size_t i=0; i<itBanks.size(); i++)
+ {
+ size_t groupSize = 1000;
+ bool deleteSynchro = true;
+
+ /** We fill the partitions. Each thread will read synchronously and will call FillPartitions
+ * in a synchronous way (in order to have global BanksStats correctly computed). */
+
+ getDispatcher()->iterate(
+ itBanks[i],
+ FillPartitions<span, false>(
+ model, _config._nb_passes, pass, _config._nb_partitions,
+ _config._nb_cached_items_per_core_per_part, _progress, _bankStats,
+ _tmpPartitions, *_repartitor, pInfo, _superKstorage),
+ groupSize, deleteSynchro);
+
+ /** We flush the partitions in order to be sure to have the exact number of items per partition. */
+ _tmpPartitions->flush();
+
+ /** We get a snapshot of items number in each partition. */
+ vector<size_t> nbItems;
+ for (size_t p=0; p<_config._nb_partitions; p++)
+ {
+ nbItems.push_back ((*_tmpPartitions)[p].getNbItems()); //todo for multi count
+ }
+
+ /** We add the current number of kmers in each partition for the reached ith bank. */
+ _nbKmersPerPartitionPerBank.push_back (nbItems);
+
+ //GR: close the input bank here with call to finalize
+ itBanks[i]->finalize();
+ }
}
-
-
}
/*********************************************************************
@@ -1424,8 +1432,8 @@ void SortingCountAlgorithm<span>::fillSolidKmers_aux (ICountProcessor<span>* pro
/** If we have several input banks, we may have to compute kmer solidity for each bank, which
* can be currently done only with sorted vector. */
- bool forceVector = _nbKmersPerPartitionPerBank.size() > 1 && \
- ( _config._solidityKind != KMER_SOLIDITY_SUM);
+ bool forceVector = _nbKmersPerPartitionPerBank.size() > 1 && ( _config._solidityKind != KMER_SOLIDITY_SUM);
+
ICommand* cmd = 0;
=====================================
gatb-core/src/gatb/template/TemplateSpecialization10.cpp.in
=====================================
@@ -31,9 +31,9 @@ template void bglue<${KSIZE}>(Storage* storage,
template class graph3<${KSIZE}>; // graph3<span> switch
template void link_tigs<${KSIZE}>
- (std::string unitigs_filename, int kmerSize, int nb_threads, uint64_t &nb_unitigs, bool verbose);
+ (std::string unitigs_filename, int kmerSize, int nb_threads, uint64_t &nb_unitigs, bool verbose, bool renumber_unitigs = false);
-template void link_unitigs_pass<${KSIZE}>(const std::string unitigs_filename, bool verbose, const int pass, const int kmerSize);
+template void link_unitigs_pass<${KSIZE}>(const std::string unitigs_filename, bool verbose, const int pass, const int kmerSize, const bool renumber_unitigs);
/********************************************************************************/
=====================================
gatb-core/src/gatb/tools/collections/impl/MapMPHF.hpp
=====================================
@@ -141,8 +141,8 @@ namespace gatb {
- _abundanceDiscretization[256] = UINT_MAX;
-
+ _abundanceDiscretization[256] = total;
+
=====================================
gatb-core/src/gatb/tools/misc/impl/Histogram.hpp
=====================================
@@ -85,6 +85,7 @@ public:
{
FREE (_histogram);
FREE (_histogram_smoothed);
+ FREE (_histogram2D);
}
/** \copydoc IHistogram::inc */
=====================================
gatb-core/src/gatb/tools/misc/impl/Pool.hpp
=====================================
@@ -347,11 +347,9 @@ public:
//clear all previous allocs, and alloc pool capacity
void reserve(u_int64_t size)
{
- if(size ==0 && mainbuffer !=NULL)
+ if(mainbuffer !=NULL)
{
FREE (mainbuffer);
- capacity = used_space = 0;
- mainbuffer = NULL ;
}
/** We add a little bit of memory in case "align" method is called often.
=====================================
gatb-core/test/jenkins/test-suite-fedora20-gcc-4.8.sh
=====================================
@@ -37,7 +37,7 @@ JOB_NAME : ${JOB_NAME}
# Make sure, we use the appropriate cmake
-export PATH=/home/ci-gatb/cmake-3.7.2-Linux-x86_64/bin:$PATH
+export PATH=/home/ci-gatb/cmake-3.10.0-Linux-x86_64/bin:$PATH
error_code () { [ "$DO_NOT_STOP_AT_ERROR" = "true" ] && { return 0 ; } }
View it on GitLab: https://salsa.debian.org/med-team/gatb-core/commit/8febe06b31d34f7554d7dc689e1c4ecfa95d3b17
--
View it on GitLab: https://salsa.debian.org/med-team/gatb-core/commit/8febe06b31d34f7554d7dc689e1c4ecfa95d3b17
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20190820/70b14e62/attachment-0001.html>
More information about the debian-med-commit
mailing list