Merge pull request #272 from waveygang/smooth-chain

Smooth chain
waveygang · Sep 20, 2024 · b731e41 · b731e41
2 parents 11c7a80 + 240a0d3
commit b731e41
Show file tree

Hide file tree

Showing 6 changed files with 352 additions and 287 deletions.
diff --git a/.github/workflows/test_on_push.yml b/.github/workflows/test_on_push.yml
@@ -53,8 +53,12 @@ jobs:
           override: true
       - name: Install wgatools
         run: cargo install --git https://github.com/wjwei-handsome/wgatools.git
+      - name: Install wgatools
+        run: cargo install --git https://github.com/ekg/pafcheck.git
       - name: Run wfmash and generate PAF
-        run: build/bin/wfmash -t 8 -n 1 -k 19 -s 5000 -p 90 -c 30k -P 50k -T SGDref -Q S288C -Y '#' data/scerevisiae8.fa.gz data/scerevisiae8.fa.gz > test.paf
+        run: build/bin/wfmash -t 8 -T SGDref -Q S288C -Y '#' data/scerevisiae8.fa.gz > test.paf
+      - name: check PAF coordinates and extended CIGAR validity
+        run: pafcheck --query-fasta data/scerevisiae8.fa.gz --paf test.paf
       - name: Convert PAF to MAF using wgatools
         run: wgatools paf2maf --target data/scerevisiae8.fa.gz --query data/scerevisiae8.fa.gz test.paf > test.maf
       - name: Check if MAF file is not empty

diff --git a/src/common/wflign/src/wflign_patch.cpp b/src/common/wflign/src/wflign_patch.cpp
@@ -1397,17 +1397,17 @@ void write_merged_alignment(
                     query_pos = query_length;
                     target_pos = target_length;
 
-                    // Add safety checks first
-                    uint64_t new_query_end = query_offset + query_length;
-                    uint64_t new_target_end = target_offset + target_length + actual_extension;
+                    // Calculate new ends relative to the segment being aligned
+                    uint64_t new_query_end = query_length;
+                    uint64_t new_target_end = target_length + actual_extension;
 
-                    if (new_query_end > query_total_length || new_target_end > target_total_length) {
+                    if (query_offset + new_query_end > query_total_length || target_offset + new_target_end > target_total_length) {
                         std::cerr << "[wfmash::patch] Warning: Alignment extends beyond sequence bounds. Truncating." << std::endl;
                     }
 
-                    // Adjust query_end and target_end, ensuring we don't exceed the total lengths
-                    query_end = std::min(new_query_end, query_total_length);
-                    target_end = std::min(new_target_end, target_total_length);
+                    // Adjust query_end and target_end, ensuring we don't exceed the segment lengths
+                    query_end = std::min(new_query_end, query_length);
+                    target_end = std::min(new_target_end, target_length + actual_extension);
                 }
             }
 

diff --git a/src/interface/parse_args.hpp b/src/interface/parse_args.hpp
@@ -237,8 +237,7 @@ void parse_args(int argc,
 
     // If there are no queries, go in all-vs-all mode with the sequences specified in `target_sequence_file`
     if (target_sequence_file && map_parameters.querySequences.empty()) {
-        map_parameters.skip_self = true;
-        std::cerr << "[mashmap] Skipping self mappings for single file all-vs-all mapping." << std::endl;
+        std::cerr << "[mashmap] Performing all-vs-all mapping including self mappings." << std::endl;
         map_parameters.querySequences.push_back(map_parameters.refSequences.back());
         align_parameters.querySequences.push_back(align_parameters.refSequences.back());
     }

diff --git a/src/map/include/base_types.hpp b/src/map/include/base_types.hpp
@@ -157,10 +157,6 @@ namespace skch
     seqno_t refSeqId;                                   //internal sequence id of the reference contig
     seqno_t querySeqId;                                 //internal sequence id of the query sequence
     offset_t blockLength;                                    //the block length of the mapping
-    offset_t blockRefStartPos;
-    offset_t blockRefEndPos;
-    offset_t blockQueryStartPos;
-    offset_t blockQueryEndPos;
     float blockNucIdentity;
 
     float nucIdentity;                                  //calculated identity
@@ -178,6 +174,8 @@ namespace skch
     uint8_t discard;                                    // set to 1 for deletion
     bool overlapped;                                    // set to true if this mapping is overlapped with another mapping
     bool selfMapFilter;                                 // set to true if a long-to-short mapping in all-vs-all mode (we report short as the query)
+    double chainPairScore;                              // best score for potential chain pair
+    int64_t chainPairId;                                // best partner mapping for potential chain pair
 
     offset_t qlen() {                                   //length of this mapping on query axis
       return queryEndPos - queryStartPos + 1;