More work on integrating gaf.

kcleal · kcleal · commit 3c2fb51b1232 · 2025-02-19T16:19:58.000Z
diff --git a/src/alignment_format.cpp b/src/alignment_format.cpp
@@ -45,7 +45,7 @@ namespace AlignFormat {
                         cur = cur->duplicate();
                         other_segments.push_back(cur);
                     }
-                    cur->core.flag = (line[i] == '>') ? 0 : 16;
+                    cur->flag = (line[i] == '>') ? 0 : 16;
                     ++match_count;
                     break;
                 case 1:
@@ -62,16 +62,16 @@ namespace AlignFormat {
                     break;
                 case 2:
                     j = line.find("-", i);
-                    std::from_chars(line.data() + i, line.data() + j, cur->core.pos);
+                    std::from_chars(line.data() + i, line.data() + j, cur->pos);
                     ++match_count;
                     len = j - i;
                     i += len;
                     break;
                 default:
                     j = line.find_first_of("><", i);
-                    std::from_chars(line.data() + i, line.data() + j, cur->core.end);
+                    std::from_chars(line.data() + i, line.data() + j, cur->end);
                     if (!have_cigar) {
-                        cur->blocks.emplace_back() = {(uint32_t)cur->core.pos, (uint32_t)cur->core.end};
+                        cur->blocks.emplace_back() = {(uint32_t)cur->pos, (uint32_t)cur->end};
                     }
                     if (j == std::string::npos) {
                         i = pe;
@@ -89,15 +89,15 @@ namespace AlignFormat {
             cur = g;
             size_t next_segment = 0;
             uint32_t num = 0;
-            uint32_t current_pos = cur->core.pos;
+            uint32_t current_pos = cur->pos;
             uint32_t block_start = current_pos;
             bool in_match = false;
 
             auto flush_block = [&]() {
                 if (in_match && cur) {
                     // Ensure block doesn't extend beyond segment boundary and check for overflow
-                    uint32_t end_pos = (current_pos < (uint32_t)cur->core.end) ?
-                                       std::min(current_pos, (uint32_t)cur->core.end) : (uint32_t)cur->core.end;
+                    uint32_t end_pos = (current_pos < (uint32_t)cur->end) ?
+                                       std::min(current_pos, (uint32_t)cur->end) : (uint32_t)cur->end;
                     if (block_start < end_pos) {
                         cur->blocks.emplace_back() = {block_start, end_pos};
                     }
@@ -111,7 +111,7 @@ namespace AlignFormat {
                     if (next) {
                         cur = next;
                         next_segment++;
-                        current_pos = cur->core.pos;
+                        current_pos = cur->pos;
                         block_start = current_pos;
                         in_match = false;
                         return true;
@@ -136,8 +136,8 @@ namespace AlignFormat {
                         uint32_t remaining = num;
                         while (remaining > 0 && cur != nullptr) {
                             // Calculate how much of the operation fits in current segment
-                            uint32_t segment_remaining = (current_pos < (uint32_t)cur->core.end) ?
-                                                         ((uint32_t)cur->core.end - current_pos) : 0;
+                            uint32_t segment_remaining = (current_pos < (uint32_t)cur->end) ?
+                                                         ((uint32_t)cur->end - current_pos) : 0;
 
                             // Protect against overflow in advance calculation
                             uint32_t max_advance = UINT32_MAX - current_pos;
@@ -152,7 +152,7 @@ namespace AlignFormat {
                             remaining -= advance;
 
                             // If we've reached the end of the segment
-                            if (current_pos >= (uint32_t)cur->core.end) {
+                            if (current_pos >= (uint32_t)cur->end) {
                                 flush_block();
                                 if (!advance_to_next_segment()) {
                                     break;  // No more segments
@@ -167,16 +167,16 @@ namespace AlignFormat {
                         flush_block();
                         uint32_t remaining = num;
                         while (remaining > 0 && cur != nullptr) {
-                            uint32_t segment_remaining = (current_pos < (uint32_t)cur->core.end) ?
-                                                         ((uint32_t)cur->core.end - current_pos) : 0;
+                            uint32_t segment_remaining = (current_pos < (uint32_t)cur->end) ?
+                                                         ((uint32_t)cur->end - current_pos) : 0;
 
                             uint32_t max_advance = UINT32_MAX - current_pos;
                             uint32_t advance = std::min({remaining, segment_remaining, max_advance});
 
                             current_pos += advance;
                             remaining -= advance;
 
-                            if (current_pos >= (uint32_t)cur->core.end) {
+                            if (current_pos >= (uint32_t)cur->end) {
                                 if (!advance_to_next_segment()) {
                                     break;  // No more segments
                                 }
@@ -199,12 +199,12 @@ namespace AlignFormat {
             flush_block();
         }
 
-//        std::cout << "g is " << g->core.pos << " - " << g->core.end << std::endl;
+//        std::cout << "g is " << g->pos << " - " << g->end << std::endl;
 //        for (auto item : g->blocks) {
 //            std::cout << item.start << "-" << item.end << " ";
 //        } std::cout << std::endl;
 //        for (auto gg : other_segments) {
-//            std::cout << "gg is " << gg->core.pos << " - " << gg->core.end << std::endl;
+//            std::cout << "gg is " << gg->pos << " - " << gg->end << std::endl;
 //            for (auto item : gg->blocks) {
 //                std::cout << item.start << "-" << item.end << " ";
 //            } std::cout << std::endl;
@@ -231,8 +231,7 @@ namespace AlignFormat {
 // ds:Z::36*ct:129*tg*tc:2+[ag]gcgcag:1*ag:228+[ag]:19+ca:4*ta:11*ga-[caggcgcagaga]ggcgcgccgcgcc[ggcg]:21*cg:4*tc:7*ta:29*ag:3*ct:2*ct:1*gc:32+[g]:33*ac:14*gt:1*tg:15*ga:28+[g]:11*ag:3*tg+[g]:33*tc:20*ag:19*ac:12*ca:223*ta:17*ag:47*cg:50*ag:33*at:22*gt:206*ga:400-[g]:24*gc:184*ga:11+[g]:388*ga:211-[c]:120*tg:1*ag:183*ct:24*gc:329-ag:237*ca:61-[c]:21*tc:365-[c]:109+[aaga]g:217*gc:25*cg:42+[g]:165*ag:22*ag:187*ag:92*tg:62*at:542*gt:2*gt:85*ac+[g]:15*ag:33*ga:111*tc:34*tg:194*ct:79*tc:155*ct:36*ga:74*ct:160*gc:1*cg:46*ag:1307*tc:684*cg:322*ag:68*ct:534*ag:351*ga:92*gc:5*tc:796*ct:89-[c]:465*ct:577*ag:23*tc:586+[g]:302*tc:83*ga:18*gt:25*ac:102+[g]:56*cg:35*ag:312*at:1*gc:170*ca:93*ga:156*ct:136*ct:453*ag:108*ag:67*cg:11*ga:3*gc:17*gt:4*gc:1*ca:15*cg:347+[g]:8*ag:10*cg*cg:21*ac:93+[a]:292*ag:142*ct:107*tc:281*ct:444*ga:284*ag:57*ga:77*ga:309*ga:166+[t]:267*ta:24
 
     void gafParser(std::string& line,
-                   ankerl::unordered_dense::map< std::string, SuperIntervals<int, GAF_t *>>& cached_alignments,
-                   ankerl::unordered_dense::map< std::string, uint32_t>& rids) {
+                   ankerl::unordered_dense::map< std::string, SuperIntervals<int, GAF_t *>>& cached_alignments) {
         GAF_t *g = new GAF_t();
 
         size_t pos = 0;
@@ -247,7 +246,8 @@ namespace AlignFormat {
         };
 
         next = line.find('\t', pos);
-        g->qname = line.substr(0, next);
+        std::string qname = line.substr(0, next);
+
         pos = next + 1;
 
         parse_int(g->qlen);
@@ -287,7 +287,7 @@ namespace AlignFormat {
         pos = next + 1;
 
         //12 	int 	Mapping quality (0-255; 255 for missing)
-        parse_int(g->core.qual);
+        parse_int(g->qual);
 
         size_t cigar_start = line.find("cg:Z:");
         size_t cigar_end = std::string::npos;
@@ -299,13 +299,13 @@ namespace AlignFormat {
         std::vector<GAF_t*> other_segments;
         extractAlignmentPath(line, g, ps, pe, other_segments, cigar_start, cigar_end);
 
-        cached_alignments[g->chrom].add(g->core.pos, g->core.end, g);
+        cached_alignments[g->chrom].add(g->pos, g->end, g);
 
-//        std::cout << g->core.flag << " " << g->chrom << " " << g->core.pos << " " << g->core.end << std::endl;
-        for (auto & v : other_segments) {
-            cached_alignments[v->chrom].add(v->core.pos, v->core.end, g);
-//            std::cout << v->core.flag << " " << v->chrom << " " << v->core.pos << " " << v->core.end << std::endl;
-        }
+//        std::cout << g->flag << " " << g->chrom << " " << g->pos << " " << g->end << std::endl;
+//        for (auto & v : other_segments) {
+//            cached_alignments[v->chrom].add(v->pos, v->end, g);
+//            std::cout << v->flag << " " << v->chrom << " " << v->pos << " " << v->end << std::endl;
+//        }
 //        std::exit(0);
 
     }
@@ -314,26 +314,26 @@ namespace AlignFormat {
         int start, end;
     };
 
-    void gafFindY(std::vector<GAF_t *>& gafAlignments) {
+    void gafFindY(std::vector<AlignFormat::GAF_t *>& gafAlignments) {
         std::vector<TrackRange> trackLevels;
         for (const auto &b : gafAlignments) {
             size_t memLen = trackLevels.size();
             size_t i = 0;
             for (; i < memLen; ++i) {
-                if (b->core.pos > trackLevels[i].end) {
-                    trackLevels[i].end = b->core.end + 2;
+                if (b->pos > trackLevels[i].end) {
+                    trackLevels[i].end = b->end;
                     b->y = (int)i;
                     break;
                 }
             }
             if (i == memLen) {
-                trackLevels.emplace_back() = {b->core.pos, b->core.end + 2};
+                trackLevels.emplace_back() = {b->pos, b->end};
                 b->y = (int)memLen;
             }
         }
     }
 
-    void GwAlignment::open(const std::string& file_path, const std::string& reference, int threads, faidx_t* fai) {
+    void AlignFormat::GwAlignment::open(const std::string& file_path, const std::string& reference, int threads) {
         path = file_path;
         if (Utils::endsWith(path, "bam") || Utils::endsWith(path, "cram")) {
             type = AlignmentType::HTSLIB_t;
@@ -371,24 +371,16 @@ namespace AlignFormat {
                 throw std::exception();
             }
 #endif
-            // Reference id's
-            ankerl::unordered_dense::map< std::string, uint32_t> rids;
-            int num_sequences = faidx_nseq(fai);
-            for (int i = 0; i < num_sequences; ++i) {
-                const char* seq_name = faidx_iseq(fai, i);
-                rids[seq_name] = i;
-            }
 
             std::string tp;
             while (true) {
-                auto got_line = (bool)getline(*fpu, tp);
-                if (!got_line) {
+                if (!(bool)getline(*fpu, tp)) {
                     break;
                 }
                 if (tp[0] == '#') {
                     continue;
                 }
-                gafParser(tp, this->cached_alignments, rids);
+                gafParser(tp, this->cached_alignments);
             }
             for (auto& item : cached_alignments) {
                 item.second.index();
@@ -397,10 +389,48 @@ namespace AlignFormat {
         }
     }
 
-    GwAlignment::~GwAlignment() {
+    AlignFormat::GwAlignment::~GwAlignment() {
         if (index) hts_idx_destroy(index);
         if (header) sam_hdr_destroy(header);
         if (bam) sam_close(bam);
+        if (!cached_alignments.empty()) {
+            for (auto& item : cached_alignments) {
+                for (auto &g : item.second.data) {
+                    delete g;
+                }
+            }
+        }
+    }
+
+    void gafToAlign(AlignFormat::GAF_t* gaf, AlignFormat::Align* align) {
+
+        align->cov_start = gaf->pos;
+        align->cov_end = gaf->end;
+        align->orient_pattern = AlignFormat::Pattern::NORMAL;
+        align->left_soft_clip = 0;
+        align->right_soft_clip = 0;
+        align->y = gaf->y;
+        align->edge_type = 1;
+        align->pos = gaf->pos;
+        align->reference_end = gaf->end;
+        align->has_SA = false;
+        align->blocks = gaf->blocks;
+        std::cout << gaf->pos << " " << gaf->end << " " << gaf->y << std::endl;
+        align->delegate = bam_init1();
+        const char* qname = gaf->qname.c_str();
+        size_t l_qname = gaf->qname.size();
+        // int bam_set1(bam1_t *bam,
+        //             size_t l_qname, const char *qname,
+        //             uint16_t flag, int32_t tid, hts_pos_t pos, uint8_t mapq,
+        //             size_t n_cigar, const uint32_t *cigar,
+        //             int32_t mtid, hts_pos_t mpos, hts_pos_t isize,
+        //             size_t l_seq, const char *seq, const char *qual,
+        //             size_t l_aux);
+        int res = bam_set1(align->delegate, l_qname, qname, (uint64_t)gaf->flag, 0,
+                     (hts_pos_t)align->pos, (uint8_t)gaf->qual, 0, nullptr, 0, 0, 0, 0, nullptr, nullptr, 0);
+        if (res < 0) {
+            return;
+        }
     }
 
 }
diff --git a/src/alignment_format.h b/src/alignment_format.h
@@ -82,6 +82,7 @@ namespace AlignFormat {
         std::vector<AlignFormat::ModItem> any_mods;
 
         // Constructor
+        Align() {}
         Align(bam1_t *src) { delegate = src; }
 
         // Destructor
@@ -146,16 +147,11 @@ namespace AlignFormat {
         GAF_t
     };
 
-    class Core {
-    public:
-        int tid, pos, end, qual, flag, n_cigar, l_qseq, mtid, mpos;
-    };
-
-    class GAF_t {
+    EXPORT class GAF_t {
     public:
         std::string qname, chrom;
-        Core core;
         int qlen, qstart, qend;
+        int pos{0}, end{0}, qual, flag{0};
         char strand;  // this is strand of the query, not reference
         int y{-1};
 
@@ -168,9 +164,6 @@ namespace AlignFormat {
         }
     };
 
-    void gafParser(std::string& line,
-                  ankerl::unordered_dense::map< std::string, SuperIntervals<int, GAF_t *>>& cached_alignments);
-
     EXPORT class GwAlignment {
     public:
         AlignmentType type;
@@ -186,8 +179,9 @@ namespace AlignFormat {
         GwAlignment() = default;
         ~GwAlignment();
 
-        void open(const std::string& path, const std::string& reference, int threads, faidx_t* fai);
+        void open(const std::string& path, const std::string& reference, int threads);
     };
 
+    void gafToAlign(AlignFormat::GAF_t* gaf, AlignFormat::Align* align);
 
 }
diff --git a/src/hts_funcs.cpp b/src/hts_funcs.cpp
@@ -169,24 +169,25 @@ namespace HGW {
                                  bool coverage, std::vector<Parse::Parser> &filters, BS::thread_pool &pool,
                                  const int parse_mods_threshold) {
 
-        std::vector<AlignFormat::GAF_t*>& readQueueGAF = col.readQueueGAF;
+        std::vector<AlignFormat::GAF_t*> &readQueueGAF = col.readQueueGAF;
         if (!readQueueGAF.empty()) {
             readQueueGAF.clear();
         }
         if (!col.alignmentFile->cached_alignments.contains(region->chrom)) {
             return;
         }
         col.alignmentFile->cached_alignments[region->chrom].findOverlaps(region->start, region->end, readQueueGAF);
+        std::reverse(readQueueGAF.begin(), readQueueGAF.end());
         if (coverage) {
             std::fill(col.covArr.begin(), col.covArr.end(), 0);
             for (const auto &i : readQueueGAF) {
-//                std::cerr << i->core.pos << "-" << i->core.end << std::endl << "  ";
-//                for (auto &v : i->blocks) {
-//                    std::cerr << v.start << "," << v.end << "  ";
-//                } std::cerr << std::endl << std::endl;
                 Segs::addToCovArray(col.covArr, i->blocks, region->start, region->end);
             }
         }
+        col.readQueue.resize(readQueueGAF.size());
+        for (size_t i=0; i < readQueueGAF.size(); ++i) {
+            AlignFormat::gafToAlign(readQueueGAF[i], &col.readQueue[i]);
+        }
         col.collection_processed = false;
     }
 
diff --git a/src/plot_manager.cpp b/src/plot_manager.cpp
@@ -96,7 +96,7 @@ namespace Manager {
             hts_idx_t *idx = sam_index_load(f, fn.c_str());
             indexes.push_back(idx);
 
-            inputAlignments[i].open(fn, reference, opt.threads, fai);
+            inputAlignments[i].open(fn, reference, opt.threads);
             ++i;
         }
 
diff --git a/src/segments.cpp b/src/segments.cpp
@@ -734,7 +734,7 @@ namespace Segs {
                     pos += l;
                     break;
                 case BAM_CREF_SKIP:
-                    op = BAM_CDEL;
+//                    op = BAM_CDEL;
                     pos += l;
                     break;
                 case BAM_CSOFT_CLIP:
@@ -763,17 +763,6 @@ namespace Segs {
             self->has_SA = false;
         }
 
-//        if (sort_state > 0) {
-//            if (sort_state == Utils::SortType::HP) {
-//                uint8_t *HP_tag = bam_aux_get(self->delegate, "HP");
-//                self->sort_tag = (HP_tag != nullptr) ? (int) bam_aux2i(HP_tag) : 0;
-//            } else if (sort_state == Utils::SortType::STRAND) {
-//                self->sort_tag = (int)(src->core.flag & BAM_FREVERSE) ? 1 : 0;
-//            } else if (sort_by_pos >= 0 && ref_base != '\0') {
-//                self->sort_tag = getSortCode(src, sort_by_pos, ref_base);
-//            }
-//        }
-
         if (parse_mods_threshold > 0) {
             hts_base_mod_state* mod_state = new hts_base_mod_state;
             int res = bam_parse_basemod_gw(src, mod_state, 0);
@@ -855,6 +844,9 @@ namespace Segs {
         }
     }
 
+
+
+
     constexpr int POS_MASK = 0b0000111;
 
     void setAlignSortCode(AlignFormat::Align &a, Utils::SortType sort_state, int target_pos, char ref_base) {

Original file line number	Diff line number	Diff line change
`@@ -169,24 +169,25 @@ namespace HGW {`
`169`	`169`	`bool coverage, std::vector<Parse::Parser> &filters, BS::thread_pool &pool,`
`170`	`170`	`const int parse_mods_threshold) {`
`171`	`171`
`172`		`- std::vector<AlignFormat::GAF_t*>& readQueueGAF = col.readQueueGAF;`
	`172`	`+ std::vector<AlignFormat::GAF_t*> &readQueueGAF = col.readQueueGAF;`
`173`	`173`	`if (!readQueueGAF.empty()) {`
`174`	`174`	`readQueueGAF.clear();`
`175`	`175`	`}`
`176`	`176`	`if (!col.alignmentFile->cached_alignments.contains(region->chrom)) {`
`177`	`177`	`return;`
`178`	`178`	`}`
`179`	`179`	`col.alignmentFile->cached_alignments[region->chrom].findOverlaps(region->start, region->end, readQueueGAF);`
	`180`	`+ std::reverse(readQueueGAF.begin(), readQueueGAF.end());`
`180`	`181`	`if (coverage) {`
`181`	`182`	`std::fill(col.covArr.begin(), col.covArr.end(), 0);`
`182`	`183`	`for (const auto &i : readQueueGAF) {`
`183`		`-// std::cerr << i->core.pos << "-" << i->core.end << std::endl << " ";`
`184`		`-// for (auto &v : i->blocks) {`
`185`		`-// std::cerr << v.start << "," << v.end << " ";`
`186`		`-// } std::cerr << std::endl << std::endl;`
`187`	`184`	`Segs::addToCovArray(col.covArr, i->blocks, region->start, region->end);`
`188`	`185`	`}`
`189`	`186`	`}`
	`187`	`+ col.readQueue.resize(readQueueGAF.size());`
	`188`	`+ for (size_t i=0; i < readQueueGAF.size(); ++i) {`
	`189`	`+ AlignFormat::gafToAlign(readQueueGAF[i], &col.readQueue[i]);`
	`190`	`+ }`
`190`	`191`	`col.collection_processed = false;`
`191`	`192`	`}`
`192`	`193`
Original file line number	Diff line number	Diff line change
`@@ -96,7 +96,7 @@ namespace Manager {`
`96`	`96`	`hts_idx_t *idx = sam_index_load(f, fn.c_str());`
`97`	`97`	`indexes.push_back(idx);`
`98`	`98`
`99`		`- inputAlignments[i].open(fn, reference, opt.threads, fai);`
	`99`	`+ inputAlignments[i].open(fn, reference, opt.threads);`
`100`	`100`	`++i;`
`101`	`101`	`}`
`102`	`102`