v2.6.2

* Fix for a rare segfault when freeing cram idx

CMakeLists.txt

@@ -6,7 +6,7 @@ include(ExternalProject)
 # The version number
 set (graphtyper_VERSION_MAJOR 2)
 set (graphtyper_VERSION_MINOR 6)
-set (graphtyper_VERSION_PATCH 1)
+set (graphtyper_VERSION_PATCH 2)
 set(STATIC_DIR "" CACHE STRING "If set, GraphTyper will be built as a static binary using libraries from the given STATIC_DIR.")
 
 # Get the current working branch

include/graphtyper/typer/bucket.hpp

@@ -1,7 +1,10 @@
 #pragma once
 
+#include <set>
+
 #include <parallel_hashmap/phmap_fwd_decl.h>
 
+#include <graphtyper/graph/snp_event.hpp>
 #include <graphtyper/typer/event.hpp>
 #include <graphtyper/typer/read.hpp>
 
@@ -19,6 +22,17 @@ class BucketFirstPass
   //Tindel_events indel_events; // type is std::map<Event, EventInfo>
 };
 
+class SnpEvent;
+
+class BucketLR
+{
+public:
+  std::vector<BaseCount> pileup;
+};
+
+
+void merge_bucket_lr(std::vector<BucketLR> & old_bucket, std::vector<BucketLR> const & new_bucket);
+
 
 class Bucket
 {
@@ -56,12 +70,12 @@ add_indel_event_to_bucket(std::vector<TBucket> & buckets,
                           std::vector<char> const & reference_sequence,
                           long ref_offset);
 
-//void
-//add_base_to_bucket(std::vector<Bucket> & buckets,
-//                   int32_t pos,
-//                   char seq,
-//                   char qual,
-//                   long const region_begin,
-//                   long const BUCKET_SIZE);
+void
+add_base_to_bucket(std::vector<BucketLR> & buckets,
+                   int32_t pos,
+                   char seq,
+                   char qual,
+                   long const region_begin,
+                   long const BUCKET_SIZE);
 
 } // namespace gyper

include/graphtyper/typer/caller.hpp

@@ -44,4 +44,10 @@ streamlined_discovery(std::vector<std::string> const & hts_paths,
                       std::string const & region_str,
                       gyper::Vcf & vcf);
 
+void
+streamlined_lr_genotyping(std::vector<std::string> const & hts_paths,
+                          std::string const & reference_fn,
+                          std::string const & region_str,
+                          gyper::Vcf & vcf);
+
 } // namespace gyper

include/graphtyper/typer/event.hpp

@@ -1,5 +1,6 @@
 #pragma once
 
+#include <array>
 #include <cstdint>
 #include <string>
 #include <sstream>
@@ -8,9 +9,6 @@
 #include <vector>
 
 
-#include <graphtyper/graph/snp_event.hpp>
-
-
 namespace gyper
 {
 
@@ -20,6 +18,22 @@ long
 base2index(char const base);
 
 
+struct BaseCount
+{
+  std::array<int32_t, 4> acgt;
+  std::array<int64_t, 4> acgt_qualsum;
+  int32_t deleted{0};
+  int32_t unknown{0};
+
+  long get_depth_without_deleted() const;
+  long get_depth_with_deleted() const;
+  std::string to_string() const;
+
+  void add_base(char seq, char qual);
+
+};
+
+
 class Event
 {
 public:

include/graphtyper/utilities/genotype_lr.hpp

@@ -0,0 +1,14 @@
+#include <string>
+
+
+namespace gyper
+{
+
+void
+genotype_lr_regions(std::string ref_path,
+                    std::vector<std::string> const & sams,
+                    std::vector<GenomicRegion> const & regions,
+                    std::string const & output_path,
+                    bool const is_copy_reference);
+
+} // namespace gyper

include/graphtyper/utilities/hts_parallel_reader.hpp

@@ -22,7 +22,6 @@ class VariantMap;
 class HtsParallelReader
 {
 private:
-  std::vector<HtsReader> hts_files;
   std::vector<HtsRecord> heap;
   HtsStore store;
   std::vector<std::string> samples; // List of sample names
@@ -36,6 +35,9 @@ class HtsParallelReader
   HtsParallelReader & operator=(HtsParallelReader &&) = delete;
   ~HtsParallelReader();
 
+  // all hts files the parallel reader has
+  std::vector<HtsReader> hts_files;
+
   // Opens multiple hts paths
   void open(std::vector<std::string> const & hts_file_paths,
             std::string const & reference = "",
@@ -59,6 +61,9 @@ class HtsParallelReader
   // get the number of read groups
   long get_num_rg() const;
 
+  // get the number of samples
+  long get_num_samples() const;
+
   // get pointer to a header for the opened hts_files
   bam_hdr_t * get_header() const;
 

include/graphtyper/utilities/hts_reader.hpp

@@ -24,6 +24,7 @@ class HtsReader
   std::vector<std::string> samples;
   int sample_index_offset{0};
   int rg_index_offset{0};
+  std::vector<int> rg2sample_i; // uses the rg index to determine the sample index
 
 private:
   int ret{0}; // return value of the most recently read record
@@ -32,7 +33,6 @@ class HtsReader
   // TODO check if it is faster to first check if the records are sorted
   HtsStore & store;
   std::unordered_map<std::string, long> rg2index; // associates read groups with indices of that rg
-  std::vector<int> rg2sample_i; // uses the rg index to determine the sample index
 
   void set_reference(std::string const & reference_path);
 

include/graphtyper/utilities/options.hpp

@@ -42,6 +42,7 @@ class Options
   bool filter_on_strand_bias{true};
   bool no_filter_on_begin_pos{false};
   bool no_filter_on_coverage{false};
+  int lr_mapq_filter{5};
 
   /*******************
    * GENERAL OPTIONS *

src/CMakeLists.txt

@@ -46,6 +46,7 @@ set(graphtyper_SOURCE_FILES
   utilities/bamshrink.cpp
   utilities/genotype.cpp
   utilities/genotype_camou.cpp
+  utilities/genotype_lr.cpp
   utilities/genotype_sv.cpp
   utilities/hash_seqan.cpp
   utilities/hts_parallel_reader.cpp

src/graph/constructor.cpp

@@ -1878,9 +1878,8 @@ construct_graph(std::string const & reference_filename,
   if (!is_sorted)
   {
     std::sort(var_records.begin(), var_records.end());
-
-    BOOST_LOG_TRIVIAL(info) << __HERE__ << " Input VCF records are not sorted. "
-                                           "This is normal if the input VCF was not generated by graphtyper.";
+    //BOOST_LOG_TRIVIAL(info) << __HERE__ << " Input VCF records are not sorted. "
+    //                                       "This is normal if the input VCF was not generated by graphtyper.";
   }
 
   graph.add_genomic_region(std::move(reference_sequence),

src/graph/reference_depth.cpp

@@ -3,6 +3,8 @@
 #include <sstream> // std::stringstream
 #include <vector> // std::vector
 
+#include <boost/log/trivial.hpp>
+
 #include <graphtyper/graph/graph.hpp>
 #include <graphtyper/graph/reference_depth.hpp>
 #include <graphtyper/typer/genotype_paths.hpp>
@@ -126,6 +128,16 @@ ReferenceDepth::get_total_read_depth_of_samples(VariantCandidate const & var,
 void
 ReferenceDepth::add_genotype_paths(GenotypePaths const & geno, long const sample_index)
 {
+  assert(depths.size() > 0);
+
+  if (sample_index >= static_cast<long>(depths.size()))
+  {
+    BOOST_LOG_TRIVIAL(error) << __HERE__ << " Odd.. sample_index >= expected ("
+                             << sample_index << " >= " << depths.size() << ")";
+    return;
+    //std::exit(1);
+  }
+
   assert(sample_index < static_cast<long>(depths.size()));
 
   if (geno.paths.size() == 0 || geno.paths[0].size() < 63)

src/graph/sv.cpp

@@ -297,8 +297,7 @@ reformat_sv_vcf_records(std::vector<Variant> & variants, ReferenceDepth const &
               combined_call.filter = 0; // PASS
             }
             else if (max_gq > 40 &&
-                     (var1.calls[i].phred[index] + var2.calls[i].phred[index]) <=
-                     20)
+                     (var1.calls[i].phred[index] + var2.calls[i].phred[index]) <= 20)
             {
               combined_call.filter = 0; // PASS
             }
@@ -338,7 +337,10 @@ reformat_sv_vcf_records(std::vector<Variant> & variants, ReferenceDepth const &
 
         var.infos["SVTYPE"] = sv.get_type();
 
-        if (sv.end != 0)
+        // Make sure we do not write an END which is smaller then the begin position - otherwise bcftools is not happy
+        if (sv.end < sv.begin)
+          var.infos["END"] = std::to_string(sv.begin);
+        else
           var.infos["END"] = std::to_string(sv.end);
 
         if (sv.length != 0)