Skip to content

Commit 2224cd0

Browse files
h-2h-2
committed
[fix] pairs counting
1 parent 40b49a0 commit 2224cd0

2 files changed

Lines changed: 51 additions & 27 deletions

File tree

‎src/search_algo.hpp‎

Lines changed: 39 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -860,7 +860,7 @@ inline void _writeRecord(TBlastRecord & record, TLocalHolder & lH)
860860
return std::tie(m1._n_sId, m1.qStart, m1.qEnd, m1.sStart, m1.sEnd, m1.qFrameShift, m1.sFrameShift) ==
861861
std::tie(m2._n_sId, m2.qStart, m2.qEnd, m2.sStart, m2.sEnd, m2.qFrameShift, m2.sFrameShift);
862862
});
863-
lH.stats.hitsDuplicate += before - record.matches.size();
863+
lH.stats.hitsDuplicate2 += before - record.matches.size();
864864

865865
// sort by evalue before writing
866866
record.matches.sort([](auto const & m1, auto const & m2) { return m1.bitScore > m2.bitScore; });
@@ -873,6 +873,14 @@ inline void _writeRecord(TBlastRecord & record, TLocalHolder & lH)
873873
}
874874
lH.stats.hitsFinal += record.matches.size();
875875

876+
/* count uniq qry-subj-pairs */
877+
lH.uniqSubjIds.clear();
878+
lH.uniqSubjIds.reserve(record.matches.size());
879+
for (auto const & bm : record.matches)
880+
lH.uniqSubjIds.insert(bm._n_sId);
881+
882+
lH.stats.pairs += lH.uniqSubjIds.size();
883+
876884
// compute LCA
877885
if (lH.options.computeLCA)
878886
{
@@ -922,7 +930,8 @@ inline void _widenMatch(Match & m, TLocalHolder const & lH)
922930
uint64_t band = _bandSize(lH.transQrySeqs[m.qryId].size());
923931

924932
// end on subject is beginning plus full query length plus band
925-
m.subjEnd = std::min<size_t>(m.subjStart + lH.transQrySeqs[m.qryId].size() + band, lH.gH.transSbjSeqs[m.subjId].size());
933+
m.subjEnd =
934+
std::min<size_t>(m.subjStart + lH.transQrySeqs[m.qryId].size() + band, lH.gH.transSbjSeqs[m.subjId].size());
926935

927936
// account for band in subj start
928937
m.subjStart = (band < m.subjStart) ? m.subjStart - band : 0;
@@ -1125,19 +1134,19 @@ inline void _performAlignment(TDepSetH & depSetH,
11251134
}
11261135

11271136
template <typename TLocalHolder>
1128-
inline void _widenAndPreprocessMatches(TLocalHolder & lH)
1137+
inline void _widenAndPreprocessMatches(std::span<Match> & matches, TLocalHolder & lH)
11291138
{
1130-
auto before = lH.matches.size();
1139+
auto before = matches.size();
11311140

1132-
for (Match & m : lH.matches)
1141+
for (Match & m : matches)
11331142
_widenMatch<TLocalHolder>(m, lH);
11341143

1135-
std::ranges::sort(lH.matches);
1144+
std::ranges::sort(matches);
11361145

1137-
if (lH.matches.size() > 1)
1146+
if (matches.size() > 1)
11381147
{
11391148
// pairwise merge from left to right
1140-
for (auto it = lH.matches.begin(); it < lH.matches.end() - 1; ++it)
1149+
for (auto it = matches.begin(); it < matches.end() - 1; ++it)
11411150
{
11421151
Match & l = *it;
11431152
Match & r = *(it + 1);
@@ -1149,7 +1158,7 @@ inline void _widenAndPreprocessMatches(TLocalHolder & lH)
11491158
}
11501159

11511160
// pairwise "swallow" from right to left
1152-
for (auto it = lH.matches.rbegin(); it < lH.matches.rend() - 1; ++it)
1161+
for (auto it = matches.rbegin(); it < matches.rend() - 1; ++it)
11531162
{
11541163
Match & r = *it;
11551164
Match & l = *(it + 1);
@@ -1159,14 +1168,14 @@ inline void _widenAndPreprocessMatches(TLocalHolder & lH)
11591168
}
11601169
}
11611170

1162-
auto const ret = std::ranges::unique(lH.matches);
1163-
lH.matches.erase(ret.begin(), ret.end());
1164-
lH.stats.hitsDuplicate += (before - lH.matches.size());
1171+
auto [new_end, old_end] = std::ranges::unique(matches); // move non-uniq to the end
1172+
matches = std::span<Match>{matches.begin(), new_end}; // "resize" of the span
1173+
lH.stats.hitsDuplicate += (before - matches.size());
11651174
}
11661175
}
11671176

11681177
template <typename TLocalHolder>
1169-
inline void iterateMatchesFullSimd(TLocalHolder & lH, bsDirection const dir = bsDirection::fwd)
1178+
inline void iterateMatchesFullSimd(std::span<Match> lambdaMatches, TLocalHolder & lH, bsDirection const dir)
11701179
{
11711180
using TGlobalHolder = typename TLocalHolder::TGlobalHolder;
11721181
using TBlastMatch = typename TLocalHolder::TBlastMatch;
@@ -1176,7 +1185,7 @@ inline void iterateMatchesFullSimd(TLocalHolder & lH, bsDirection const dir = bs
11761185
// statistics
11771186
#ifdef LAMBDA_MICRO_STATS
11781187
++lH.stats.numQueryWithExt;
1179-
lH.stats.numExtScore += seqan::length(lH.matches);
1188+
lH.stats.numExtScore += seqan::length(lambdaMatches);
11801189

11811190
double start = sysTime();
11821191
#endif
@@ -1186,19 +1195,12 @@ inline void iterateMatchesFullSimd(TLocalHolder & lH, bsDirection const dir = bs
11861195
seqan::StringSet<typename seqan::Source<typename TLocalHolder::TAlignRow1>::Type> depSetV;
11871196

11881197
// pre-sort and filter
1189-
_widenAndPreprocessMatches(lH);
1198+
_widenAndPreprocessMatches(lambdaMatches, lH);
11901199

11911200
// create blast matches from Lambda matches
11921201
std::list<TBlastMatch> blastMatches;
1193-
for (Match const & m : lH.matches)
1202+
for (Match const & m : lambdaMatches)
11941203
{
1195-
/* TODO we need to externalise this because right now we sort and merge hits that are not need in BS mode */
1196-
// In BS-mode, skip those results that have wrong orientation
1197-
if constexpr (TLocalHolder::TGlobalHolder::c_redAlph == AlphabetEnum::DNA3BS)
1198-
{
1199-
if ((dir == bsDirection::fwd && (m.subjId % 2)) || (dir == bsDirection::rev && !(m.subjId % 2)))
1200-
continue;
1201-
}
12021204
// create blastmatch in list without copy or move
12031205
blastMatches.emplace_back(lH.qryIds[m.qryId / TGlobalHolder::qryNumFrames],
12041206
const_gH.indexFile.ids[m.subjId / TGlobalHolder::sbjNumFrames]);
@@ -1362,12 +1364,24 @@ inline void writeRecords(TLocalHolder & lH)
13621364
template <typename TLocalHolder>
13631365
inline void iterateMatches(TLocalHolder & lH)
13641366
{
1365-
iterateMatchesFullSimd(lH, bsDirection::fwd);
13661367
if constexpr (TLocalHolder::TGlobalHolder::c_redAlph == AlphabetEnum::DNA3BS)
13671368
{
1368-
iterateMatchesFullSimd(lH, bsDirection::rev);
1369+
std::ranges::sort(lH.matches,
1370+
[](Match const & l, Match const & r) {
1371+
return std::tuple<bool, Match const &>{l.subjId % 2, l} <
1372+
std::tuple<bool, Match const &>{r.subjId % 2, r};
1373+
});
1374+
1375+
auto it = std::ranges::find_if(lH.matches, [](Match const & m) { return m.subjId % 2; });
1376+
1377+
iterateMatchesFullSimd(std::span{lH.matches.begin(), it}, lH, bsDirection::fwd);
1378+
iterateMatchesFullSimd(std::span{it, lH.matches.end()}, lH, bsDirection::rev);
13691379
lH.blastMatches.sort([](auto const & lhs, auto const & rhs) { return lhs._n_qId < rhs._n_qId; });
13701380
}
1381+
else
1382+
{
1383+
iterateMatchesFullSimd(lH.matches, lH, bsDirection::fwd);
1384+
}
13711385
}
13721386

13731387
//-----------------------------------------------------------------------

‎src/search_datastructures.hpp‎

Lines changed: 12 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -108,10 +108,12 @@ struct StatsHolder
108108
uint64_t hitsFailedExtendEValueTest;
109109
uint64_t hitsAbundant;
110110
uint64_t hitsDuplicate;
111+
uint64_t hitsDuplicate2;
111112

112113
// final
113114
uint64_t hitsFinal;
114115
uint64_t qrysWithHit;
116+
uint64_t pairs;
115117

116118
#ifdef LAMBDA_MICRO_STATS
117119
// times
@@ -146,9 +148,11 @@ struct StatsHolder
146148
hitsFailedExtendEValueTest = 0;
147149
hitsAbundant = 0;
148150
hitsDuplicate = 0;
151+
hitsDuplicate2 = 0;
149152

150153
hitsFinal = 0;
151154
qrysWithHit = 0;
155+
pairs = 0;
152156

153157
#ifdef LAMBDA_MICRO_STATS
154158
seedLengths.clear();
@@ -178,9 +182,11 @@ struct StatsHolder
178182
hitsFailedExtendEValueTest += rhs.hitsFailedExtendEValueTest;
179183
hitsAbundant += rhs.hitsAbundant;
180184
hitsDuplicate += rhs.hitsDuplicate;
185+
hitsDuplicate2 += rhs.hitsDuplicate2;
181186

182187
hitsFinal += rhs.hitsFinal;
183188
qrysWithHit += rhs.qrysWithHit;
189+
pairs += rhs.pairs;
184190

185191
#ifdef LAMBDA_MICRO_STATS
186192
seqan::append(seedLengths, rhs.seedLengths);
@@ -239,6 +245,8 @@ void printStats(StatsHolder const & stats, LambdaOptions const & options)
239245
std::cout << "\n - failed %-identity test " << R << stats.hitsFailedExtendPercentIdentTest << RR
240246
<< (rem -= stats.hitsFailedExtendPercentIdentTest);
241247
std::cout << "\n - duplicates " << R << stats.hitsDuplicate << RR << (rem -= stats.hitsDuplicate);
248+
std::cout << "\n - late duplicates " << R << stats.hitsDuplicate2 << RR
249+
<< (rem -= stats.hitsDuplicate2);
242250
std::cout << "\n - abundant " << R << stats.hitsAbundant << "\033[1m" << RR
243251
<< (rem -= stats.hitsAbundant) << "\033[0m\n\n";
244252

@@ -280,7 +288,8 @@ void printStats(StatsHolder const & stats, LambdaOptions const & options)
280288
if (options.verbosity >= 1)
281289
{
282290
auto const w = seqan::_numberOfDigits(stats.hitsFinal);
283-
std::cout << "Number of valid hits: " << std::setw(w) << stats.hitsFinal
291+
std::cout << "Number of total hits: " << std::setw(w) << stats.hitsFinal
292+
<< "\nNumber of Query-Subject pairs: " << std::setw(w) << stats.pairs
284293
<< "\nNumber of Queries with at least one valid hit: " << std::setw(w) << stats.qrysWithHit << "\n";
285294
}
286295
}
@@ -472,7 +481,8 @@ class LocalDataHolder
472481
std::vector<std::string>, // not used
473482
std::string_view,
474483
uint32_t>;
475-
std::list<TBlastMatch> blastMatches;
484+
std::list<TBlastMatch> blastMatches;
485+
std::unordered_set<uint64_t> uniqSubjIds;
476486

477487
// regarding the gathering of stats
478488
StatsHolder stats{};

0 commit comments

Comments
 (0)