include original sequence and colored nodes

include/kreeq.h

@@ -123,19 +123,29 @@ struct DBGkmer32 {
 
 };
 
-using ParallelMap = phmap::parallel_flat_hash_map<uint64_t, DBGkmer,
+struct DBGkmer32color : public DBGkmer32 {
+    uint8_t color = 0;
+
+    DBGkmer32color() {}
+
+    DBGkmer32color(const DBGkmer32& dbgkmer32) {
+        std::copy(std::begin(dbgkmer32.fw), std::end(dbgkmer32.fw), std::begin(fw));
+        std::copy(std::begin(dbgkmer32.bw), std::end(dbgkmer32.bw), std::begin(bw));
+        cov = dbgkmer32.cov;
+    }
+};
+
+template<typename T>
+class PM : public phmap::parallel_flat_hash_map<uint64_t, T,
                                           std::hash<uint64_t>,
                                           std::equal_to<uint64_t>,
-                                          std::allocator<std::pair<const uint64_t, DBGkmer>>,
+                                          std::allocator<std::pair<const uint64_t, T>>,
                                           8,
-                                          phmap::NullMutex>;
+phmap::NullMutex> {};
 
-using ParallelMap32 = phmap::parallel_flat_hash_map<uint64_t, DBGkmer32,
-std::hash<uint64_t>,
-                                          std::equal_to<uint64_t>,
-                                          std::allocator<std::pair<const uint64_t, DBGkmer32>>,
-                                          8,
-                                          phmap::NullMutex>;
+using ParallelMap = PM<DBGkmer>;
+using ParallelMap32 = PM<DBGkmer32>;
+using ParallelMap32color = PM<DBGkmer32color>;
 
 class DBG : public Kmap<DBG, UserInputKreeq, uint64_t, DBGkmer, DBGkmer32> { // CRTP
 
@@ -145,8 +155,8 @@ class DBG : public Kmap<DBG, UserInputKreeq, uint64_t, DBGkmer, DBGkmer32> { //
     InSequencesDBG *genome;
 
     // subgraph objects
-    ParallelMap32 *DBGsubgraph = new ParallelMap32;
-    std::vector<ParallelMap32*> DBGTmpSubgraphs;
+    ParallelMap32color *DBGsubgraph = new ParallelMap32color;
+    std::vector<ParallelMap32color*> DBGTmpSubgraphs;
     InSequences GFAsubgraph;
 
     uint64_t totEdgeCount = 0;
@@ -185,7 +195,8 @@ class DBG : public Kmap<DBG, UserInputKreeq, uint64_t, DBGkmer, DBGkmer32> { //
 
     bool mergeSubMaps(ParallelMap32* map1, ParallelMap32* map2, uint8_t subMapIndex);
 
-    bool unionSum(ParallelMap32* map1, ParallelMap32* map2);
+    template<typename MAPTYPE>
+    bool unionSum(MAPTYPE* map1, MAPTYPE* map2);
 
     void kunion();
 
@@ -229,9 +240,9 @@ class DBG : public Kmap<DBG, UserInputKreeq, uint64_t, DBGkmer, DBGkmer32> { //
 
     void DFS();
 
-    void summary(ParallelMap32& DBGsubgraph);
+    void summary(ParallelMap32color& DBGsubgraph);
 
-    ParallelMap32 DFSpass(ParallelMap32* candidates, std::array<uint16_t, 2> mapRange);
+    ParallelMap32color DFSpass(ParallelMap32color* candidates, std::array<uint16_t, 2> mapRange);
 
     void mergeSubgraphs();
 

src/graph-builder.cpp

@@ -475,8 +475,8 @@ bool DBG::mergeSubMaps(ParallelMap32* map1, ParallelMap32* map2, uint8_t subMapI
 
 }
 
-
-bool DBG::unionSum(ParallelMap32* map1, ParallelMap32* map2) {
+template<typename MAPTYPE>
+bool DBG::unionSum(MAPTYPE* map1, MAPTYPE* map2) {
 
     std::vector<std::function<bool()>> jobs;
 

src/kreeq.cpp

@@ -913,7 +913,7 @@ void DBG::subgraph() {
 
 }
 
-void DBG::summary(ParallelMap32& DBGsubgraph) {
+void DBG::summary(ParallelMap32color& DBGsubgraph) {
 
     uint64_t tot = 0, kmersUnique = 0, kmersDistinct = DBGsubgraph.size(), edgeCount = 0;
     phmap::parallel_flat_hash_map<uint64_t, uint64_t> hist;
@@ -948,7 +948,7 @@ bool DBG::DBGsubgraphFromSegment(InSegment *inSegment, std::array<uint16_t, 2> m
     std::string sHeader = inSegment->getSeqHeader();
     ParallelMap *map;
     ParallelMap32 *map32;
-    ParallelMap32 *segmentSubmap = new ParallelMap32;
+    ParallelMap32color *segmentSubmap = new ParallelMap32color;
     uint64_t key, i;
     bool isFw = false;
     std::vector<uint64_t> segmentCoordinates;
@@ -990,12 +990,43 @@ bool DBG::DBGsubgraphFromSegment(InSegment *inSegment, std::array<uint16_t, 2> m
 
                 if (got != map->end()) {
                     if (got->second.cov != 255) {
-                        segmentSubmap->insert(*got);
+                        DBGkmer32color dbgKmer32color(got->second);
+                        dbgKmer32color.color = 1;
+                        segmentSubmap->insert(std::make_pair(got->first,dbgKmer32color));
                     }else{
                         map32 = maps32[i];
                         auto got = map32->find(key);
-                        segmentSubmap->insert(*got);
+                        DBGkmer32color dbgKmer32color(got->second);
+                        dbgKmer32color.color = 1;
+                        segmentSubmap->insert(std::make_pair(got->first,dbgKmer32color));
                     }
+                }else{ // construct the kmer
+
+                    DBGkmer32color dbgKmer32color;
+                    dbgKmer32color.color = 2;
+                    edgeBit edges;
+
+                    if (isFw){
+                        if (ctoi[*(first+p+k)] <= 3)
+                            edges.assign(ctoi[*(first+p+k)]);
+                        if (p > 0 && *(str+p-1) <= 3)
+                            edges.assign(4+*(str+p-1));
+                    }else{
+                        if (p > 0 && *(str+p-1) <= 3)
+                            edges.assign(3-*(str+p-1));
+                        if (ctoi[*(first+p+k)] <= 3)
+                            edges.assign(4+3-ctoi[*(first+p+k)]);
+                    }
+
+                    for (uint64_t w = 0; w<4; ++w) { // update weights
+
+                        dbgKmer32color.fw[w] += edges.read(w);
+                        dbgKmer32color.bw[w] += edges.read(4+w);
+                    }
+                    if (dbgKmer32color.cov < LARGEST)
+                        ++dbgKmer32color.cov; // increase kmer coverage
+
+                    segmentSubmap->insert(std::make_pair(key,dbgKmer32color));
                 }
             }
         }
@@ -1011,8 +1042,8 @@ bool DBG::DBGsubgraphFromSegment(InSegment *inSegment, std::array<uint16_t, 2> m
 
 void DBG::DFS() {
 
-    ParallelMap32 candidates, newCandidates;
-    ParallelMap32* subgraph = DBGsubgraph;
+    ParallelMap32color candidates, newCandidates;
+    ParallelMap32color* subgraph = DBGsubgraph;
 
     std::array<uint16_t, 2> mapRange = {0,0};
     for (uint8_t i = 0; i < userInput.kmerDepth; ++i) {
@@ -1032,9 +1063,9 @@ void DBG::DFS() {
     DBGsubgraph->insert(candidates.begin(), candidates.end());
 }
 
-ParallelMap32 DBG::DFSpass(ParallelMap32* subgraph, std::array<uint16_t, 2> mapRange) {
+ParallelMap32color DBG::DFSpass(ParallelMap32color* subgraph, std::array<uint16_t, 2> mapRange) {
 
-    ParallelMap32 newCandidates;
+    ParallelMap32color newCandidates;
 
     for (auto pair : *subgraph) {
 
@@ -1066,11 +1097,13 @@ ParallelMap32 DBG::DFSpass(ParallelMap32* subgraph, std::array<uint16_t, 2> mapR
                         if (got != map->end()) {
 
                             if (got->second.cov != 255) {
-                                newCandidates.insert(*got);
+                                DBGkmer32color dbgKmer32color(got->second);
+                                newCandidates.insert(std::make_pair(got->first,dbgKmer32color));
                             }else{
                                 map32 = maps32[m];
                                 auto got = map32->find(key);
-                                newCandidates.insert(*got);
+                                DBGkmer32color dbgKmer32color(got->second);
+                                newCandidates.insert(std::make_pair(got->first,dbgKmer32color));
                             }
                         }
                     }
@@ -1106,11 +1139,13 @@ ParallelMap32 DBG::DFSpass(ParallelMap32* subgraph, std::array<uint16_t, 2> mapR
 
                         if (got != map->end()) {
                             if (got->second.cov != 255) {
-                                newCandidates.insert(*got);
+                                DBGkmer32color dbgKmer32color(got->second);
+                                newCandidates.insert(std::make_pair(got->first,dbgKmer32color));
                             }else{
                                 map32 = maps32[m];
                                 auto got = map32->find(key);
-                                newCandidates.insert(*got);
+                                DBGkmer32color dbgKmer32color(got->second);
+                                newCandidates.insert(std::make_pair(got->first,dbgKmer32color));
                             }
                         }
                     }
@@ -1123,11 +1158,10 @@ ParallelMap32 DBG::DFSpass(ParallelMap32* subgraph, std::array<uint16_t, 2> mapR
 
 void DBG::mergeSubgraphs() {
 
-    for (ParallelMap32 *map1 : DBGTmpSubgraphs) {
+    for (ParallelMap32color *map1 : DBGTmpSubgraphs) {
         unionSum(map1, DBGsubgraph);
         delete map1;
     }
-
 }
 
 void DBG::DBGgraphToGFA() {

test.gfa

@@ -0,0 +1,17 @@
+H	VN:Z:1.2
+S	0	GAAATGATAGCATGACTCAGACTGATCAGATCGA	RC:i:42
+S	1	ATCGATCGATGCATGCATTCGTACATCACTGCATGTACG	RC:i:42
+S	2	CTGAGTCATGCTATCATTTCAATCGATCGATGCATGCATTC	RC:i:21
+S	3	CTACGATCAGATCGACTGACTCGTACATGCAGTGATGTACG	RC:i:21
+S	4	GAATGCATGCATCGATCGATCGAAATGATAGCATGACTCAG	RC:i:21
+S	5	CTACGATCAGATCGACTGACACGTACATGCAGTGATGTACG	RC:i:21
+S	6	GTCAGTCGATCTGATCGTAGTATG	RC:i:42
+L	0	-	2	+	1N20M	KC:i:1
+L	2	+	1	+	1N20M	KC:i:1
+L	4	-	1	+	1N20M	KC:i:1
+L	1	+	3	-	1N20M	KC:i:1
+L	3	-	6	+	1N20M	KC:i:1
+L	5	+	1	-	1N20M	KC:i:1
+L	4	+	0	+	1N20M	KC:i:1
+L	6	-	5	+	1N20M	KC:i:1
+P	path1	0-,2+,1+,3-,6+	1N20M,1N20M,1N20M,1N20M

testFiles/random11.fasta

@@ -0,0 +1,2 @@
+>sequence7 - two point errors overlapping (25:T>A, 35:C>T)
+CATACTACGATCAGATCGACTGACaCGTACATGCtGTGATGTACGAATGCATGCATCGATCGATCGAAATGATAGCATGACTCAGACTGATCAGATCGA

testFiles/random11.fastq

@@ -0,0 +1,4 @@
+@sequence1 - no errors
+CATACTACGATCAGATCGACTGACTCGTACATGCAGTGATGTACGAATGCATGCATCGATCGATCGAAATGATAGCATGACTCAGACTGATCAGATCGA
++
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

testFiles/random11.kreeq/.index

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