Implementation of soft-clipping and CIGAR computation

include/ProgOpts.hpp

@@ -123,8 +123,13 @@ class AlignmentOpts{
   bool recoverOrphans{false};
   bool mimicBt2Default{false};
   bool mimicBt2Strict{false};
-  bool allowOverhangSoftclip{false};
   bool allowSoftclip{false};
+  bool computeCIGAR{false};
+  bool debug{false};
+  bool end2end{true};
+  double maxSoftclipFractionGeneral{0.2};
+  double maxSoftclipFractionOverhang{0.2};
+  uint32_t endBonus{5};
   bool useAlignmentCache{true};
   uint32_t alignmentStreamLimit{10000};
   double preMergeChainSubThresh{0.9};

include/PuffAligner.hpp

@@ -83,10 +83,10 @@ class PuffAligner {
   PuffAligner(PuffAligner&& other) = delete;
   PuffAligner& operator=(PuffAligner&& other) = delete;
 
-  int32_t calculateAlignments(std::string& rl, std::string& rr, pufferfish::util::JointMems& jointHit, HitCounters& hctr, bool isMultimapping, bool verbose);
-  int32_t calculateAlignments(std::string& read, pufferfish::util::JointMems& jointHit, HitCounters& hctr, bool isMultimapping, bool verbose);
+  int32_t calculateAlignments(std::string& rl, std::string& rr, pufferfish::util::JointMems& jointHit, HitCounters& hctr, bool isMultimapping, std::shared_ptr<spdlog::logger>& logger_, bool verbose);
+  int32_t calculateAlignments(std::string& read, pufferfish::util::JointMems& jointHit, HitCounters& hctr, bool isMultimapping, std::shared_ptr<spdlog::logger>& logger_, bool verbose);
 
-  bool alignRead(std::string& read, std::string& read_rc, const std::vector<pufferfish::util::MemInfo>& mems, uint64_t queryChainHash, bool perfectChain, bool isFw, size_t tid, AlnCacheMap& alnCache, HitCounters& hctr, AlignmentResult& arOut, bool verbose);
+  bool alignRead(std::string& read, std::string& read_rc, const std::vector<pufferfish::util::MemInfo>& mems, uint64_t queryChainHash, bool perfectChain, bool isFw, size_t tid, AlnCacheMap& alnCache, HitCounters& hctr, AlignmentResult& arOut, std::shared_ptr<spdlog::logger>& logger_, bool verbose);
 
   int32_t align_ungapped(const char* const query, const char* const target, int len);
   bool recoverSingleOrphan(std::string& rl, std::string& rr, pufferfish::util::MemCluster& clust, std::vector<pufferfish::util::MemCluster> &recoveredMemClusters, uint32_t tid, bool anchorIsLeft, bool verbose);

include/Util.hpp

@@ -397,89 +397,42 @@ Compile-time selection between list-like and map-like printing.
       struct CIGARGenerator {
         // TODO: @fataltes --- think about just replacing this
         // with CIGAR Op class or some such.
-        std::vector<uint32_t> cigar_counts;
+        std::vector<int32_t> cigar_counts;
         std::string cigar_types;
         int32_t begin_softclip_len{0}, end_softclip_len{0};
-        bool beginOverhang{false}, endOverhang{false};
-        bool allowOverhangSoftClip{false};
 
         void clear() {
           cigar_counts.clear(); cigar_types.clear();
           begin_softclip_len = end_softclip_len = 0;
-          beginOverhang = endOverhang = false;
         }
 
-        void add_item(uint32_t count, char type) {
+        void add_item(int32_t count, char type) {
           cigar_counts.push_back(count);
           cigar_types.push_back(type);
         }
 
-        void get_approx_cigar(int32_t readLen, std::string& cigar) {
-          if (begin_softclip_len > 0) {
-            cigar += std::to_string(begin_softclip_len);
-            cigar += beginOverhang and allowOverhangSoftClip ? "I" : "S";
-          }
-          cigar += std::to_string(readLen - (begin_softclip_len + end_softclip_len));
-          cigar += "M";
-          if (end_softclip_len > 0) {
-            cigar += std::to_string(end_softclip_len);
-            cigar += endOverhang and allowOverhangSoftClip ? "I" : "S";
-          }
-        }
+        std::string get_cigar() {
+          if (cigar_counts.size() == 0) return "*";
+          if (cigar_counts.size() != cigar_types.size()) return "!";
 
-        std::string get_cigar(uint32_t readLen, bool &cigar_fixed) {
-          cigar_fixed = false;
           std::string cigar = "";
-          if (cigar_counts.size() != cigar_types.size() or cigar_counts.size() == 0) {
-            return "!";
-          }
-          if (cigar_counts.size() == 0) {
-            return cigar;
-          }
-
-          uint32_t cigar_length = 0;
-          uint32_t count = cigar_counts[0];
-
-          if (cigar_counts.size() == 1) {
-            if (count != readLen) {
-              count = readLen;
-              cigar_fixed = true;
-            }
-            cigar += std::to_string(count);
-            cigar += cigar_types[0];
-            return cigar;
-          }
-
-          char type = cigar_types[0];
-          if (type == 'I' or type == 'M')
-            cigar_length += count;
+          int32_t last_count = cigar_counts[0];
+          char last_type = cigar_types[0];
           for (size_t i = 1; i < cigar_counts.size(); i++) {
-            if (cigar_types[i] == 'I' or cigar_types[i] == 'M')
-              cigar_length += cigar_counts[i];
-            if (type == cigar_types[i]) {
-              count += cigar_counts[i];
-            } else {
-              cigar += std::to_string(count);
-              cigar += type;
-              count = cigar_counts[i];
-              type = cigar_types[i];
+            if (cigar_types[i] == last_type) {
+              last_count += cigar_counts[i];
             }
-            if (i == cigar_counts.size() - 1) {
-              cigar += std::to_string(count);
-              cigar += type;
-              if (cigar_length < readLen) {
-                cigar_fixed = true;
-                count = readLen - cigar_length;
-                cigar += std::to_string(count);
-                cigar += 'I';
-              } else if (cigar_length > readLen) {
-                cigar_fixed = true;
-                count = cigar_length - readLen;
-                cigar += std::to_string(count);
-                cigar += 'I';
-              }
+            else {
+              cigar += std::to_string(last_count);
+              cigar += last_type;
+              // reset
+              last_count = cigar_counts[i];
+              last_type = cigar_types[i];
             }
           }
+          // add last
+          cigar += std::to_string(last_count);
+          cigar += last_type;
           return cigar;
         }
       };
@@ -857,8 +810,6 @@ Compile-time selection between list-like and map-like printing.
     };
     */
 
-      enum class PuffAlignmentMode : uint8_t { SCORE_ONLY, APPROXIMATE_CIGAR,  EXACT_CIGAR};
-
       struct AlignmentConfig {
         int32_t refExtendLength{20};
         bool fullAlignment{false};
@@ -869,12 +820,14 @@ Compile-time selection between list-like and map-like printing.
         double minScoreFraction{0.0};
         bool mimicBT2{false};
         bool mimicBT2Strict{false};
-        bool allowOverhangSoftclip{false};
         bool allowSoftclip{false};
+        bool computeCIGAR{false};
+        bool end2end{true};
+        double maxSoftclipFractionGeneral{0.2};
+        double maxSoftclipFractionOverhang{0.2};
         bool useAlignmentCache{true};
         bool noDovetail{false};
         uint32_t maxFragmentLength{1000};
-        PuffAlignmentMode alignmentMode{PuffAlignmentMode::SCORE_ONLY};
         bool bestStrata{false};
         bool decoyPresent{false};
       };
@@ -1342,15 +1295,17 @@ Compile-time selection between list-like and map-like printing.
 
         struct AlignmentResult {
             AlignmentResult(bool isFwIn, int32_t scoreIn, std::string cigarIn, uint32_t openGapLenIn, 
-                            uint16_t softclip_start_in ) :
-                    isFw(isFwIn), score(scoreIn), cigar(cigarIn), openGapLen(openGapLenIn), softclip_start(softclip_start_in) {}
+                            uint16_t softclip_start_in, uint16_t softclip_end_in ) :
+                    isFw(isFwIn), score(scoreIn), cigar(cigarIn), openGapLen(openGapLenIn), softclip_start(softclip_start_in),
+                    softclip_end(softclip_end_in) {}
 
-            AlignmentResult() : isFw(true), score(0), cigar(""), openGapLen(0), softclip_start(0) {}
+            AlignmentResult() : isFw(true), score(0), cigar(""), openGapLen(0), softclip_start(0), softclip_end(0) {}
             bool isFw;
             int32_t score;
             std::string cigar;
             uint32_t openGapLen;
             uint16_t softclip_start{0};
+            uint16_t softclip_end{0};
         };
 
       void joinReadsAndFilterSingle( pufferfish::util::CachedVectorMap<size_t, std::vector<pufferfish::util::MemCluster>, std::hash<size_t>>& leftMemClusters,

include/ksw2pp/KSW2Aligner.hpp

@@ -38,7 +38,7 @@ struct KSW2Config {
   int dropoff = -1;
   int flag = 0;
   int alphabetSize = 5;
-  int end_bonus = 10;
+  int end_bonus = 5;
   KSW2AlignmentType atype = KSW2AlignmentType::GLOBAL;
 };
 
@@ -66,15 +66,15 @@ class KSW2Aligner {
 
   int operator()(const char* const queryOriginal, const int queryLength,
                  const char* const targetOriginal, const int targetLength,
-                 ksw_extz_t* ez, int cutoff, EnumToType<KSW2AlignmentType::EXTENSION>);
+                 ksw_extz_t* ez, int cutoff, int maxSoftclip, EnumToType<KSW2AlignmentType::EXTENSION>);
 
   int operator()(const uint8_t* const queryOriginal, const int queryLength,
                  const uint8_t* const targetOriginal, const int targetLength,
                  ksw_extz_t* ez, EnumToType<KSW2AlignmentType::GLOBAL>);
 
   int operator()(const uint8_t* const queryOriginal, const int queryLength,
                  const uint8_t* const targetOriginal, const int targetLength,
-                 ksw_extz_t* ez, int cutoff, EnumToType<KSW2AlignmentType::EXTENSION>);
+                 ksw_extz_t* ez, int cutoff, int maxSoftclip, EnumToType<KSW2AlignmentType::EXTENSION>);
 
   /**
    * Variants of the operator that do not require an output

include/ksw2pp/ksw2.h

@@ -57,11 +57,11 @@ void ksw_extz(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t
 			  int8_t q, int8_t e, int w, int zdrop, int flag, ksw_extz_t *ez);
 
 void ksw_extz2_sse(/*unsigned int simd, */void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez, int cutoff);
+				   int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez, int cutoff, int max_sc_len);
 void ksw_extz2_sse41(/*unsigned int simd, */void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-           int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez, int cutoff);
+           int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez, int cutoff, int max_sc_len);
 void ksw_extz2_sse2(/*unsigned int simd, */void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-           int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez, int cutoff);
+           int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez, int cutoff, int max_sc_len);
 
 
 void ksw_extd(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,