Gather scaffolds (#154)

* Gather short scaffolds into a single bin

* Count all minimisers by default

Author: eaasna

include/valik/shared.hpp

@@ -93,8 +93,8 @@ struct build_arguments final : public split_arguments
     bool manual_parameters{false};
     bool input_is_minimiser{false};
 
-    uint8_t kmer_count_min_cutoff{2};
-    uint8_t kmer_count_max_cutoff{64};
+    uint8_t kmer_count_min_cutoff{0};
+    uint8_t kmer_count_max_cutoff{254};
     bool use_filesize_dependent_cutoff{false};
 
     std::filesystem::path ref_meta_path{};

include/valik/split/metadata.hpp

@@ -145,6 +145,12 @@ SEQAN3_WORKAROUND_GCC_BOGUS_MEMCPY_STOP
             seq_vec = std::move(ind_vec);
         }
 
+        void add_sequence(sequence_stats const & seq)
+        {
+            seq_vec.push_back(seq.ind);
+            len += seq.len;
+        }
+
         std::string unique_id()
         {
             std::string str_id{};
@@ -175,8 +181,6 @@ SEQAN3_WORKAROUND_GCC_BOGUS_MEMCPY_STOP
 
         inline bool operator() (segment_stats const & left, segment_stats const & right)
         {
-            if (left.seq_vec.size() > 1 || right.seq_vec.size() > 1)
-                throw std::runtime_error("Can't order sets of sets of sequences.");
             return (left.seq_vec[0] < right.seq_vec[0]);
         }
     };
@@ -265,15 +269,35 @@ SEQAN3_WORKAROUND_GCC_BOGUS_MEMCPY_STOP
          * @param overlap Length of overlap between adjacent segments.
          * @param seq_it Iterator to first sequence of sufficient length.
          */
-        template <typename it_t>
-        void make_exactly_n_segments(size_t const & n, size_t const & overlap, it_t & seq_it)
+        void make_exactly_n_segments(uint32_t const & n, size_t const & overlap)
         {
-            size_t remaining_db_len = total_len;
-            for (auto it = seq_it; it != sequences.end(); it++)
+            uint64_t remaining_db_len = total_len;
+            size_t len_lower_bound = default_seg_len / 10;
+            auto first_long_seq = std::find_if(sequences.begin(), sequences.end(), [&](auto const & seq){return (seq.len > len_lower_bound);});
+            
+            segment_stats short_sequences{};
+            short_sequences.start = 0;
+            short_sequences.id = segments.size();
+
+            uint64_t short_len{0};
+            for (auto it = sequences.begin(); it < sequences.end(); it++)
             {
                 auto seq = *it;
                 size_t start = 0;
-                if (seq.len <= default_seg_len * 1.5)
+                if (seq.len < len_lower_bound)
+                {
+                    short_len += seq.len;
+                    short_sequences.add_sequence(seq);
+                    if ((short_len >= default_seg_len) || 
+                        ((it + 1) == first_long_seq))
+                    {
+                        segments.push_back(short_sequences);
+                        short_sequences.seq_vec.clear();
+                        remaining_db_len -= short_len;
+                        short_len = 0;
+                    }
+                }
+                else if (seq.len <= default_seg_len * 1.5)
                 {
                     // database sequence is single segment
                     add_segment(seq.ind, start, seq.len);
@@ -319,10 +343,9 @@ SEQAN3_WORKAROUND_GCC_BOGUS_MEMCPY_STOP
          * @param overlap Length of overlap between adjacent segments.
          * @param seq_it Iterator to first sequence of sufficient length.
          */
-        template <typename it_t>
-        void make_equal_length_segments(size_t const & overlap, it_t & seq_it)
+        void make_equal_length_segments(size_t const & overlap)
         {
-            for (auto it = seq_it; it != sequences.end(); it++)
+            for (auto it = sequences.begin(); it != sequences.end(); it++)
             {
                 auto seq = *it;
                 size_t start = 0;
@@ -371,28 +394,16 @@ SEQAN3_WORKAROUND_GCC_BOGUS_MEMCPY_STOP
                                          std::to_string(arguments.pattern_size) + "bp.\nDecrease the overlap or the number of segments.");
             }
 
-            size_t len_lower_bound = default_seg_len / 10;
-            // Check how many sequences are discarded for being too short
-            //!TODO: gather short sequences
-            auto first_long_seq = std::find_if(sequences.begin(), sequences.end(), [&](auto const & seq){return (seq.len > len_lower_bound);});
-            size_t discarded_short_sequences = first_long_seq - sequences.begin();
-
-            for (auto it = sequences.begin(); it < first_long_seq; it++)
-            {
-                seqan3::debug_stream << "Sequence: " << (*it).id << " is too short and will be skipped.\n";
-                total_len -= (*it).len;
-            }
-
-            if (arguments.seg_count < (sequences.size() - discarded_short_sequences))
+            if (arguments.seg_count < sequences.size())
             {
                 throw std::runtime_error("Can not split " + std::to_string(sequences.size()) +
                                          " sequences into " + std::to_string(arguments.seg_count) + " segments.");
             }
 
             if constexpr (std::is_same<arg_t, build_arguments>::value)
-                make_exactly_n_segments(arguments.seg_count, arguments.pattern_size, first_long_seq);
+                make_exactly_n_segments(arguments.seg_count, arguments.pattern_size);
             else
-                make_equal_length_segments(arguments.pattern_size, first_long_seq);
+                make_equal_length_segments(arguments.pattern_size);
 
             std::stable_sort(sequences.begin(), sequences.end(), fasta_order());
             std::stable_sort(segments.begin(), segments.end(), fasta_order());
@@ -413,8 +424,7 @@ SEQAN3_WORKAROUND_GCC_BOGUS_MEMCPY_STOP
                                          std::to_string(arguments.pattern_size) + "bp.\nDecrease the overlap or the number of segments.");
             }
 
-            auto seq_begin = sequences.begin();
-            make_equal_length_segments(arguments.pattern_size, seq_begin);
+            make_equal_length_segments(arguments.pattern_size);
             seg_count = segments.size();
 
             std::stable_sort(sequences.begin(), sequences.end(), fasta_order());
@@ -561,7 +571,7 @@ SEQAN3_WORKAROUND_GCC_BOGUS_MEMCPY_STOP
                 segment_stats seg = segments[seg_id];
                 out_str << seg_id << '\t';
                 for (size_t ind : seg.seq_vec) 
-                    out_str << ind << '\t';
+                    out_str << ind << ';';
                 out_str << seg.start << '\t' << seg.len << '\n';
             }
 

test/api/valik/split/metadata_test.cpp

@@ -94,7 +94,7 @@ TEST_F(split_options, split_various_length)
 
     try
     {
-        for (uint32_t b : std::vector<uint32_t>{4, 16})
+        for (uint32_t b : std::vector<uint32_t>{5, 16})
         {
             arguments.seg_count = b;
             auto expected_segments = valik::metadata(segment_metadata_path(arguments.pattern_size, arguments.seg_count));
@@ -119,18 +119,23 @@ TEST_F(split_options, split_few_long)
     arguments.manual_parameters = true;
     arguments.pattern_size = 20;
     arguments.bin_path.emplace_back(data("various_chromosome_lengths.fasta"));
+    const uint8_t record_count{5};
 
     try
     {
-        for (uint32_t b : std::vector<uint32_t>{4, 12, 19})
+        for (uint32_t b : std::vector<uint32_t>{5, 12, 19})
         {
             arguments.seg_count = b;
             valik::metadata meta(arguments);
-
-            if (arguments.seg_count > meta.seq_count) // one-to-many pairing of sequences and segments
+            std::unordered_set<size_t> sequenes_ids{};
+            for (auto & seg : meta.segments)
             {
-                EXPECT_GE(0.2f, meta.segment_length_cv());  // create segments of roughly equal length
+                for (auto & seq_id : seg.seq_vec)
+                {
+                    sequenes_ids.insert(seq_id);
+                }
             }
+            EXPECT_EQ(sequenes_ids.size(), record_count); // all sequences should be used
         }
     }
     catch( const std::runtime_error& e )

test/data/split/api_test_output.sh

@@ -7,7 +7,7 @@ set -Eeuo pipefail
 seg_input="various_chromosome_lengths.fasta"
 for p in 20
 do
-    for b in 4 16
+    for b in 5 16
     do
         echo "Splitting the genome into $b segments that overlap by $p"
         seg_meta="multi/"$p"overlap"$b"bins.index"