Version 2.5 (#34)

* Added option for inputting custom BAM/CRAM index files in the genotype subcommand.

* Added advanved option in the 'genotype' called '--normal_and_no_variant_overlapping' which will output both of those things in the VCF output.

* Made copying tabix index a warning instead of error. Tweaked how merging of vcf files worked. Fixed a rare bug in graph construction.

* Improved performance of creating bgzipped VCF.

* Minor bugfix in failed regression test

* Removed haplotype impurty filter as it seems to be causing some issues in large cohorts

Author: hannespetur

CMakeLists.txt

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

include/graphtyper/graph/graph.hpp

@@ -172,10 +172,9 @@ class Graph
   /**********************
    * GRAPH MODIFICATION *
    **********************/
-  void add_reference(unsigned end_pos,
+  bool add_reference(unsigned end_pos,
                      unsigned const & num_var,
-                     std::vector<char> const & reference_sequence
-                     );
+                     std::vector<char> const & reference_sequence);
 
   void add_variants(VarRecord && record);
 

include/graphtyper/graph/haplotype.hpp

@@ -101,7 +101,7 @@ class Haplotype
   /*********************
    * CLASS INFORMATION *
    *********************/
-  std::vector<uint16_t> get_haplotype_calls(std::vector<double> & haplotype_impurity) const;
+  std::vector<uint16_t> get_haplotype_calls() const;
   uint32_t get_genotype_num() const;
   bool has_too_many_genotypes() const;
   std::vector<uint32_t> get_genotype_ids() const;

include/graphtyper/graph/haplotype_calls.hpp

@@ -1,4 +1,5 @@
 #pragma once
+
 #include <iostream>
 #include <vector>
 
@@ -7,6 +8,7 @@
 #include <graphtyper/graph/haplotype.hpp>
 #include <graphtyper/typer/vcf_writer.hpp>
 
+
 namespace gyper
 {
 
@@ -18,7 +20,6 @@ class HaplotypeCall
   std::vector<uint16_t> calls;
   std::vector<Genotype> gts;
   std::vector<std::vector<ReadStrand> > read_strand;
-  std::vector<double> haplotype_impurity;
   long num_samples{0};
 
   HaplotypeCall() = default;
@@ -58,6 +59,6 @@ class HaplotypeCalls
 
 
 void save_calls(HaplotypeCalls & calls, std::string const & filename);
-std::vector<HaplotypeCall> load_calls(std::string filename);
+std::vector<HaplotypeCall> load_calls(std::string const & filename);
 
 }

include/graphtyper/typer/sample_call.hpp

@@ -3,6 +3,9 @@
 #include <cstdint> // uint8_t, uint16_t
 #include <vector> // std::vector
 
+#include <boost/serialization/access.hpp>
+
+
 namespace gyper
 {
 
@@ -12,6 +15,8 @@ class ReferenceDepth;
 
 class SampleCall
 {
+  friend class boost::serialization::access;
+
 public:
   SampleCall() noexcept;
   SampleCall(std::vector<uint8_t> && phred,
@@ -36,6 +41,9 @@ class SampleCall
   uint8_t get_gq() const;
   int8_t check_filter(long gq) const;
 
+private:
+  template <class Archive>
+  void serialize(Archive & ar, unsigned int version);
 
 };
 

include/graphtyper/typer/variant.hpp

@@ -4,6 +4,8 @@
 #include <string> // std::string
 #include <vector> // std::vector
 
+#include <boost/serialization/access.hpp>
+
 #include <graphtyper/graph/genotype.hpp> // gyper::Genotype
 #include <graphtyper/graph/haplotype.hpp> // gyper::Haplotype
 #include <graphtyper/typer/sample_call.hpp> // gyper::SampleCall
@@ -16,6 +18,8 @@ class VariantCandidate;
 
 class Variant
 {
+  friend class boost::serialization::access;
+
 public:
   uint32_t abs_pos;
   std::vector<std::vector<char> > seqs;
@@ -71,6 +75,10 @@ class Variant
   bool operator==(Variant const & b) const;
   bool operator!=(Variant const & b) const;
   bool operator<(Variant const & b) const;
+
+private:
+  template <class Archive>
+  void serialize(Archive & ar, unsigned int version);
 };
 
 // Hash function for Variant

include/graphtyper/typer/vcf.hpp

@@ -5,6 +5,8 @@
 #include <unordered_map> // std::unordered_map
 #include <vector> // std::vector
 
+#include <boost/serialization/access.hpp>
+
 #include <graphtyper/constants.hpp>
 #include <graphtyper/graph/genotype.hpp>
 #include <graphtyper/graph/haplotype.hpp>
@@ -32,6 +34,8 @@ enum VCF_FILE_MODE
 
 class Vcf
 {
+  friend class boost::serialization::access;
+
 public:
   Vcf();
   explicit Vcf(VCF_FILE_MODE _filemode, std::string const & filename);
@@ -52,22 +56,21 @@ class Vcf
   void read_samples();
   bool read_record(bool SITES_ONLY = false);
   void read(bool SITES_ONLY = false); /** \brief Reads the VCF file. */
-  void open_for_writing();
+  void open_for_writing(long const n_threads = 1);
   void write_header();
 
   void write_record(Variant const & var,
                     std::string const & suffix = "",
-                    const bool FILTER_ZERO_QUAL = false
-                    );
+                    const bool FILTER_ZERO_QUAL = false);
 
   void write_tbi_index() const;
   void write_segments();
-  void write(std::string const & region = "."); /** \brief Writes the VCF file. */
+  void write(std::string const & region = ".", long const n_threads = 1); /** \brief Writes the VCF file. */
+
   void write_records(uint32_t region_begin,
                      uint32_t region_end,
                      bool FILTER_ZERO_QUAL,
-                     std::vector<Variant> const & vars
-                     );
+                     std::vector<Variant> const & vars);
 
   void write_records(std::string const & region = ".", bool FILTER_ZERO_QUAL = false);
   void close_vcf_file();
@@ -83,23 +86,24 @@ class Vcf
 
   void add_haplotypes_for_extraction(std::vector<HaplotypeCall> const & hap_calls, bool const is_splitting_vars);
 
-  // Modify data
-  //void post_process_camou_variants(long const ploidy);
-
-  //void post_process_variants(bool NORMALIZE = true,
-  //                           bool TRIM_SEQUENCES = true
-  //                           );
-
   VCF_FILE_MODE filemode;
   std::string filename;
   std::vector<std::string> sample_names;
   std::vector<Variant> variants;
   std::vector<Segment> segments;
+
+private:
+  template <class Archive>
+  void serialize(Archive & ar, unsigned int version);
 };
 
 std::vector<std::size_t> get_all_pos(std::string const & line, char const delim = '\t');
 std::string get_string_at_tab_index(std::string const & line,
                                     std::vector<std::size_t> const & tabs,
                                     int index);
 
+void save_vcf(Vcf const & vcf, std::string const & filename);
+void load_vcf(Vcf & vcf, std::string const & filename, long n_batch);
+bool append_vcf(Vcf & vcf, std::string const & filename, long n_batch);
+
 } // namespace gyper

include/graphtyper/typer/vcf_operations.hpp

@@ -17,10 +17,11 @@ void vcf_concatenate(std::vector<std::string> const & vcfs,
 
 void vcf_break_down(std::string const & vcf, std::string const & output, std::string const & region);
 
-void vcf_merge_and_break(std::vector<std::string> & vcfs,
+void vcf_merge_and_break(std::vector<std::string> const & vcfs,
                          std::string const & output,
                          std::string const & region,
-                         bool const FILTER_ZERO_QUAL);
+                         bool const FILTER_ZERO_QUAL,
+                         bool const force_no_variant_overlap);
 
 void vcf_update_info(std::string const & vcf, std::string const & output);
 

include/graphtyper/utilities/bamshrink.hpp

@@ -33,20 +33,21 @@ namespace gyper
 {
 
 void
-bamshrink(std::string const & chrom,
+bamshrink(std::string const chrom,
           int begin,
           int end,
-          std::string const & path_in,
-          std::string const & path_out,
+          std::string const path_in,
+          std::string const sam_index_in,
+          std::string const path_out,
           double const avg_cov_by_readlen,
-          std::string const & ref_fn);
+          std::string const ref_fn);
 
 
 void
-bamshrink_multi(std::string const & interval_fn,
-                std::string const & path_in,
-                std::string const & path_out,
+bamshrink_multi(std::string const interval_fn,
+                std::string const path_in,
+                std::string const path_out,
                 double const avg_cov_by_readlen,
-                std::string const & ref_fn);
+                std::string const ref_fn);
 
 } // namespace gyper

include/graphtyper/utilities/bgzf_stream.hpp

@@ -19,45 +19,46 @@ class BGZF_stream
 {
 private:
   BGZF * fp = nullptr;
-  std::ostringstream ss;
-
 
 public:
+  std::ostringstream ss;
+
   BGZF_stream() = default;
-  BGZF_stream(std::string const & filename, std::string const & filemode);
-  ~BGZF_stream();
+  ~BGZF_stream(); // custom
+  BGZF_stream(BGZF_stream const &) = delete;
+  BGZF_stream(BGZF_stream &&) = delete;
+  BGZF_stream & operator=(BGZF_stream const &) = delete;
+  BGZF_stream & operator=(BGZF_stream &&) = delete;
 
   template <class T>
   BGZF_stream & operator<<(T const & x);
+  void check_cache();
   void flush();
-  void open(std::string const & filename, std::string const & filemode);
+  void open(std::string const & filename, std::string const & filemode, long const n_threads);
   bool is_open() const;
   void close(); // Close BGZF file
 
-  long MAX_CACHE_SIZE = 10000000ll;
+  long MAX_CACHE_SIZE{1000000ll};
 };
 
 
-inline
-BGZF_stream::BGZF_stream(std::string const & filename, std::string const & filemode)
-{
-  open(filename, filemode);
-}
-
-
 template <class T>
 inline
 BGZF_stream &
 BGZF_stream::operator<<(T const & x)
 {
   // Add to stringstream
   ss << x;
+  return *this;
+}
+
 
-  // Check if we should flush
+inline
+void
+BGZF_stream::check_cache()
+{
   if (static_cast<long>(ss.tellp()) > this->MAX_CACHE_SIZE)
     flush();
-
-  return *this;
 }
 
 
@@ -72,7 +73,7 @@ BGZF_stream::flush()
   }
   else
   {
-    std::string const & str = ss.str();
+    std::string str = ss.str();
     int ret = bgzf_write(fp, str.data(), str.size());
 
     if (ret < 0)
@@ -89,14 +90,17 @@ BGZF_stream::flush()
 
 inline
 void
-BGZF_stream::open(std::string const & filename, std::string const & filemode)
+BGZF_stream::open(std::string const & filename, std::string const & filemode, long const n_threads)
 {
   if (fp)
     close();
 
   if (filename.size() > 0 && filename != "-")
   {
     fp = bgzf_open(filename.c_str(), filemode.c_str());
+
+    if (n_threads > 1)
+      bgzf_mt(fp, n_threads, 256);
   }
 }
 

include/graphtyper/utilities/genotype.hpp

@@ -17,24 +17,27 @@ namespace gyper
 
 std::vector<std::string>
 run_bamshrink(std::vector<std::string> const & sams,
+              std::vector<std::string> const & sams_index,
               std::string const & ref_fn,
               GenomicRegion const & region,
               std::vector<double> const & avg_cov_by_readlen,
               std::string const & tmp);
 
+// bamshrink variant for multiple regions
 std::vector<std::string>
-run_bamshrink(std::vector<std::string> const & sams,
-              std::string const & ref_fn,
-              std::string const & interval_fn,
-              std::vector<double> const & avg_cov_by_readlen,
-              std::string const & tmp);
+run_bamshrink_multi(std::vector<std::string> const & sams,
+                    std::string const & ref_fn,
+                    std::string const & interval_fn,
+                    std::vector<double> const & avg_cov_by_readlen,
+                    std::string const & tmp);
 
 void
 run_samtools_merge(std::vector<std::string> & shrinked_sams, std::string const & tmp);
 
 void
 genotype(std::string ref_path,
          std::vector<std::string> const & sams,
+         std::vector<std::string> const & sams_index,
          GenomicRegion const & region,
          std::string const & output_path,
          std::vector<double> const & avg_cov_by_readlen,
@@ -44,6 +47,7 @@ genotype(std::string ref_path,
 void
 genotype_regions(std::string const & ref_path,
                  std::vector<std::string> const & sams,
+                 std::vector<std::string> const & sams_index,
                  std::vector<GenomicRegion> const & regions,
                  std::string const & output_path,
                  std::vector<double> const & avg_cov_by_readlen,

include/graphtyper/utilities/options.hpp

@@ -29,6 +29,7 @@ class Options
   bool no_decompose{false};
   bool no_bamshrink{false};
   bool no_variant_overlapping{false};
+  bool normal_and_no_variant_overlapping{false};
   bool is_all_biallelic{false};
   bool is_only_cigar_discovery{false};
   bool is_discovery_only_for_paired_reads{false};
@@ -94,6 +95,8 @@ class Options
   long genotype_dis_min_support{8};
   double genotype_dis_min_support_ratio{0.30};
 
+  // internal vcf options
+  long num_alleles_in_batch{500};
 
   /********************************
    * HAPLOTYPE EXTRACTION OPTIONS *