Merge remote-tracking branch 'gitlab/master'

Conflicts:
	paw

Author: hannespetur

CMakeLists.txt

@@ -7,7 +7,7 @@ include(CheckIncludeFileCXX)
 # The version number
 set (graphtyper_VERSION_MAJOR 2)
 set (graphtyper_VERSION_MINOR 7)
-set (graphtyper_VERSION_PATCH 5)
+set (graphtyper_VERSION_PATCH 6)
 
 # Graphtyper's headers
 include_directories("${CMAKE_CURRENT_SOURCE_DIR}/include")

include/graphtyper/graph/sv.hpp

@@ -44,8 +44,8 @@ class SV
   int32_t length = 0; // Length of alt. allele minus ref. allele
   int32_t size = 0;   // Total size, i.e. length of larger allele excluding padding base
   int32_t end = 0;
-  int16_t n_clusters = 0;
-  int16_t num_merged_svs = -1;
+  int32_t n_clusters = 0;
+  int32_t num_merged_svs = -1;
   int32_t or_start = -1; // Start coordinate of the sequence origin
   int32_t or_end = -1;   // End coordinate of the sequence origin
   int32_t related_sv = -1;

paw

@@ -1283,6 +1283,9 @@ void add_var_record(std::vector<VarRecord> & var_records,
     seqan::StringSet<seqan::CharString> infos;
     seqan::strSplit(infos, vcf_record.info, seqan::EqualsChar<';'>());
 
+    // for a special case for Dragen 3.7+ files.
+    bool is_a_dup = false;
+
     for (auto const & info : infos)
     {
       long const EQ_SIGN_POS = std::distance(seqan::begin(info), std::find(seqan::begin(info), seqan::end(info), '='));
@@ -1292,6 +1295,9 @@ void add_var_record(std::vector<VarRecord> & var_records,
                               ? std::string(seqan::begin(info) + EQ_SIGN_POS + 1, seqan::end(info))
                               : std::string("");
 
+      if (key == "DUPSVLEN")
+        is_a_dup = true;
+
       parse_info_sv_type("SVTYPE", key, val, sv.type) ||                 // SVTYPE
         parse_info_int("END", key, val, sv.end) ||                       // END
         parse_info_int("SVSIZE", key, val, sv.size) ||                   // SVSIZE
@@ -1304,6 +1310,7 @@ void add_var_record(std::vector<VarRecord> & var_records,
         parse_info_str("SVINSSEQ", key, val, sv.ins_seq) ||              // SVINSSEQ
         parse_info_str("LEFT_SVINSSEQ", key, val, sv.ins_seq_left) ||    // LEFT_SVINSSEQ
         parse_info_str("RIGHT_SVINSSEQ", key, val, sv.ins_seq_right) ||  // RIGHT_SVINSSEQ
+        parse_info_str("DUPSVINSSEQ", key, val, sv.ins_seq) ||           // DUPSVINSSEQ
         parse_info_inv_type(key, sv.inv_type);                           // INV3 and INV5
     }
 
@@ -1318,6 +1325,9 @@ void add_var_record(std::vector<VarRecord> & var_records,
       std::exit(1);
     }
 
+    if (sv.type == INS && is_a_dup)
+      sv.type = DUP; // special case for Dragen 3.7 files, there small duplications have SVTYPE=INS
+
     if (sv.length < 0)
       sv.length = -sv.length; // Make SVLEN positive
 

src/graph/constructor.cpp

@@ -890,6 +890,7 @@ int subcmd_genotype_sv(paw::Parser & parser)
   std::string sam{};
   std::string sams{};
   std::string sv_vcf{};
+  std::string encoding{"vcf"};
   bool force_copy_reference{false};
   bool force_no_copy_reference{false};
   bool see_advanced_options{false};
@@ -1006,6 +1007,12 @@ int subcmd_genotype_sv(paw::Parser & parser)
                                "(0-level) in the header. The byte range of these blocks will also be printed in "
                                "${prefix}.samples_byte_range.");
 
+  parser.parse_advanced_option(encoding, ' ', "encoding", "Select output encoding. Available are: vcf, popvcf");
+
+  // set default compression level as 9 when popvcf, since level 9 is already very fast anyway in that encoding mode
+  if (encoding == "popvcf")
+    opts.bgzf_compression_level = 9;
+
   // Changed behaviour such that zero qual SVs are not filtered out by default
   if (!force_filter_zero_qual)
     opts.force_no_filter_zero_qual = true;
@@ -1016,6 +1023,7 @@ int subcmd_genotype_sv(paw::Parser & parser)
   parser.finalize();
   setup_logger();
 
+  opts.encoding = (encoding == "popvcf") ? 'p' : 'v';
   print_log(gyper::log_severity::info, "Running the 'genotype_sv' subcommand.");
 
 #ifndef NDEBUG
@@ -1302,20 +1310,31 @@ int subcmd_vcf_concatenate(paw::Parser & parser)
 
 int subcmd_vcf_merge(paw::Parser & parser)
 {
+  gyper::Options & opts = *(gyper::Options::instance());
+
   std::vector<std::string> vcfs;
   std::string output_fn;
   std::string file_list;
+  std::string encoding{"vcf"};
   bool is_sv_vcf{false};
 
   parser.parse_option(output_fn, 'o', "output", "Output VCF file name.");
   parser.parse_option(file_list, ' ', "file_list", "File containing VCFs to merge.");
   parser.parse_option(is_sv_vcf, ' ', "sv", "Set if the input VCFs were generated from genotype_sv.");
+  parser.parse_option(encoding, ' ', "encoding", "Select output encoding. Available are: vcf, popvcf");
+
+  // set default compression level as 9 when popvcf, since level 9 is already very fast anyway in that encoding mode
+  if (encoding == "popvcf")
+    opts.bgzf_compression_level = 9;
 
   parser.parse_remaining_positional_arguments(vcfs, "vcfs...", "VCFs to merge");
 
   parser.finalize();
   setup_logger();
 
+  opts.encoding = (encoding == "popvcf") ? 'p' : 'v';
+  opts.is_on_final_output = true;
+
   if (is_sv_vcf)
     gyper::graph.is_sv_graph = true;
 

src/main.cpp

@@ -549,8 +549,14 @@ void run_first_pass(bam1_t * hts_rec,
 
     if (ref_offset >= REF_SIZE)
     {
-      print_log(log_severity::error, __HERE__, " Unexpected ref_offset = ", ref_offset);
-      std::exit(1);
+      if (ref_offset > REF_SIZE)
+      {
+        print_log(log_severity::warning, __HERE__, " Unexpected ref_offset=", ref_offset, " > REF_SIZE=", REF_SIZE);
+      }
+
+      // ref_offset==REF_SIZE occurs if the first read happens to be on the very last base. In this case is safe to just
+      // ignore it.
+      break;
     }
 
     Read read;

src/typer/caller.cpp

@@ -1632,9 +1632,9 @@ std::vector<Variant> break_down_variant(Variant && var,
   else if (!is_no_variant_overlapping)
   {
     // Use the skyr
-    print_log(log_severity::debug, "Using the skyr");
+    print_debug("Using the skyr");
     std::vector<Variant> new_broken_down_vars = break_down_skyr(std::move(var), reach);
-    print_log(log_severity::debug, "skyr finished.");
+    print_debug("skyr finished.");
 
     std::move(new_broken_down_vars.begin(), new_broken_down_vars.end(), std::back_inserter(broken_down_vars));
   }

src/typer/variant.cpp

@@ -1119,8 +1119,9 @@ void Vcf::write_record(Variant const & var,
         assert(new_line_size > old_line_size);
         double const bytes_per_call = static_cast<double>(new_line_size - old_line_size) / static_cast<double>(i + 1);
         double total_bytes_expected = old_line_size + bytes_per_call * static_cast<double>(var.calls.size());
+        double constexpr MAX_BYTES{static_cast<double>(std::numeric_limits<int32_t>::max()) * 0.9};
 
-        if (total_bytes_expected >= 2000000000.0)
+        if (total_bytes_expected >= MAX_BYTES)
         {
           print_log(log_severity::warning,
                     " Skipping variant with extreme expected line size=",

src/typer/vcf.cpp

@@ -1,7 +1,9 @@
-#include <cmath>   // sqrt
-#include <sstream> // std::ostringstream
-#include <string>  // std::string
-#include <vector>  // std::vector
+#include <cmath>      // sqrt
+#include <libgen.h>   // dirname
+#include <sstream>    // std::ostringstream
+#include <string>     // std::string
+#include <sys/stat.h> // mkdir
+#include <vector>     // std::vector
 
 #include <graphtyper/graph/absolute_position.hpp>
 #include <graphtyper/graph/genomic_region.hpp>
@@ -149,6 +151,16 @@ void vcf_merge(std::vector<std::string> & vcfs, std::string const & output)
   if (vcfs.size() == 0)
     return;
 
+  // Create output directory
+  {
+    int const l = output.size();
+    char * c = new char[l + 1];
+    std::copy(output.begin(), output.end(), c);
+    c[l] = '\0';
+    mkdir(dirname(c), 0755);
+    delete[] c;
+  }
+
   gyper::Vcf vcf;
   vcf.open(READ_MODE, vcfs.at(0));
   vcf.read();                   // Read the entire file