Merge pull request #94 from smithlabcode/format-bsmap

Format bsmap

Author: andrewdavidsmith

src/analysis/bsrate.cpp

@@ -332,12 +332,10 @@ main_bsrate(int argc, const char **argv) {
         const int32_t the_tid = get_tid(aln);
 
         // make sure all reads from same chrom are contiguous in the file
-        if (chroms_seen.find(the_tid) != end(chroms_seen)) {
-          cerr << the_tid << '\t' << current_tid << endl;
+        if (chroms_seen.find(the_tid) != end(chroms_seen))
           throw runtime_error("chroms out of order in mapped reads file");
-        }
+
         current_tid = the_tid;
-        if (VERBOSE) cerr << "processing " << the_tid << endl;
 
         chroms_seen.insert(the_tid);
 
@@ -346,6 +344,9 @@ main_bsrate(int argc, const char **argv) {
           throw runtime_error("could not find chrom: " + the_tid);
 
         chrom_idx = chrom_itr->second;
+
+        if (VERBOSE) cerr << "processing " << names[chrom_idx] << endl;
+
         use_this_chrom = seq_to_use.empty() || chrom_idx == chrom_idx_to_use;
       }
 

src/common/bam_record_utils.cpp

@@ -527,15 +527,12 @@ merge_by_byte(const bam1_t *a, const bam1_t *b, bam1_t *c) {
 
 static inline void
 flip_conversion(bam1_t *aln) {
-  if (aln->core.flag & BAM_FREVERSE)
-    aln->core.flag = aln->core.flag & (~BAM_FREVERSE);
-  else
-    aln->core.flag = aln->core.flag | BAM_FREVERSE;
+  aln->core.flag ^= BAM_FREVERSE;  // ADS: flip the "reverse" bit
 
   revcomp_seq_by_byte(aln);
 
   // ADS: don't like *(cv + 1) below, but no HTSlib function for it?
-  uint8_t *cv = bam_aux_get(aln, "CV");
+  auto cv = bam_aux_get(aln, "CV");
   if (!cv) throw dnmt_error("bam_aux_get failed for CV");
   *(cv + 1) = 'T';
 }
@@ -812,31 +809,40 @@ standardize_format(const string &input_format, bam1_t *aln) {
   if (input_format == "abismal" || input_format == "walt") return;
 
   if (input_format == "bsmap") {
-    // A/T rich; get the ZS tag value
-    auto zs_tag = bam_aux_get(aln, "ZS");
+    // A/T rich: get the ZS tag value
+    const auto zs_tag = bam_aux_get(aln, "ZS");
     if (!zs_tag) throw dnmt_error("bam_aux_get for ZS (invalid bsmap)");
     // ADS: test for errors on the line below
-    const uint8_t cv = string(bam_aux2Z(zs_tag))[1] == '-' ? 'A' : 'T';
+    const auto zs_tag_value = string(bam_aux2Z(zs_tag));
+    if (zs_tag_value.empty()) throw dnmt_error("empty ZS tag in bsmap format");
+    if (zs_tag_value[0] != '-' && zs_tag_value[0] != '+')
+      throw dnmt_error("invalid ZS tag in bsmap format");
+    const uint8_t cv = zs_tag_value[1] == '-' ? 'A' : 'T';
     // get the "mismatches" tag
-    auto nm_tag = bam_aux_get(aln, "NM");
-    if (!nm_tag) throw dnmt_error("bam_aux_get for NM (invalid bsmap)");
+    const auto nm_tag = bam_aux_get(aln, "NM");
+    if (!nm_tag) throw dnmt_error("invalid NM tag in bsmap format");
     const int64_t nm = bam_aux2i(nm_tag);
 
-    aln->l_data = bam_get_aux(aln) - aln->data;  // del aux (no data resize)
+    // ADS: this should delete the aux data by truncating the used
+    // range within the bam1_t while avoiding resizing memory
+    aln->l_data = bam_get_aux(aln) - aln->data;
 
     /* add the tags we want */
     // "udpate" for "int" because it determines the right size; even
-    // though we just deleted all tags, it will add it back here.
+    // though we just deleted all tags, it will add it back here
     err_code = bam_aux_update_int(aln, "NM", nm);
-    if (err_code < 0) throw dnmt_error(err_code, "bam_aux_update_int");
+    if (err_code < 0)
+      throw dnmt_error(err_code, "error setting NM in bsmap format");
+
     // "append" for "char" because there is no corresponding update
     err_code = bam_aux_append(aln, "CV", 'A', 1, &cv);
-    if (err_code < 0) throw dnmt_error(err_code, "bam_aux_append");
+    if (err_code < 0)
+      throw dnmt_error(err_code, "error setting conversion in bsmap format");
 
-    if (bam_is_rev(aln))
-      revcomp_seq_by_byte(aln);  // reverse complement if needed
+    // reverse complement if needed
+    if (bam_is_rev(aln)) revcomp_seq_by_byte(aln);
   }
-  if (input_format == "bismark") {
+  else if (input_format == "bismark") {
     // ADS: Previously we modified the read names at the first
     // underscore. Even if the names are still that way, it should no
     // longer be needed since we compare names up to a learned suffix.
@@ -864,6 +870,9 @@ standardize_format(const string &input_format, bam1_t *aln) {
     if (bam_is_rev(aln))
       revcomp_seq_by_byte(aln);  // reverse complement if needed
   }
+  // ADS: the condition below should be checked much earlier, ideally
+  // before the output file is created
+  else throw runtime_error("incorrect format specified: " + input_format);
 
   // Be sure this doesn't depend on mapper! Removes the "qual" part of
   // the data in a bam1_t struct but does not change its uncompressed

src/utils/format-reads.cpp

@@ -57,9 +57,9 @@ using std::vector;
 
 using bamxx::bam_rec;
 
-static size_t
+static int32_t
 merge_mates(const size_t range, bam_rec &one, bam_rec &two, bam_rec &merged) {
-  if (!are_mates(one, two)) return -std::numeric_limits<int>::max();
+  if (!are_mates(one, two)) return -std::numeric_limits<int32_t>::max();
 
   // arithmetic easier using base 0 so subtracting 1 from pos
   const int one_s = get_pos(one);
@@ -309,7 +309,7 @@ swap(bam_rec &a, bam_rec &b) {
 static void
 format(const string &cmd, const size_t n_threads, const string &inputfile,
        const string &outfile, const bool bam_format, const string &input_format,
-       const size_t suff_len, const size_t max_frag_len) {
+       const size_t suff_len, const int32_t max_frag_len) {
   static const dnmt_error bam_write_err{"error writing bam"};
 
   bamxx::bam_tpool tp(n_threads);
@@ -345,8 +345,7 @@ format(const string &cmd, const size_t n_threads, const string &inputfile,
       if (same_name(prev_aln, aln, suff_len)) {
         // below: essentially check for dovetail
         if (!bam_is_rev(aln)) swap(prev_aln, aln);
-        const size_t frag_len =
-          merge_mates(max_frag_len, prev_aln, aln, merged);
+        const auto frag_len = merge_mates(max_frag_len, prev_aln, aln, merged);
         if (frag_len > 0 && frag_len < max_frag_len) {
           if (is_a_rich(merged)) flip_conversion(merged);
           if (!out.write(hdr, merged)) throw bam_write_err;
@@ -385,7 +384,7 @@ main_format(int argc, const char **argv) {
 
     string input_format;
     string outfile;
-    int max_frag_len = std::numeric_limits<int>::max();
+    int32_t max_frag_len = 10000;
     size_t suff_len = 0;
     bool single_end = false;
     bool VERBOSE = false;
@@ -409,7 +408,7 @@ main_format(int argc, const char **argv) {
     opt_parse.add_opt("single-end", '\0',
                       "assume single-end [do not use with -suff]", false,
                       single_end);
-    opt_parse.add_opt("max-frag", 'L', "maximum allowed insert size", false,
+    opt_parse.add_opt("max-frag", 'L', "max allowed insert size", false,
                       max_frag_len);
     opt_parse.add_opt("check", 'c',
                       "check this many reads to validate read name suffix",
@@ -443,6 +442,12 @@ main_format(int argc, const char **argv) {
            << opt_parse.about_message() << endl;
       return EXIT_FAILURE;
     }
+    if (max_frag_len <= 0) {
+      cerr << "specified maximum fragment size: " << max_frag_len << endl
+           << opt_parse.help_message() << endl
+           << opt_parse.about_message() << endl;
+      return EXIT_FAILURE;
+    }
     const string infile(leftover_args.front());
     if (leftover_args.size() == 2 && !use_stdout)
       outfile = leftover_args.back();