update v0.10

Assembly.cpp

@@ -12,6 +12,7 @@
 #include "kthread.h"
 
 void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain);
+void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list);
 
 All_reads R_INF;
 
@@ -343,6 +344,64 @@ long long push_final_overlaps(ma_hit_t_alloc* paf, ma_hit_t_alloc* reverse_paf_l
 
 
 
+            ///for overlap_list, the x_strand of all overlaps are 0, so the tmp.rev is the same as the y_strand
+            tmp.rev = overlap_list->list[i].y_pos_strand;
+
+
+            /**********************target***************************/
+            tmp.tn = overlap_list->list[i].y_id;
+            if(tmp.rev == 1)
+            {
+                long long y_readLen = R_INF.read_length[overlap_list->list[i].y_id];
+                tmp.ts = y_readLen - overlap_list->list[i].y_pos_e - 1;
+                tmp.te = y_readLen - overlap_list->list[i].y_pos_s - 1;
+            }
+            else
+            {
+                tmp.ts = overlap_list->list[i].y_pos_s;
+                tmp.te = overlap_list->list[i].y_pos_e;
+            }
+            ///the end pos is open
+            tmp.te++;
+            /**********************target***************************/
+
+            tmp.bl = R_INF.read_length[overlap_list->list[i].y_id];
+            tmp.ml = overlap_list->list[i].strong;
+            tmp.no_l_indel = overlap_list->list[i].without_large_indel;
+
+            tmp.el = overlap_list->list[i].shared_seed;
+
+            add_ma_hit_t_alloc(paf, &tmp);
+        }
+    }
+
+    return available_overlaps;
+}
+
+
+
+long long push_final_overlaps_increment(ma_hit_t_alloc* paf, ma_hit_t_alloc* reverse_paf_list, overlap_region_alloc* overlap_list, int flag)
+{
+    long long i = 0;
+    long long available_overlaps = paf->length;
+    ma_hit_t tmp;
+    ///clear_ma_hit_t_alloc(paf); // paf has been preallocated, so we don't need preallocation
+    for (i = 0; i < (long long)overlap_list->length; i++)
+    {
+        if (overlap_list->list[i].is_match == flag)
+        {
+            available_overlaps++;
+            /**********************query***************************/
+            //the interval of overlap is half-open [start, end) 
+            tmp.qns = overlap_list->list[i].x_id;
+            tmp.qns = tmp.qns << 32;
+            tmp.qns = tmp.qns | (uint64_t)(overlap_list->list[i].x_pos_s);
+            ///the end pos is open
+            tmp.qe = overlap_list->list[i].x_pos_e + 1;
+            /**********************query***************************/
+
+
+
             ///for overlap_list, the x_strand of all overlaps are 0, so the tmp.rev is the same as the y_strand
             tmp.rev = overlap_list->list[i].y_pos_strand;
 
@@ -505,6 +564,7 @@ static void worker_ovec(void *data, long i, int tid)
 	}
 }
 
+
 static void worker_ovec_related_reads(void *data, long i, int tid)
 {
 	ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
@@ -529,6 +589,7 @@ static void worker_ovec_related_reads(void *data, long i, int tid)
 	}
 }
 
+
 static inline long long get_N_occ(char* seq, long long length)
 {
 	long long j, N_occ = 0;
@@ -681,6 +742,7 @@ void ha_overlap_and_correct(int round)
 	free(e);
 }
 
+
 void update_overlaps(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf, 
 UC_Read* g_read, UC_Read* overlap_read, int is_match, int is_exact)
 {
@@ -762,6 +824,116 @@ UC_Read* g_read, UC_Read* overlap_read, int is_match, int is_exact)
     }
 }
 
+
+int check_chain_indels(Fake_Cigar* chain, long long xBeg, long long xEnd, float indel_rate)
+{
+    uint64_t i = 0;
+    long long indels = 0, xOffset;
+    if(chain->length != 0)
+    {
+        indels += abs(get_fake_gap_shift(chain, 0));
+        xOffset = get_fake_gap_pos(chain, 0);
+        if(indels > (xOffset - xBeg + 1) * indel_rate) return 0;
+
+        for (i = 1; i < chain->length; i++)
+        {
+            indels += abs((get_fake_gap_shift(chain, i) - get_fake_gap_shift(chain, i-1)));
+            xOffset = get_fake_gap_pos(chain, i);
+            if(indels > (xOffset - xBeg + 1) * indel_rate) return 0;
+        }
+    }
+
+    if(indels > (xEnd - xBeg + 1) * indel_rate) return 0;
+    return 1;
+}
+
+void update_overlaps_chain_width(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf, 
+UC_Read* g_read, UC_Read* overlap_read, int is_match, int is_exact, float indel_rate)
+{
+
+    uint64_t inner_j = 0;
+    uint64_t j = 0;
+    long long x_overlapLen, y_overlapLen;
+    while (j < overlap_list->length && inner_j < paf->length)
+    {
+        if(overlap_list->list[j].y_id < paf->buffer[inner_j].tn)
+        {
+            j++;
+        }
+        else if(overlap_list->list[j].y_id > paf->buffer[inner_j].tn)
+        {
+            inner_j++;
+        }
+        else
+        {
+            if(check_chain_indels(&(overlap_list->list[j].f_cigar), overlap_list->list[j].x_pos_s, 
+            overlap_list->list[j].x_pos_e, indel_rate) == 1)
+            {
+                if(overlap_list->list[j].y_pos_strand == paf->buffer[inner_j].rev)
+                {
+                    x_overlapLen = Get_qe(paf->buffer[inner_j]) - Get_qs(paf->buffer[inner_j]) + 1;
+                    y_overlapLen = Get_te(paf->buffer[inner_j]) - Get_ts(paf->buffer[inner_j]) + 1;
+                    if(x_overlapLen < y_overlapLen) x_overlapLen = y_overlapLen;
+                    x_overlapLen = x_overlapLen * 0.1;
+
+                    // if(
+                    // ((DIFF(overlap_list->list[j].x_pos_s, Get_qs(paf->buffer[inner_j])) < x_overlapLen)
+                    // && (DIFF(overlap_list->list[j].x_pos_e, Get_qe(paf->buffer[inner_j])) < x_overlapLen))
+                    // ||
+                    // ((DIFF(overlap_list->list[j].y_pos_s, Get_ts(paf->buffer[inner_j])) < x_overlapLen)
+                    // && (DIFF(overlap_list->list[j].y_pos_e, Get_te(paf->buffer[inner_j])) < x_overlapLen)))
+                    if(
+                    ((DIFF(overlap_list->list[j].x_pos_s, Get_qs(paf->buffer[inner_j])) < (uint64_t)x_overlapLen)
+                    && (DIFF(overlap_list->list[j].x_pos_e, Get_qe(paf->buffer[inner_j])) < (uint64_t)x_overlapLen))
+                    || 
+                    ((DIFF(overlap_list->list[j].y_pos_s, Get_ts(paf->buffer[inner_j])) < (uint64_t)x_overlapLen)
+                    && (DIFF(overlap_list->list[j].y_pos_e, Get_te(paf->buffer[inner_j])) < (uint64_t)x_overlapLen))
+                    )
+                    {
+                        overlap_list->list[j].is_match = is_match;
+                        overlap_list->list[j].strong = paf->buffer[inner_j].ml;
+                        overlap_list->list[j].without_large_indel = paf->buffer[inner_j].no_l_indel;
+                        if(is_exact == 1)
+                        {
+                            if(overlap_list->list[j].y_pos_strand == 0)
+                            {
+                                recover_UC_Read(overlap_read, &R_INF, overlap_list->list[j].y_id);
+                            }
+                            else
+                            {
+                                recover_UC_Read_RC(overlap_read, &R_INF, overlap_list->list[j].y_id);
+                            }
+                            if(if_exact_match(g_read->seq, g_read->length, overlap_read->seq, overlap_read->length, 
+                            overlap_list->list[j].x_pos_s, overlap_list->list[j].x_pos_e, 
+                            overlap_list->list[j].y_pos_s, overlap_list->list[j].y_pos_e))
+                            {
+                                overlap_list->list[j].shared_seed = 1;
+                            }
+                            else
+                            {
+                                overlap_list->list[j].shared_seed = 0;
+                            }
+                        }
+                    }
+                    else
+                    {
+                        overlap_list->list[j].is_match = 3;
+                    }
+                }
+                else
+                {
+                    overlap_list->list[j].is_match = 3;
+                }
+            }
+
+            j++;
+            inner_j++;
+        }
+    }
+}
+
+
+
 void update_exact_overlaps(overlap_region_alloc* overlap_list, UC_Read* g_read, UC_Read* overlap_read)
 {
     uint64_t j;
@@ -769,6 +941,12 @@ void update_exact_overlaps(overlap_region_alloc* overlap_list, UC_Read* g_read,
     {
         if (overlap_list->list[j].is_match != 1)
         {
+            if((overlap_list->list[j].x_pos_e + 1 - overlap_list->list[j].x_pos_s) != 
+               (overlap_list->list[j].y_pos_e + 1 - overlap_list->list[j].y_pos_s))
+            {
+                continue;
+            }
+
             if(overlap_list->list[j].y_pos_strand == 0)
             {
                 recover_UC_Read(overlap_read, &R_INF, overlap_list->list[j].y_id);
@@ -819,9 +997,11 @@ void ha_print_ovlp_stat(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, long long
 void fill_chain(Fake_Cigar* chain, char* x_string, char* y_string, long long xBeg, long long yBeg,
 long long x_readLen, long long y_readLen, Cigar_record* cigar, uint8_t* c2n)
 {
-    long long i, xOffset, yOffset, xRegionLen, yRegionLen, /**bandLen,**/ maxXpos, maxYpos, mapScore, zdroped;
+    /**
+    long long i, xOffset, yOffset, xRegionLen, yRegionLen, maxXpos, maxYpos, mapGlobalScore, mapExtentScore, zdroped;
+    long long xBuoundaryScore, yBuoundaryScore;
     ///float band_rate = 0.08;
-    int endbouns;
+    int endbouns,mode;
     if(chain->length <= 0) return;
 
     kvec_t(uint8_t) x_num;
@@ -850,7 +1030,9 @@ long long x_readLen, long long y_readLen, Cigar_record* cigar, uint8_t* c2n)
         ///text is x, query is y
         afine_gap_alignment(x_string, x_num.a, xRegionLen, y_string, y_num.a, yRegionLen, 
         c2n, BACKWARD_KSW, MATCH_SCORE_KSW, MISMATCH_SCORE_KSW, GAP_OPEN_KSW, GAP_EXT_KSW,
-        /**bandLen,**/BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapScore, &zdroped);
+        BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapGlobalScore, 
+        &mapExtentScore, &xBuoundaryScore, &yBuoundaryScore, &zdroped);
+
         // fprintf(stderr, "* xOffset: %lld, yOffset: %lld, xRegionLen: %lld, yRegionLen: %lld, bandLen: %lld, maxXpos: %lld, maxYpos: %lld, zdroped: %lld\n",
         // xOffset, yOffset, xRegionLen, yRegionLen, BAND_KSW, maxXpos, maxYpos, zdroped);
     }
@@ -891,7 +1073,8 @@ long long x_readLen, long long y_readLen, Cigar_record* cigar, uint8_t* c2n)
         ///text is x, query is y
         afine_gap_alignment(x_string+xOffset, x_num.a, xRegionLen, y_string+yOffset, y_num.a, yRegionLen, 
         c2n, FORWARD_KSW, MATCH_SCORE_KSW, MISMATCH_SCORE_KSW, GAP_OPEN_KSW, GAP_EXT_KSW,
-        /**bandLen,**/BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapScore, &zdroped);
+        BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapGlobalScore, 
+        &mapExtentScore, &xBuoundaryScore, &yBuoundaryScore, &zdroped);
         // fprintf(stderr, "# xOffset: %lld, yOffset: %lld, xRegionLen: %lld, yRegionLen: %lld, bandLen: %lld, maxXpos: %lld, maxYpos: %lld, zdroped: %lld\n",
         // xOffset, yOffset, xRegionLen, yRegionLen, BAND_KSW, maxXpos, maxYpos, zdroped);
     }
@@ -899,6 +1082,7 @@ long long x_readLen, long long y_readLen, Cigar_record* cigar, uint8_t* c2n)
 
     kv_destroy(x_num);
     kv_destroy(y_num);
+    **/
 }
 void Final_phasing(overlap_region_alloc* overlap_list, Cigar_record_alloc* cigarline,
 UC_Read* g_read, UC_Read* overlap_read, uint8_t* c2n)
@@ -952,10 +1136,6 @@ UC_Read* g_read, UC_Read* overlap_read, uint8_t* c2n)
 static void worker_ov_final(void *data, long i, int tid)
 {
 	ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
-    uint8_t c2n[256]; // this may be moved to ha_ovec_buf_t, but it should be fast to populate anyway
-    memset(c2n, 4, 256);
-    c2n[(uint8_t)'A'] = c2n[(uint8_t)'a'] = 0; c2n[(uint8_t)'C'] = c2n[(uint8_t)'c'] = 1;
-	c2n[(uint8_t)'G'] = c2n[(uint8_t)'g'] = 2; c2n[(uint8_t)'T'] = c2n[(uint8_t)'t'] = 3; // build the encoding table
 
 	//get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.001, 0);
 	ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.001, asm_opt.max_n_chain, 0);
@@ -982,6 +1162,80 @@ static void worker_ov_final(void *data, long i, int tid)
 	push_final_overlaps(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, &b->olist, 2);
 }
 
+
+
+
+void reset_final_overlaps(overlap_region_alloc *overlap_list)
+{
+    
+    uint64_t i;
+    for (i = 0; i < overlap_list->length; i++)
+    {
+        if (overlap_list->list[i].is_match == 1 || overlap_list->list[i].is_match == 2)
+        {
+            overlap_list->list[i].x_pos_s = overlap_list->list[i].x_pos_e = (uint32_t)-1;
+            overlap_list->list[i].y_pos_s = overlap_list->list[i].y_pos_e = (uint32_t)-1;
+            overlap_list->list[i].is_match = 0;
+        }
+    }
+
+    ha_sort_list_by_anchor(overlap_list);
+}
+
+void debug_affine_gap_alignment(overlap_region_alloc *overlap_list, UC_Read* g_read, UC_Read* overlap_read)
+{
+    uint64_t i;
+    kvec_t(uint8_t) x_num;
+    kvec_t(uint8_t) y_num;
+    kv_init(x_num);
+    kv_init(y_num);
+    for (i = 0; i < overlap_list->length; i++)
+    {
+        if (overlap_list->list[i].is_match == 1 && overlap_list->list[i].shared_seed == 1)
+        {
+
+            kv_resize(uint8_t, x_num, (uint64_t)(Get_READ_LENGTH(R_INF, overlap_list->list[i].x_id)));
+            kv_resize(uint8_t, y_num, (uint64_t)(Get_READ_LENGTH(R_INF, overlap_list->list[i].y_id)));
+
+            get_affine_gap_score(&(overlap_list->list[i]), g_read, overlap_read, x_num.a, y_num.a,
+            overlap_list->list[i].x_pos_e + 1 - overlap_list->list[i].x_pos_s, 
+            overlap_list->list[i].y_pos_e + 1 - overlap_list->list[i].y_pos_s);
+        }
+    }
+
+    kv_destroy(x_num);
+    kv_destroy(y_num);
+}
+
+static void worker_ov_final_high_het(void *data, long i, int tid)
+{
+    ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
+
+    ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, HIGH_HET_ERROR_RATE, asm_opt.max_n_chain, 1);
+
+    overlap_region_sort_y_id(b->olist.list, b->olist.length);
+    ma_hit_sort_tn(R_INF.paf[i].buffer, R_INF.paf[i].length);
+    ma_hit_sort_tn(R_INF.reverse_paf[i].buffer, R_INF.reverse_paf[i].length);
+
+
+    ///update_overlaps(&b->olist, &(R_INF.paf[i]), &b->self_read, &b->ovlp_read, 1, 1);
+    update_overlaps_chain_width(&b->olist, &(R_INF.paf[i]), &b->self_read, &b->ovlp_read, 1, 1, 0.002);
+    update_overlaps(&b->olist, &(R_INF.reverse_paf[i]), &b->self_read, &b->ovlp_read, 2, 0);
+    ///recover missing exact overlaps
+    update_exact_overlaps(&b->olist, &b->self_read, &b->ovlp_read);
+
+    
+    ///Final_phasing(&overlap_list, &cigarline, &g_read, &overlap_read, c2n);
+    push_final_overlaps(&(R_INF.paf[i]), R_INF.reverse_paf, &b->olist, 1);
+    push_final_overlaps(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, &b->olist, 2);
+
+    ///debug_affine_gap_alignment(&b->olist, &b->self_read, &b->ovlp_read);
+
+    reset_final_overlaps(&b->olist);
+    correct_overlap_high_het(&b->olist, &R_INF, &b->self_read, &b->correct, &b->ovlp_read);
+    push_final_overlaps_increment(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, &b->olist, 2);
+}
+
 void Output_PAF()
 {
     fprintf(stderr, "Writing PAF to disk ...... \n");
@@ -1183,9 +1437,17 @@ void ha_overlap_final(void)
 	ha_ovec_buf_t **b;
 	CALLOC(b, asm_opt.thread_num);
 	for (i = 0; i < asm_opt.thread_num; ++i)
-		b[i] = ha_ovec_init(1, 1);
+		b[i] = ha_ovec_init(asm_opt.flag & HA_F_HIGH_HET, 1);///b[i] = ha_ovec_init(1, 1);
 	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF, &hom_cov, &het_cov); // build the index
-	kt_for(asm_opt.thread_num, worker_ov_final, b, R_INF.total_reads);
+    if(asm_opt.flag & HA_F_HIGH_HET)
+    {
+        kt_for(asm_opt.thread_num, worker_ov_final_high_het, b, R_INF.total_reads);
+    }
+    else
+    {
+        kt_for(asm_opt.thread_num, worker_ov_final, b, R_INF.total_reads);
+    }
+    
 	ha_pt_destroy(ha_idx);
 	ha_idx = 0;
 	for (i = 0; i < asm_opt.thread_num; ++i)

CommandLines.cpp

@@ -23,6 +23,7 @@ static ko_longopt_t long_options[] = {
 	{ "ex-iter",       ko_required_argument, 308 },
     { "purge-cov",     ko_required_argument, 309 },
     { "pri-range",     ko_required_argument, 310 },
+    { "high-het",      ko_no_argument, 311 },
 	{ 0, 0, 0 }
 };
 
@@ -77,6 +78,8 @@ void Print_H(hifiasm_opt_t* asm_opt)
                                      asm_opt->purge_overlap_len);
     fprintf(stderr, "    --purge-cov INT\n");
     fprintf(stderr, "                coverage upper bound of Purge-dups [auto]\n");
+    fprintf(stderr, "    --high-het  enable this mode for high heterozygosity sample\n");
+    
 
     fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n");
     fprintf(stderr, "See `man ./hifiasm.1' for detailed description of these command-line options.\n");
@@ -433,7 +436,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
             {
                 asm_opt->recover_atg_cov_min = asm_opt->recover_atg_cov_max = -1;
             }
-        } 
+        }
+        else if (c == 311) asm_opt->flag |= HA_F_HIGH_HET;
         else if (c == 'l')
         {   ///0: disable purge_dup; 1: purge containment; 2: purge overlap
             asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg);

CommandLines.h

@@ -3,7 +3,7 @@
 
 #include <pthread.h>
 
-#define HA_VERSION "0.9-r289"
+#define HA_VERSION "0.9-r297"
 
 #define VERBOSE 0
 
@@ -17,6 +17,7 @@
 #define HA_F_PURGE_JOIN      0x80
 #define HA_F_BAN_POST_JOIN   0x100
 #define HA_F_BAN_ASSEMBLY    0x200
+#define HA_F_HIGH_HET        0x400
 
 #define HA_MIN_OV_DIFF       0.02 // min sequence divergence in an overlap
 

Correct.cpp

@@ -6,7 +6,7 @@
 #include "Levenshtein_distance.h"
 #include "Assembly.h"
 #include "CommandLines.h"
-///#include "ksw2.h"
+#include "ksw2.h"
 
 
 
@@ -7483,39 +7483,46 @@ void add_existing_cell_to_cigar_record_with_different_base(Cigar_record* dummy,
     }
 }
 
-
-
-void afine_gap_alignment(const char *tseq, uint8_t* tnum, const int tl, 
-const char *qseq, uint8_t* qnum, const int ql, const uint8_t *c2n, const int strand,
+void afine_gap_alignment(const char *qseq, uint8_t* qnum, const int ql, 
+const char *tseq, uint8_t* tnum, const int tl, const uint8_t *c2n, const int strand,
 int sc_mch, int sc_mis, int gapo, int gape, int bandLen, int zdrop, int end_bonus,
-long long* max_t_pos, long long* max_q_pos, long long* score, long long* droped)
+long long* max_q_pos, long long* max_t_pos, long long* global_score, 
+long long* extention_score, long long* q_boundary_score, long long* q_boundary_t_coordinate,
+long long* t_boundary_score, long long* t_boundary_q_coordinate,
+long long* droped, int mode)
 {
     /**************for ksw2**************/
-    // (*max_t_pos) = (*max_q_pos) = -1;
-	// int i, a = sc_mch, b = sc_mis < 0? sc_mis : -sc_mis; // a>0 and b<0
+    (*max_t_pos) = (*max_q_pos) = -1;
+	int i, a = sc_mch, b = sc_mis < 0? sc_mis : -sc_mis; // a>0 and b<0
     
-	// int8_t mat[25] = {(int8_t)a,(int8_t)b,(int8_t)b,(int8_t)b,0, 
-    // (int8_t)b,(int8_t)a,(int8_t)b,(int8_t)b,0, (int8_t)b,(int8_t)b,(int8_t)a,(int8_t)b,0, 
-    // (int8_t)b,(int8_t)b,(int8_t)b,(int8_t)a,0, 0,0,0,0,0};
-	// ksw_extz_t ez;
-	// memset(&ez, 0, sizeof(ksw_extz_t));
+	int8_t mat[25] = {(int8_t)a,(int8_t)b,(int8_t)b,(int8_t)b,0, 
+    (int8_t)b,(int8_t)a,(int8_t)b,(int8_t)b,0, (int8_t)b,(int8_t)b,(int8_t)a,(int8_t)b,0, 
+    (int8_t)b,(int8_t)b,(int8_t)b,(int8_t)a,0, 0,0,0,0,0};
+	ksw_extz_t ez;
+	memset(&ez, 0, sizeof(ksw_extz_t));
 
-    // if(strand == FORWARD_KSW)
-    // {
-    //     for (i = 0; i < tl; ++i) tnum[i] = c2n[(uint8_t)tseq[i]]; // encode to 0/1/2/3
-	//     for (i = 0; i < ql; ++i) qnum[i] = c2n[(uint8_t)qseq[i]];
-    // }
-    // else if(strand == BACKWARD_KSW)
-    // {
-    //     for (i = 0; i < tl; ++i) tnum[i] = c2n[(uint8_t)tseq[tl - i - 1]]; // encode to 0/1/2/3
-	//     for (i = 0; i < ql; ++i) qnum[i] = c2n[(uint8_t)qseq[ql - i - 1]];
-    // }
-    // ksw_extz2_sse(0, ql, qnum, tl, tnum, 5, mat, gapo, gape, bandLen, zdrop, end_bonus, 0, &ez);
-    // (*score) = ez.max;
-    // (*max_t_pos) = ez.max_t;
-    // (*max_q_pos) = ez.max_q;
-    // (*droped) = ez.zdropped;
-    // free(ez.cigar);
+    if(strand == FORWARD_KSW)
+    {
+         for (i = 0; i < tl; ++i) tnum[i] = c2n[(uint8_t)tseq[i]]; // encode to 0/1/2/3
+	     for (i = 0; i < ql; ++i) qnum[i] = c2n[(uint8_t)qseq[i]];
+    }
+    else if(strand == BACKWARD_KSW)
+    {
+         for (i = 0; i < tl; ++i) tnum[i] = c2n[(uint8_t)tseq[tl - i - 1]]; // encode to 0/1/2/3
+	     for (i = 0; i < ql; ++i) qnum[i] = c2n[(uint8_t)qseq[ql - i - 1]];
+    }
+    ksw_extz2_sse(0, ql, qnum, tl, tnum, 5, mat, gapo, gape, bandLen, zdrop, end_bonus, 
+    mode, &ez);
+    (*global_score) = ez.score;
+    (*extention_score) = ez.max;
+    (*q_boundary_score) = ez.mqe;
+    (*q_boundary_t_coordinate) = ez.mqe_t;
+    (*t_boundary_score) = ez.mte;
+    (*t_boundary_q_coordinate) = ez.mte_q;
+    (*max_t_pos) = ez.max_t;
+    (*max_q_pos) = ez.max_q;
+    (*droped) = ez.zdropped;
+    free(ez.cigar);
     
     /**
 	for (i = 0; i < ez.n_cigar; ++i) // print CIGAR
@@ -7524,3 +7531,791 @@ long long* max_t_pos, long long* max_q_pos, long long* score, long long* droped)
     **/
    /**************for ksw2**************/
 }
+
+int fill_chain_by_affine_gap_debug(Fake_Cigar* chain, char* x_string, char* y_string, overlap_region* ovc,
+long long x_readLen, long long y_readLen, Cigar_record* cigar, uint8_t* c2n, uint8_t* x_num, uint8_t* y_num,
+long long* minus_score_thres, long long* final_scores)
+{
+    long long i, xOffset, yOffset, xRegionLen, yRegionLen, /**bandLen,**/ maxXpos, maxYpos, zdroped;
+    long long mapGlobalScore, mapExtentScore;
+    long long xBuoundaryScore, xBuoundaryYcoordinate, yBuoundaryScore, yBuoundaryXcoordinate;
+    ///float band_rate = 0.08;
+    int endbouns, mode;
+    long long xBeg, yBeg;
+    xBeg = ovc->x_pos_s;
+    yBeg = ovc->y_pos_s;
+    if(chain->length <= 0) return 0;
+    // long long minus_score_thres = (EstimateOlen*HIGH_HET_ERROR_RATE*(MATCH_SCORE_KSW+(MAX(MISMATCH_SCORE_KSW,GAP_EXT_KSW))));
+    // long long total_score_thres = EstimateOlen*MATCH_SCORE_KSW - minus_score_thres;
+    long long sum_score = 0, current_ovlp = 0, zdrop_occ = 0;
+    long long new_xBeg, new_yBeg, new_xEnd, new_yEnd;
+    new_xBeg = ovc->x_pos_s;
+    new_yBeg = ovc->y_pos_s;
+    new_xEnd = ovc->x_pos_e;
+    new_yEnd = ovc->y_pos_e;
+    ///long long sub_score_sum;
+    ///deal with region 0 backward
+    i = 0;
+    endbouns = 0;
+
+    xOffset = get_fake_gap_pos(chain, 0);
+    xOffset = xOffset - 1;
+    yOffset = (xOffset - xBeg) + yBeg + get_fake_gap_shift(chain, 0);
+    if(xOffset >= 0 && yOffset >= 0)
+    {
+        xRegionLen = xOffset + 1;
+        yRegionLen = yOffset + 1;
+        //note here cannot use DIFF(xRegionLen, yRegionLen)
+        // bandLen = (MIN(xRegionLen, yRegionLen))*band_rate;
+        // if(bandLen == 0) bandLen = MIN(xRegionLen, yRegionLen);
+
+
+        ///do alignment backward
+        ///for beginning part and end part, must use exact mode
+        mode = KSW_EZ_SCORE_ONLY;
+        afine_gap_alignment(x_string, x_num, xRegionLen, y_string, y_num, yRegionLen, 
+        c2n, BACKWARD_KSW, MATCH_SCORE_KSW, MISMATCH_SCORE_KSW, GAP_OPEN_KSW, GAP_EXT_KSW,
+        /**bandLen,**/BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapGlobalScore, 
+        &mapExtentScore, &xBuoundaryScore, &xBuoundaryYcoordinate, 
+        &yBuoundaryScore, &yBuoundaryXcoordinate, &zdroped, mode);
+
+        
+        
+
+        if(!zdroped)
+        {
+            if(xRegionLen <= yRegionLen)
+            {
+                sum_score += xBuoundaryScore;
+                new_yBeg = yRegionLen - xBuoundaryYcoordinate - 1;
+            }
+            else
+            {
+                sum_score += yBuoundaryScore;
+                new_xBeg = xRegionLen - yBuoundaryXcoordinate - 1;
+            }
+        }
+        else
+        {   ///return 0;
+            sum_score += mapExtentScore;
+            if(xRegionLen <= yRegionLen)
+            {
+                
+                sum_score -= (GAP_OPEN_KSW + (xRegionLen - maxXpos)*GAP_EXT_KSW);
+            }
+            else
+            {
+                sum_score -= (GAP_OPEN_KSW + (yRegionLen - maxYpos)*GAP_EXT_KSW);
+            }
+        }   
+    }
+
+    ///align forward
+    for (i = 0; i < (long long)chain->length; i++)
+    {
+        // xOffset = get_fake_gap_pos(chain, i);
+        // yOffset = xOffset + get_fake_gap_shift(chain, i);
+        xOffset = get_fake_gap_pos(chain, i);
+        yOffset = (xOffset - xBeg) + yBeg + get_fake_gap_shift(chain, i);
+        ///last region
+        if(i == (long long)(chain->length - 1))
+        {
+            endbouns = 0;
+            xRegionLen = x_readLen - xOffset;
+            yRegionLen = y_readLen - yOffset;
+            ///for beginning part and end part, must use exact mode
+            mode = KSW_EZ_SCORE_ONLY;
+            //note here cannot use DIFF(xRegionLen, yRegionLen)
+            // bandLen = (MIN(xRegionLen, yRegionLen))*band_rate;
+            // if(bandLen == 0) bandLen = MIN(xRegionLen, yRegionLen);
+        }
+        else
+        {
+            ///higher endbouns for middle regions
+            endbouns = MATCH_SCORE_KSW;
+            xRegionLen = get_fake_gap_pos(chain, i+1) - xOffset;
+            yRegionLen = (get_fake_gap_pos(chain, i+1) + get_fake_gap_shift(chain, i+1)) - 
+                         (get_fake_gap_pos(chain, i) + get_fake_gap_shift(chain, i));
+            mode = KSW_EZ_SCORE_ONLY | KSW_EZ_APPROX_MAX | KSW_EZ_APPROX_DROP;
+            // bandLen = MAX((MIN(xRegionLen, yRegionLen))*band_rate, DIFF(xRegionLen, yRegionLen));
+            // if(bandLen == 0) bandLen = MIN(xRegionLen, yRegionLen);
+        }
+
+
+        if(minus_score_thres)
+        {
+            current_ovlp = MIN((xOffset + 1 - xBeg), (yOffset + 1 - yBeg));
+            current_ovlp = current_ovlp*MATCH_SCORE_KSW;
+            if(current_ovlp - sum_score > (*minus_score_thres))
+            {
+                return 0;
+            }
+        }
+
+        if(xOffset < 0) xOffset = 0;
+        if(yOffset < 0) yOffset = 0;
+        if(xRegionLen < 0) xRegionLen = 0;
+        if(yRegionLen < 0) yRegionLen = 0;
+
+        ///do alignment forward
+        ///text is x, query is y
+        afine_gap_alignment(x_string+xOffset, x_num, xRegionLen, y_string+yOffset, y_num, yRegionLen, 
+        c2n, FORWARD_KSW, MATCH_SCORE_KSW, MISMATCH_SCORE_KSW, GAP_OPEN_KSW, GAP_EXT_KSW,
+        /**bandLen,**/BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapGlobalScore, 
+        &mapExtentScore, &xBuoundaryScore, &xBuoundaryYcoordinate, 
+        &yBuoundaryScore, &yBuoundaryXcoordinate, &zdroped, mode);
+        // fprintf(stderr, "# xOffset: %lld, yOffset: %lld, xRegionLen: %lld, yRegionLen: %lld, bandLen: %lld, maxXpos: %lld, maxYpos: %lld, zdroped: %lld\n",
+        // xOffset, yOffset, xRegionLen, yRegionLen, BAND_KSW, maxXpos, maxYpos, zdroped);
+
+        if(!zdroped)
+        {
+            if(i != (long long)(chain->length - 1))
+            {
+                sum_score += mapGlobalScore;
+            }
+            else
+            {
+                if(xRegionLen <= yRegionLen)
+                {
+                    sum_score += xBuoundaryScore;
+                    new_yEnd = yOffset + xBuoundaryYcoordinate;
+                }
+                else
+                {
+                    sum_score += yBuoundaryScore;
+                    new_xEnd = xOffset + yBuoundaryXcoordinate;
+                    // if(new_xEnd != (long long)ovc->x_pos_e)
+                    // {
+                    //     fprintf(stderr, "\n******direction: %u, new_xBeg: %lld, new_xEnd: %lld, new_yBeg: %lld, new_yEnd: %lld, old_xBeg: %u, old_xEnd: %u, old_yBeg: %u, old_yEnd: %u\n", 
+                    //     ovc->y_pos_strand, new_xBeg, new_xEnd, new_yBeg, new_yEnd, ovc->x_pos_s, ovc->x_pos_e, ovc->y_pos_s, ovc->y_pos_e);
+                    //     fprintf(stderr, "x_readLen: %lld, y_readLen: %lld\n", x_readLen, y_readLen);
+                    //     fprintf(stderr, "xID: %lld, yID: %lld\n", ovc->x_id, ovc->y_id);
+                    //     fprintf(stderr, "xRegionLen: %lld, yRegionLen: %lld\n", xRegionLen, yRegionLen);
+                    //     fprintf(stderr, "xOffset: %lld, yOffset: %lld\n", xOffset, yOffset);
+                    //     fprintf(stderr, "yBuoundaryXcoordinate: %lld\n", yBuoundaryXcoordinate);
+                    // }
+                }
+            }
+        }
+        else
+        {
+            ///return 0;
+            if(i != (long long)(chain->length - 1)) zdrop_occ++;
+            if(zdrop_occ > 1) return 0;
+            sum_score += mapExtentScore;
+            if(xRegionLen <= yRegionLen)
+            {
+                sum_score -= (GAP_OPEN_KSW + (xRegionLen - maxXpos)*GAP_EXT_KSW);
+            }
+            else
+            {
+                sum_score -= (GAP_OPEN_KSW + (yRegionLen - maxYpos)*GAP_EXT_KSW);
+            }
+        }
+    }
+
+    (*final_scores) = sum_score;
+    if(new_xBeg != (long long)ovc->x_pos_s || new_xEnd != (long long)ovc->x_pos_e || 
+       new_yBeg != (long long)ovc->y_pos_s || new_yEnd != (long long)ovc->y_pos_e)
+    {   
+        // fprintf(stderr, "\ntttdirection: %u, new_xBeg: %lld, new_xEnd: %lld, new_yBeg: %lld, new_yEnd: %lld, old_xBeg: %u, old_xEnd: %u, old_yBeg: %u, old_yEnd: %u\n", 
+        // ovc->y_pos_strand, new_xBeg, new_xEnd, new_yBeg, new_yEnd, ovc->x_pos_s, ovc->x_pos_e, ovc->y_pos_s, ovc->y_pos_e);
+        // fprintf(stderr, "x_readLen: %lld, y_readLen: %lld\n", x_readLen, y_readLen);
+        // fprintf(stderr, "xID: %lld, yID: %lld\n", ovc->x_id, ovc->y_id);
+        // for (i = 0; i < (long long)chain->length; i++)
+        // {
+        //     fprintf(stderr,"i: %lld, x_pos: %d, offset: %d\n", 
+        //     i, get_fake_gap_pos(chain, i), get_fake_gap_shift(chain, i));
+        // }
+    }
+    else
+    {
+        // fprintf(stderr, "\nkkkdirection: %u, new_xBeg: %lld, new_xEnd: %lld, new_yBeg: %lld, new_yEnd: %lld, old_xBeg: %u, old_xEnd: %u, old_yBeg: %u, old_yEnd: %u\n", 
+        // ovc->y_pos_strand, new_xBeg, new_xEnd, new_yBeg, new_yEnd, ovc->x_pos_s, ovc->x_pos_e, ovc->y_pos_s, ovc->y_pos_e);
+        // fprintf(stderr, "x_readLen: %lld, y_readLen: %lld\n", x_readLen, y_readLen);
+        // fprintf(stderr, "xID: %lld, yID: %lld\n", ovc->x_id, ovc->y_id);
+        // for (i = 0; i < (long long)chain->length; i++)
+        // {
+        //     fprintf(stderr,"i: %lld, x_pos: %d, offset: %d\n", 
+        //     i, get_fake_gap_pos(chain, i), get_fake_gap_shift(chain, i));
+        // }
+    }
+    
+    return 1;
+}
+
+
+
+int fill_chain_by_affine_gap(Fake_Cigar* chain, char* x_string, char* y_string, overlap_region* ovc,
+long long x_readLen, long long y_readLen, Cigar_record* cigar, uint8_t* c2n, uint8_t* x_num, uint8_t* y_num,
+long long* minus_score_thres, long long* final_scores)
+{
+    long long i, xOffset, yOffset, xRegionLen, yRegionLen, /**bandLen,**/ maxXpos, maxYpos, zdroped;
+    long long mapGlobalScore, mapExtentScore;
+    long long xBuoundaryScore, xBuoundaryYcoordinate, yBuoundaryScore, yBuoundaryXcoordinate;
+    ///float band_rate = 0.08;
+    int endbouns, mode;
+    long long xBeg, yBeg;
+    xBeg = ovc->x_pos_s;
+    yBeg = ovc->y_pos_s;
+    if(chain->length <= 0) return 0;
+    long long sum_score = 0, current_ovlp = 0, zdrop_occ = 0;
+    long long chain_num = (long long)chain->length - 1;
+    
+    ///align forward
+    for (i = 0; i < chain_num; i++)
+    {
+        xOffset = get_fake_gap_pos(chain, i);
+        yOffset = (xOffset - xBeg) + yBeg + get_fake_gap_shift(chain, i);
+        ///last region
+        
+        ///higher endbouns for middle regions
+        endbouns = MATCH_SCORE_KSW;
+        xRegionLen = get_fake_gap_pos(chain, i+1) - xOffset;
+        yRegionLen = (get_fake_gap_pos(chain, i+1) + get_fake_gap_shift(chain, i+1)) - 
+                        (get_fake_gap_pos(chain, i) + get_fake_gap_shift(chain, i));
+        ///last region
+        if(i == chain_num - 1)
+        {
+            xRegionLen++;
+            yRegionLen++;
+        }
+
+        mode = KSW_EZ_SCORE_ONLY | KSW_EZ_APPROX_MAX | KSW_EZ_APPROX_DROP;
+            
+
+        if(minus_score_thres)
+        {
+            current_ovlp = MIN((xOffset - xBeg), (yOffset - yBeg));
+            current_ovlp = current_ovlp*MATCH_SCORE_KSW;
+            if(current_ovlp - sum_score > (*minus_score_thres))
+            {
+                return 0;
+            }
+        }
+
+        if(xOffset < 0) xOffset = 0;
+        if(yOffset < 0) yOffset = 0;
+        if(xRegionLen < 0) xRegionLen = 0;
+        if(yRegionLen < 0) yRegionLen = 0;
+
+        ///do alignment forward
+        ///text is x, query is y
+        afine_gap_alignment(x_string+xOffset, x_num, xRegionLen, y_string+yOffset, y_num, yRegionLen, 
+        c2n, FORWARD_KSW, MATCH_SCORE_KSW, MISMATCH_SCORE_KSW, GAP_OPEN_KSW, GAP_EXT_KSW,
+        /**bandLen,**/BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapGlobalScore, 
+        &mapExtentScore, &xBuoundaryScore, &xBuoundaryYcoordinate, 
+        &yBuoundaryScore, &yBuoundaryXcoordinate, &zdroped, mode);
+        
+
+        if(!zdroped)
+        {
+            sum_score += mapGlobalScore;
+        }
+        else
+        {
+            ///return 0;
+            zdrop_occ++;
+            ///if(zdrop_occ > 1) return 0;
+            sum_score += mapExtentScore;
+            if(xRegionLen <= yRegionLen)
+            {
+                sum_score -= (GAP_OPEN_KSW + (xRegionLen - maxXpos)*GAP_EXT_KSW);
+            }
+            else
+            {
+                sum_score -= (GAP_OPEN_KSW + (yRegionLen - maxYpos)*GAP_EXT_KSW);
+            }
+        }
+    }
+
+    (*final_scores) = sum_score;
+    return 1;
+}
+
+
+
+long long get_affine_gap_score(overlap_region* ovc, UC_Read* g_read, UC_Read* overlap_read, uint8_t* x_num,
+uint8_t* y_num, uint64_t EstimateXOlen, uint64_t EstimateYOlen)
+{
+    char* x_string;
+    char* y_string;
+    uint64_t yStrand;
+    long long minus_score_thres = (MAX(EstimateXOlen, EstimateYOlen)*HIGH_HET_ERROR_RATE*(MATCH_SCORE_KSW+(MAX(MISMATCH_SCORE_KSW,GAP_EXT_KSW))));
+    long long total_score_thres = MAX(EstimateXOlen, EstimateYOlen)*MATCH_SCORE_KSW - minus_score_thres;
+
+    yStrand = ovc->y_pos_strand;
+
+    if(yStrand == 0)
+    {
+        recover_UC_Read(overlap_read, &R_INF, ovc->y_id);
+    }
+    else
+    {
+        recover_UC_Read_RC(overlap_read, &R_INF, ovc->y_id);
+    }
+    x_string = g_read->seq;
+    y_string = overlap_read->seq;
+
+    long long sum;
+    
+    if(fill_chain_by_affine_gap(&(ovc->f_cigar), x_string, y_string, ovc, Get_READ_LENGTH(R_INF, ovc->x_id), 
+    Get_READ_LENGTH(R_INF, ovc->y_id), NULL, seq_nt6_table, x_num, y_num, &minus_score_thres, &sum) == 0)
+    {
+        return 0;
+    }
+
+    if(sum >= total_score_thres) return 1;
+    return 0;
+}
+
+void recalcate_high_het_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, 
+                        UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read)
+{
+    long long j, k, i;
+    int threshold;
+    long long y_id;
+    int y_strand;
+    long long y_readLen;
+    long long x_start;
+    long long x_end;
+    long long x_len;
+    long long total_y_start;
+    long long total_y_end;
+    long long y_start;
+    long long Window_Len;
+    char* x_string;
+    char* y_string;
+    int end_site;
+    unsigned int error;
+    int real_y_start;
+    long long overlap_length;
+    int extra_begin, extra_end;
+    long long o_len;
+    kvec_t(uint8_t) x_num;
+    kvec_t(uint8_t) y_num;
+    kv_init(x_num);
+    kv_init(y_num);
+
+    for (j = 0; j < (long long)overlap_list->length; j++)
+    {
+        
+        if(overlap_list->list[j].w_list_length == 0) continue;
+        y_id = overlap_list->list[j].y_id;
+        y_strand = overlap_list->list[j].y_pos_strand;
+        y_readLen = Get_READ_LENGTH((*R_INF), y_id);
+
+        //i corresponding to each window of a overlap
+        //utilize the the end pos of pre-window in backwards
+        for (i = overlap_list->list[j].w_list_length - 1; i >= 0; i--)
+        {
+            ///the first matched window
+            if(overlap_list->list[j].w_list[i].y_end != -1)
+            {
+                ///note!!! need notification
+                ///this is the actual end postion in ystring
+                total_y_start = overlap_list->list[j].w_list[i].y_end 
+                            - overlap_list->list[j].w_list[i].extra_begin + 1;
+
+                ///k corresponding to all unmatched windows at the right side of overlap_list->list[j].w_list[i]
+                ///so k starts from i + 1, and end to the first matched window
+                for (k = i + 1; k < (long long)overlap_list->list[j].w_list_length && overlap_list->list[j].w_list[k].y_end == -1; k++)
+                {
+                    extra_begin = extra_end = 0;
+
+                    ///if y_start > y_readLen, direct terminate
+                    if (total_y_start >= y_readLen)
+                    {
+                        break;
+                    }
+                    
+                    ///there is no problem for x
+                    x_start = overlap_list->list[j].w_list[k].x_start;
+                    x_end = overlap_list->list[j].w_list[k].x_end;
+                    x_len = x_end - x_start + 1;
+                    ///there are two potiential reasons for unmatched window:
+                    ///1. this window has a large number of differences
+                    ///2. DP does not start from the right offset
+                    threshold = double_error_threshold(overlap_list->list[j].w_list[k].error_threshold, x_len);
+                    
+                    y_start = total_y_start;
+                    Window_Len = x_len + (threshold << 1);
+
+                    if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, R_INF, 
+                    &extra_begin, &extra_end, &y_start, &o_len))
+                    {
+                        break;
+                    }
+
+                    if(o_len + threshold < x_len)
+                    {
+                        break;
+                    }
+                    
+                    fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, 
+                    R_INF, y_id, extra_begin, extra_end);
+
+                    x_string = g_read->seq + x_start;
+                    y_string = dumy->overlap_region;
+
+                    ///note!!! need notification
+                    end_site = Reserve_Banded_BPM(y_string, Window_Len, x_string, x_len, threshold, &error);
+                    
+                    ///if error==-1, unmatched
+                    if (error!=(unsigned int)-1)
+                    {
+                        overlap_list->list[j].w_list[k].cigar.length = -1;
+                        overlap_list->list[j].w_list[k].y_start = y_start;
+                        overlap_list->list[j].w_list[k].y_end = y_start + end_site;
+                        overlap_list->list[j].w_list[k].error = (int)error;
+                        ///note!!! need notification
+                        overlap_list->list[j].w_list[k].extra_begin = extra_begin;
+                        overlap_list->list[j].w_list[k].extra_end = extra_end;
+                        overlap_list->list[j].w_list[k].error_threshold = threshold;
+                        
+                        overlap_list->list[j].align_length += x_len;
+                    }
+                    else
+                    {
+                        break;
+                    }
+
+                    ///note!!! need notification
+                    total_y_start = y_start + end_site - extra_begin + 1;
+                }
+                
+            }
+            
+        }
+        
+
+
+        //i corresponding to each window of a overlap
+        //utilize the the start pos of next window in forward
+        for (i = 0; i < (long long)overlap_list->list[j].w_list_length; i++)
+        {
+            ///find the first matched window, which should not be the first window
+            ///the pre-window of this matched window must be unmatched
+            if(overlap_list->list[j].w_list[i].y_end != -1 && i != 0 && overlap_list->list[j].w_list[i - 1].y_end == -1)
+            {
+                ///check if the start pos of this matched window has been calculated
+                if(overlap_list->list[j].w_list[i].cigar.length == -1)
+                {
+                    ///there is no problem for x
+                    x_start = overlap_list->list[j].w_list[i].x_start;
+                    x_end = overlap_list->list[j].w_list[i].x_end;
+                    x_len = x_end - x_start + 1;
+                    /****************************may have bugs********************************/
+                    threshold = overlap_list->list[j].w_list[i].error_threshold;
+                    /****************************may have bugs********************************/
+                    /****************************may have bugs********************************/
+                    ///should not adjust threshold, since this window can be matched by the old threshold
+                    ///threshold = Adjust_Threshold(threshold, x_len);
+                    /****************************may have bugs********************************/
+                    Window_Len = x_len + (threshold << 1);
+
+
+                    ///y_start is the real y_start
+                    y_start = overlap_list->list[j].w_list[i].y_start;
+                    extra_begin = overlap_list->list[j].w_list[i].extra_begin;
+                    extra_end = overlap_list->list[j].w_list[i].extra_end;
+                    o_len = Window_Len - extra_end - extra_begin;
+                    fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, 
+                    R_INF, y_id, extra_begin, extra_end);
+                    x_string = g_read->seq + x_start;
+                    y_string = dumy->overlap_region;
+
+                    ///note!!! need notification
+                    end_site = Reserve_Banded_BPM_PATH(y_string, Window_Len, x_string, x_len, threshold, &error, &real_y_start,
+                    &(dumy->path_length), dumy->matrix_bit, dumy->path, 
+                    overlap_list->list[j].w_list[i].error, overlap_list->list[j].w_list[i].y_end - y_start);
+
+
+                    ///y_start has already been calculated
+                    if (error != (unsigned int)-1)
+                    {
+                        ///this condition is always wrong
+                        ///in best case, real_y_start = threshold, end_site = Window_Len - threshold - 1
+                        if (end_site == Window_Len - 1 || real_y_start == 0)
+                        {
+                            if(fix_boundary(x_string, x_len, threshold, y_start, real_y_start, 
+                            end_site, extra_begin, extra_end, y_id, Window_Len, R_INF, dumy, 
+                            y_strand, error, &y_start, &real_y_start, &end_site, &extra_begin, 
+                            &extra_end, &error))
+                            {
+                                overlap_list->list[j].w_list[i].error = error;
+                                overlap_list->list[j].w_list[i].extra_begin = extra_begin;
+                                overlap_list->list[j].w_list[i].extra_end = extra_end;
+                            }
+                        }
+                                                 
+                        generate_cigar(dumy->path, dumy->path_length, &(overlap_list->list[j].w_list[i]),
+                        &real_y_start, &end_site, &error, x_string, x_len, y_string);   
+
+                        ///note!!! need notification
+                        real_y_start = y_start + real_y_start - extra_begin;
+                        overlap_list->list[j].w_list[i].y_start = real_y_start; 
+                        ///I forget why don't reduce the extra_begin for y_end
+                        ///it seems extra_begin will be reduced at the end of this function 
+                        overlap_list->list[j].w_list[i].y_end = y_start + end_site;
+                        overlap_list->list[j].w_list[i].error = error;                         
+                    }
+                    else
+                    {
+                        fprintf(stderr, "error\n");
+                    }
+                }
+                else
+                {
+                    real_y_start = overlap_list->list[j].w_list[i].y_start;
+                }
+
+                
+                ///the end pos for pre window is real_y_start - 1
+                total_y_end = real_y_start - 1;
+                ///find the unmatched window on the left of current matched window
+                ///k starts from i - 1
+                for (k = i - 1; k >= 0 && overlap_list->list[j].w_list[k].y_end == -1; k--)
+                {  
+                    ///there is no problem in x
+                    x_start = overlap_list->list[j].w_list[k].x_start;
+                    x_end = overlap_list->list[j].w_list[k].x_end;
+                    x_len = x_end - x_start + 1;
+                    ///there are two potiential reasons for unmatched window:
+                    ///1. this window has a large number of differences
+                    ///2. DP does not start from the right offset
+                    threshold = double_error_threshold(overlap_list->list[j].w_list[k].error_threshold, x_len);
+
+                    Window_Len = x_len + (threshold << 1);
+
+                    if(total_y_end <= 0)
+                    {
+                        break;
+                    }
+
+                    ///y_start might be less than 0
+                    y_start = total_y_end - x_len + 1;
+                    if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, R_INF,
+                    &extra_begin, &extra_end, &y_start, &o_len))
+                    {
+                        break;
+                    }
+
+                    if(o_len + threshold < x_len)
+                    {
+                        break;
+                    }
+
+                    fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, 
+                    R_INF, y_id, extra_begin, extra_end);
+                    x_string = g_read->seq + x_start;
+                    y_string = dumy->overlap_region;
+
+                    ///note!!! need notification
+                    end_site = Reserve_Banded_BPM_PATH(y_string, Window_Len, x_string, x_len, threshold, &error, &real_y_start,
+                    &(dumy->path_length), dumy->matrix_bit, dumy->path, -1, -1);
+
+                    if (error!=(unsigned int)-1)
+                    { 
+                        ///this condition is always wrong
+                        ///in best case, real_y_start = threshold, end_site = Window_Len - threshold - 1
+                        if (end_site == Window_Len - 1 || real_y_start == 0)
+                        {
+                            fix_boundary(x_string, x_len, threshold, y_start, real_y_start, end_site,
+                            extra_begin, extra_end, y_id, Window_Len, R_INF, dumy, y_strand, error,
+                            &y_start, &real_y_start, &end_site,
+                            &extra_begin, &extra_end, &error);
+                        }
+
+                        generate_cigar(dumy->path, dumy->path_length, &(overlap_list->list[j].w_list[k]),
+                        &real_y_start, &end_site, &error, x_string, x_len, y_string);  
+
+                        ///y_start has no shift, but y_end has shift           
+                        overlap_list->list[j].w_list[k].y_start = y_start + real_y_start - extra_begin;
+                        overlap_list->list[j].w_list[k].y_end = y_start + end_site;
+                        overlap_list->list[j].w_list[k].error = error;
+                        overlap_list->list[j].align_length += x_len;
+                        overlap_list->list[j].w_list[k].extra_begin = extra_begin;
+                        overlap_list->list[j].w_list[k].extra_end = extra_end;
+                        overlap_list->list[j].w_list[k].error_threshold = threshold;
+                    }
+                    else
+                    {
+                        break;
+                    }
+
+                    total_y_end = y_start + real_y_start - 1 - extra_begin;
+                }
+            }
+        }
+    }
+
+
+
+    overlap_list->mapped_overlaps_length = 0;
+    
+    double error_rate;
+    int is_update = 0;
+    for (j = 0; j < (long long)overlap_list->length; j++)
+    {
+        y_id = overlap_list->list[j].y_id;
+        y_strand = overlap_list->list[j].y_pos_strand;
+        y_readLen = Get_READ_LENGTH((*R_INF), y_id);
+        overlap_length = overlap_list->list[j].x_pos_e - overlap_list->list[j].x_pos_s + 1;
+        overlap_list->list[j].is_match = 0;
+        is_update = 0;
+        ///debug_scan_cigar(&(overlap_list->list[j]));
+        if(overlap_list->list[j].w_list_length == 0 || overlap_length == 0 || overlap_list->list[j].align_length == 0) continue;
+        ///only calculate cigar for high quality overlaps
+        if (overlap_length * OVERLAP_THRESHOLD_FILTER <=  overlap_list->list[j].align_length)
+        {
+            
+            for (i = 0; i < (long long)overlap_list->list[j].w_list_length; i++)
+            {
+                ///first we need to check if this window is matched
+                if(overlap_list->list[j].w_list[i].y_end != -1)
+                {
+                    ///second check if the cigar of this window has been got 
+                    if(overlap_list->list[j].w_list[i].cigar.length == -1)
+                    {
+                        ///there is no problem for x
+                        x_start = overlap_list->list[j].w_list[i].x_start;
+                        x_end = overlap_list->list[j].w_list[i].x_end;
+                        x_len = x_end - x_start + 1;
+                        /****************************may have bugs********************************/
+                        ///threshold = x_len * asm_opt.max_ov_diff_ec;
+                        threshold = overlap_list->list[j].w_list[i].error_threshold;
+                        /****************************may have bugs********************************/
+                        /****************************may have bugs********************************/
+                        ///should not adjust threshold, since this window can be matched by the old threshold
+                        ///threshold = Adjust_Threshold(threshold, x_len);
+                        /****************************may have bugs********************************/
+                        Window_Len = x_len + (threshold << 1);
+
+
+                        ///y_start is the real y_start
+                        ///for the window with cigar, y_start has already reduced extra_begin
+                        y_start = overlap_list->list[j].w_list[i].y_start;
+                        extra_begin = overlap_list->list[j].w_list[i].extra_begin;
+                        extra_end = overlap_list->list[j].w_list[i].extra_end;
+                        o_len = Window_Len - extra_end - extra_begin;
+                        fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, 
+                        R_INF, y_id, extra_begin, extra_end);
+                        x_string = g_read->seq + x_start;
+                        y_string = dumy->overlap_region;
+
+
+                        ///note!!! need notification
+                        end_site = Reserve_Banded_BPM_PATH(y_string, Window_Len, x_string, x_len, threshold, &error, &real_y_start,
+                        &(dumy->path_length), dumy->matrix_bit, dumy->path, 
+                        overlap_list->list[j].w_list[i].error, overlap_list->list[j].w_list[i].y_end - y_start);
+
+                        if (error != (unsigned int)-1)
+                        {
+                            if (end_site == Window_Len - 1 || real_y_start == 0)
+                            {
+                                
+                                if(fix_boundary(x_string, x_len, threshold, y_start, real_y_start, end_site,
+                                extra_begin, extra_end, y_id, Window_Len, R_INF, dumy, y_strand, error,
+                                &y_start, &real_y_start, &end_site,
+                                &extra_begin, &extra_end, &error))
+                                {
+                                    overlap_list->list[j].w_list[i].error = error;
+                                    overlap_list->list[j].w_list[i].extra_begin = extra_begin;
+                                    overlap_list->list[j].w_list[i].extra_end = extra_end;
+                                }
+                                
+                            }
+
+                            generate_cigar(dumy->path, dumy->path_length, &(overlap_list->list[j].w_list[i]),
+                            &real_y_start, &end_site, &error, x_string, x_len, y_string);    
+
+                            ///note!!! need notification
+                            real_y_start = y_start + real_y_start - extra_begin;
+                            overlap_list->list[j].w_list[i].y_start = real_y_start;  
+                            overlap_list->list[j].w_list[i].y_end = y_start + end_site - extra_begin;
+                            overlap_list->list[j].w_list[i].error = error;                              
+                        }
+                        else
+                        {
+                            fprintf(stderr, "error\n");
+                        }
+                    }
+                    else
+                    {
+                        overlap_list->list[j].w_list[i].y_end -= overlap_list->list[j].w_list[i].extra_begin;
+                    }
+
+
+                }
+            }
+
+            error_rate = non_trim_error_rate(overlap_list, j, R_INF, dumy, g_read);
+            
+            
+            if (error_rate <= HIGH_HET_ERROR_RATE)
+            {
+                is_update = 1;
+            }
+        }
+
+        if((is_update == 0) && (overlap_list->list[j].align_length >= WINDOW) && 
+           (overlap_length * HIGH_HET_OVERLAP_THRESHOLD_FILTER <=  overlap_list->list[j].align_length))
+        {
+            kv_resize(uint8_t, x_num, (uint64_t)(Get_READ_LENGTH((*R_INF), overlap_list->list[j].x_id)));
+            kv_resize(uint8_t, y_num, (uint64_t)(Get_READ_LENGTH((*R_INF), overlap_list->list[j].y_id)));
+            is_update = get_affine_gap_score(&(overlap_list->list[j]), g_read, overlap_read, x_num.a, y_num.a,
+            overlap_list->list[j].x_pos_e + 1 - overlap_list->list[j].x_pos_s, 
+            overlap_list->list[j].y_pos_e + 1 - overlap_list->list[j].y_pos_s);
+        }
+
+        if(is_update)
+        {
+            overlap_list->list[j].is_match = 2;
+        }
+    }
+
+    kv_destroy(x_num);
+    kv_destroy(y_num);
+}
+
+
+
+void correct_overlap_high_het(overlap_region_alloc* overlap_list, All_reads* R_INF, 
+                        UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read)
+{
+    clear_Correct_dumy(dumy, overlap_list);
+
+    long long window_start, window_end;
+
+    Window_Pool w_inf;
+
+    init_Window_Pool(&w_inf, g_read->length, WINDOW, (int)(1.0/asm_opt.max_ov_diff_ec));
+
+    int flag = 0;
+
+    while(get_Window(&w_inf, &window_start, &window_end) && flag != -2)
+    {
+        dumy->length = 0;
+        dumy->lengthNT = 0;
+        flag = get_interval(window_start, window_end, overlap_list, dumy);
+
+        switch (flag)
+        {
+            case 1:    ///no match here
+                break;
+            case 0:    ///no match here
+                break;
+            case -2: ///if flag == -2, loop would be terminated
+                break;
+        }
+
+        ///dumy->lengthNT represent how many overlaps that the length of them is not equal to WINDOW; may larger or less than WINDOW
+        ///dumy->length represent how many overlaps that the length of them is WINDOW
+        /****************************may improve**************************/
+        ///now the windows which are larger than WINDOW are verified one-by-one, to improve it, we can do it group-bygroup
+        verify_window(window_start, window_end, overlap_list, dumy, R_INF, g_read->seq);
+    }
+
+    recalcate_high_het_overlap(overlap_list, R_INF, g_read, dumy, overlap_read);
+}

Correct.h

@@ -1162,10 +1162,17 @@ void init_Cigar_record_alloc(Cigar_record_alloc* x);
 void resize_Cigar_record_alloc(Cigar_record_alloc* x, long long new_size);
 void destory_Cigar_record_alloc(Cigar_record_alloc* x);
 
-void afine_gap_alignment(const char *tseq, uint8_t* tnum, const int tl, 
-const char *qseq, uint8_t* qnum, const int ql, const uint8_t *c2n, const int strand,
+void afine_gap_alignment(const char *qseq, uint8_t* qnum, const int ql, 
+const char *tseq, uint8_t* tnum, const int tl, const uint8_t *c2n, const int strand,
 int sc_mch, int sc_mis, int gapo, int gape, int bandLen, int zdrop, int end_bonus,
-long long* max_t_pos, long long* max_q_pos, long long* score, long long* droped);
+long long* max_q_pos, long long* max_t_pos, long long* global_score, 
+long long* extention_score, long long* q_boundary_score, long long* q_boundary_t_coordinate,
+long long* t_boundary_score, long long* t_boundary_q_coordinate,
+long long* droped, int mode);
+void correct_overlap_high_het(overlap_region_alloc* overlap_list, All_reads* R_INF, 
+                        UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read);
+long long get_affine_gap_score(overlap_region* ovc, UC_Read* g_read, UC_Read* overlap_read, uint8_t* x_num,
+uint8_t* y_num, uint64_t EstimateXOlen, uint64_t EstimateYOlen);
 
 #define FORWARD_KSW 0
 #define BACKWARD_KSW 1
@@ -1174,5 +1181,5 @@ long long* max_t_pos, long long* max_q_pos, long long* score, long long* droped)
 #define GAP_OPEN_KSW 4
 #define GAP_EXT_KSW 2
 #define Z_DROP_KSW 400
-#define BAND_KSW 50
+#define BAND_KSW 500
 #endif

Hash_Table.h

@@ -12,6 +12,8 @@
 #define WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY 25
 #define THRESHOLD  15
 #define OVERLAP_THRESHOLD_FILTER 0.9
+#define HIGH_HET_OVERLAP_THRESHOLD_FILTER 0.3
+#define HIGH_HET_ERROR_RATE 0.08
 #define THRESHOLD_MAX_SIZE  31
 
 #define GROUP_SIZE 4

Makefile

@@ -4,7 +4,7 @@ CPPFLAGS=
 INCLUDES=
 OBJS=		CommandLines.o Process_Read.o Assembly.o Hash_Table.o \
 			POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o Purge_Dups.o \
-			htab.o hist.o sketch.o anchor.o extract.o sys.o
+			htab.o hist.o sketch.o anchor.o extract.o sys.o ksw2_extz2_sse.o
 EXE=		hifiasm
 LIBS=		-lz -lpthread -lm
 
@@ -36,7 +36,7 @@ Assembly.o: Assembly.h CommandLines.h Process_Read.h Overlaps.h kvec.h kdq.h
 Assembly.o: Hash_Table.h htab.h POA.h Correct.h Levenshtein_distance.h
 Assembly.o: kthread.h
 CommandLines.o: CommandLines.h ketopt.h
-Correct.o: Correct.h Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h
+Correct.o: Correct.h Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h ksw2.h
 Correct.o: kdq.h CommandLines.h Levenshtein_distance.h POA.h Assembly.h
 Hash_Table.o: Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h kdq.h
 Hash_Table.o: CommandLines.h ksort.h
@@ -65,3 +65,4 @@ main.o: CommandLines.h Process_Read.h Overlaps.h kvec.h kdq.h Assembly.h
 main.o: Levenshtein_distance.h htab.h
 sketch.o: kvec.h htab.h Process_Read.h Overlaps.h kdq.h CommandLines.h
 sys.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
+ksw2_extz2_sse.o: ksw2.h

Overlaps.cpp

@@ -9384,13 +9384,19 @@ ma_hit_t_alloc* sources, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp)
 
 
 uint32_t get_ug_coverage(ma_utg_t* u, asg_t* read_g, const ma_sub_t* coverage_cut, 
-ma_hit_t_alloc* sources, R_to_U* ruIndex)
+ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag)
 {
     uint32_t k, j, rId, tn, is_Unitig;
     long long R_bases = 0, C_bases = 0;
     ma_hit_t *h;
     if(u->m == 0) return 0;
 
+    for (k = 0; k < u->n; k++)
+    {
+        rId = u->a[k]>>33;
+        r_flag[rId] = 1;
+    }
+
     for (k = 0; k < u->n; k++)
     {
         rId = u->a[k]>>33;
@@ -9407,16 +9413,25 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex)
                 if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue;
             }
             if(read_g->seq[tn].del == 1) continue;
+            if(r_flag[tn] != 1) continue;
             C_bases += (Get_qe((*h)) - Get_qs((*h)));
         }
     }
 
+
+    for (k = 0; k < u->n; k++)
+    {
+        rId = u->a[k]>>33;
+        r_flag[rId] = 0;
+    }
+
     return C_bases/R_bases;
 }
 
 void ma_ug_print2(const ma_ug_t *ug, All_reads *RNF, asg_t* read_g, const ma_sub_t *coverage_cut, 
 ma_hit_t_alloc* sources, R_to_U* ruIndex, int print_seq, const char* prefix, FILE *fp)
 {
+    uint8_t* primary_flag = (uint8_t*)calloc(read_g->n_seq, sizeof(uint8_t));
 	uint32_t i, j, l;
 	char name[32];
 	for (i = 0; i < ug->u.n; ++i) { // the Segment lines in GFA
@@ -9424,9 +9439,9 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int print_seq, const char* prefix, FIL
         if(p->m == 0) continue;
         sprintf(name, "%s%.6d%c", prefix, i + 1, "lc"[p->circ]);
 		if (print_seq) fprintf(fp, "S\t%s\t%s\tLN:i:%d\trd:i:%u\n", name, p->s? p->s : "*", p->len, 
-        get_ug_coverage(p, read_g, coverage_cut, sources, ruIndex));
+        get_ug_coverage(p, read_g, coverage_cut, sources, ruIndex, primary_flag));
 		else fprintf(fp, "S\t%s\t*\tLN:i:%d\trd:i:%u\n", name, p->len, 
-        get_ug_coverage(p, read_g, coverage_cut, sources, ruIndex));
+        get_ug_coverage(p, read_g, coverage_cut, sources, ruIndex, primary_flag));
         // if (print_seq) fprintf(fp, "S\t%s\t%s\tLN:i:%d\n", name, p->s? p->s : "*", p->len);
 		// else fprintf(fp, "S\t%s\t*\tLN:i:%d\n", name, p->len);
         
@@ -9452,7 +9467,8 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int print_seq, const char* prefix, FIL
 		fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\n", 
         prefix, (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1],
 		prefix,	(v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], ug->g->arc[i].ol, asg_arc_len(ug->g->arc[i]));
-	}    
+	}
+    free(primary_flag);    
 }
 
 void ma_ug_print(const ma_ug_t *ug, All_reads *RNF, asg_t* read_g, const ma_sub_t *coverage_cut, 
@@ -11334,6 +11350,45 @@ const char* command)
 }
 
 
+uint32_t print_debug_gfa(asg_t *read_g, ma_ug_t *ug, ma_sub_t* coverage_cut, const char* output_file_name, 
+ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp)
+{
+    kvec_asg_arc_t_warp new_rtg_edges;
+    kv_init(new_rtg_edges.a);
+
+
+    uint32_t i;
+    for (i = 0; i < ug->u.n; ++i) 
+    {
+        ma_utg_t *u = &ug->u.a[i];
+        if(u->m == 0 || ug->g->seq[i].c == ALTER_LABLE)
+        {
+            asg_seq_del(ug->g, i);
+            if(ug->u.a[i].m!=0)
+            {
+                ug->u.a[i].m = ug->u.a[i].n = 0;
+                free(ug->u.a[i].a);
+                ug->u.a[i].a = NULL;
+            }
+        }
+    }
+
+    ma_ug_seq(ug, read_g, &R_INF, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp);
+
+    fprintf(stderr, "Writing raw unitig GFA to disk... \n");
+    char* gfa_name = (char*)malloc(strlen(output_file_name)+25);
+    sprintf(gfa_name, "%s.r_utg.noseq.gfa", output_file_name);
+    FILE* output_file = fopen(gfa_name, "w");
+    ma_ug_print_simple(ug, &R_INF, read_g, coverage_cut, sources, ruIndex, "utg", output_file);
+    fclose(output_file);
+
+    free(gfa_name);
+    ma_ug_destroy(ug);
+    kv_destroy(new_rtg_edges.a);
+    exit(0);
+}
+
+
 uint32_t print_untig_by_read(ma_ug_t *g, const char* name, uint32_t in, ma_hit_t_alloc* sources, 
 ma_hit_t_alloc* reverse_sources, const char* info)
 {
@@ -11346,6 +11401,7 @@ ma_hit_t_alloc* reverse_sources, const char* info)
     {
         for (i = 0; i < R_INF.total_reads; ++i) 
         {
+            if(Get_NAME_LENGTH(R_INF, i) != strlen(name)) continue;
             if(memcmp(name, Get_NAME(R_INF, i), Get_NAME_LENGTH(R_INF, i)) == 0)
             {
                 fprintf(stderr, "%s: i: %u, >%.*s\n", info, i, 
@@ -11364,7 +11420,7 @@ ma_hit_t_alloc* reverse_sources, const char* info)
         return (uint32_t)-1;
     }
 
-    print_specfic_read_ovlp(rId, sources, reverse_sources, info); 
+    if(sources && reverse_sources) print_specfic_read_ovlp(rId, sources, reverse_sources, info); 
 
     if(g != NULL)
     {
@@ -11376,7 +11432,8 @@ ma_hit_t_alloc* reverse_sources, const char* info)
             {
                 if(rId == (u->a[k]>>33))
                 {
-                    fprintf(stderr, "%s: %s is the %u-th read at %u-th unitig\n", info, name, k, i);
+                    fprintf(stderr, "%s: %s is the %u-th read at %u-th unitig (label: %u)\n", 
+                    info, name, k, i, g->g->seq[i].c);
                     return i;
                 }
             }
@@ -12165,10 +12222,11 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag)
         return 0;
     }
     long long R_bases = 0, C_bases = 0, C_bases_primary = 0, C_bases_alter = 0;
+    long long total_C_bases = 0, total_C_bases_primary = 0;
     uint32_t available_reads = 0, k, j, rId, tn, is_Unitig;
     ma_hit_t *h;
     if(u->m == 0) return 0;
-    available_reads = 0;
+    total_C_bases = total_C_bases_primary = available_reads = 0;
     for (k = 0; k < u->n; k++)
     {
         rId = u->a[k]>>33;
@@ -12202,15 +12260,15 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag)
             }
         }
 
-        ///fprintf(stderr, "C_bases_primary: %lld, C_bases_alter: %lld\n", C_bases_primary, C_bases_alter);
-
         C_bases = C_bases_primary + C_bases_alter;
-        if(C_bases_primary < C_bases * ALTER_COV_THRES) continue;
+        total_C_bases += C_bases; 
+        ///if(C_bases_primary < C_bases * ALTER_COV_THRES) continue;
 
         C_bases = C_bases/R_bases;
         if(C_bases >= asm_opt.recover_atg_cov_min && C_bases <= asm_opt.recover_atg_cov_max)
         {
-            available_reads++;
+            if(C_bases_primary >= (C_bases_primary + C_bases_alter) * ALTER_COV_THRES) available_reads++;
+            total_C_bases_primary += C_bases_primary;
         }
     }
 
@@ -12222,13 +12280,14 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag)
 
     ///fprintf(stderr, "available_reads: %u, u->n: %u\n", available_reads, u->n);
 
-    if(available_reads < (u->n * 0.8) || available_reads == 0)
+    //if(available_reads < (u->n * 0.8) || available_reads == 0)
+    if(((available_reads < (u->n * 0.8)) && (total_C_bases_primary < (total_C_bases * 0.8))) 
+       || available_reads == 0)
     {
         return 0;
     }
     else
     {
-        ///fprintf(stderr, "*****************\n");
         return 1;
     }
 }
@@ -12260,12 +12319,39 @@ R_to_U* ruIndex, uint32_t positive_flag, float drop_rate)
             nsu = &(ug->u.a[av[i].v>>1]);
             get_unitig_trio_flag(nsu, flag, &flag_occ, &non_flag_occ, &ambigious);
             ///we may need it or not
-            if((flag_occ <= ((non_flag_occ+flag_occ)*drop_rate))||
-                                                (flag_occ+non_flag_occ) == 0)
+            if((flag_occ <= ((non_flag_occ+flag_occ)*drop_rate))||((flag_occ+non_flag_occ) == 0))
             {
                 keep_node++;
                 continue;
-            } 
+            }
+
+            if(nsu->n >= 100)
+            {
+                if(flag_occ < nsu->n*DOUBLE_CHECK_THRES)
+                {
+                    keep_node++;
+                    continue;
+                }
+            }
+            else if(nsu->n >= 50)
+            {
+                if(flag_occ < nsu->n*DOUBLE_CHECK_THRES*0.5)
+                {
+                    keep_node++;
+                    continue;
+                }
+            }
+            else
+            {
+                if(flag_occ < nsu->n*DOUBLE_CHECK_THRES*0.25)
+                {
+                    keep_node++;
+                    continue;
+                }
+            }
+            
+            
+
             del_node++; 
         }
 
@@ -12278,12 +12364,33 @@ R_to_U* ruIndex, uint32_t positive_flag, float drop_rate)
 
             get_unitig_trio_flag(nsu, flag, &flag_occ, &non_flag_occ, &ambigious);
             ///we may need it or not
-            if((flag_occ <= ((non_flag_occ+flag_occ)*drop_rate))||
-                                                (flag_occ+non_flag_occ) == 0)
+            if((flag_occ <= ((non_flag_occ+flag_occ)*drop_rate))||((flag_occ+non_flag_occ) == 0))
             {
                 continue;
             }
 
+            if(nsu->n >= 100)
+            {
+                if(flag_occ < nsu->n*DOUBLE_CHECK_THRES)
+                {
+                    continue;
+                }
+            }
+            else if(nsu->n >= 50)
+            {
+                if(flag_occ < nsu->n*DOUBLE_CHECK_THRES*0.5)
+                {
+                    continue;
+                }
+            }
+            else
+            {
+                if(flag_occ < nsu->n*DOUBLE_CHECK_THRES*0.25)
+                {
+                    continue;
+                }
+            }
+
             if(if_primary_unitig(nsu, read_sg, coverage_cut, sources, ruIndex, primary_flag))
             {
                 continue;
@@ -12435,7 +12542,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, uint32_t stops_threshold
 R_to_U* ruIndex)
 {
     double startTime = Get_T();
-    uint32_t i, k, v_i, v_beg, v_end, selfLen, w1, w2, wv, nw, n_vtx = g->n_seq * 2, n_reduced = 0, convex, convex_T;
+    uint32_t i, k, v_i, v_beg, v_end, selfLen, w1, w2, wv, nw, n_vtx = g->n_seq * 2, n_reduced = 0, convex, convex_T, read_num;
     asg_arc_t *aw;
     long long ll, tmp, max_stop_nodeLen, max_stop_baseLen;
 
@@ -12508,7 +12615,7 @@ R_to_U* ruIndex)
             ///if(convex == convex_T) break;
             for (k = 0; k < b_0.b.n; k++)
             {
-                if(b_0.b.a[k] == (convex_T>>1)) break;
+                if((b_0.b.a[k]>>1) == (convex_T>>1)) break;
             }
             if(k != b_0.b.n) break;
             
@@ -12554,7 +12661,7 @@ R_to_U* ruIndex)
             ///if(convex == convex_T) break;
             for (k = 0; k < b_0.b.n; k++)
             {
-                if(b_0.b.a[k] == (convex_T>>1)) break;
+                if((b_0.b.a[k]>>1) == (convex_T>>1)) break;
             }
             if(k != b_0.b.n) break;
 
@@ -12566,17 +12673,15 @@ R_to_U* ruIndex)
         }
         if(i!=nw) continue;
 
-
-
-
-
         n_reduced++;
         b_0.b.n = 0;
+        read_num = 0;
         ///get_long_tip_length(g, &(ug->u), v_beg^1, &v_end, &b_0);
         get_unitig(g, ug, v_beg^1, &v_end, &ll, &tmp, &max_stop_nodeLen, &max_stop_baseLen, 1, &b_0);
         for (k = 0; k < b_0.b.n; k++)
         {
             g->seq[b_0.b.a[k]>>1].c = ALTER_LABLE;
+            read_num += ug->u.a[b_0.b.a[k]>>1].n;
         }
 
         for (k = 0; k < b_0.b.n; k++)
@@ -12584,6 +12689,14 @@ R_to_U* ruIndex)
             asg_seq_drop(g, b_0.b.a[k]>>1);
         }
 
+        if(read_num <= CHIMERIC_TRIM_THRES)
+        {
+            for (k = 0; k < b_0.b.n; k++)
+            {
+                g->seq[b_0.b.a[k]>>1].c = CUT;
+            }
+        }
+
     }
 
     asg_cleanup(g);
@@ -12756,12 +12869,12 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate)
     uint32_t is_first = 1;
 
     redo:
-    ///if(trio_flag == MOTHER) fprintf(stderr, "(0) c: %u, del: %u, n: %u\n", ug->g->seq[28141].c, ug->g->seq[28141].del, ug->u.a[28141].n);
+    ///print_untig((ug), 61955, "i-0:", 0);
+
     asg_pop_bubble_primary_trio(ug, bubble_dist, trio_flag, DROP);
     untig_asg_arc_simple_large_bubbles_trio(ug, read_g, reverse_sources, 2, ruIndex, trio_flag, DROP); 
-    magic_trio_phasing(g, ug, read_g, coverage_cut, sources, reverse_sources, 2, ruIndex, trio_flag, trio_drop_rate);
+    magic_trio_phasing(g, ug, read_g, coverage_cut, sources, reverse_sources, 2, ruIndex, trio_flag, trio_drop_rate);    
     ///drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex);
-    
     /**********debug**********/
     if(just_bubble_pop == 0)
     {
@@ -12769,11 +12882,6 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate)
         cut_trio_tip_primary(g, ug, tipsLen, trio_flag, 0, read_g, reverse_sources, ruIndex, 
         2);
     }
-
-    // print_read_all(ug, "m64011_190901_095311/1442493/ccs", "**convex**");
-    // print_read_all(ug, "m64011_190830_220126/173672768/ccs", "**kept**");
-    // print_read_all(ug, "m64011_190901_095311/6621017/ccs", "**missed**");
-
     /**********debug**********/
     long long pre_cons = get_graph_statistic(g);
     long long cur_cons = 0;
@@ -12787,17 +12895,13 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate)
         if(just_bubble_pop == 0)
         {
             ///need consider tangles
-            asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 
-            2, tip_drop_ratio);
-
-            asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, trio_flag);
-            
+            asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio);            
+            asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, trio_flag);            
             asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex,
             2, tip_drop_ratio, stops_threshold);
-
             asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 
             2, ruIndex, stops_threshold);
-
+            ///print_debug_gfa(read_g, ug, coverage_cut, "debug_chimeric", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len);
             detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate,
             ruIndex);
             ///need consider tangles
@@ -12815,9 +12919,9 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate)
         2);
     }
 
-    resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, trio_flag, drop_ratio);
+    resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, trio_flag, drop_ratio);    
     drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex);
-    all_to_all_deduplicate(ug, read_g, coverage_cut, sources, trio_flag, trio_drop_rate, reverse_sources, ruIndex);
+    all_to_all_deduplicate(ug, read_g, coverage_cut, sources, trio_flag, trio_drop_rate, reverse_sources, ruIndex, DOUBLE_CHECK_THRES);
 
     if(is_first)
     {
@@ -12928,7 +13032,7 @@ void set_drop_trio_flag(ma_ug_t *ug)
 
 void update_unitig_graph(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut, 
 ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, 
-uint8_t is_double_check, uint8_t flag, float drop_rate)
+uint8_t is_final_check, float double_check_rate, uint8_t flag, float drop_rate)
 {
     asg_t* nsg = ug->g;
     uint32_t v, n_vtx = nsg->n_seq, k, rId, flag_occ, non_flag_occ, hap_label_occ, n_reduce = 1;
@@ -12961,14 +13065,44 @@ uint8_t is_double_check, uint8_t flag, float drop_rate)
             }
 
             if(hap_label_occ == u->n) continue;
-            if(is_double_check && non_flag_occ < u->n*DOUBLE_CHECK_THRES) continue;
+            ///if(is_double_check && non_flag_occ < u->n*DOUBLE_CHECK_THRES) continue;
+            ///if(is_double_check && non_flag_occ < u->n*double_check_rate) continue;
+            if(is_final_check)
+            {
+                if(non_flag_occ < u->n*double_check_rate) continue;
+            }
+            else
+            {
+                if(u->n >= 100)
+                {
+                    if(non_flag_occ < u->n*double_check_rate)
+                    {
+                        continue;
+                    }
+                }
+                else if(u->n >= 50)
+                {
+                    if(non_flag_occ < u->n*double_check_rate*0.5)
+                    {
+                        continue;
+                    }
+                }
+                else
+                {
+                    if(non_flag_occ < u->n*double_check_rate*0.25)
+                    {
+                        continue;
+                    }
+                }
+            }
+            
 
             if(non_flag_occ > ((non_flag_occ+flag_occ)*drop_rate))
             {
                 if(if_primary_unitig(u, read_g, coverage_cut, sources, ruIndex, primary_flag))
                 {
                     continue;
-                }
+                }                         
                 if(u->m != 0)
                 {
                     u->circ = u->end = u->len = u->m = u->n = u->start = 0;
@@ -13109,7 +13243,7 @@ uint32_t* non_require, uint32_t* ambigious)
 ///note: to use this function, don't renew unitig graph!!!!!!!!!
 void all_to_all_deduplicate(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut, 
 ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, 
-ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex)
+ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate)
 {
     
 
@@ -13170,6 +13304,27 @@ ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex)
                 ///we may need it or not
                 if(flag_occ <= ((non_flag_occ+flag_occ)*drop_rate)) continue;
                 if((flag_occ+non_flag_occ) == 0) continue;
+                if(nsu->n >= 100)
+                {
+                    if(flag_occ < nsu->n*double_check_rate)
+                    {
+                        continue;
+                    }
+                }
+                else if(nsu->n >= 50)
+                {
+                    if(flag_occ < nsu->n*double_check_rate*0.5)
+                    {
+                        continue;
+                    }
+                }
+                else
+                {
+                    if(flag_occ < nsu->n*double_check_rate*0.25)
+                    {
+                        continue;
+                    }
+                }
             }
             
 
@@ -13493,9 +13648,6 @@ kvec_asg_arc_t_warp* new_rtg_edges)
 {
     asg_t* nsg = (*ug)->g;
     uint32_t v, n_vtx = nsg->n_seq;    
-    //if(flag == MOTHER) print_untig_by_read(*ug, "m54329U_190617_231905/65340614/ccs", (uint32_t)-1, sources, reverse_sources, "beg1");
-    //if(flag == MOTHER) print_untig_by_read(*ug, "m54329U_190827_173812/67108993/ccs", (uint32_t)-1, sources, reverse_sources, "beg2");
-    
     /**
     kvec_t_u32_warp new_rtg_nodes;
     kv_init(new_rtg_nodes.a);
@@ -13524,10 +13676,11 @@ kvec_asg_arc_t_warp* new_rtg_edges)
     
 
     ///primary_flag = get_utg_attributes(*ug, read_g, coverage_cut, sources, ruIndex);
+    update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, 0, 
+    DOUBLE_CHECK_THRES, flag, drop_rate);
 
-
-    update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, 0, flag, drop_rate);
     adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex);
+
     nsg = (*ug)->g;
     n_vtx = nsg->n_seq;
     for (v = 0; v < n_vtx; ++v) 
@@ -13536,7 +13689,6 @@ kvec_asg_arc_t_warp* new_rtg_edges)
         nsg->seq[v].c = PRIMARY_LABLE;
         EvaluateLen((*ug)->u, v) = (*ug)->u.a[v].n;
     }
-
     clean_trio_untig_graph(*ug, read_g, coverage_cut, sources, reverse_sources, bubble_dist, 
     tipsLen, tip_drop_ratio, stops_threshold, ruIndex, NULL, NULL, 0, 0, 0, 
     chimeric_rate, 0, 0, drop_ratio, flag, drop_rate);
@@ -13548,7 +13700,8 @@ kvec_asg_arc_t_warp* new_rtg_edges)
 
     update_hap_label(*ug, read_g);
 
-    update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, 0, flag, drop_rate);
+    update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, 0, 
+    DOUBLE_CHECK_THRES, flag, drop_rate);
 
     renew_utg(ug, read_g, new_rtg_edges);
 
@@ -13565,7 +13718,8 @@ kvec_asg_arc_t_warp* new_rtg_edges)
         renew_utg(ug, read_g, new_rtg_edges);
     }
     
-    update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, 1, flag, drop_rate);
+    update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, 1, 
+    FINAL_DOUBLE_CHECK_THRES, flag, drop_rate);
 
     update_hap_label(NULL, read_g);
 
@@ -22432,14 +22586,18 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex)
     uint32_t v, n_vtx = nsg->n_seq, k, j, rId, available_reads = 0, keep_atg = 0, tn, is_Unitig;
     ma_utg_t* u = NULL;
     ma_hit_t *h;
+
+    
+    ///print_untig_by_read(atg, "SRR11606870.634978", -1, NULL, NULL, "debug");
     long long R_bases = 0, C_bases = 0, C_bases_primary = 0, C_bases_alter = 0;
+    long long total_C_bases = 0, total_C_bases_alter = 0;
     for (v = 0; v < n_vtx; ++v) 
     {
         if(nsg->seq[v].del) continue;
         u = &(atg->u.a[v]);
         if(u->m == 0) continue;
-        available_reads = 0;
-
+        total_C_bases = total_C_bases_alter = available_reads = 0;
+        
         for (k = 0; k < u->n; k++)
         {
             rId = u->a[k]>>33;
@@ -22468,16 +22626,19 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex)
             }
 
             C_bases = C_bases_primary + C_bases_alter;
-            if(C_bases_alter < C_bases * ALTER_COV_THRES) continue;
+            total_C_bases += C_bases; 
 
             C_bases = C_bases/R_bases;
             if(C_bases >= asm_opt.recover_atg_cov_min && C_bases <= asm_opt.recover_atg_cov_max)
             {
-                available_reads++;
+                if(C_bases_alter >= (C_bases_primary + C_bases_alter) * ALTER_COV_THRES) available_reads++;
+                total_C_bases_alter += C_bases_alter;
             }
         }
 
-        if(available_reads < (u->n * 0.8) || available_reads == 0)
+
+        if(((available_reads < (u->n * 0.8)) && (total_C_bases_alter < (total_C_bases * 0.8))) 
+           || available_reads == 0)
         {
             asg_seq_del(nsg, v);
 
@@ -26642,9 +26803,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
     normalize_ma_hit_t_single_side_advance(sources, n_read);
     normalize_ma_hit_t_single_side_advance(reverse_sources, n_read);
     
-    // debug_info_of_specfic_read("m64062_190803_042216/122882911/ccs", sources, reverse_sources, 
-    // -1, "clean");
-
 
     if (ha_opt_triobin(&asm_opt))
     {
@@ -26670,6 +26828,9 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
     ///just need to deal with trio here
     ma_hit_contained_advance(sources, n_read, coverage_cut, ruIndex, max_hang_length, mini_overlap_length);
 
+    // debug_info_of_specfic_read("m54329U_190617_231905/176226420/ccs", sources, reverse_sources, 
+    // -1, "clean");
+
     sg = ma_sg_gen(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length);
     asg_arc_del_trans(sg, gap_fuzz);
     asm_opt.coverage = get_coverage(sources, coverage_cut, n_read);

Overlaps.h

@@ -22,7 +22,9 @@
 #define ALTER_LABLE 1
 #define HAP_LABLE 2
 #define TRIO_THRES 0.9
-#define DOUBLE_CHECK_THRES 0.2
+#define DOUBLE_CHECK_THRES 0.1
+#define FINAL_DOUBLE_CHECK_THRES 0.2
+#define CHIMERIC_TRIM_THRES 4
 // #define PRIMARY_LABLE 1
 // #define ALTER_LABLE 2
 // #define HAP_LABLE 4
@@ -1048,7 +1050,7 @@ uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges);
 void deduplicate(ma_ug_t *src, asg_t *read_g, ma_hit_t_alloc* reverse_sources, long long minLongUntig, 
 long long maxShortUntig, float l_untig_rate, float max_node_threshold, R_to_U* ruIndex, uint32_t resolve_tangle);
 void all_to_all_deduplicate(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut, 
-ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
+ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate);
 void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
 void rescue_wrong_overlaps_to_unitigs(ma_ug_t *i_ug, asg_t *r_g,  ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
 ma_sub_t *coverage_cut, R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, kvec_asg_arc_t_warp* keep_edges);

anchor.cpp

@@ -174,3 +174,10 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 
 	ks_introsort_or_xs(overlap_list->length, overlap_list->list);
 }
+
+
+
+void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list)
+{
+	ks_introsort_or_xs(overlap_list->length, overlap_list->list);
+}

ksw2.h

@@ -0,0 +1,177 @@
+#ifndef KSW2_H_
+#define KSW2_H_
+
+#include <stdint.h>
+
+#define KSW_NEG_INF -0x40000000
+
+#define KSW_EZ_SCORE_ONLY  0x01 // don't record alignment path/cigar
+#define KSW_EZ_RIGHT       0x02 // right-align gaps
+#define KSW_EZ_GENERIC_SC  0x04 // without this flag: match/mismatch only; last symbol is a wildcard
+#define KSW_EZ_APPROX_MAX  0x08 // approximate max; this is faster with sse
+#define KSW_EZ_APPROX_DROP 0x10 // approximate Z-drop; faster with sse
+#define KSW_EZ_EXTZ_ONLY   0x40 // only perform extension
+#define KSW_EZ_REV_CIGAR   0x80 // reverse CIGAR in the output
+#define KSW_EZ_SPLICE_FOR  0x100
+#define KSW_EZ_SPLICE_REV  0x200
+#define KSW_EZ_SPLICE_FLANK 0x400
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef struct {
+	uint32_t max:31, zdropped:1;
+	int max_q, max_t;      // max extension coordinate
+	int mqe, mqe_t;        // max score when reaching the end of query
+	int mte, mte_q;        // max score when reaching the end of target
+	int score;             // max score reaching both ends; may be KSW_NEG_INF
+	int m_cigar, n_cigar;
+	int reach_end;
+	uint32_t *cigar;
+} ksw_extz_t;
+
+/**
+ * NW-like extension
+ *
+ * @param km        memory pool, when used with kalloc
+ * @param qlen      query length
+ * @param query     query sequence with 0 <= query[i] < m
+ * @param tlen      target length
+ * @param target    target sequence with 0 <= target[i] < m
+ * @param m         number of residue types
+ * @param mat       m*m scoring mattrix in one-dimension array
+ * @param gapo      gap open penalty; a gap of length l cost "-(gapo+l*gape)"
+ * @param gape      gap extension penalty
+ * @param w         band width (<0 to disable)
+ * @param zdrop     off-diagonal drop-off to stop extension (positive; <0 to disable)
+ * @param flag      flag (see KSW_EZ_* macros)
+ * @param ez        (out) scores and cigar
+ */
+void ksw_extz(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 flag, ksw_extz_t *ez);
+
+void ksw_extz2_sse(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);
+
+void ksw_extd(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
+			  int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int flag, ksw_extz_t *ez);
+
+void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
+				   int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
+
+void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
+				   int8_t gapo, int8_t gape, int8_t gapo2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez);
+
+void ksw_extf2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t mch, int8_t mis, int8_t e, int w, int xdrop, ksw_extz_t *ez);
+
+/**
+ * Global alignment
+ *
+ * (first 10 parameters identical to ksw_extz_sse())
+ * @param m_cigar   (modified) max CIGAR length; feed 0 if cigar==0
+ * @param n_cigar   (out) number of CIGAR elements
+ * @param cigar     (out) BAM-encoded CIGAR; caller need to deallocate with kfree(km, )
+ *
+ * @return          score of the alignment
+ */
+int ksw_gg(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_);
+int ksw_gg2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_);
+int ksw_gg2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_);
+
+void *ksw_ll_qinit(void *km, int size, int qlen, const uint8_t *query, int m, const int8_t *mat);
+int ksw_ll_i16(void *q, int tlen, const uint8_t *target, int gapo, int gape, int *qe, int *te);
+
+#ifdef __cplusplus
+}
+#endif
+
+/************************************
+ *** Private macros and functions ***
+ ************************************/
+
+#ifdef HAVE_KALLOC
+#include "kalloc.h"
+#else
+#include <stdlib.h>
+#define kmalloc(km, size) malloc((size))
+#define kcalloc(km, count, size) calloc((count), (size))
+#define krealloc(km, ptr, size) realloc((ptr), (size))
+#define kfree(km, ptr) free((ptr))
+#endif
+
+static inline uint32_t *ksw_push_cigar(void *km, int *n_cigar, int *m_cigar, uint32_t *cigar, uint32_t op, int len)
+{
+	if (*n_cigar == 0 || op != (cigar[(*n_cigar) - 1]&0xf)) {
+		if (*n_cigar == *m_cigar) {
+			*m_cigar = *m_cigar? (*m_cigar)<<1 : 4;
+			cigar = (uint32_t*)krealloc(km, cigar, (*m_cigar) << 2);
+		}
+		cigar[(*n_cigar)++] = len<<4 | op;
+	} else cigar[(*n_cigar)-1] += len<<4;
+	return cigar;
+}
+
+// In the backtrack matrix, value p[] has the following structure:
+//   bit 0-2: which type gets the max - 0 for H, 1 for E, 2 for F, 3 for \tilde{E} and 4 for \tilde{F}
+//   bit 3/0x08: 1 if a continuation on the E state (bit 5/0x20 for a continuation on \tilde{E})
+//   bit 4/0x10: 1 if a continuation on the F state (bit 6/0x40 for a continuation on \tilde{F})
+static inline void ksw_backtrack(void *km, int is_rot, int is_rev, int min_intron_len, const uint8_t *p, const int *off, const int *off_end, int n_col, int i0, int j0,
+								 int *m_cigar_, int *n_cigar_, uint32_t **cigar_)
+{ // p[] - lower 3 bits: which type gets the max; bit
+	int n_cigar = 0, m_cigar = *m_cigar_, i = i0, j = j0, r, state = 0;
+	uint32_t *cigar = *cigar_, tmp;
+	while (i >= 0 && j >= 0) { // at the beginning of the loop, _state_ tells us which state to check
+		int force_state = -1;
+		if (is_rot) {
+			r = i + j;
+			if (i < off[r]) force_state = 2;
+			if (off_end && i > off_end[r]) force_state = 1;
+			tmp = force_state < 0? p[(size_t)r * n_col + i - off[r]] : 0;
+		} else {
+			if (j < off[i]) force_state = 2;
+			if (off_end && j > off_end[i]) force_state = 1;
+			tmp = force_state < 0? p[(size_t)i * n_col + j - off[i]] : 0;
+		}
+		if (state == 0) state = tmp & 7; // if requesting the H state, find state one maximizes it.
+		else if (!(tmp >> (state + 2) & 1)) state = 0; // if requesting other states, _state_ stays the same if it is a continuation; otherwise, set to H
+		if (state == 0) state = tmp & 7; // TODO: probably this line can be merged into the "else if" line right above; not 100% sure
+		if (force_state >= 0) state = force_state;
+		if (state == 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 0, 1), --i, --j; // match
+		else if (state == 1 || (state == 3 && min_intron_len <= 0)) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 2, 1), --i; // deletion
+		else if (state == 3 && min_intron_len > 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 3, 1), --i; // intron
+		else cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, 1), --j; // insertion
+	}
+	if (i >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, min_intron_len > 0 && i >= min_intron_len? 3 : 2, i + 1); // first deletion
+	if (j >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, j + 1); // first insertion
+	if (!is_rev)
+		for (i = 0; i < n_cigar>>1; ++i) // reverse CIGAR
+			tmp = cigar[i], cigar[i] = cigar[n_cigar-1-i], cigar[n_cigar-1-i] = tmp;
+	*m_cigar_ = m_cigar, *n_cigar_ = n_cigar, *cigar_ = cigar;
+}
+
+static inline void ksw_reset_extz(ksw_extz_t *ez)
+{
+	ez->max_q = ez->max_t = ez->mqe_t = ez->mte_q = -1;
+	ez->max = 0, ez->score = ez->mqe = ez->mte = KSW_NEG_INF;
+	ez->n_cigar = 0, ez->zdropped = 0, ez->reach_end = 0;
+}
+
+static inline int ksw_apply_zdrop(ksw_extz_t *ez, int is_rot, int32_t H, int a, int b, int zdrop, int8_t e)
+{
+	int r, t;
+	if (is_rot) r = a, t = b;
+	else r = a + b, t = a;
+	if (H > (int32_t)ez->max) {
+		ez->max = H, ez->max_t = t, ez->max_q = r - t;
+	} else if (t >= ez->max_t && r - t >= ez->max_q) {
+		int tl = t - ez->max_t, ql = (r - t) - ez->max_q, l;
+		l = tl > ql? tl - ql : ql - tl;
+		if (zdrop >= 0 && ez->max - H > zdrop + l * e) {
+			ez->zdropped = 1;
+			return 1;
+		}
+	}
+	return 0;
+}
+#endif

ksw2_extz2_sse.c

@@ -0,0 +1,305 @@
+#include <string.h>
+#include <assert.h>
+#include "ksw2.h"
+
+#ifdef __SSE2__
+#include <emmintrin.h>
+
+#ifdef KSW_SSE2_ONLY
+#undef __SSE4_1__
+#endif
+
+#ifdef __SSE4_1__
+#include <smmintrin.h>
+#endif
+
+#ifdef KSW_CPU_DISPATCH
+#ifdef __SSE4_1__
+void ksw_extz2_sse41(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)
+#else
+void ksw_extz2_sse2(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)
+#endif
+#else
+void ksw_extz2_sse(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)
+#endif // ~KSW_CPU_DISPATCH
+{
+#define __dp_code_block1 \
+	z = _mm_add_epi8(_mm_load_si128(&s[t]), qe2_); \
+	xt1 = _mm_load_si128(&x[t]);                     /* xt1 <- x[r-1][t..t+15] */ \
+	tmp = _mm_srli_si128(xt1, 15);                   /* tmp <- x[r-1][t+15] */ \
+	xt1 = _mm_or_si128(_mm_slli_si128(xt1, 1), x1_); /* xt1 <- x[r-1][t-1..t+14] */ \
+	x1_ = tmp; \
+	vt1 = _mm_load_si128(&v[t]);                     /* vt1 <- v[r-1][t..t+15] */ \
+	tmp = _mm_srli_si128(vt1, 15);                   /* tmp <- v[r-1][t+15] */ \
+	vt1 = _mm_or_si128(_mm_slli_si128(vt1, 1), v1_); /* vt1 <- v[r-1][t-1..t+14] */ \
+	v1_ = tmp; \
+	a = _mm_add_epi8(xt1, vt1);                      /* a <- x[r-1][t-1..t+14] + v[r-1][t-1..t+14] */ \
+	ut = _mm_load_si128(&u[t]);                      /* ut <- u[t..t+15] */ \
+	b = _mm_add_epi8(_mm_load_si128(&y[t]), ut);     /* b <- y[r-1][t..t+15] + u[r-1][t..t+15] */
+
+#define __dp_code_block2 \
+	z = _mm_max_epu8(z, b);                          /* z = max(z, b); this works because both are non-negative */ \
+	z = _mm_min_epu8(z, max_sc_); \
+	_mm_store_si128(&u[t], _mm_sub_epi8(z, vt1));    /* u[r][t..t+15] <- z - v[r-1][t-1..t+14] */ \
+	_mm_store_si128(&v[t], _mm_sub_epi8(z, ut));     /* v[r][t..t+15] <- z - u[r-1][t..t+15] */ \
+	z = _mm_sub_epi8(z, q_); \
+	a = _mm_sub_epi8(a, z); \
+	b = _mm_sub_epi8(b, z);
+
+	int r, t, qe = q + e, n_col_, *off = 0, *off_end = 0, tlen_, qlen_, last_st, last_en, wl, wr, max_sc, min_sc;
+	int with_cigar = !(flag&KSW_EZ_SCORE_ONLY), approx_max = !!(flag&KSW_EZ_APPROX_MAX);
+	int32_t *H = 0, H0 = 0, last_H0_t = 0;
+	uint8_t *qr, *sf, *mem, *mem2 = 0;
+	__m128i q_, qe2_, zero_, flag1_, flag2_, flag8_, flag16_, sc_mch_, sc_mis_, sc_N_, m1_, max_sc_;
+	__m128i *u, *v, *x, *y, *s, *p = 0;
+
+	ksw_reset_extz(ez);
+	if (m <= 0 || qlen <= 0 || tlen <= 0) return;
+
+	zero_   = _mm_set1_epi8(0);
+	q_      = _mm_set1_epi8(q);
+	qe2_    = _mm_set1_epi8((q + e) * 2);
+	flag1_  = _mm_set1_epi8(1);
+	flag2_  = _mm_set1_epi8(2);
+	flag8_  = _mm_set1_epi8(0x08);
+	flag16_ = _mm_set1_epi8(0x10);
+	sc_mch_ = _mm_set1_epi8(mat[0]);
+	sc_mis_ = _mm_set1_epi8(mat[1]);
+	sc_N_   = mat[m*m-1] == 0? _mm_set1_epi8(-e) : _mm_set1_epi8(mat[m*m-1]);
+	m1_     = _mm_set1_epi8(m - 1); // wildcard
+	max_sc_ = _mm_set1_epi8(mat[0] + (q + e) * 2);
+
+	if (w < 0) w = tlen > qlen? tlen : qlen;
+	wl = wr = w;
+	tlen_ = (tlen + 15) / 16;
+	n_col_ = qlen < tlen? qlen : tlen;
+	n_col_ = ((n_col_ < w + 1? n_col_ : w + 1) + 15) / 16 + 1;
+	qlen_ = (qlen + 15) / 16;
+	for (t = 1, max_sc = mat[0], min_sc = mat[1]; t < m * m; ++t) {
+		max_sc = max_sc > mat[t]? max_sc : mat[t];
+		min_sc = min_sc < mat[t]? min_sc : mat[t];
+	}
+	if (-min_sc > 2 * (q + e)) return; // otherwise, we won't see any mismatches
+
+	mem = (uint8_t*)kcalloc(km, tlen_ * 6 + qlen_ + 1, 16);
+	u = (__m128i*)(((size_t)mem + 15) >> 4 << 4); // 16-byte aligned
+	v = u + tlen_, x = v + tlen_, y = x + tlen_, s = y + tlen_, sf = (uint8_t*)(s + tlen_), qr = sf + tlen_ * 16;
+	if (!approx_max) {
+		H = (int32_t*)kmalloc(km, tlen_ * 16 * 4);
+		for (t = 0; t < tlen_ * 16; ++t) H[t] = KSW_NEG_INF;
+	}
+	if (with_cigar) {
+		mem2 = (uint8_t*)kmalloc(km, ((size_t)(qlen + tlen - 1) * n_col_ + 1) * 16);
+		p = (__m128i*)(((size_t)mem2 + 15) >> 4 << 4);
+		off = (int*)kmalloc(km, (qlen + tlen - 1) * sizeof(int) * 2);
+		off_end = off + qlen + tlen - 1;
+	}
+
+	for (t = 0; t < qlen; ++t) qr[t] = query[qlen - 1 - t];
+	memcpy(sf, target, tlen);
+
+	for (r = 0, last_st = last_en = -1; r < qlen + tlen - 1; ++r) {
+		int st = 0, en = tlen - 1, st0, en0, st_, en_;
+		int8_t x1, v1;
+		uint8_t *qrr = qr + (qlen - 1 - r), *u8 = (uint8_t*)u, *v8 = (uint8_t*)v;
+		__m128i x1_, v1_;
+		// find the boundaries
+		if (st < r - qlen + 1) st = r - qlen + 1;
+		if (en > r) en = r;
+		if (st < (r-wr+1)>>1) st = (r-wr+1)>>1; // take the ceil
+		if (en > (r+wl)>>1) en = (r+wl)>>1; // take the floor
+		if (st > en) {
+			ez->zdropped = 1;
+			break;
+		}
+		st0 = st, en0 = en;
+		st = st / 16 * 16, en = (en + 16) / 16 * 16 - 1;
+		// set boundary conditions
+		if (st > 0) {
+			if (st - 1 >= last_st && st - 1 <= last_en)
+				x1 = ((uint8_t*)x)[st - 1], v1 = v8[st - 1]; // (r-1,s-1) calculated in the last round
+			else x1 = v1 = 0; // not calculated; set to zeros
+		} else x1 = 0, v1 = r? q : 0;
+		if (en >= r) ((uint8_t*)y)[r] = 0, u8[r] = r? q : 0;
+		// loop fission: set scores first
+		if (!(flag & KSW_EZ_GENERIC_SC)) {
+			for (t = st0; t <= en0; t += 16) {
+				__m128i sq, st, tmp, mask;
+				sq = _mm_loadu_si128((__m128i*)&sf[t]);
+				st = _mm_loadu_si128((__m128i*)&qrr[t]);
+				mask = _mm_or_si128(_mm_cmpeq_epi8(sq, m1_), _mm_cmpeq_epi8(st, m1_));
+				tmp = _mm_cmpeq_epi8(sq, st);
+#ifdef __SSE4_1__
+				tmp = _mm_blendv_epi8(sc_mis_, sc_mch_, tmp);
+				tmp = _mm_blendv_epi8(tmp,     sc_N_,   mask);
+#else
+				tmp = _mm_or_si128(_mm_andnot_si128(tmp,  sc_mis_), _mm_and_si128(tmp,  sc_mch_));
+				tmp = _mm_or_si128(_mm_andnot_si128(mask, tmp),     _mm_and_si128(mask, sc_N_));
+#endif
+				_mm_storeu_si128((__m128i*)((uint8_t*)s + t), tmp);
+			}
+		} else {
+			for (t = st0; t <= en0; ++t)
+				((uint8_t*)s)[t] = mat[sf[t] * m + qrr[t]];
+		}
+		// core loop
+		x1_ = _mm_cvtsi32_si128(x1);
+		v1_ = _mm_cvtsi32_si128(v1);
+		st_ = st / 16, en_ = en / 16;
+		assert(en_ - st_ + 1 <= n_col_);
+		if (!with_cigar) { // score only
+			for (t = st_; t <= en_; ++t) {
+				__m128i z, a, b, xt1, vt1, ut, tmp;
+				__dp_code_block1;
+#ifdef __SSE4_1__
+				z = _mm_max_epi8(z, a);                          // z = z > a? z : a (signed)
+#else // we need to emulate SSE4.1 intrinsics _mm_max_epi8()
+				z = _mm_and_si128(z, _mm_cmpgt_epi8(z, zero_));  // z = z > 0? z : 0;
+				z = _mm_max_epu8(z, a);                          // z = max(z, a); this works because both are non-negative
+#endif
+				__dp_code_block2;
+#ifdef __SSE4_1__
+				_mm_store_si128(&x[t], _mm_max_epi8(a, zero_));
+				_mm_store_si128(&y[t], _mm_max_epi8(b, zero_));
+#else
+				tmp = _mm_cmpgt_epi8(a, zero_);
+				_mm_store_si128(&x[t], _mm_and_si128(a, tmp));
+				tmp = _mm_cmpgt_epi8(b, zero_);
+				_mm_store_si128(&y[t], _mm_and_si128(b, tmp));
+#endif
+			}
+		} else if (!(flag&KSW_EZ_RIGHT)) { // gap left-alignment
+			__m128i *pr = p + (size_t)r * n_col_ - st_;
+			off[r] = st, off_end[r] = en;
+			for (t = st_; t <= en_; ++t) {
+				__m128i d, z, a, b, xt1, vt1, ut, tmp;
+				__dp_code_block1;
+				d = _mm_and_si128(_mm_cmpgt_epi8(a, z), flag1_); // d = a > z? 1 : 0
+#ifdef __SSE4_1__
+				z = _mm_max_epi8(z, a);                          // z = z > a? z : a (signed)
+				tmp = _mm_cmpgt_epi8(b, z);
+				d = _mm_blendv_epi8(d, flag2_, tmp);             // d = b > z? 2 : d
+#else // we need to emulate SSE4.1 intrinsics _mm_max_epi8() and _mm_blendv_epi8()
+				z = _mm_and_si128(z, _mm_cmpgt_epi8(z, zero_));  // z = z > 0? z : 0;
+				z = _mm_max_epu8(z, a);                          // z = max(z, a); this works because both are non-negative
+				tmp = _mm_cmpgt_epi8(b, z);
+				d = _mm_or_si128(_mm_andnot_si128(tmp, d), _mm_and_si128(tmp, flag2_)); // d = b > z? 2 : d; emulating blendv
+#endif
+				__dp_code_block2;
+				tmp = _mm_cmpgt_epi8(a, zero_);
+				_mm_store_si128(&x[t], _mm_and_si128(tmp, a));
+				d = _mm_or_si128(d, _mm_and_si128(tmp, flag8_));  // d = a > 0? 0x08 : 0
+				tmp = _mm_cmpgt_epi8(b, zero_);
+				_mm_store_si128(&y[t], _mm_and_si128(tmp, b));
+				d = _mm_or_si128(d, _mm_and_si128(tmp, flag16_)); // d = b > 0? 0x10 : 0
+				_mm_store_si128(&pr[t], d);
+			}
+		} else { // gap right-alignment
+			__m128i *pr = p + (size_t)r * n_col_ - st_;
+			off[r] = st, off_end[r] = en;
+			for (t = st_; t <= en_; ++t) {
+				__m128i d, z, a, b, xt1, vt1, ut, tmp;
+				__dp_code_block1;
+				d = _mm_andnot_si128(_mm_cmpgt_epi8(z, a), flag1_); // d = z > a? 0 : 1
+#ifdef __SSE4_1__
+				z = _mm_max_epi8(z, a);                          // z = z > a? z : a (signed)
+				tmp = _mm_cmpgt_epi8(z, b);
+				d = _mm_blendv_epi8(flag2_, d, tmp);             // d = z > b? d : 2
+#else // we need to emulate SSE4.1 intrinsics _mm_max_epi8() and _mm_blendv_epi8()
+				z = _mm_and_si128(z, _mm_cmpgt_epi8(z, zero_));  // z = z > 0? z : 0;
+				z = _mm_max_epu8(z, a);                          // z = max(z, a); this works because both are non-negative
+				tmp = _mm_cmpgt_epi8(z, b);
+				d = _mm_or_si128(_mm_andnot_si128(tmp, flag2_), _mm_and_si128(tmp, d)); // d = z > b? d : 2; emulating blendv
+#endif
+				__dp_code_block2;
+				tmp = _mm_cmpgt_epi8(zero_, a);
+				_mm_store_si128(&x[t], _mm_andnot_si128(tmp, a));
+				d = _mm_or_si128(d, _mm_andnot_si128(tmp, flag8_));  // d = 0 > a? 0 : 0x08
+				tmp = _mm_cmpgt_epi8(zero_, b);
+				_mm_store_si128(&y[t], _mm_andnot_si128(tmp, b));
+				d = _mm_or_si128(d, _mm_andnot_si128(tmp, flag16_)); // d = 0 > b? 0 : 0x10
+				_mm_store_si128(&pr[t], d);
+			}
+		}
+		if (!approx_max) { // find the exact max with a 32-bit score array
+			int32_t max_H, max_t;
+			// compute H[], max_H and max_t
+			if (r > 0) {
+				int32_t HH[4], tt[4], en1 = st0 + (en0 - st0) / 4 * 4, i;
+				__m128i max_H_, max_t_, qe_;
+				max_H = H[en0] = en0 > 0? H[en0-1] + u8[en0] - qe : H[en0] + v8[en0] - qe; // special casing the last element
+				max_t = en0;
+				max_H_ = _mm_set1_epi32(max_H);
+				max_t_ = _mm_set1_epi32(max_t);
+				qe_    = _mm_set1_epi32(q + e);
+				for (t = st0; t < en1; t += 4) { // this implements: H[t]+=v8[t]-qe; if(H[t]>max_H) max_H=H[t],max_t=t;
+					__m128i H1, tmp, t_;
+					H1 = _mm_loadu_si128((__m128i*)&H[t]);
+					t_ = _mm_setr_epi32(v8[t], v8[t+1], v8[t+2], v8[t+3]);
+					H1 = _mm_add_epi32(H1, t_);
+					H1 = _mm_sub_epi32(H1, qe_);
+					_mm_storeu_si128((__m128i*)&H[t], H1);
+					t_ = _mm_set1_epi32(t);
+					tmp = _mm_cmpgt_epi32(H1, max_H_);
+#ifdef __SSE4_1__
+					max_H_ = _mm_blendv_epi8(max_H_, H1, tmp);
+					max_t_ = _mm_blendv_epi8(max_t_, t_, tmp);
+#else
+					max_H_ = _mm_or_si128(_mm_and_si128(tmp, H1), _mm_andnot_si128(tmp, max_H_));
+					max_t_ = _mm_or_si128(_mm_and_si128(tmp, t_), _mm_andnot_si128(tmp, max_t_));
+#endif
+				}
+				_mm_storeu_si128((__m128i*)HH, max_H_);
+				_mm_storeu_si128((__m128i*)tt, max_t_);
+				for (i = 0; i < 4; ++i)
+					if (max_H < HH[i]) max_H = HH[i], max_t = tt[i] + i;
+				for (; t < en0; ++t) { // for the rest of values that haven't been computed with SSE
+					H[t] += (int32_t)v8[t] - qe;
+					if (H[t] > max_H)
+						max_H = H[t], max_t = t;
+				}
+			} else H[0] = v8[0] - qe - qe, max_H = H[0], max_t = 0; // special casing r==0
+			// update ez
+			if (en0 == tlen - 1 && H[en0] > ez->mte)
+				ez->mte = H[en0], ez->mte_q = r - en;
+			if (r - st0 == qlen - 1 && H[st0] > ez->mqe)
+				ez->mqe = H[st0], ez->mqe_t = st0;
+			if (ksw_apply_zdrop(ez, 1, max_H, r, max_t, zdrop, e)) break;
+			if (r == qlen + tlen - 2 && en0 == tlen - 1)
+				ez->score = H[tlen - 1];
+		} else { // find approximate max; Z-drop might be inaccurate, too.
+			if (r > 0) {
+				if (last_H0_t >= st0 && last_H0_t <= en0 && last_H0_t + 1 >= st0 && last_H0_t + 1 <= en0) {
+					int32_t d0 = v8[last_H0_t] - qe;
+					int32_t d1 = u8[last_H0_t + 1] - qe;
+					if (d0 > d1) H0 += d0;
+					else H0 += d1, ++last_H0_t;
+				} else if (last_H0_t >= st0 && last_H0_t <= en0) {
+					H0 += v8[last_H0_t] - qe;
+				} else {
+					++last_H0_t, H0 += u8[last_H0_t] - qe;
+				}
+				if ((flag & KSW_EZ_APPROX_DROP) && ksw_apply_zdrop(ez, 1, H0, r, last_H0_t, zdrop, e)) break;
+			} else H0 = v8[0] - qe - qe, last_H0_t = 0;
+			if (r == qlen + tlen - 2 && en0 == tlen - 1)
+				ez->score = H0;
+		}
+		last_st = st, last_en = en;
+		//for (t = st0; t <= en0; ++t) printf("(%d,%d)\t(%d,%d,%d,%d)\t%d\n", r, t, ((int8_t*)u)[t], ((int8_t*)v)[t], ((int8_t*)x)[t], ((int8_t*)y)[t], H[t]); // for debugging
+	}
+	kfree(km, mem);
+	if (!approx_max) kfree(km, H);
+	if (with_cigar) { // backtrack
+		int rev_cigar = !!(flag & KSW_EZ_REV_CIGAR);
+		if (!ez->zdropped && !(flag&KSW_EZ_EXTZ_ONLY)) {
+			ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*16, tlen-1, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
+		} else if (!ez->zdropped && (flag&KSW_EZ_EXTZ_ONLY) && ez->mqe + end_bonus > (int)ez->max) {
+			ez->reach_end = 1;
+			ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*16, ez->mqe_t, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
+		} else if (ez->max_t >= 0 && ez->max_q >= 0) {
+			ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*16, ez->max_t, ez->max_q, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
+		}
+		kfree(km, mem2); kfree(km, off);
+	}
+}
+#endif // __SSE2__