From 251d39f057bc0a03892a0f39d7f6d89aef91070a Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sat, 19 Oct 2024 12:52:04 -0700 Subject: [PATCH] src/load_data_for_complexity.cpp: fixing a bug in checking whether a read is mapped -- the condition was inverted --- src/load_data_for_complexity.cpp | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/src/load_data_for_complexity.cpp b/src/load_data_for_complexity.cpp index d8aa7a0..ed44001 100644 --- a/src/load_data_for_complexity.cpp +++ b/src/load_data_for_complexity.cpp @@ -520,6 +520,11 @@ load_coverage_counts_GR(const string &input_file_name, const uint64_t seed, #ifdef HAVE_HTSLIB // Deal with SAM/BAM format only if we have htslib +static inline void +is_unmapped(const bamxx::bam_rec &aln) { + return get_tid(aln) == -1; +} + static inline void swap(bamxx::bam_rec &a, bamxx::bam_rec &b) { std::swap(a.b, b.b); @@ -593,12 +598,12 @@ load_counts_BAM(const uint32_t n_threads, const string &inputfile, // find first mapped read to start bamxx::bam_rec aln; - while (hts.read(hdr, aln) && get_tid(aln) == -1) + while (hts.read(hdr, aln) && is_unmapped(aln)) ; size_t n_reads{}; // if all reads unmapped, must return - if (get_tid(aln) == -1) + if (is_unmapped(aln)) return n_reads; // to check that reads are sorted properly @@ -611,7 +616,7 @@ load_counts_BAM(const uint32_t n_threads, const string &inputfile, size_t current_count = 1; while (hts.read(hdr, aln)) { - if (get_tid(aln) != -1) + if (is_unmapped(aln)) continue; // skip unmapped reads const aln_pos_t curr{aln}; @@ -739,11 +744,11 @@ load_coverage_counts_BAM(const uint32_t n_threads, const string &inputfile, // find first mapped read to start bamxx::bam_rec aln; - while (hts.read(hdr, aln) && get_tid(aln) == -1) + while (hts.read(hdr, aln) && is_unmapped(aln)) ; size_t n_reads{}; - if (get_tid(aln) == -1) // no reads unmapped + if (is_unmapped(aln)) // no reads unmapped return 0; // to check reads are sorted properly @@ -761,7 +766,7 @@ load_coverage_counts_BAM(const uint32_t n_threads, const string &inputfile, const hts_pos_t max_dist = bin_size + max_width; while (hts.read(hdr, aln)) { - if (get_tid(aln) == -1) + if (is_unmapped(aln)) continue; // check that read is mapped const hts_pos_t rlen = rlen_from_cigar(aln);