diff --git a/src/aln.cpp b/src/aln.cpp index 05bd4f29..099ba248 100644 --- a/src/aln.cpp +++ b/src/aln.cpp @@ -56,7 +56,7 @@ bool reverse_nam_if_needed(Nam& nam, const Read& read, const References& referen std::string ref_end_kmer = references.sequences[nam.ref_id].substr(nam.ref_end-k, k); std::string seq, seq_rc; - if (nam.is_rc) { + if (nam.is_revcomp) { seq = read.rc; seq_rc = read.seq; } else { @@ -77,7 +77,7 @@ bool reverse_nam_if_needed(Nam& nam, const Read& read, const References& referen read_start_kmer = seq_rc.substr(q_start_tmp, k); read_end_kmer = seq_rc.substr(q_end_tmp - k, k); if (ref_start_kmer == read_start_kmer && ref_end_kmer == read_end_kmer) { - nam.is_rc = !nam.is_rc; + nam.is_revcomp = !nam.is_revcomp; nam.query_start = q_start_tmp; nam.query_end = q_end_tmp; return true; @@ -220,7 +220,7 @@ inline Alignment extend_seed( const Read& read, bool consistent_nam ) { - const std::string query = nam.is_rc ? read.rc : read.seq; + const std::string query = nam.is_revcomp ? read.rc : read.seq; const std::string& ref = references.sequences[nam.ref_id]; const auto projected_ref_start = nam.projected_ref_start(); @@ -269,7 +269,7 @@ inline Alignment extend_seed( alignment.score = info.sw_score; alignment.ref_start = result_ref_start; alignment.length = info.ref_span(); - alignment.is_rc = nam.is_rc; + alignment.is_revcomp = nam.is_revcomp; alignment.is_unaligned = false; alignment.ref_id = nam.ref_id; alignment.gapped = gapped; @@ -335,7 +335,7 @@ inline std::vector get_best_scoring_pairs( for (auto &a2 : alignments2) { float dist = std::abs(a1.ref_start - a2.ref_start); double score = a1.score + a2.score; - if ((a1.is_rc ^ a2.is_rc) && (dist < mu + 4 * sigma)) { + if ((a1.is_revcomp ^ a2.is_revcomp) && (dist < mu + 4 * sigma)) { score += log(normal_pdf(dist, mu, sigma)); } else { // individual score @@ -350,17 +350,17 @@ inline std::vector get_best_scoring_pairs( } bool is_proper_nam_pair(const Nam nam1, const Nam nam2, float mu, float sigma) { - if (nam1.ref_id != nam2.ref_id || nam1.is_rc == nam2.is_rc) { + if (nam1.ref_id != nam2.ref_id || nam1.is_revcomp == nam2.is_revcomp) { return false; } int r1_ref_start = nam1.projected_ref_start(); int r2_ref_start = nam2.projected_ref_start(); // r1 ---> <---- r2 - bool r1_r2 = nam2.is_rc && (r1_ref_start <= r2_ref_start) && (r2_ref_start - r1_ref_start < mu + 10*sigma); + bool r1_r2 = nam2.is_revcomp && (r1_ref_start <= r2_ref_start) && (r2_ref_start - r1_ref_start < mu + 10*sigma); // r2 ---> <---- r1 - bool r2_r1 = nam1.is_rc && (r2_ref_start <= r1_ref_start) && (r1_ref_start - r2_ref_start < mu + 10*sigma); + bool r2_r1 = nam1.is_revcomp && (r2_ref_start <= r1_ref_start) && (r1_ref_start - r2_ref_start < mu + 10*sigma); return r1_r2 || r2_r1; } @@ -458,7 +458,7 @@ inline Alignment rescue_align( std::string r_tmp; auto read_len = read.size(); - if (mate_nam.is_rc) { + if (mate_nam.is_revcomp) { r_tmp = read.seq; a = mate_nam.projected_ref_start() - (mu+5*sigma); b = mate_nam.projected_ref_start() + read_len/2; // at most half read overlap @@ -477,7 +477,7 @@ inline Alignment rescue_align( alignment.edit_distance = read_len; alignment.score = 0; alignment.ref_start = 0; - alignment.is_rc = mate_nam.is_rc; + alignment.is_revcomp = mate_nam.is_revcomp; alignment.ref_id = mate_nam.ref_id; alignment.is_unaligned = true; // std::cerr << "RESCUE: Caught Bug3! ref start: " << ref_start << " ref end: " << ref_end << " ref len: " << ref_len << std::endl; @@ -490,7 +490,7 @@ inline Alignment rescue_align( alignment.edit_distance = read_len; alignment.score = 0; alignment.ref_start = 0; - alignment.is_rc = mate_nam.is_rc; + alignment.is_revcomp = mate_nam.is_revcomp; alignment.ref_id = mate_nam.ref_id; alignment.is_unaligned = true; return alignment; @@ -502,7 +502,7 @@ inline Alignment rescue_align( alignment.edit_distance = info.edit_distance; alignment.score = info.sw_score; alignment.ref_start = ref_start + info.ref_start; - alignment.is_rc = !mate_nam.is_rc; + alignment.is_revcomp = !mate_nam.is_revcomp; alignment.ref_id = mate_nam.ref_id; alignment.is_unaligned = info.cigar.empty(); alignment.length = info.ref_span(); @@ -859,8 +859,8 @@ std::vector align_paired( a2_indv_max = a2; } - bool r1_r2 = a2.is_rc && (a1.ref_start <= a2.ref_start) && ((a2.ref_start - a1.ref_start) < mu + 10*sigma); // r1 ---> <---- r2 - bool r2_r1 = a1.is_rc && (a2.ref_start <= a1.ref_start) && ((a1.ref_start - a2.ref_start) < mu + 10*sigma); // r2 ---> <---- r1 + bool r1_r2 = a2.is_revcomp && (a1.ref_start <= a2.ref_start) && ((a2.ref_start - a1.ref_start) < mu + 10*sigma); // r1 ---> <---- r2 + bool r2_r1 = a1.is_revcomp && (a2.ref_start <= a1.ref_start) && ((a1.ref_start - a2.ref_start) < mu + 10*sigma); // r2 ---> <---- r1 double combined_score; if (r1_r2 || r2_r1) { @@ -1019,6 +1019,61 @@ bool has_shared_substring(const std::string& read_seq, const std::string& ref_se return false; } + +/* + * Obtain NAMs for a sequence record, doing rescue if needed. + * Return NAMs sorted by decreasing score. + */ +std::vector get_nams( + const KSeq& record, + const StrobemerIndex& index, + AlignmentStatistics& statistics, + Details& details, + const MappingParameters &map_param, + const IndexParameters& index_parameters, + std::minstd_rand& random_engine +) { + // Compute randstrobes + Timer strobe_timer; + auto query_randstrobes = randstrobes_query(record.seq, index_parameters); + statistics.n_randstrobes += query_randstrobes.size(); + statistics.tot_construct_strobemers += strobe_timer.duration(); + + // Find NAMs + Timer nam_timer; + auto [nonrepetitive_fraction, n_hits, nams] = find_nams(query_randstrobes, index, map_param.use_mcs); + statistics.tot_find_nams += nam_timer.duration(); + statistics.n_hits += n_hits; + details.nams = nams.size(); + + // Rescue if requested and needed + if (map_param.rescue_level > 1 && (nams.empty() || nonrepetitive_fraction < 0.7)) { + Timer rescue_timer; + int n_rescue_hits; + std::tie(n_rescue_hits, nams) = find_nams_rescue(query_randstrobes, index, map_param.rescue_cutoff, map_param.use_mcs); + statistics.n_rescue_hits += n_rescue_hits; + details.rescue_nams = nams.size(); + details.nam_rescue = true; + statistics.tot_time_rescue += rescue_timer.duration(); + } + + // Sort by score + Timer nam_sort_timer; + std::sort(nams.begin(), nams.end(), by_score); + shuffle_top_nams(nams, random_engine); + statistics.tot_sort_nams += nam_sort_timer.duration(); + +#ifdef TRACE + std::cerr << "Query: " << record.name << '\n'; + std::cerr << "Found " << nams.size() << " NAMs\n"; + for (const auto& nam : nams) { + std::cerr << "- " << nam << '\n'; + } +#endif + + return nams; +} + void align_or_map_paired( const KSeq &record1, const KSeq &record2, @@ -1037,50 +1092,13 @@ void align_or_map_paired( std::array details; std::array, 2> nams_pair; - for (size_t is_revcomp : {0, 1}) { - const auto& record = is_revcomp == 0 ? record1 : record2; - - Timer strobe_timer; - auto query_randstrobes = randstrobes_query(record.seq, index_parameters); - statistics.n_randstrobes += query_randstrobes.size(); - statistics.tot_construct_strobemers += strobe_timer.duration(); - - // Find NAMs - Timer nam_timer; - auto [nonrepetitive_fraction, n_hits, nams] = find_nams(query_randstrobes, index, map_param.use_mcs); - statistics.tot_find_nams += nam_timer.duration(); - statistics.n_hits += n_hits; - details[is_revcomp].nams = nams.size(); - - if (map_param.rescue_level > 1) { - Timer rescue_timer; - if (nams.empty() || nonrepetitive_fraction < 0.7) { - int n_rescue_hits; - std::tie(n_rescue_hits, nams) = find_nams_rescue(query_randstrobes, index, map_param.rescue_cutoff, map_param.use_mcs); - details[is_revcomp].nam_rescue = true; - details[is_revcomp].rescue_nams = nams.size(); - statistics.n_rescue_hits += n_rescue_hits; - } - statistics.tot_time_rescue += rescue_timer.duration(); - } - Timer nam_sort_timer; - std::sort(nams.begin(), nams.end(), by_score); - shuffle_top_nams(nams, random_engine); - statistics.tot_sort_nams += nam_sort_timer.duration(); - nams_pair[is_revcomp] = nams; - } - + for (size_t is_r1 : {0, 1}) { + const auto& record = is_r1 == 0 ? record1 : record2; #ifdef TRACE - for (int is_revcomp : {0, 1}) { - const auto& record = is_revcomp == 0 ? record1 : record2; - std::cerr << "R" << is_revcomp + 1 << " name: " << record.name << '\n'; - const auto& nams = nams_pair[is_revcomp]; - std::cerr << "Found " << nams.size() << " NAMs\n"; - for (auto& nam : nams) { - std::cerr << "- " << nam << '\n'; - } - } + std::cerr << "R" << is_r1 + 1 << '\n'; #endif + nams_pair[is_r1] = get_nams(record, index, statistics, details[is_r1], map_param, index_parameters, random_engine); + } Timer extend_timer; if (map_param.output_format != OutputFormat::SAM) { // PAF or abundance @@ -1184,42 +1202,7 @@ void align_or_map_single( std::vector &abundances ) { Details details; - Timer strobe_timer; - auto query_randstrobes = randstrobes_query(record.seq, index_parameters); - statistics.n_randstrobes += query_randstrobes.size(); - statistics.tot_construct_strobemers += strobe_timer.duration(); - - // Find NAMs - Timer nam_timer; - auto [nonrepetitive_fraction, n_hits, nams] = find_nams(query_randstrobes, index, map_param.use_mcs); - statistics.tot_find_nams += nam_timer.duration(); - statistics.n_hits += n_hits; - details.nams = nams.size(); - - if (map_param.rescue_level > 1) { - Timer rescue_timer; - if (nams.empty() || nonrepetitive_fraction < 0.7) { - int n_rescue_hits; - std::tie(n_rescue_hits, nams) = find_nams_rescue(query_randstrobes, index, map_param.rescue_cutoff, map_param.use_mcs); - statistics.n_rescue_hits += n_rescue_hits; - details.rescue_nams = nams.size(); - details.nam_rescue = true; - } - statistics.tot_time_rescue += rescue_timer.duration(); - } - - Timer nam_sort_timer; - std::sort(nams.begin(), nams.end(), by_score); - shuffle_top_nams(nams, random_engine); - statistics.tot_sort_nams += nam_sort_timer.duration(); - -#ifdef TRACE - std::cerr << "Query: " << record.name << '\n'; - std::cerr << "Found " << nams.size() << " NAMs\n"; - for (auto& nam : nams) { - std::cerr << "- " << nam << '\n'; - } -#endif + std::vector nams = get_nams(record, index, statistics, details, map_param, index_parameters, random_engine); Timer extend_timer; size_t n_best = 0; diff --git a/src/nam.cpp b/src/nam.cpp index fbd81e94..e484842d 100644 --- a/src/nam.cpp +++ b/src/nam.cpp @@ -120,7 +120,7 @@ void merge_matches_into_nams( n.query_prev_match_startpos = m.query_start; n.ref_prev_match_startpos = m.ref_start; n.n_matches = 1; - n.is_rc = is_revcomp; + n.is_revcomp = is_revcomp; // n.score += (float)1 / (float)h.count; open_nams.push_back(n); } @@ -187,7 +187,7 @@ std::vector merge_matches_into_nams_forward_and_reverse( * Return the fraction of nonrepetitive hits (those not above the filter_cutoff threshold) */ std::tuple> find_nams( - const QueryRandstrobeVector &query_randstrobes, + const std::vector& query_randstrobes, const StrobemerIndex& index, bool use_mcs ) { @@ -204,7 +204,7 @@ std::tuple> find_nams( continue; } nr_good_hits++; - add_to_matches_map_full(matches_map[q.is_reverse], q.start, q.end, index, position); + add_to_matches_map_full(matches_map[q.is_revcomp], q.start, q.end, index, position); } else if (use_mcs) { size_t partial_pos = index.find_partial(q.hash); @@ -214,7 +214,7 @@ std::tuple> find_nams( continue; } nr_good_hits++; - add_to_matches_map_partial(matches_map[q.is_reverse], q.start, q.start + index.k(), index, partial_pos); + add_to_matches_map_partial(matches_map[q.is_revcomp], q.start, q.start + index.k(), index, partial_pos); } } } @@ -229,7 +229,7 @@ std::tuple> find_nams( continue; } nr_good_hits++; - add_to_matches_map_partial(matches_map[q.is_reverse], q.start, q.start + index.k(), index, partial_pos); + add_to_matches_map_partial(matches_map[q.is_revcomp], q.start, q.start + index.k(), index, partial_pos); } } } @@ -246,7 +246,7 @@ std::tuple> find_nams( * Return the number of hits and the vector of NAMs. */ std::pair> find_nams_rescue( - const QueryRandstrobeVector &query_randstrobes, + const std::vector& query_randstrobes, const StrobemerIndex& index, unsigned int rescue_cutoff, bool use_mcs @@ -275,14 +275,14 @@ std::pair> find_nams_rescue( if (position != index.end()) { unsigned int count = index.get_count_full(position); RescueHit rh{position, count, qr.start, qr.end, false}; - hits[qr.is_reverse].push_back(rh); + hits[qr.is_revcomp].push_back(rh); } else if (use_mcs) { size_t partial_pos = index.find_partial(qr.hash); if (partial_pos != index.end()) { unsigned int partial_count = index.get_count_partial(partial_pos); RescueHit rh{partial_pos, partial_count, qr.start, qr.start + index.k(), true}; - hits[qr.is_reverse].push_back(rh); + hits[qr.is_revcomp].push_back(rh); } } } @@ -312,6 +312,6 @@ std::pair> find_nams_rescue( } std::ostream& operator<<(std::ostream& os, const Nam& n) { - os << "Nam(ref_id=" << n.ref_id << ", query: " << n.query_start << ".." << n.query_end << ", ref: " << n.ref_start << ".." << n.ref_end << ", rc=" << static_cast(n.is_rc) << ", score=" << n.score << ")"; + os << "Nam(ref_id=" << n.ref_id << ", query: " << n.query_start << ".." << n.query_end << ", ref: " << n.ref_start << ".." << n.ref_end << ", rc=" << static_cast(n.is_revcomp) << ", score=" << n.score << ")"; return os; } diff --git a/src/nam.hpp b/src/nam.hpp index 0ca8c02a..a2d4b3b2 100644 --- a/src/nam.hpp +++ b/src/nam.hpp @@ -20,7 +20,7 @@ struct Nam { float score; // unsigned int previous_query_start; // unsigned int previous_ref_start; - bool is_rc = false; + bool is_revcomp = false; int ref_span() const { return ref_end - ref_start; @@ -36,13 +36,13 @@ struct Nam { }; std::tuple> find_nams( - const QueryRandstrobeVector &query_randstrobes, + const std::vector &query_randstrobes, const StrobemerIndex& index, bool use_mcs ); std::pair> find_nams_rescue( - const QueryRandstrobeVector &query_randstrobes, + const std::vector &query_randstrobes, const StrobemerIndex& index, unsigned int rescue_cutoff, bool use_mcs diff --git a/src/paf.cpp b/src/paf.cpp index ede81e25..4e9690bf 100644 --- a/src/paf.cpp +++ b/src/paf.cpp @@ -26,7 +26,7 @@ void output_hits_paf_PE(std::string &paf_output, const Nam &n, const std::string paf_output.append("\t"); paf_output.append(std::to_string(n.query_end)); paf_output.append("\t"); - paf_output.append(n.is_rc ? "-" : "+"); + paf_output.append(n.is_revcomp ? "-" : "+"); paf_output.append("\t"); paf_output.append(references.names[n.ref_id]); paf_output.append("\t"); diff --git a/src/python/strobealign.cpp b/src/python/strobealign.cpp index 3a2d9b26..729cb248 100644 --- a/src/python/strobealign.cpp +++ b/src/python/strobealign.cpp @@ -135,12 +135,12 @@ NB_MODULE(strobealign_extension, m_) { .def_ro("hash", &QueryRandstrobe::hash) .def_ro("start", &QueryRandstrobe::start) .def_ro("end", &QueryRandstrobe::end) - .def_ro("is_reverse", &QueryRandstrobe::is_reverse) + .def_ro("is_revcomp", &QueryRandstrobe::is_revcomp) .def("__repr__", [](const QueryRandstrobe& qr) { std::stringstream s; s << qr; return s.str(); }) ; - nb::bind_vector(m, "QueryRandstrobeVector"); + nb::bind_vector>(m, "QueryRandstrobeVector"); nb::class_(m, "Nam") .def_ro("query_start", &Nam::query_start) @@ -150,7 +150,7 @@ NB_MODULE(strobealign_extension, m_) { .def_ro("score", &Nam::score) .def_ro("n_hits", &Nam::n_matches) .def_ro("reference_index", &Nam::ref_id) - .def_rw("is_rc", &Nam::is_rc) + .def_rw("is_revcomp", &Nam::is_revcomp) .def_prop_ro("ref_span", &Nam::ref_span) .def_prop_ro("query_span", &Nam::query_span) .def("__repr__", [](const Nam& nam) { @@ -160,7 +160,7 @@ NB_MODULE(strobealign_extension, m_) { nb::bind_vector>(m, "NamVector"); m.def("randstrobes_query", &randstrobes_query); - m.def("find_nams", [](const QueryRandstrobeVector &query_randstrobes, const StrobemerIndex& index, bool use_mcs) -> std::vector { + m.def("find_nams", [](const std::vector& query_randstrobes, const StrobemerIndex& index, bool use_mcs) -> std::vector { auto [nonrepetitive_fraction, n_hits, nams] = find_nams(query_randstrobes, index, use_mcs); return nams; }, nb::arg("query_randstrobes"), nb::arg("index"), nb::arg("use_mcs")); diff --git a/src/randstrobes.cpp b/src/randstrobes.cpp index 9170437b..4ba37a19 100644 --- a/src/randstrobes.cpp +++ b/src/randstrobes.cpp @@ -149,7 +149,7 @@ std::ostream& operator<<(std::ostream& os, const QueryRandstrobe& randstrobe) { os << "QueryRandstrobe(hash=" << randstrobe.hash << ", start=" << randstrobe.start << ", end=" << randstrobe.end - << ", is_reverse=" << randstrobe.is_reverse + << ", is_revcomp=" << randstrobe.is_revcomp << ")"; return os; @@ -223,8 +223,8 @@ Randstrobe RandstrobeGenerator::next() { /* * Generate randstrobes for a query sequence and its reverse complement. */ -QueryRandstrobeVector randstrobes_query(const std::string_view seq, const IndexParameters& parameters) { - QueryRandstrobeVector randstrobes; +std::vector randstrobes_query(const std::string_view seq, const IndexParameters& parameters) { + std::vector randstrobes; if (seq.length() < parameters.randstrobe.w_max) { return randstrobes; } diff --git a/src/randstrobes.hpp b/src/randstrobes.hpp index ae308966..af281696 100644 --- a/src/randstrobes.hpp +++ b/src/randstrobes.hpp @@ -68,14 +68,12 @@ struct QueryRandstrobe { randstrobe_hash_t hash; unsigned int start; unsigned int end; - bool is_reverse; + bool is_revcomp; }; std::ostream& operator<<(std::ostream& os, const QueryRandstrobe& randstrobe); -using QueryRandstrobeVector = std::vector; - -QueryRandstrobeVector randstrobes_query(const std::string_view seq, const IndexParameters& parameters); +std::vector randstrobes_query(const std::string_view seq, const IndexParameters& parameters); struct Randstrobe { randstrobe_hash_t hash; diff --git a/src/sam.cpp b/src/sam.cpp index 254ce812..92d75d67 100644 --- a/src/sam.cpp +++ b/src/sam.cpp @@ -138,7 +138,7 @@ void Sam::add( assert(!alignment.is_unaligned); int flags = 0; - if (!alignment.is_unaligned && alignment.is_rc) { + if (!alignment.is_unaligned && alignment.is_revcomp) { flags |= REVERSE; } if (!is_primary) { @@ -274,7 +274,7 @@ void Sam::add_pair( pos1 = -1; reference_name1 = "*"; } else { - if (alignment1.is_rc) { + if (alignment1.is_revcomp) { f1 |= REVERSE; f2 |= MREVERSE; } @@ -290,7 +290,7 @@ void Sam::add_pair( pos2 = -1; reference_name2 = "*"; } else { - if (alignment2.is_rc) { + if (alignment2.is_revcomp) { f1 |= MREVERSE; f2 |= REVERSE; } @@ -333,8 +333,8 @@ bool is_proper_pair(const Alignment& alignment1, const Alignment& alignment2, fl const int dist = alignment2.ref_start - alignment1.ref_start; const bool same_reference = alignment1.ref_id == alignment2.ref_id; const bool both_aligned = same_reference && !alignment1.is_unaligned && !alignment2.is_unaligned; - const bool r1_r2 = !alignment1.is_rc && alignment2.is_rc && dist >= 0; // r1 ---> <---- r2 - const bool r2_r1 = !alignment2.is_rc && alignment1.is_rc && dist <= 0; // r2 ---> <---- r1 + const bool r1_r2 = !alignment1.is_revcomp && alignment2.is_revcomp && dist >= 0; // r1 ---> <---- r2 + const bool r2_r1 = !alignment2.is_revcomp && alignment1.is_revcomp && dist <= 0; // r2 ---> <---- r1 const bool rel_orientation_good = r1_r2 || r2_r1; const bool insert_good = std::abs(dist) <= mu + 6 * sigma; @@ -350,7 +350,7 @@ std::ostream& operator<<(std::ostream& os, const Alignment& alignment) { << ", global_ed=" << alignment.global_ed << ", score=" << alignment.score << ", length=" << alignment.length - << ", is_rc=" << alignment.is_rc + << ", is_revcomp=" << alignment.is_revcomp << ", is_unaligned=" << alignment.is_unaligned << ", gapped=" << alignment.gapped << ")"; diff --git a/src/sam.hpp b/src/sam.hpp index 7d3fbdda..6c97d270 100644 --- a/src/sam.hpp +++ b/src/sam.hpp @@ -17,7 +17,7 @@ struct Alignment { int global_ed; int score; int length; - bool is_rc; + bool is_revcomp; bool is_unaligned{false}; // Whether a gapped alignment function was used to obtain this alignment // (even if true, the alignment can still be without gaps) diff --git a/src/unused.cpp b/src/unused.cpp index 3282f5f0..05b01d74 100644 --- a/src/unused.cpp +++ b/src/unused.cpp @@ -1306,7 +1306,7 @@ static inline bool sort_hits(const hit &a, const hit &b) static inline void find_nams_rescue( std::vector &final_nams, robin_hood::unordered_map> &hits_per_ref, - const QueryRandstrobeVector &query_mers, + const std::vector &query_mers, const RefRandstrobeVector &ref_mers, RandstrobeMap &mers_index, int k, @@ -1603,7 +1603,7 @@ static inline void find_nams_rescue( static inline std::pair find_nams( std::vector &final_nams, robin_hood::unordered_map> &hits_per_ref, - const QueryRandstrobeVector &query_mers, + const std::vector &query_mers, const RefRandstrobeVector &ref_mers, RandstrobeMap &mers_index, int k, @@ -1863,8 +1863,8 @@ static inline Alignment get_alignment_unused( const auto read_len = read.size(); const int diff = std::abs(nam.ref_span() - nam.query_span()); - const std::string r_tmp = nam.is_rc ? read.rc : read.seq; - const bool is_rc = nam.is_rc; + const std::string r_tmp = nam.is_revcomp ? read.rc : read.seq; + const bool is_rc = nam.is_revcomp; int ext_left; int ext_right; @@ -1889,7 +1889,7 @@ static inline Alignment get_alignment_unused( alignment.edit_distance = info.edit_distance; alignment.score = info.sw_score; // aln_params.match*(read_len-hamming_dist) - aln_params.mismatch*hamming_dist; alignment.ref_start = ref_start + ext_left + info.query_start; - alignment.is_rc = is_rc; + alignment.is_revcomp = is_rc; alignment.is_unaligned = false; return alignment; } @@ -1960,7 +1960,7 @@ static inline Alignment get_alignment_unused( alignment.edit_distance = alignment_segm_left.edit_distance + alignment_segm_right.edit_distance; alignment.score = alignment_segm_left.score + alignment_segm_right.score; alignment.ref_start = alignment_segm_left.ref_start; - alignment.is_rc = nam.is_rc; + alignment.is_revcomp = nam.is_revcomp; alignment.is_unaligned = false; } else { // std::cerr << "NOOOO MAX BREAKPOINT " << break_point << " candidates: " << n.query_s << " " << n.query_e - k << std::endl; diff --git a/tests/test_sam.cpp b/tests/test_sam.cpp index d7e0dc9e..01f15fe3 100644 --- a/tests/test_sam.cpp +++ b/tests/test_sam.cpp @@ -60,7 +60,7 @@ TEST_CASE("Sam::add") { Alignment aln; aln.ref_id = 0; aln.is_unaligned = false; - aln.is_rc = true; + aln.is_revcomp = true; aln.ref_start = 2; aln.edit_distance = 3; aln.score = 9; @@ -97,7 +97,7 @@ TEST_CASE("Pair with one unmapped SAM record") { aln1.ref_id = 0; aln1.is_unaligned = false; aln1.ref_start = 2; - aln1.is_rc = true; + aln1.is_revcomp = true; aln1.edit_distance = 17; aln1.score = 9; aln1.cigar = Cigar("2M"); @@ -153,7 +153,7 @@ TEST_CASE("TLEN zero when reads map to different contigs") { aln1.ref_id = 0; aln1.is_unaligned = false; aln1.ref_start = 2; - aln1.is_rc = false; + aln1.is_revcomp = false; aln1.edit_distance = 17; aln1.score = 9; aln1.cigar = Cigar("2M"); @@ -162,7 +162,7 @@ TEST_CASE("TLEN zero when reads map to different contigs") { aln2.is_unaligned = false; aln2.ref_id = 1; aln2.ref_start = 3; - aln2.is_rc = false; + aln2.is_revcomp = false; aln2.edit_distance = 2; aln2.score = 4; aln2.cigar = Cigar("3M");