Allowing unmapped reads in counts, bsrate, format and uniq

src/analysis/bsrate.cpp

@@ -580,12 +580,14 @@ main_bsrate(int argc, const char **argv) {
     unordered_set<int32_t> chroms_seen;
 
     while (hts.read(hdr, aln)) {
+      const int32_t the_tid = get_tid(aln);
+      if (the_tid == -1)
+        continue;
+
       if (reads_are_a_rich) flip_conversion(aln);
 
       // get the correct chrom if it has changed
-      if (get_tid(aln) != current_tid) {
-        const int32_t the_tid = get_tid(aln);
-
+      if (the_tid != current_tid) {
         // make sure all reads from same chrom are contiguous in the file
         if (chroms_seen.find(the_tid) != end(chroms_seen))
           throw runtime_error("chroms out of order in mapped reads file");

src/analysis/methcounts.cpp

@@ -446,22 +446,30 @@ process_reads(const bool VERBOSE, const bool show_progress,
   vector<string>::const_iterator chrom_itr;
 
   while (hts.read(hdr, aln)) {
-    if (get_tid(aln) == prev_tid) {
+    const int32_t tid = get_tid(aln);
+    if (tid == -1)  // ADS: skip reads that have no tid -- they are not mapped
+      continue;
+    if (tid == prev_tid) {
       // do the work for this mapped read, depending on strand
       if (bam_is_rev(aln))
         count_states_neg(aln, counts);
       else
         count_states_pos(aln, counts);
     }
-    else {  // chrom changes, output results, get the next one
+    else {  // chrom has changed, so output results and get the next chrom
 
-      // write output if there is any; should fail only once
+      // write output if there is any; counts is empty only for first chrom
       if (!counts.empty())
         write_output<require_covered>(hdr, out, prev_tid,
                                       *chrom_itr, counts, CPG_ONLY);
-
-      const int32_t tid = get_tid(aln);
-      if (tid < prev_tid) throw dnmt_error("SAM file is not sorted");
+      // make sure reads are sorted chrom tid number in header
+      if (tid < prev_tid) {
+        const std::string message = "SAM file is not sorted "
+          "previous tid: " +
+          std::to_string(prev_tid) +
+          " current tid: " + std::to_string(tid);
+        throw dnmt_error(message);
+      }
 
       if (!require_covered)
         for (auto i = prev_tid + 1; i < tid; ++i)

src/utils/format-reads.cpp

@@ -341,6 +341,10 @@ format(const string &cmd, const size_t n_threads, const string &inputfile,
     swap(aln, prev_aln);  // start with prev_aln being first read
 
     while (hts.read(hdr, aln)) {
+      // ADS: skip reads that have no tid -- they are not mapped
+      if (get_tid(aln) == -1)
+        continue;
+
       standardize_format(input_format, aln);
       if (same_name(prev_aln, aln, suff_len)) {
         // below: essentially check for dovetail

src/utils/uniq.cpp

@@ -227,8 +227,14 @@ uniq(const bool add_dup_count, const uint32_t max_buffer_size,
   }
 
   bam_rec aln;
-  if (hts.read(hdr, aln)) {  // valid SAM/BAM can have 0 reads
+  bool found_mapped_read{false};  // valid SAM/BAM can have 0 reads
+  while (!found_mapped_read && hts.read(hdr, aln)) {
+    // ADS: skip reads that have no tid -- they are not mapped
+    if (get_tid(aln) != -1)
+      found_mapped_read = true;
+  }
 
+  if (found_mapped_read) {
     rs_in.update(aln);  // update stats for input we just got
 
     vector<bam_rec> buffer(1, aln);  // select output from this buffer
@@ -238,6 +244,9 @@ uniq(const bool add_dup_count, const uint32_t max_buffer_size,
     int32_t cur_chrom = get_tid(aln);
 
     while (hts.read(hdr, aln)) {
+      // ADS: skip reads that have no tid -- they are not mapped
+      if (get_tid(aln) == -1)
+        continue;
       rs_in.update(aln);
 
       // below works because buffer reset at every new chrom