DFS subgraph working prototype with test files

include/input.h

@@ -24,7 +24,7 @@ class InSequencesDBG : public InSequences {
 
 struct UserInputKreeq : UserInput {
 
-    uint32_t covCutOff = 0;
+    uint32_t covCutOff = 0, kmerDepth = 11; // kmer search is in both directions
     uint8_t depth = 3, backtrackingSpan = 5;
     uint64_t maxMem = 0;
 

include/kreeq.h

@@ -139,6 +139,10 @@ class DBG : public Kmap<DBG, UserInputKreeq, uint64_t, DBGkmer, DBGkmer32> { //
 
     void subgraph();
 
+    void DFS();
+
+    ParallelMap32 DFSpass(ParallelMap32* candidates, std::array<uint16_t, 2> mapRange);
+
     void mergeSubgraphs();
 
     void DBGgraphToGFA();

src/generate-tests.cpp

@@ -80,7 +80,7 @@ int main(void) {
     }
     // test subgraph
     file_args = {
-        {{"-d testFiles/test1.kreeq -f testFiles/random1.fasta"}, {""}}
+        {{"-d testFiles/test1.kreeq -f testFiles/random1.fasta --search-depth 0"}, {""}}
     //  {{set of test inputs}, {list of command line args to run with}}
     };
     for(const auto &pair : file_args) {

src/kreeq.cpp

@@ -906,8 +906,8 @@ void DBG::subgraph() {
         deleteMapRange(mapRange);
 
     }
-
     mergeSubgraphs();
+    DFS();
     DBGgraphToGFA();
 
 }
@@ -981,6 +981,106 @@ bool DBG::DBGsubgraphFromSegment(InSegment *inSegment, std::array<uint16_t, 2> m
     return true;
 }
 
+void DBG::DFS() {
+
+    ParallelMap32 candidates, newCandidates;
+    ParallelMap32* subgraph = DBGsubgraph;
+
+    std::array<uint16_t, 2> mapRange = {0,0};
+    for (uint8_t i = 0; i < userInput.kmerDepth; ++i) {
+
+        mapRange = {0,0};
+
+        while (mapRange[1] < mapCount) {
+
+            mapRange = computeMapRange(mapRange);
+            loadMapRange(mapRange);
+            newCandidates = DFSpass(subgraph, mapRange);
+            deleteMapRange(mapRange);
+            candidates.insert(newCandidates.begin(), newCandidates.end());
+            subgraph = &newCandidates;
+        }
+    }
+    DBGsubgraph->insert(candidates.begin(), candidates.end());
+}
+
+ParallelMap32 DBG::DFSpass(ParallelMap32* subgraph, std::array<uint16_t, 2> mapRange) {
+
+    ParallelMap32 newCandidates;
+
+    for (auto pair : *subgraph) {
+
+        for (uint8_t i = 0; i<4; ++i) { // forward edges
+            if (pair.second.fw[i] != 0) {
+
+                uint8_t nextKmer[k];
+                std::string firstKmer = reverseHash(pair.first);
+                firstKmer.push_back(itoc[i]);
+                for (uint8_t e = 0; e<k; ++e)
+                    nextKmer[e] = ctoi[(unsigned char)firstKmer[e+1]];
+
+                ParallelMap *map;
+                ParallelMap32 *map32;
+                bool isFw = false;
+                uint64_t key = hash(nextKmer, &isFw);
+                uint64_t m = key % mapCount;
+
+                if (m >= mapRange[0] && m < mapRange[1]) {
+
+                    map = maps[m];
+                    auto got = map->find(key);
+
+                    if (got != map->end()) {
+
+                        if (got->second.cov != 255) {
+                            newCandidates.insert(*got);
+                        }else{
+                            map32 = maps32[m];
+                            auto got = map32->find(key);
+                            newCandidates.insert(*got);
+                        }
+                    }
+                }
+            }
+        }
+        for (uint8_t i = 0; i<4; ++i) { // reverse edges
+            if (pair.second.bw[i] != 0) {
+
+                uint8_t nextKmer[k];
+                std::string firstKmer;
+                firstKmer.push_back(itoc[i]);
+                firstKmer.append(reverseHash(pair.first));
+
+                for (uint8_t e = 0; e<k; ++e)
+                    nextKmer[e] = ctoi[(unsigned char)firstKmer[e]];
+
+                ParallelMap *map;
+                ParallelMap32 *map32;
+                bool isFw = false;
+                uint64_t key = hash(nextKmer, &isFw);
+                uint64_t m = key % mapCount;
+
+                if (m >= mapRange[0] && m < mapRange[1]) {
+
+                    map = maps[m];
+                    auto got = map->find(key);
+
+                    if (got != map->end()) {
+                        if (got->second.cov != 255) {
+                            newCandidates.insert(*got);
+                        }else{
+                            map32 = maps32[m];
+                            auto got = map32->find(key);
+                            newCandidates.insert(*got);
+                        }
+                    }
+                }
+            }
+        }
+    }
+    return newCandidates;
+}
+
 void DBG::mergeSubgraphs() {
 
     for (ParallelMap32 *map1 : DBGTmpSubgraphs) {

src/main.cpp

@@ -196,6 +196,7 @@ int main(int argc, char **argv) {
                         printf("\t-m --max-memory use at most this amount of memory (in Gb, default: 0.9 of max).\n");
                         printf("\t-j --threads <n> numbers of threads (default: max).\n");
                         printf("\t-v --version software version.\n");
+                        printf("\t--search-depth the max depth for graph traversal (default: 3).\n");
                         printf("\t--cmd print $0 to stdout.\n");
                         exit(0);
                 }
@@ -290,7 +291,9 @@ int main(int argc, char **argv) {
             static struct option long_options[] = { // struct mapping long options
                 {"database", required_argument, 0, 'd'},
                 {"input-sequence", required_argument, 0, 'f'},
+                {"search-depth", required_argument, 0, 0},
                 {"out-format", required_argument, 0, 'o'},
+                {"input-positions", required_argument, 0, 'p'},
 
                 {"threads", required_argument, 0, 'j'},
                 {"verbose", no_argument, &verbose_flag, 1},
@@ -302,7 +305,7 @@ int main(int argc, char **argv) {
             while (true) { // loop through argv
 
                 int option_index = 1;
-                c = getopt_long(argc, argv, "-:d:f:j:o:h",
+                c = getopt_long(argc, argv, "-:d:f:j:o:p:h",
                                 long_options, &option_index);
 
                 if (c == -1) // exit the loop if run out of options
@@ -321,8 +324,8 @@ int main(int argc, char **argv) {
                         break;
                     default: // handle positional arguments
                     case 0: // case for long options without short options
-                        //                if (strcmp(long_options[option_index].name,"line-length") == 0)
-                        //                  splitLength = atoi(optarg);
+                        if(strcmp(long_options[option_index].name,"search-depth") == 0)
+                            userInput.kmerDepth = atoi(optarg);
                         break;
                     case 'd': // input sequence
 
@@ -350,13 +353,19 @@ int main(int argc, char **argv) {
                     case 'o': // handle output (file or stdout)
                         userInput.outFile = optarg;
                         break;
+                    case 'p': // input coordinates
+                        ifFileExists(optarg);
+                        userInput.inBedInclude = optarg;
+                        break;
                     case 'h': // help
                         printf("kreeq subgraph [options]\n");
                         printf("\nOptions:\n");
                         printf("\t-d --database DBG database.\n");
                         printf("\t-f --input-sequence sequence input file (fasta).\n");
+                        printf("\t--search-depth the max depth for graph traversal (default: 3).\n");
                         printf("\t-j --threads <n> numbers of threads (default: max).\n");
                         printf("\t-o --out-format generates various kinds of outputs (currently supported: .gfa1/2).\n");
+                        printf("\t-p --input-positions BED coordinates of positions to extract kmers from.\n");
                         printf("\t--cmd print $0 to stdout.\n");
                         exit(0);
                 }

testFiles/random5.fasta

@@ -0,0 +1,2 @@
+>sequence1 - no errors
+CATACTACGATCAGATCGACTGACTCGTACATGCAGTGATGTACGAATGCATGCATCGATCGATCGAAATGATAGCATGACTCAGACTGATCAGATCGA

testFiles/random5.fastq

@@ -0,0 +1,8 @@
+@sequence1 - no errors
+CATACTACGATCAGATCGACTGACTCGTACATGCAGTGATGTACGAATGCATGCATCGATCGATCGAAATGATAGCATGACTCAGACTGATCAGATCGA
++
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+@sequence2 - two point errors apart (25:T>A, 65:C>T)
+CATACTACGATCAGATCGACTGACaCGTACATGCAGTGATGTACGAATGCATGCATCGATCGATtGAAATGATAGCATGACTCAGACTGATCAGATCGA
++
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

testFiles/random5.kreeq/.index

@@ -0,0 +1,2 @@
+21
+128