improve performance and reduce false positive

OpenGene · Feb 18, 2022 · 69ddeba · 69ddeba
1 parent a52d31b
commit 69ddeba
Show file tree

Hide file tree

Showing 2 changed files with 23 additions and 8 deletions.
diff --git a/src/common.h b/src/common.h
@@ -38,5 +38,9 @@ static const int PACK_SIZE = 1000;
 // if the number of in memory packs is full, the producer thread should sleep
 static const int PACK_IN_MEM_LIMIT = 100;
 
+// the key dup in normal level will be kept, in high level will be skipped
+static const int DUPE_NORMAL_LEVEL = -1;
+static const int DUPE_HIGH_LEVEL = -2;
+
 
 #endif /* COMMON_H */
diff --git a/src/indexer.cpp b/src/indexer.cpp
@@ -85,16 +85,24 @@ void Indexer::indexContig(int ctg, string seq, int start) {
         if(mKmerPos.count(kmer) >0 ){
             GenePos gp = mKmerPos[kmer];
             // already marked as a dupe
-            if(gp.contig < 0) {
-                mDupeList[gp.position].push_back(site);
+            if(gp.contig == DUPE_HIGH_LEVEL){
+                // skip it since it's already a high level dupe
+                continue;
+            } else if(gp.contig == DUPE_NORMAL_LEVEL) {
+                if(mDupeList[gp.position].size() >= GlobalSettings::skipKeyDupThreshold) {
+                    mKmerPos[kmer].contig = DUPE_HIGH_LEVEL;
+                    mDupeList[gp.position]=vector<GenePos>();
+                }
+                else
+                    mDupeList[gp.position].push_back(site);
             } else {
                 // else make a new dupe entry
                 vector<GenePos> gps;
                 gps.push_back(gp);
                 gps.push_back(site);
                 mDupeList.push_back(gps);
                 // and mark it as a dupe
-                mKmerPos[kmer].contig = -1;
+                mKmerPos[kmer].contig = DUPE_NORMAL_LEVEL;
                 mKmerPos[kmer].position = mDupeList.size() -1;
                 mUniquePos--;
                 mDupePos++;
@@ -126,18 +134,18 @@ vector<SeqMatch> Indexer::mapRead(Read* r) {
         }
         GenePos gp = mKmerPos[kmer];
         // is a dupe
-        if(gp.contig < 0) {
+        if(gp.contig == DUPE_HIGH_LEVEL) {
             // too much keys in this dupe, then skip it
-            if(mDupeList[gp.position].size() > GlobalSettings::skipKeyDupThreshold)
-                continue;
+            continue;
+        } else if(gp.contig == DUPE_NORMAL_LEVEL) {
             for(int g=0; g<mDupeList[gp.position].size();g++) {
                 long gplong = gp2long(shift(mDupeList[gp.position][g], i));
                 if(kmerStat.count(gplong)==0)
                     kmerStat[gplong] = 1;
                 else
                     kmerStat[gplong] += 1;
             }
-        } else {
+        } else { // not a dupe
             long gplong = gp2long(shift(gp, i));
             if(kmerStat.count(gplong)==0)
                 kmerStat[gplong] = 1;
@@ -183,7 +191,10 @@ vector<SeqMatch> Indexer::mapRead(Read* r) {
             continue;
         GenePos gp = mKmerPos[kmer];
         // is a dupe
-        if(gp.contig < 0) {
+        if(gp.contig == DUPE_HIGH_LEVEL) {
+            // too much keys in this dupe, then skip it
+            continue;
+        } else if(gp.contig == DUPE_NORMAL_LEVEL) {
             for(int g=0; g<mDupeList[gp.position].size();g++) {
                 long gplong = gp2long(shift(mDupeList[gp.position][g], i));
                 if(abs(gplong - gp1) <= 1)