Merge pull request #77 from bacpop/reads_in_hdf5

Add reads attribute to databases

Author: johnlees

CMakeLists.txt

@@ -45,8 +45,6 @@ if(DEFINED ENV{CONDA_PREFIX})
     include_directories($ENV{CONDA_PREFIX}/include)
     link_directories($ENV{CONDA_PREFIX}/lib)
     link_directories($ENV{CONDA_PREFIX}/lib/intel64)
-else()
-    find_package(OpenMP)
 endif()
 
 # Add libraries
@@ -148,6 +146,7 @@ if(DEFINED ENV{CONDA_PREFIX} AND (NOT APPLE OR CMAKE_COMPILER_IS_GNUCC OR ENV{SK
     target_link_libraries("${TARGET_NAME}" PRIVATE gomp z)
 else()
     target_link_libraries("${TARGET_NAME}" PRIVATE ZLIB::ZLIB)
+    find_package(OpenMP)
     if(OpenMP_CXX_FOUND)
         target_link_libraries("${TARGET_NAME}" PRIVATE OpenMP::OpenMP_CXX)
     endif()

src/database/database.cpp

@@ -68,6 +68,8 @@ void Database::add_sketch(const Reference &ref)
   length_a.write(ref.seq_length());
   HighFive::Attribute missing_a = sketch_group.createAttribute<unsigned long int>("missing_bases", HighFive::DataSpace::From(ref.missing_bases()));
   missing_a.write(ref.missing_bases());
+  HighFive::Attribute reads_a = sketch_group.createAttribute<bool>("reads", HighFive::DataSpace::From(ref.is_reads()));
+  reads_a.write(ref.is_reads());
 
   // Write base composition and k-mer length vectors as further group attributes
   const std::vector<double> bases = ref.base_composition();
@@ -117,6 +119,7 @@ Reference Database::load_sketch(const std::string &name)
   size_t seq_size = DEFAULT_LENGTH;
   std::vector<double> bases{0.25, 0.25, 0.25, 0.25};
   unsigned long int missing_bases = 0;
+  bool reads = false;
   std::vector<std::string> attributes_keys = sketch_group.listAttributeNames();
   for (const auto &attr : attributes_keys)
   {
@@ -132,9 +135,13 @@ Reference Database::load_sketch(const std::string &name)
     {
       sketch_group.getAttribute("missing_bases").read(missing_bases);
     }
+    else if (attr == "reads")
+    {
+      sketch_group.getAttribute("reads").read(reads);
+    }
   }
 
-  Reference new_ref(name, bbits, sketchsize64, seq_size, bases, missing_bases);
+  Reference new_ref(name, bbits, sketchsize64, seq_size, bases, missing_bases, reads);
   for (auto kmer_it = kmer_lengths.cbegin(); kmer_it != kmer_lengths.cend(); kmer_it++)
   {
     std::vector<uint64_t> usigs;

src/gpu/gpu_api.cpp

@@ -138,7 +138,8 @@ std::vector<Reference> create_sketches_cuda(
           sketches[i + j] =
               Reference(names[i + j], usigs, def_bbits, sketchsize64,
                         seq_length, seq_in_batch[j]->get_composition(),
-                        seq_in_batch[j]->missing_bases(), use_rc, densified);
+                        seq_in_batch[j]->missing_bases(), use_rc, densified,
+                        seq_in_batch[j]->is_reads());
           sketch_db.add_sketch(sketches[i + j]);
           if (densified) {
             std::cerr << "NOTE: " << names[i + j] << " required densification"

src/reference.cpp

@@ -23,6 +23,7 @@ Reference::Reference()
     : _bbits(def_bbits),
       _sketchsize64(def_sketchsize64),
       _use_rc(true),
+      _reads(false),
       _seq_size(0),
       _densified(false) {}
 
@@ -38,6 +39,7 @@ Reference::Reference(const std::string &name,
       _bbits(def_bbits),
       _sketchsize64(sketchsize64),
       _use_rc(use_rc),
+      _reads(sequence.is_reads()),
       _seq_size(0),
       _densified(false)
 {
@@ -78,9 +80,11 @@ Reference::Reference(const std::string &name,
                      const size_t sketchsize64,
                      const size_t seq_size,
                      const std::vector<double> bases,
-                     const unsigned long int missing_bases)
+                     const unsigned long int missing_bases,
+                     const bool reads)
     : _name(name), _bbits(bbits), _sketchsize64(sketchsize64), _use_rc(true),
-      _seq_size(seq_size), _missing_bases(missing_bases), _densified(false)
+      _reads(reads), _seq_size(seq_size), _missing_bases(missing_bases),
+      _densified(false)
 {
   _bases.a = bases[0];
   _bases.c = bases[1];
@@ -97,10 +101,11 @@ Reference::Reference(const std::string &name,
                      const BaseComp<double> &bases,
                      const unsigned long int missing_bases,
                      const bool use_rc,
-                     const bool densified)
+                     const bool densified,
+                     const bool reads)
     : _name(name), _bbits(bbits), _sketchsize64(sketchsize64), _use_rc(use_rc),
-      _seq_size(seq_size), _missing_bases(missing_bases), _densified(densified),
-      _bases(bases), usigs(sketch) {}
+      _reads(reads), _seq_size(seq_size), _missing_bases(missing_bases),
+      _densified(densified), _bases(bases), usigs(sketch) {}
 
 double Reference::jaccard_dist(Reference &query, const int kmer_len, const double random_jaccard)
 {

src/reference.hpp

@@ -44,7 +44,8 @@ class Reference
             const size_t sketchsize64,
             const size_t seq_size,
             const std::vector<double> bases,
-            const unsigned long int missing_bases);
+            const unsigned long int missing_bases,
+            const bool reads);
 
   // Initialise from GPU sketch
   Reference(const std::string &name,
@@ -55,7 +56,8 @@ class Reference
             const BaseComp<double> &bases,
             const unsigned long int missing_bases,
             const bool use_rc,
-            const bool densified);
+            const bool densified,
+            const bool reads);
 
   const std::vector<uint64_t> &get_sketch(const int kmer_len) const;
   void add_kmer_sketch(const std::vector<uint64_t> &sketch, const int kmer_len);
@@ -74,6 +76,7 @@ class Reference
   size_t seq_length() const { return _seq_size; }
   bool densified() const { return _densified; }
   bool rc() const { return _use_rc; }
+  bool is_reads() const { return _reads; }
   std::vector<double> base_composition() const { return {_bases.a, _bases.c, _bases.g, _bases.t}; }
   unsigned long int missing_bases() const { return _missing_bases; }
 
@@ -93,6 +96,7 @@ class Reference
   size_t _bbits;
   size_t _sketchsize64;
   bool _use_rc;
+  bool _reads;
 
   // Sequence statistics
   size_t _seq_size;