fixing merge conflict

Author: guilhermesena1

Makefile.am

@@ -25,14 +25,17 @@ lib_LIBRARIES = libsmithlab_cpp.a
 libsmithlab_cpp_a_SOURCES = GenomicRegion.cpp MappedRead.cpp		\
 OptionParser.cpp QualityScore.cpp bisulfite_utils.cpp			\
 chromosome_utils.cpp sim_utils.cpp smithlab_os.cpp smithlab_utils.cpp	\
-zlib_wrapper.cpp dna_four_bit.cpp cigar_utils.cpp
+zlib_wrapper.cpp dna_four_bit.cpp cigar_utils.cpp sam_record.cpp
 
 if ENABLE_HTS
-libsmithlab_cpp_a_SOURCES += htslib_wrapper_deprecated.cpp sam_record.cpp
+libsmithlab_cpp_a_SOURCES += htslib_wrapper_deprecated.cpp htslib_wrapper.cpp
 endif
 
 include_HEADERS = GenomicRegion.hpp MappedRead.hpp OptionParser.hpp	\
 QualityScore.hpp bisulfite_utils.hpp chromosome_utils.hpp		\
 sim_utils.hpp smithlab_os.hpp smithlab_utils.hpp zlib_wrapper.hpp	\
-dna_four_bit.hpp cigar_utils.hpp htslib_wrapper_deprecated.hpp \
-sam_record.hpp
+dna_four_bit.hpp cigar_utils.hpp sam_record.hpp
+
+if ENABLE_HTS
+include_HEADERS += htslib_wrapper.hpp htslib_wrapper_deprecated.hpp
+endif

cigar_utils.hpp

@@ -64,8 +64,7 @@ cigar_qseq_ops(InputItr itr, const InputItr last) {
   size_t op_count = 0;
   while (itr != last) {
     const size_t curr_val = extract_op_count(itr, last);
-    if (consumes_query(*itr++))
-      op_count += curr_val;
+    op_count += curr_val*consumes_query(*itr++);
   }
   return op_count;
 }
@@ -82,8 +81,7 @@ cigar_rseq_ops(InputItr itr, const InputItr last) {
   size_t op_count = 0;
   while (itr != last) {
     const size_t curr_val = extract_op_count(itr, last);
-    if (consumes_reference(*itr++))
-      op_count += curr_val;
+    op_count += curr_val*consumes_reference(*itr++);
   }
   return op_count;
 }

dna_four_bit.hpp

@@ -27,40 +27,24 @@ extern char dna_four_bit_decoding[16];
 
 template <typename uint_type> constexpr
 uint_type
-get_low_nibble(const uint_type x) {return x & 15u;}
-
-template <typename uint_type> constexpr
-uint_type
-get_high_nibble(const uint_type x) {return (x >> 4) & 15u;}
-
-template <typename uint_type> constexpr
-char
-decode_dna_four_bit_low(const uint_type x) {
-  return dna_four_bit_decoding[get_low_nibble(x)];
-}
-
-template <typename uint_type> constexpr
-char
-decode_dna_four_bit_high(const uint_type x) {
-  return dna_four_bit_decoding[get_high_nibble(x)];
+get_nibble(const uint_type x, const size_t offset) {
+  return (x >> (4*offset)) & 15ul;
 }
 
 template <typename uint_type> constexpr
 char
 decode_dna_four_bit(const uint_type x,
-                    const base_in_byte b = base_in_byte::left) {
-  return b == base_in_byte::left ?
-    decode_dna_four_bit_low(x) :
-    decode_dna_four_bit_high(x);
+                    const size_t offset) {
+  return dna_four_bit_decoding[get_nibble(x, offset)];
 }
 
 template<class InputItr, class OutputIt>
 OutputIt
 decode_dna_four_bit(InputItr first, InputItr last, OutputIt d_first) {
   // ADS: assume destination has enough space
   while (first != last) {
-    *d_first++ = decode_dna_four_bit(*first, base_in_byte::left);
-    *d_first++ = decode_dna_four_bit(*first, base_in_byte::right);
+    for (size_t offset = 0; offset < 16; ++offset)
+      *d_first++ = decode_dna_four_bit(*first, offset);
     ++first;
   }
   // if original sequence length is odd and encoding not padded at the front,
@@ -73,95 +57,88 @@ void
 decode_dna_four_bit(const InCtr &source, OutCtr &dest) {
   // expand out the bytes as pairs (do this backwards in case source == dest)
   const size_t source_size = source.size();
-  dest.resize(2*source_size);
+  dest.resize(16*source_size);
   size_t i = source_size;
   size_t j = dest.size();
   while (i > 0) {
     dest[--j] = source[--i];
     dest[--j] = source[i];
   }
-  for (i = 0; i < dest.size(); i += 2) {
-    dest[i] = decode_dna_four_bit(dest[i], base_in_byte::left);
-    dest[i+1] = decode_dna_four_bit(dest[i+1], base_in_byte::right);
+  for (i = 0; i < dest.size(); i += 16) {
+    for (size_t offset = 0; offset < 16; ++offset)
+      dest[i + offset] = decode_dna_four_bit(dest[i], offset);
   }
 }
 
 extern uint8_t dna_four_bit_encoding[128];
-
-template <typename uint_type> constexpr
-uint8_t
-encode_dna_four_bit_low(const uint_type x) {
-  return dna_four_bit_encoding[static_cast<unsigned>(x)];
-}
-
 template <typename uint_type> constexpr
-uint8_t
-encode_dna_four_bit_high(const uint_type x) {
-  return dna_four_bit_encoding[static_cast<unsigned>(x)] << 4;
-}
-
-template <typename uint_type> constexpr
-uint8_t
+size_t
 encode_dna_four_bit(const uint_type x,
-                    const base_in_byte b = base_in_byte::left) {
-  return b == base_in_byte::left ?
-    encode_dna_four_bit_low(x) :
-    encode_dna_four_bit_high(x);
+                    const size_t offset) {
+  return (static_cast<size_t>(
+           dna_four_bit_encoding[static_cast<unsigned>(x)])
+         ) << (4*offset);
 }
 
 template<class InputItr, class OutputIt>
 OutputIt
 encode_dna_four_bit(InputItr first, InputItr last, OutputIt d_first) {
   while (first != last) {
-    *d_first  = encode_dna_four_bit(*first++, base_in_byte::left);
-    *d_first |= (first == last ? 0 :
-                 encode_dna_four_bit(*first++, base_in_byte::right));
+    *d_first = 0;
+    for (size_t i = 0; i < 16 && first != last; ++i)
+      *d_first |= encode_dna_four_bit(*first++, i);
     ++d_first;
   }
   return d_first;
 }
 
-// ADS: indented to be used as pointer to 4-bit encoding of DNA within a vector
-// of uint8_t values
+// GS: intended to be used as pointer to 4-bit encoding of DNA within a vector
+// of size_t values
 struct genome_four_bit_itr {
-  genome_four_bit_itr(const std::vector<uint8_t>::const_iterator itr_,
-                      const bool odd_ = false) : itr(itr_), itr_odd(odd_) {}
+  genome_four_bit_itr(const std::vector<size_t>::const_iterator itr_,
+                      const int off_ = 0) : itr(itr_), offset(off_) {}
 
-  uint8_t operator*() const {
-    return (!itr_odd ? *itr : (*itr >> 4)) & 15;
+  size_t operator*() const {
+    return (*itr >> (offset << 2)) & 15ul;
   }
   genome_four_bit_itr& operator++() {
-    itr += itr_odd;
-    itr_odd ^= 1ul;
+    offset = (offset + 1) & 15ul;
+    itr += (offset == 0);
     return *this;
   }
   genome_four_bit_itr operator++(int) {
     genome_four_bit_itr tmp(*this);
-    itr += itr_odd;
-    itr_odd ^= 1ul;
+    offset = (offset + 1) & 15ul;
+    itr += (offset == 0);
     return tmp;
   }
   genome_four_bit_itr& operator--() {
-    itr -= !itr_odd;
-    itr_odd ^= 1ul;
+    itr -= (offset == 0);
+
+    // GS: will underflow on 0 but it's ok?
+    offset = (offset - 1) & 15ul;
     return *this;
   }
   genome_four_bit_itr operator--(int) {
     genome_four_bit_itr tmp(*this);
-    itr -= !itr_odd;
-    itr_odd ^= 1ul;
+    itr -= (offset == 0);
+    offset = (offset - 1) & 15ul;
     return tmp;
   }
-  genome_four_bit_itr operator+(const size_t offset) const {
-    const size_t offset_odd = offset & 1ul;
-    return genome_four_bit_itr(itr + offset/2 + (itr_odd & offset_odd),
-                               itr_odd != offset_odd);
+  genome_four_bit_itr operator+(const size_t step) const {
+    // whether the sum of offsets is >= 16 
+    const bool shift_one_pos =
+      (((offset + (static_cast<int>(step) & 15)) & 16) >> 4);
+
+    const int new_offset = (offset + step) & 15;
+    return genome_four_bit_itr(itr + step/16 + shift_one_pos,
+                               new_offset);
   }
   bool operator!=(const genome_four_bit_itr &rhs) const {
-    return itr != rhs.itr || itr_odd != rhs.itr_odd;
+    return itr != rhs.itr || offset != rhs.offset;
   }
-  std::vector<uint8_t>::const_iterator itr;
-  size_t itr_odd;
+  std::vector<size_t>::const_iterator itr;
+  int offset;
 };
 
 #endif

sam_record.cpp

@@ -34,6 +34,7 @@ using std::ostream;
 using std::begin;
 using std::end;
 using std::ostringstream;
+using std::to_string;
 
 // ADS: this is for debugging purposes
 string
@@ -73,49 +74,44 @@ format_sam_flags(const uint16_t the_flags) {
 //   return regex_match(qual, regex("[!-~]+"));
 // }
 
+size_t
+sam_rec::estimate_line_size() const {
+  static const size_t all_field_estimates = 100;
+  return qname.size() + rname.size() + qual.size() + all_field_estimates;
+}
+
 string
 sam_rec::tostring() const {
-  ostringstream oss;
-  oss << qname << '\t'
-      << flags << '\t'
-      << rname << '\t'
-      << pos << '\t'
-      << static_cast<unsigned>(mapq) << '\t'
-      << cigar << '\t'
-      << rnext << '\t'
-      << pnext << '\t'
-      << tlen << '\t'
-      << seq << '\t'
-      << qual;
-
+  string out;
+  out.reserve(estimate_line_size());
+  out.append(qname + "\t" +
+               to_string(flags) + "\t" +
+               rname + "\t" +
+               to_string(pos) + "\t" +
+               to_string(static_cast<unsigned>(mapq)) + "\t" +
+               cigar + "\t" +
+               rnext + "\t" +
+               to_string(pnext) + "\t" +
+               to_string(tlen) + "\t" +
+               seq + "\t" +
+               qual);
   for (auto it(begin(tags)); it != end(tags); ++it)
-    oss << '\t' << *it;
-  return oss.str() + "\n";
+    out.append("\t" + *it);
+
+  return out;
 }
 
 ostream &
 operator<<(std::ostream &the_stream, const sam_rec &r) {
-  the_stream << r.qname << '\t'
-             << r.flags << '\t'
-             << r.rname << '\t'
-             << r.pos << '\t'
-             << static_cast<unsigned>(r.mapq) << '\t'
-             << r.cigar << '\t'
-             << r.rnext << '\t'
-             << r.pnext << '\t'
-             << r.tlen << '\t'
-             << r.seq << '\t'
-             << r.qual;
-
-  for (auto it(begin(r.tags)); it != end(r.tags); ++it)
-    the_stream << '\t' << *it;
+  the_stream << r.tostring();
   return the_stream;
 }
 
 sam_rec::sam_rec(const string &line) {
-
+  /*
   istringstream iss; // ADS: change to set the buffer from "line"
-  iss.rdbuf()->pubsetbuf(const_cast<char*>(line.c_str()), line.size());
+  iss.rdbuf()->pubsetbuf(const_cast<char*>(line.c_str()), line.size());*/
+  istringstream iss(line);
   uint32_t will_become_mapq = 0; // to not read mapq as character
                                  // since it's uint8_t
   if (!(iss >>

sam_record.hpp

@@ -22,6 +22,8 @@
 #include <string>
 #include <vector>
 #include <iostream>
+#include <sstream>
+#include <iterator>
 
 // from 30 April 2020 SAM documentation
 // 1    0x1   template having multiple segments in sequencing
@@ -118,6 +120,7 @@ class sam_rec {
     qual(_qual) {}
   void add_tag(const std::string &the_tag) {tags.push_back(the_tag);}
   bool empty() const;
+  size_t estimate_line_size() const;
   std::string tostring() const;
 };
 
@@ -146,4 +149,36 @@ void
 inflate_with_cigar(const sam_rec &sr, std::string &to_inflate,
                    const char inflation_symbol = 'N');
 
+template<typename T>
+static void
+write_sam_header(const std::vector<std::string> &chrom_names,
+                 const std::vector<T> &chrom_starts,
+                 const std::string program_name,
+                 const std::string program_version,
+                 const int argc, const char **argv,
+                 std::ostream &out) {
+  static const std::string SAM_VERSION = "1.0";
+
+  // sam version
+  out <<"@HD" << '\t' << "VN:" << SAM_VERSION << '\n'; // sam version
+
+  // chromosome sizes
+  const size_t n_chroms = chrom_names.size() - 1;
+  for (size_t i = 1; i < n_chroms; ++i) {
+    out << "@SQ" << '\t'
+        << "SN:" << chrom_names[i] << '\t'
+        << "LN:" << chrom_starts[i+1] - chrom_starts[i] << '\n';
+  }
+
+  // program details
+  out << "@PG" << '\t'
+      << "ID:" << program_name << '\t'
+      << "VN:" << program_version << '\t';
+
+  // how the program was run
+  std::ostringstream the_command;
+  copy(argv, argv + argc, std::ostream_iterator<const char*>(the_command, " "));
+  out << "CL:\"" << the_command.str() << "\"" << std::endl;
+}
+
 #endif

smithlab_os.cpp

@@ -390,7 +390,11 @@ read_fasta_file_short_names(const string &filename,
   string line;
   while (getline(in, line)) {
     if (line[0] == '>') {
-      names.push_back(string(begin(line) + 1,
+      const auto first_space = line.find_first_of(" \t", 1);
+      if (first_space == string::npos)
+        names.push_back(line.substr(1));
+      else
+        names.push_back(string(begin(line) + 1,
                              begin(line) + line.find_first_of(" \t", 1)));
       sequences.push_back(string());
     }