Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Match analysis process became 200X faster in Windows OS #84

Merged
merged 5 commits into from
Aug 31, 2024
Merged

Match analysis process became 200X faster in Windows OS #84

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
8a0a6e1
reduce using push_back
jaebeom-kim Aug 28, 2024
d94db1e
find optimal path and build match connectivity graph at the same time
jaebeom-kim Aug 29, 2024
77b75af
find optimal path and build match connectivity graph at the same time
jaebeom-kim Aug 29, 2024
dedce1d
find optimal path and build match connectivity graph at the same time
jaebeom-kim Aug 29, 2024
ec4391e
Match analysis process became 200X faster in Windows OS
jaebeom-kim Aug 29, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
67 changes: 38 additions & 29 deletions src/commons/KmerMatcher.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -170,9 +170,9 @@ bool KmerMatcher::matchKmers(QueryKmerBuffer * queryKmerBuffer,
remainder--;
}
bool needLastTargetBlock = true;
queryAA = AminoAcidPart(queryKmerList[startIdx].ADkmer);
queryAA = AMINO_ACID_PART(queryKmerList[startIdx].ADkmer);
for (size_t j = 0; j < numOfDiffIdxSplits_use; j ++) {
if (queryAA <= AminoAcidPart(diffIdxSplits.data[j].ADkmer)) {
if (queryAA <= AMINO_ACID_PART(diffIdxSplits.data[j].ADkmer)) {
j = j - (j != 0);
querySplits.emplace_back(startIdx, endIdx, endIdx - startIdx + 1, diffIdxSplits.data[j]);
needLastTargetBlock = false;
Expand Down Expand Up @@ -231,9 +231,14 @@ querySplits, queryKmerList, matchBuffer, cout, targetDiffIdxFileName, numOfDiffI
std::vector<size_t> selectedMatches;
std::vector<uint16_t> selectedHammings;
std::vector<uint32_t> selectedDnaEncodings;
selectedHammingSum.resize(1024);
selectedMatches.resize(1024);
selectedHammings.resize(1024);
selectedDnaEncodings.resize(1024);
size_t selectedMatchCnt = 0;

size_t posToWrite;

int currMatchNum;
size_t idx;
#pragma omp for schedule(dynamic, 1)
for (size_t i = 0; i < querySplits.size(); i++) {
Expand Down Expand Up @@ -268,9 +273,8 @@ querySplits, queryKmerList, matchBuffer, cout, targetDiffIdxFileName, numOfDiffI
// Reuse the comparison data if queries are exactly identical
if (currentQuery == queryKmerList[j].ADkmer
&& (currentQueryInfo.frame/3 == queryKmerList[j].info.frame/3)) {
currMatchNum = selectedMatches.size();
// If local buffer is full, copy them to the shared buffer.
if (matchCnt + currMatchNum > localBufferSize) {
if (matchCnt + selectedMatchCnt > localBufferSize) {
// Check if the shared buffer is full.
posToWrite = matchBuffer->reserveMemory(matchCnt);
if (posToWrite + matchCnt >= matchBuffer->bufferSize) {
Expand All @@ -283,7 +287,7 @@ querySplits, queryKmerList, matchBuffer, cout, targetDiffIdxFileName, numOfDiffI
lastMovedQueryIdx = j;
}
}
for (int k = 0; k < currMatchNum; k++) {
for (int k = 0; k < selectedMatchCnt; k++) {
idx = selectedMatches[k];
// Check if candidateKmerInfos[idx].sequenceID is valid
// if (taxId2genusId.find(candidateKmerInfos[idx].sequenceID) == taxId2genusId.end() ||
Expand All @@ -301,19 +305,15 @@ querySplits, queryKmerList, matchBuffer, cout, targetDiffIdxFileName, numOfDiffI
}
continue;
}
selectedMatches.clear();
selectedHammingSum.clear();
selectedHammings.clear();
selectedDnaEncodings.clear();
selectedMatchCnt = 0;

// Reuse the candidate target k-mers to compare in DNA level if queries are the same at amino acid level but not at DNA level
if (currentQueryAA == AminoAcidPart(queryKmerList[j].ADkmer)) {
if (currentQueryAA == AMINO_ACID_PART(queryKmerList[j].ADkmer)) {
compareDna(queryKmerList[j].ADkmer, candidateTargetKmers, selectedMatches,
selectedHammingSum, selectedHammings, selectedDnaEncodings, queryKmerList[j].info.frame);
currMatchNum = selectedMatches.size();
selectedHammingSum, selectedHammings, selectedDnaEncodings, selectedMatchCnt, queryKmerList[j].info.frame);

// If local buffer is full, copy them to the shared buffer.
if (matchCnt + currMatchNum > localBufferSize) {
if (matchCnt + selectedMatchCnt > localBufferSize) {
// Check if the shared buffer is full.
posToWrite = matchBuffer->reserveMemory(matchCnt);
if (posToWrite + matchCnt >= matchBuffer->bufferSize) {
Expand All @@ -326,7 +326,7 @@ querySplits, queryKmerList, matchBuffer, cout, targetDiffIdxFileName, numOfDiffI
lastMovedQueryIdx = j;
}
}
for (int k = 0; k < currMatchNum; k++) {
for (int k = 0; k < selectedMatchCnt; k++) {
idx = selectedMatches[k];
matches[matchCnt] = {queryKmerList[j].info,
candidateKmerInfos[idx].sequenceID,
Expand All @@ -338,7 +338,7 @@ querySplits, queryKmerList, matchBuffer, cout, targetDiffIdxFileName, numOfDiffI
matchCnt++;
}
currentQuery = queryKmerList[j].ADkmer;
currentQueryAA = AminoAcidPart(currentQuery);
currentQueryAA = AMINO_ACID_PART(currentQuery);
currentQueryInfo = queryKmerList[j].info;
continue;
}
Expand All @@ -347,12 +347,12 @@ querySplits, queryKmerList, matchBuffer, cout, targetDiffIdxFileName, numOfDiffI

// Get next query, and start to find
currentQuery = queryKmerList[j].ADkmer;
currentQueryAA = AminoAcidPart(currentQuery);
currentQueryAA = AMINO_ACID_PART(currentQuery);
currentQueryInfo = queryKmerList[j].info;

// Skip target k-mers that are not matched in amino acid level
while (diffIdxPos != numOfDiffIdx
&& (currentQueryAA > AminoAcidPart(currentTargetKmer))) {
&& (currentQueryAA > AMINO_ACID_PART(currentTargetKmer))) {
if (unlikely(BufferSize < diffIdxBufferIdx + 7)){
loadBuffer(diffIdxFp, diffIdxBuffer, diffIdxBufferIdx, BufferSize, ((int)(BufferSize - diffIdxBufferIdx)) * -1 );
}
Expand All @@ -361,12 +361,12 @@ querySplits, queryKmerList, matchBuffer, cout, targetDiffIdxFileName, numOfDiffI
kmerInfoBufferIdx ++;
}

if (currentQueryAA != AminoAcidPart(currentTargetKmer)) // Move to next query k-mer if there isn't any match.
if (currentQueryAA != AMINO_ACID_PART(currentTargetKmer)) // Move to next query k-mer if there isn't any match.
continue;

// Load target k-mers that are matched in amino acid level
while (diffIdxPos != numOfDiffIdx &&
currentQueryAA == AminoAcidPart(currentTargetKmer)) {
currentQueryAA == AMINO_ACID_PART(currentTargetKmer)) {
candidateTargetKmers.push_back(currentTargetKmer);
candidateKmerInfos.push_back(getKmerInfo(BufferSize, kmerInfoFp, kmerInfoBuffer, kmerInfoBufferIdx));
// Print the target k-mer
Expand Down Expand Up @@ -397,13 +397,19 @@ querySplits, queryKmerList, matchBuffer, cout, targetDiffIdxFileName, numOfDiffI
kmerInfoBufferIdx ++;
}

if (candidateTargetKmers.size() > selectedMatches.size()) {
selectedMatches.resize(candidateTargetKmers.size());
selectedHammingSum.resize(candidateTargetKmers.size());
selectedHammings.resize(candidateTargetKmers.size());
selectedDnaEncodings.resize(candidateTargetKmers.size());
}

// Compare the current query and the loaded target k-mers and select
compareDna(currentQuery, candidateTargetKmers, selectedMatches, selectedHammingSum,
selectedHammings, selectedDnaEncodings, queryKmerList[j].info.frame);
selectedHammings, selectedDnaEncodings, selectedMatchCnt, queryKmerList[j].info.frame);

// If local buffer is full, copy them to the shared buffer.
currMatchNum = selectedMatches.size();
if (matchCnt + currMatchNum > localBufferSize) {
if (matchCnt + selectedMatchCnt > localBufferSize) {
// Check if the shared buffer is full.
posToWrite = matchBuffer->reserveMemory(matchCnt);
if (posToWrite + matchCnt >= matchBuffer->bufferSize) { // full -> write matches to file first
Expand All @@ -417,7 +423,7 @@ querySplits, queryKmerList, matchBuffer, cout, targetDiffIdxFileName, numOfDiffI
}
}

for (int k = 0; k < currMatchNum; k++) {
for (int k = 0; k < selectedMatchCnt; k++) {
idx = selectedMatches[k];
// Check if candidateKmerInfos[idx].sequenceID is valid
// if (taxId2genusId.find(candidateKmerInfos[idx].sequenceID) == taxId2genusId.end() ||
Expand Down Expand Up @@ -493,6 +499,7 @@ void KmerMatcher::compareDna(uint64_t query,
std::vector<uint8_t> &selectedHammingSum,
std::vector<uint16_t> &selectedHammings,
std::vector<uint32_t> &selectedDnaEncodings,
size_t & selectedMatchIdx,
uint8_t frame) {
size_t size = targetKmersToCompare.size();
auto *hammingSums = new uint8_t[size + 1];
Expand All @@ -509,18 +516,20 @@ void KmerMatcher::compareDna(uint64_t query,
}

// Select target k-mers that passed hamming criteria
selectedMatchIdx = 0;
uint8_t maxHamming = min(minHammingSum * 2, 7);
for (size_t h = 0; h < size; h++) {
if (hammingSums[h] <= maxHamming) {
selectedMatches.push_back(h);
selectedHammingSum.push_back(hammingSums[h]);
selectedMatches[selectedMatchIdx] = h;
selectedHammingSum[selectedMatchIdx] = hammingSums[h];
if (frame < 3) { // Frame of query k-mer
selectedHammings.push_back(getHammings(query, targetKmersToCompare[h]));
selectedHammings[selectedMatchIdx] = getHammings(query, targetKmersToCompare[h]);
} else {
selectedHammings.push_back(getHammings_reverse(query, targetKmersToCompare[h]));
selectedHammings[selectedMatchIdx] = getHammings_reverse(query, targetKmersToCompare[h]);
}
// Store right 24 bits of the target k-mer in selectedDnaEncodings
selectedDnaEncodings.push_back(GET_24_BITS_UINT(targetKmersToCompare[h]));
selectedDnaEncodings[selectedMatchIdx] = GET_24_BITS_UINT(targetKmersToCompare[h]);
selectedMatchIdx++;
}
}
delete[] hammingSums;
Expand Down
5 changes: 4 additions & 1 deletion src/commons/KmerMatcher.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,8 @@

#define BufferSize 16'777'216 // 16 * 1024 * 1024 // 16 M

#define AMINO_ACID_PART(kmer) ((kmer) & MARKER)

// Input
// 1. Query K-mers
// 2. Reference K-mers
Expand Down Expand Up @@ -58,7 +60,7 @@ class KmerMatcher {
// search begins.
};

size_t AminoAcidPart(size_t kmer) const { return (kmer)&MARKER; }
inline size_t AminoAcidPart(size_t kmer) const { return (kmer)&MARKER; }

template <typename T>
static void loadBuffer(FILE *fp, T *buffer, size_t &bufferIdx, size_t size) {
Expand Down Expand Up @@ -89,6 +91,7 @@ class KmerMatcher {
std::vector<uint8_t> &selectedHammingSum,
std::vector<uint16_t> &rightEndHammings,
std::vector<uint32_t> &selectedDnaEncodings,
size_t & selectedMatchIdx,
uint8_t frame);

virtual uint8_t getHammingDistanceSum(uint64_t kmer1, uint64_t kmer2);
Expand Down
Loading