@@ -111,7 +111,7 @@ int process_mismatch_bam_record_pair(
111111 fwdbase = bam_seqi (hitseq, read_pos);
112112 if (fwdbase & (unmod_int | mod_int)) {
113113 qscore_pos = bam_get_qual (bamdata)[read_pos];
114- mod_pos = fwdbase == unmod_int ? 0.0 : 1.0 ;
114+ mod_pos = fwdbase == unmod_int ? 0 : 1 ;
115115 // check if position has already been seen, and
116116 // keep the observation with the highest qscore
117117 found = false ;
@@ -164,10 +164,12 @@ int process_mismatch_bam_record_pair(
164164 int maxdist = pair_counts.nrow () - 1 , currdist = 0 ;
165165 for (i = 0 ; i < modposref.size (); i++) {
166166 for (j = i; j < modposref.size (); j++) {
167- currdist = modposref[j] - modposref[i];
168- if (currdist > maxdist)
169- break ;
170- pair_counts (currdist, 2 * modstate[i] + modstate[j])++;
167+ // use abs() to protect against situation where modposref
168+ // is not sorted
169+ currdist = std::abs (modposref[j] - modposref[i]);
170+ if (currdist <= maxdist) {
171+ pair_counts (currdist, 2 * modstate[i] + modstate[j])++;
172+ }
171173 }
172174 }
173175
@@ -337,10 +339,10 @@ int process_mismatch_bam_record(
337339 int maxdist = pair_counts.nrow () - 1 , currdist = 0 ;
338340 for (i = 0 ; i < modposref.size (); i++) {
339341 for (j = i; j < modposref.size (); j++) {
340- currdist = modposref[j] - modposref[i];
341- if (currdist > maxdist)
342- break ;
343- pair_counts (currdist, 2 * modstate[i] + modstate[j])++;
342+ currdist = std::abs ( modposref[j] - modposref[i]) ;
343+ if (currdist <= maxdist) {
344+ pair_counts (currdist, 2 * modstate[i] + modstate[j])++ ;
345+ }
344346 }
345347 }
346348 }
@@ -494,8 +496,7 @@ Rcpp::List read_mismatchbam_cpp(std::string inname_str,
494496 std::pair<std::map<std::string,bam1_t *>::iterator, bool > inserted;
495497 hts_idx_t *idx = NULL ;
496498 hts_itr_t *iter = NULL ;
497- unsigned int regcnt = 0 , alncnt = 0 ;
498- char **regions_c = NULL ;
499+ unsigned int alncnt = 0 ;
499500 int buffer_len = 2000 ;
500501 char buffer[2000 ];
501502 const char * inname = inname_str.c_str ();
@@ -568,7 +569,7 @@ Rcpp::List read_mismatchbam_cpp(std::string inname_str,
568569 // Mode 4: random-sampling-based counting of pairs of bases by distance and modification state
569570 // -------------------------------------------------------------------------------------------
570571 success = create_multi_region_iterator_for_sampling (
571- regcnt, regions_c, n_alns_to_sample, tnames_for_sampling,
572+ n_alns_to_sample, tnames_for_sampling,
572573 keep_aln_fraction, iter, idx, in_samhdr, had_error,
573574 buffer_len, buffer);
574575
@@ -580,9 +581,8 @@ Rcpp::List read_mismatchbam_cpp(std::string inname_str,
580581 } else {
581582 // Mode 3: region-based counting of pairs of bases by distance and modification state
582583 // ----------------------------------------------------------------------------------
583- success = create_multi_region_iterator (regions, regcnt, regions_c,
584- iter, idx, in_samhdr, had_error,
585- buffer_len, buffer);
584+ success = create_multi_region_iterator (regions, iter, idx, in_samhdr,
585+ had_error, buffer_len, buffer);
586586 }
587587 if (success != 0 ) {
588588 goto end;
@@ -591,8 +591,7 @@ Rcpp::List read_mismatchbam_cpp(std::string inname_str,
591591 // iterate over regions
592592 if (verbose) {
593593 snprintf (buffer, buffer_len,
594- " counting state-pairs for alignments overlapping {%u} region{?s}" ,
595- regcnt);
594+ " counting state-pairs for alignments" );
596595 cli_alert_info (buffer);
597596 bar = cli_progress_bar (n_alns_to_sample > 0 ? n_alns_to_sample : NA_REAL,
598597 Rcpp::List::create (Rcpp::_[" clear" ] = false ,
@@ -710,15 +709,14 @@ Rcpp::List read_mismatchbam_cpp(std::string inname_str,
710709 }
711710 }
712711 // process remaining (unpaired) records in curr_records
713- for (curr_records_it = curr_records.begin ();
714- curr_records_it != curr_records.end ();
715- curr_records_it++) {
716- bamdata = curr_records_it->second ;
712+ while (!curr_records.empty ()) {
713+ curr_records_it = curr_records.begin ();
714+ bamdata2 = curr_records_it->second ;
717715
718716 if ((n_alns_to_sample == 0 ) || (R::runif (0 , 1 ) < keep_aln_fraction)) {
719717 success = process_mismatch_bam_record (
720718 MISMATCHBAM_MODE_STATE, // run mode
721- bamdata, // bam record
719+ bamdata2, // bam record
722720 bam_format, // format of bam file
723721 alncnt, // alignment counter
724722 had_error, // error flag
@@ -753,10 +751,11 @@ Rcpp::List read_mismatchbam_cpp(std::string inname_str,
753751 }
754752
755753 // remove now processed record from map
756- if (bamdata ) {
757- bam_destroy1 (bamdata );
758- bamdata = NULL ;
754+ if (bamdata2 ) {
755+ bam_destroy1 (bamdata2 );
756+ bamdata2 = NULL ;
759757 }
758+ curr_records.erase (curr_records_it);
760759
761760 if (verbose && CLI_SHOULD_TICK) { // # nocov start
762761 cli_progress_set (bar, (double )alncnt);
@@ -774,7 +773,7 @@ Rcpp::List read_mismatchbam_cpp(std::string inname_str,
774773 // Mode 2: random-sampling-based alignment reading
775774 // ---------------------------------------------------------------------
776775 success = create_multi_region_iterator_for_sampling (
777- regcnt, regions_c, n_alns_to_sample, tnames_for_sampling,
776+ n_alns_to_sample, tnames_for_sampling,
778777 keep_aln_fraction, iter, idx, in_samhdr, had_error,
779778 buffer_len, buffer);
780779 if (verbose) {
@@ -784,9 +783,8 @@ Rcpp::List read_mismatchbam_cpp(std::string inname_str,
784783 } else {
785784 // Mode 1: region-based alignment reading
786785 // ---------------------------------------------------------------------
787- success = create_multi_region_iterator (regions, regcnt, regions_c,
788- iter, idx, in_samhdr, had_error,
789- buffer_len, buffer);
786+ success = create_multi_region_iterator (regions, iter, idx, in_samhdr,
787+ had_error, buffer_len, buffer);
790788 }
791789
792790 if (success != 0 ) {
@@ -796,8 +794,7 @@ Rcpp::List read_mismatchbam_cpp(std::string inname_str,
796794 // iterate over regions
797795 if (verbose) {
798796 snprintf (buffer, buffer_len,
799- " reading alignments overlapping {%u} region{?s}" ,
800- regcnt);
797+ " reading alignments" );
801798 cli_alert_info (buffer);
802799 bar = cli_progress_bar (n_alns_to_sample > 0 ? n_alns_to_sample : NA_REAL,
803800 Rcpp::List::create (Rcpp::_[" clear" ] = false ,
@@ -866,20 +863,12 @@ Rcpp::List read_mismatchbam_cpp(std::string inname_str,
866863 }
867864 if (verbose) {
868865 cli_progress_done (bar);
869- snprintf (buffer, buffer_len,
870- " removed %llu unaligned (e.g. soft-masked) of %llu called bases" ,
871- n_unaligned, n_total);
872- cli_alert_info (buffer);
873866 snprintf (buffer, buffer_len, " read %u alignments" , alncnt);
874867 cli_alert_info (buffer);
875868 }
876869
877870 end:
878871 // cleanup
879- if (regions_c) {
880- free ((void *) regions_c);
881- regions_c = NULL ;
882- }
883872 if (in_samhdr) {
884873 sam_hdr_destroy (in_samhdr);
885874 }
@@ -895,9 +884,16 @@ Rcpp::List read_mismatchbam_cpp(std::string inname_str,
895884 if (idx) {
896885 hts_idx_destroy (idx);
897886 }
898- if (curr_records.size () > 0 ) {
887+ for (curr_records_it = curr_records.begin ();
888+ curr_records_it != curr_records.end ();
889+ curr_records_it++) { // # nocov start
890+ if (curr_records_it->second ) {
891+ bam_destroy1 (curr_records_it->second );
892+ }
893+ } // # nocov end
894+ if (curr_records.size () > 0 ) { // # nocov start
899895 curr_records.clear ();
900- }
896+ } // # nocov end
901897
902898 if (had_error) {
903899 // we encountered an error (message in `buffer`) --> stop
@@ -907,7 +903,7 @@ Rcpp::List read_mismatchbam_cpp(std::string inname_str,
907903 Rcpp::List res;
908904
909905 if (windowSize > 0 ) {
910- // Mode 3
906+ // Mode 3 or 4
911907 // create return list
912908 res = Rcpp::List::create (
913909 Rcpp::_[" pair_counts" ] = pair_counts
0 commit comments