Use ksw2 instead of parasail + Remove genotyper

Author: ldenti

CMakeLists.txt

@@ -17,7 +17,7 @@ if(DEFINED CONDAPREFIX)
 endif()
 
 project(SVDSS VERSION 2.0.0)
-add_executable(SVDSS assembler.cpp bam.cpp caller.cpp chromosomes.cpp clipper.cpp clusterer.cpp config.cpp ping_pong.cpp sfs.cpp smoother.cpp sv.cpp genotyper.cpp main.cpp)
+add_executable(SVDSS assembler.cpp bam.cpp caller.cpp chromosomes.cpp clipper.cpp clusterer.cpp config.cpp ping_pong.cpp sfs.cpp smoother.cpp sv.cpp main.cpp)
 set_target_properties(SVDSS PROPERTIES RUNTIME_OUTPUT_DIRECTORY "${CMAKE_SOURCE_DIR}")
 target_compile_options(SVDSS PRIVATE -Wall -Wextra)
 
@@ -107,22 +107,19 @@ SET(ABPOA_INCLUDE_DIR ${ABPOA_SOURCE_DIR}/include)
 # add_library(ABPOA STATIC IMPORTED)
 # set_target_properties(ABPOA PROPERTIES IMPORTED_LOCATION ${ABPOA_SOURCE_DIR}/lib/libabpoa.a)
 
-# parasail
-###########
-message(STATUS "parasail will be built from source")
-ExternalProject_Add(parasail
-  GIT_REPOSITORY https://github.com/jeffdaily/parasail.git
-  GIT_TAG 600fb26151ff19899ee39a214972dcf2b9b11ed7 # v2.6.2
-  UPDATE_COMMAND ""
-  BUILD_IN_SOURCE 1
-  CMAKE_ARGS -DBUILD_SHARED_LIBS=OFF
-  INSTALL_COMMAND ""
-  )
-ExternalProject_Get_Property(parasail SOURCE_DIR)
-SET(PARASAIL_SOURCE_DIR ${SOURCE_DIR})
-SET(PARASAIL_INCLUDE_DIR ${PARASAIL_SOURCE_DIR})
-add_library(PARASAIL STATIC IMPORTED)
-set_target_properties(PARASAIL PROPERTIES IMPORTED_LOCATION ${PARASAIL_SOURCE_DIR}/libparasail.a)
+# ksw2
+message(STATUS "ksw2 will be built from source")
+ExternalProject_Add(ksw2
+    GIT_REPOSITORY https://github.com/lh3/ksw2.git
+    GIT_TAG ""
+    BUILD_IN_SOURCE 1
+    UPDATE_COMMAND ""
+    CONFIGURE_COMMAND ""
+    BUILD_COMMAND make
+    INSTALL_COMMAND ""
+)
+ExternalProject_Get_Property(ksw2 SOURCE_DIR)
+SET(KSW2_INCLUDE_DIR ${SOURCE_DIR})
 
 # rapidfuzz-cpp
 ################
@@ -182,7 +179,7 @@ FetchContent_MakeAvailable(spdlog)
 ### final setup (includes and libraries) ###
 ############################################
 
-add_dependencies(SVDSS abpoa parasail rapidfuzz intervaltree ropebwt spdlog)
+add_dependencies(SVDSS abpoa rapidfuzz intervaltree ropebwt spdlog ksw2)
 if(NOT DEFINED CONDAPREFIX)
   # we need to build them
   add_dependencies(SVDSS deflate htslib)
@@ -191,7 +188,7 @@ endif()
 target_include_directories(SVDSS
   PRIVATE ${HTS_INCLUDE_DIR}
   PRIVATE ${ABPOA_INCLUDE_DIR}
-  PRIVATE ${PARASAIL_INCLUDE_DIR}
+  PRIVATE ${KSW2_INCLUDE_DIR}
   PRIVATE ${RAPIDFUZZ_INCLUDE_DIR}
   PRIVATE ${INTERVALTREE_INCLUDE_DIR}
   PRIVATE ${ROPEBWT_INCLUDE_DIR}
@@ -220,7 +217,8 @@ endif()
 target_link_libraries(SVDSS
   PUBLIC rapidfuzz::rapidfuzz
   PUBLIC ${BINARY_DIR}/lib/libabpoa.a
-  PUBLIC PARASAIL
+  PUBLIC ${KSW2_INCLUDE_DIR}/ksw2_extz2_sse.o
+  PUBLIC ${KSW2_INCLUDE_DIR}/ksw2_extd2_sse.o
   PUBLIC ROPEBWT
   PUBLIC spdlog::spdlog
   PUBLIC z

bam.hpp

@@ -1,6 +1,7 @@
 #ifndef BAM_HPP
 #define BAM_HPP
 
+#include <cstring>
 #include <ctype.h>
 #include <iostream>
 #include <iterator>
@@ -47,21 +48,24 @@ struct CIGAR {
     }
   }
 
-  void parse_cigar(char *cigar) {
+  void parse_cigar(const char *cigar) {
     int b = 0;
     int lc = strlen(cigar);
+    char *cigar_cpy = (char *)malloc(lc + 1);
+    strncpy(cigar_cpy, cigar, lc + 1);
     for (int i = 0; i < lc; i++) {
-      if (isdigit(cigar[i]))
+      if (isdigit(cigar_cpy[i]))
         continue;
       else {
-        char type = cigar[i];
-        cigar[i] = '\0';
-        int l = stoi(string(cigar + b));
+        char type = cigar_cpy[i];
+        cigar_cpy[i] = '\0';
+        int l = stoi(string(cigar_cpy + b));
         ops.push_back(make_pair(l, type));
-        cigar[i] = type;
+        cigar_cpy[i] = type;
         b = i + 1;
       }
     }
+    free(cigar_cpy);
   }
 
   CIGAR(char *cigar, int score_, int start_ = 0) {

caller.cpp

@@ -52,6 +52,8 @@ void Caller::run() {
     for (const SV &sv : clipped_svs)
       cout << sv << endl;
   }
+
+  destroy_chromosomes();
 }
 
 void Caller::write_vcf() {
@@ -94,7 +96,7 @@ vector<Cluster> Caller::split_cluster_by_len(const Cluster &cluster) {
   return subclusters;
 }
 
-/* Split cluster in subclusters */
+// Split cluster in subclusters
 vector<Cluster> Caller::split_cluster(const Cluster &cluster) {
   // Step 1: split cluster by haplotype tag
   Cluster cluster_0 = cluster;
@@ -256,13 +258,16 @@ string Caller::run_poa(const vector<string> &seqs) {
   uint n_seqs = seqs.size();
   abpoa_t *ab = abpoa_init();
   abpoa_para_t *abpt = abpoa_init_para();
-  abpt->disable_seeding = 1;
   abpt->align_mode = 0; // global
+  abpt->disable_seeding = 1;
+  abpt->progressive_poa = 0;
+  abpt->amb_strand = 0;
   abpt->out_msa = 0;
   abpt->out_cons = 1;
   abpt->out_gfa = 0;
   // abpt->is_diploid = 1; // TODO: maybe this works now
-  abpt->progressive_poa = 0;
+  // abpt->max_n_cons = 2;
+  // abpt->min_freq = 0.25;
   abpoa_post_set_para(abpt);
 
   // abpt->match = 2;      // match score
@@ -285,7 +290,8 @@ string Caller::run_poa(const vector<string> &seqs) {
 
   abpoa_msa(ab, abpt, n_seqs, NULL, seq_lens, bseqs, NULL, NULL);
   abpoa_cons_t *abc = ab->abc;
-  string cons = "";
+  string cons = ""; // XXX: we may avoid converting to ACGT here since we need
+                    // to reconvert back for ksw2
   if (abc->n_cons > 0)
     for (int j = 0; j < abc->cons_len[0]; ++j)
       cons += "ACGTN"[abc->cons_base[0][j]];
@@ -301,10 +307,8 @@ string Caller::run_poa(const vector<string> &seqs) {
   return cons;
 }
 
-/* Call SVs by POA+realignment */
+// Call SVs by POA+realignment
 void Caller::pcall(const vector<Cluster> &clusters) {
-  // vector<Genotyper> genotypers(config->threads);
-  vector<string> gt_strings = {"0/0", "0/1", "1/0", "1/1"};
 #pragma omp parallel for num_threads(config->threads) schedule(static, 1)
   for (size_t i = 0; i < clusters.size(); i++) {
     int t = omp_get_thread_num();
@@ -313,24 +317,6 @@ void Caller::pcall(const vector<Cluster> &clusters) {
       continue;
     string chrom = cluster.chrom;
 
-    // genotypers.at(t).posterior_sv_genotype_give_reads(cluster.reads);
-    // vector<double> gts = genotypers.at(t).get_posterior_sv_genotype();
-    Genotyper gtyper;
-    gtyper.posterior_sv_genotype_give_reads(cluster.reads);
-    vector<double> gts = gtyper.get_posterior_sv_genotype();
-    auto max_gt = max_element(gts.begin(), gts.end());
-    double gtq = *max_gt;
-    if (gtq < 0 || gtq > 1) {
-      cerr << chrom << ":" << cluster.s << "-" << cluster.e << endl;
-      for (const auto &tpl : cluster.reads)
-        cerr << get<0>(tpl) << ":" << get<1>(tpl) << " ";
-      cerr << endl;
-      for (int i = 0; i < 4; ++i)
-        cerr << gt_strings[i] << " - " << gts[i] << endl;
-    }
-    string gt = gt_strings[distance(gts.begin(), max_gt)];
-    if (config->noref && gt.compare("0/0") == 0)
-      continue;
     const vector<Cluster> &subclusters = split_cluster(cluster);
 
     // Calling from one or two clusters
@@ -342,26 +328,41 @@ void Caller::pcall(const vector<Cluster> &clusters) {
 
       string ref = string(chromosome_seqs[chrom] + cl.s, cl.e - cl.s + 1);
       string consensus = run_poa(cl.get_seqs());
-      parasail_result_t *result = NULL;
-      result = parasail_nw_trace_striped_16(consensus.c_str(), consensus.size(),
-                                            ref.c_str(), ref.size(), 10, 1,
-                                            &parasail_nuc44);
-      parasail_cigar_t *pcigar =
-          parasail_result_get_cigar(result, consensus.c_str(), consensus.size(),
-                                    ref.c_str(), ref.size(), NULL);
-      char *cigar_str = parasail_cigar_decode(pcigar);
-      int score = result->score;
-      parasail_cigar_free(pcigar);
-      parasail_result_free(result);
+
+      // ksw2 stuff - TODO: move to a separate function
+      int sc_mch = 1, sc_mis = -9, gapo = 16, gape = 2, gapo2 = 41, gape2 = 1;
+      int8_t a = (int8_t)sc_mch,
+             b = sc_mis < 0 ? (int8_t)sc_mis : -(int8_t)sc_mis; // a>0 and b<0
+      int8_t mat[25] = {a, b, b, b, 0, b, a, b, b, 0, b, b, a,
+                        b, 0, b, b, b, a, 0, 0, 0, 0, 0, 0};
+      uint tl = ref.size(), ql = consensus.size();
+      uint8_t *ts = (uint8_t *)malloc(tl);
+      uint8_t *qs = (uint8_t *)malloc(ql);
+      for (i = 0; i < tl; ++i)
+        ts[i] = _char26_table[(uint8_t)ref[i]]; // encode to 0/1/2/3
+      for (i = 0; i < ql; ++i)
+        qs[i] = _char26_table[(uint8_t)consensus[i]];
+
+      ksw_extz_t ez;
+      memset(&ez, 0, sizeof(ksw_extz_t));
+      ksw_extd2_sse(0, ql, qs, tl, ts, 5, mat, gapo, gape, gapo2, gape2, -1, -1,
+                    -1, 0, &ez);
+
+      int score = ez.score;
+      string cigar_str = "";
+      for (int i = 0; i < ez.n_cigar; ++i) {
+        cigar_str += to_string(ez.cigar[i] >> 4) + "MID"[ez.cigar[i] & 0xf];
+      }
       _p_alignments[t].push_back(
           Consensus(consensus, cigar_str, chrom, cl.s, cl.e));
+
       // -- Extracting SVs
       uint rpos = cl.s; // position on reference
       uint cpos = 0;    // position on consensus
       CIGAR cigar;
-      cigar.parse_cigar(cigar_str);
+      cigar.parse_cigar(cigar_str.c_str());
       int nv = 0;
-      for (const auto cigar_pair : cigar.ops) {
+      for (const auto &cigar_pair : cigar.ops) {
         uint l = cigar_pair.first;
         char op = cigar_pair.second;
         if (op == '=' || op == 'M') {
@@ -394,7 +395,7 @@ void Caller::pcall(const vector<Cluster> &clusters) {
       }
       for (size_t v = 0; v < _svs.size(); v++) {
         _svs[v].ngaps = nv;
-        _svs[v].set_gt(gt, max(0, (int)(gtq * 100)));
+        _svs[v].set_gt("./.", 100);
         _svs[v].set_cov(cl.cov, cl.cov0, cl.cov1, cl.cov2);
         _svs[v].set_rvec(cluster.reads);
       }
@@ -404,7 +405,7 @@ void Caller::pcall(const vector<Cluster> &clusters) {
   }
 }
 
-/* Clean same SV reported twice */
+// Clean same SV reported twice
 void Caller::clean_dups() {
   vector<SV> _svs;
   string last_chrom = "";

caller.hpp

@@ -7,16 +7,14 @@
 #include <map>
 
 #include <abpoa.h>
-#include <parasail.h>
-#include <parasail/matrices/nuc44.h>
+#include <ksw2.h>
 #include <rapidfuzz/fuzz.hpp>
 #include <spdlog/spdlog.h>
 
 #include "bam.hpp"
 #include "chromosomes.hpp"
 #include "clipper.hpp"
 #include "clusterer.hpp"
-#include "genotyper.hpp"
 #include "config.hpp"
 #include "sfs.hpp"
 #include "sv.hpp"

clipper.cpp

@@ -25,7 +25,7 @@ vector<Clip> Clipper::combine(const vector<Clip> &clips) {
   }
 // we then merge
 #pragma omp parallel for num_threads(threads) schedule(static, 1)
-  for (int i = 0; i < chromosomes.size(); i++) {
+  for (uint i = 0; i < chromosomes.size(); i++) {
     int t = i % threads;
     const string &chrom = chromosomes[i];
     for (auto it = clips_dict[chrom].begin(); it != clips_dict[chrom].end();
@@ -36,16 +36,14 @@ vector<Clip> Clipper::combine(const vector<Clip> &clips) {
           max_l = c.l;
         }
       }
-      Clip clip = Clip("", chrom, it->first, max_l,
-      it->second.front().starting,
+      Clip clip = Clip("", chrom, it->first, max_l, it->second.front().starting,
                        it->second.size());
       _p_combined_clips[t].push_back(clip);
     }
   }
   vector<Clip> combined_clips;
   for (int i = 0; i < threads; i++) {
-    combined_clips.insert(combined_clips.begin(),
-    _p_combined_clips[i].begin(),
+    combined_clips.insert(combined_clips.begin(), _p_combined_clips[i].begin(),
                           _p_combined_clips[i].end());
   }
   return combined_clips;
@@ -93,26 +91,20 @@ vector<Clip> Clipper::filter_tooclose_clips(const vector<Clip> &clips,
                                             interval_tree_t<int> &vartree) {
   vector<Clip> fclips;
   for (const Clip &c : clips) {
-    if (vartree.overlap_find({c.p, c.p + 1}) == end(vartree)) {
+    if (vartree.overlap_find({(int)c.p, (int)c.p + 1}) == end(vartree)) {
       fclips.push_back(c);
     }
   }
   return fclips;
 }
 
 // find smallest right that is larger than query
-int binary_search(const vector<Clip> &clips, int begin, int end,
+int binary_search(const vector<Clip> &clips, uint begin, uint end,
                   const Clip &query) {
-  // for (int i = 0; i < clips.size(); i++) {
-  //     if (query.p < clips[i].p) {
-  //         return i ;
-  //     }
-  // }
-  // return -1 ;
   if (begin > end || begin >= clips.size()) {
     return -1;
   }
-  int m = (begin + end) / 2;
+  uint m = (begin + end) / 2;
   if (clips[m].p == query.p) {
     if (m + 1 < clips.size()) {
       return m + 1;
@@ -130,7 +122,8 @@ int binary_search(const vector<Clip> &clips, int begin, int end,
 }
 
 void Clipper::call(int threads, interval_tree_t<int> &vartree) {
-  // lprint({"Predicting SVS from", to_string(clips.size()), "clipped SFS on", to_string(threads), "threads.."});
+  // lprint({"Predicting SVS from", to_string(clips.size()), "clipped SFS on",
+  // to_string(threads), "threads.."});
   vector<Clip> rclips;
   vector<Clip> lclips;
   for (const Clip &clip : clips) {
@@ -169,7 +162,7 @@ void Clipper::call(int threads, interval_tree_t<int> &vartree) {
   }
   // lprint({"Predicting insertions.."});
 #pragma omp parallel for num_threads(threads) schedule(static, 1)
-  for (int i = 0; i < lclips.size(); i++) {
+  for (uint i = 0; i < lclips.size(); i++) {
     const Clip &lc = lclips[i];
     int t = omp_get_thread_num();
     string chrom = lc.chrom;
@@ -194,7 +187,7 @@ void Clipper::call(int threads, interval_tree_t<int> &vartree) {
   }
   // lprint({"Predicting deletions.."});
 #pragma omp parallel for num_threads(threads) schedule(static, 1)
-  for (int i = 0; i < rclips.size(); i++) {
+  for (uint i = 0; i < rclips.size(); i++) {
     const Clip &rc = rclips[i];
     int t = omp_get_thread_num();
     string chrom = rc.chrom;

clusterer.cpp

@@ -334,7 +334,7 @@ void Clusterer::extend_alignment(bam1_t *aln, int index) {
       merged_extended_sfs.push_back(local_extended_sfs.at(i));
     }
   }
-  for (const auto mes : merged_extended_sfs)
+  for (const auto &mes : merged_extended_sfs)
     _p_extended_sfs[index].push_back(mes);
 
   if (lclip.second > 0)