correctly update edges in subgraph

Author: gf777

include/kreeq.h

@@ -243,6 +243,8 @@ class DBG : public Kmap<DBG, UserInputKreeq, uint64_t, DBGkmer, DBGkmer32> { //
 
     void bestFirst(); // best-first search traversal
 
+    void removeMissingEdges();
+    
     void summary(ParallelMap32color& DBGsubgraph);
 
     ParallelMap32color traversalPass(ParallelMap32color* candidates, std::array<uint16_t, 2> mapRange);

src/main.cpp

@@ -377,7 +377,7 @@ int main(int argc, char **argv) {
                     case 'h': // help
                         printf("kreeq subgraph [options]\n");
                         printf("\nOptions:\n");
-                        printf("\t-c --coverage-cutoff coverage cutoff.\n");
+                        printf("\t-c --coverage-cutoff exclude nodes and edges below or equal to the cutoff (default: 0).\n");
                         printf("\t-d --database DBG database.\n");
                         printf("\t-f --input-sequence sequence input file (fasta).\n");
                         printf("\t--traversal-algorithm <string> the approach used for graph search (best-first/traversal, default: best-first).\n");

src/subgraph.cpp

@@ -171,6 +171,8 @@ void DBG::subgraph() {
         fprintf(stderr, "Cannot find input algorithm (%s). Terminating.\n", userInput.travAlgorithm.c_str());
         exit(EXIT_FAILURE);
     }
+    lg.verbose("Remove missing edges");
+    removeMissingEdges();
     lg.verbose("Computing summary graph");
     summary(*DBGsubgraph);
     lg.verbose("Generating GFA");
@@ -576,6 +578,34 @@ void DBG::buildNextKmer(uint8_t* nextKmer, uint64_t hash, uint8_t nextBase, bool
     }
 }
 
+void DBG::removeMissingEdges() {
+    
+    bool isFw;
+    
+    for (auto &node : *DBGsubgraph) { // loop through all nodes
+        
+        for (uint8_t i = 0; i<4; ++i) { // forward edges
+            if (node.second.fw[i] > userInput.covCutOff) {
+                uint8_t nextKmer[k];
+                buildNextKmer(nextKmer, node.first, i, true); // compute next node
+                uint64_t key = hash(nextKmer, &isFw);
+                auto found = DBGsubgraph->find(key);
+                if (found == DBGsubgraph->end()) {
+                    node.second.fw[i] = 0;
+                }
+            }
+            if (node.second.bw[i] > userInput.covCutOff) { // backward edges
+                uint8_t nextKmer[k];
+                buildNextKmer(nextKmer, node.first, i, false); // compute next node
+                uint64_t key = hash(nextKmer, &isFw);
+                auto found = DBGsubgraph->find(key);
+                if (found == DBGsubgraph->end()) {
+                    node.second.bw[i] = 0;
+                }
+            }
+        }
+    }
+}