fix bugs for overlaps of different length

Assembly.cpp

@@ -29,6 +29,7 @@ int roundID = 0;
 long long complete_threads = 0;
 
 
+
 void* Perform_Counting(void* arg)
 {
     int thr_ID = *((int*)arg);
@@ -1237,7 +1238,7 @@ long long xBeg, long long xEnd, long long yBeg, long long yEnd)
 }
 
 long long push_final_overlaps(ma_hit_t_alloc* paf, ma_hit_t_alloc* reverse_paf_list, 
-overlap_region_alloc* overlap_list, UC_Read* x_read, UC_Read* y_read)
+overlap_region_alloc* overlap_list, UC_Read* x_read, UC_Read* y_read, int flag, int test_exact)
 {
     long long i = 0;
     long long available_overlaps = 0;
@@ -1245,7 +1246,8 @@ overlap_region_alloc* overlap_list, UC_Read* x_read, UC_Read* y_read)
     clear_ma_hit_t_alloc(paf);
     for (i = 0; i < overlap_list->length; i++)
     {
-        if (overlap_list->list[i].is_match == 1)
+        ///if (overlap_list->list[i].is_match == 1)
+        if (overlap_list->list[i].is_match == flag)
         {
             available_overlaps++;
             /**********************query***************************/
@@ -1286,8 +1288,8 @@ overlap_region_alloc* overlap_list, UC_Read* x_read, UC_Read* y_read)
 
 
 
-
+            if(test_exact == 1)
-
+            {
                 if(overlap_list->list[i].y_pos_strand == 0)
                 {
                     recover_UC_Read(y_read, &R_INF, overlap_list->list[i].y_id);
@@ -1300,10 +1302,7 @@ overlap_region_alloc* overlap_list, UC_Read* x_read, UC_Read* y_read)
                 tmp.el = if_exact_match(x_read->seq, x_read->length, y_read->seq, y_read->length, 
                 overlap_list->list[i].x_pos_s, overlap_list->list[i].x_pos_e, 
                 overlap_list->list[i].y_pos_s, overlap_list->list[i].y_pos_e);
-
+            }
-
-
-
 
             add_ma_hit_t_alloc(paf, &tmp);
         }
@@ -1911,16 +1910,23 @@ void* Overlap_calculate_heap_merge(void* arg)
     init_small_hash_table(&reverse);
 
 
+
+    uint8_t c2n[256];
+    memset(c2n, 4, 256);
+    c2n['A'] = c2n['a'] = 0; c2n['C'] = c2n['c'] = 1;
+	c2n['G'] = c2n['g'] = 2; c2n['T'] = c2n['t'] = 3; // build the encoding table
+
+
     for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num)
     {
-        get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.02);
+        get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5);
         
         clear_Cigar_record(&current_cigar);
         clear_Round2_alignment(&second_round);
 
         correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, &POA_Graph, &DAGCon,
         &matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1,
-        &current_cigar, &hap, &second_round, 0, 1, &fully_cov, &abnormal);
+        &current_cigar, &hap, &second_round, 0, 1, &fully_cov, &abnormal, c2n);
 
         num_read_base = num_read_base + g_read.length;
         num_correct_base = num_correct_base + correct.corrected_base;
@@ -2113,7 +2119,7 @@ void* Output_related_reads(void* arg)
             &&
            memcmp(required_read_name, Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0)
         {
-            get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.02);
+            get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5);
 
 
             fprintf(stderr, ">%.*s\n", Get_NAME_LENGTH((R_INF), i),
@@ -2815,6 +2821,77 @@ HeapSq* heap, Candidates_list* l)
 }
 
 
+
+
+
+void debug_print_overlap(char* y_name, overlap_region_alloc* overlap_list, All_reads* R_INF, char* func)
+{
+    fprintf(stderr, "\nafter %s\n", func);
+    long long i, j;
+    for (i = 0; i < overlap_list->length; i++)
+    {
+        if(y_name == NULL || memcmp(y_name, Get_NAME((*R_INF), overlap_list->list[i].y_id), 
+        Get_NAME_LENGTH((*R_INF), overlap_list->list[i].y_id)) == 0)
+        {
+            fprintf(stderr, "****************x_name: %.*s, x_id: %d****************\n", 
+            Get_NAME_LENGTH((*R_INF), overlap_list->list[i].x_id), 
+            Get_NAME((*R_INF), overlap_list->list[i].x_id), overlap_list->list[i].x_id);
+
+            fprintf(stderr, "y_name: %.*s, y_id: %d, is_match: %d\n", 
+            Get_NAME_LENGTH((*R_INF),overlap_list->list[i].y_id),
+            Get_NAME((*R_INF),overlap_list->list[i].y_id), overlap_list->list[i].y_id,
+            overlap_list->list[i].is_match);
+
+            fprintf(stderr, "alignLen: %d, x_s: %d, x_e: %d, y_s: %d, y_e: %d, y_dir: %d, strong: %d\n", 
+                overlap_list->list[i].align_length,
+                overlap_list->list[i].x_pos_s,
+                overlap_list->list[i].x_pos_e,
+                overlap_list->list[i].y_pos_s,
+                overlap_list->list[i].y_pos_e,
+                overlap_list->list[i].y_pos_strand,
+                overlap_list->list[i].strong);
+            
+            fprintf(stderr, "i: %d, %.*s, x_s: %d, x_e: %d, y_s: %d, y_end: %d, w_list_length: %d, dir: %d, strong: %d, is_match: %d\n", 
+                i, Get_NAME_LENGTH((*R_INF),overlap_list->list[i].y_id),
+                Get_NAME((*R_INF),overlap_list->list[i].y_id),
+                overlap_list->list[i].x_pos_s,
+                overlap_list->list[i].x_pos_e,
+                overlap_list->list[i].y_pos_s,
+                overlap_list->list[i].y_pos_e,
+                overlap_list->list[i].w_list_length,
+                overlap_list->list[i].y_pos_strand,
+                overlap_list->list[i].strong,
+                overlap_list->list[i].is_match);
+
+            for (j = 0; j < overlap_list->list[i].w_list_length; j++)
+            {
+                fprintf(stderr, "************************\ncigar_j: %d, x_s: %d, x_e: %d, y_s: %d, y_end: %d\n", 
+                j, overlap_list->list[i].w_list[j].x_start,
+                overlap_list->list[i].w_list[j].x_end, 
+                overlap_list->list[i].w_list[j].y_start,
+                overlap_list->list[i].w_list[j].y_end);
+                if(overlap_list->list[i].w_list[j].y_end == -1)
+                {
+                    fprintf(stderr, "not match\n");
+                }
+                else
+                {
+                    int cigar_i, operation, operationLen;
+                    CIGAR* cigar = &(overlap_list->list[i].w_list[j].cigar);
+                    fprintf(stderr, "length: %d\n", cigar->length);
+                    for (cigar_i = 0; cigar_i < cigar->length; cigar_i++)
+                    {
+                        operation = cigar->C_C[cigar_i];
+                        operationLen = cigar->C_L[cigar_i];
+                        fprintf(stderr, "oper: %d, Len: %d\n", operation, operationLen);
+                    }
+                }
+
+            }
+        }
+    }
+}
+
 void* Final_overlap_calculate_heap_merge(void* arg)
 {
     long long matched_overlap_0 = 0;
@@ -2892,10 +2969,36 @@ void* Final_overlap_calculate_heap_merge(void* arg)
     init_small_hash_table(&reverse);
 
 
+    long long pre_r_overlaps = 0;
+    long long cur_r_overlaps = 0;
+    long long pre_overlaps = 0;
+    long long cur_overlaps = 0;
+    if(thr_ID == 0)
+    {
+        for (i = 0; i < R_INF.total_reads; i++)
+        {
+            pre_r_overlaps += R_INF.reverse_paf[i].length;
+            pre_overlaps += R_INF.paf[i].length;
+        }
+    }
+
+    // if(thr_ID == 0)
+    // {
+    //     debug_info_of_specfic_read("m64016_190918_162737/130811282/ccs", 
+    //     R_INF.paf, R_INF.reverse_paf, -1, "xxxx");
+
+    //     debug_info_of_specfic_read("m64016_190918_162737/179635219/ccs", 
+    //     R_INF.paf, R_INF.reverse_paf, -1, "xxxx");
+    // }
+    
+
+
     for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num)
     {
-        ///0.1%
+
-        get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.001);
+        
+
+        get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5);
         /**
         correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, &POA_Graph, &DAGCon,
         &matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1,
@@ -2906,6 +3009,22 @@ void* Final_overlap_calculate_heap_merge(void* arg)
 
         overlap_region_sort_y_id(overlap_list.list, overlap_list.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);
+
+
+        // if(memcmp("m64016_190918_162737/130811282/ccs", Get_NAME((R_INF), i), 
+        // Get_NAME_LENGTH((R_INF), i)) == 0)
+        // {
+        //     fprintf(stderr, "\n1\n");
+        //     debug_print_overlap("m64016_190918_162737/179635219/ccs", &overlap_list, &R_INF, "first");
+        // }
+
+        // if(memcmp("m64016_190918_162737/179635219/ccs", Get_NAME((R_INF), i), 
+        // Get_NAME_LENGTH((R_INF), i)) == 0)
+        // {
+        //     fprintf(stderr, "\n1\n");
+        //     debug_print_overlap("m64016_190918_162737/130811282/ccs", &overlap_list, &R_INF, "first");
+        // }
 
 
         
@@ -2961,6 +3080,67 @@ void* Final_overlap_calculate_heap_merge(void* arg)
         }
 
 
+
+        // if(memcmp("m64016_190918_162737/130811282/ccs", Get_NAME((R_INF), i), 
+        // Get_NAME_LENGTH((R_INF), i)) == 0)
+        // {
+        //     fprintf(stderr, "\n2\n");
+        //     debug_print_overlap("m64016_190918_162737/179635219/ccs", &overlap_list, &R_INF, "second");
+        // }
+
+        // if(memcmp("m64016_190918_162737/179635219/ccs", Get_NAME((R_INF), i), 
+        // Get_NAME_LENGTH((R_INF), i)) == 0)
+        // {
+        //     fprintf(stderr, "\n2\n");
+        //     debug_print_overlap("m64016_190918_162737/130811282/ccs", &overlap_list, &R_INF, "second");
+        // }
+
+
+
+
+        inner_j = 0;
+        j = 0;
+        while (j < overlap_list.length && inner_j < R_INF.reverse_paf[i].length)
+        {
+            if(overlap_list.list[j].y_id < R_INF.reverse_paf[i].buffer[inner_j].tn)
+            {
+                j++;
+            }
+            else if(overlap_list.list[j].y_id > R_INF.reverse_paf[i].buffer[inner_j].tn)
+            {
+                inner_j++;
+            }
+            else
+            {
+                if(overlap_list.list[j].y_pos_strand == R_INF.reverse_paf[i].buffer[inner_j].rev)
+                {
+
+                    overlap_list.list[j].is_match = 2;
+                    overlap_list.list[j].strong = 0;
+                    overlap_list.list[j].without_large_indel = 1;
+                }
+                j++;
+                inner_j++;
+            }
+        }
+
+
+
+        // if(memcmp("m64016_190918_162737/130811282/ccs", Get_NAME((R_INF), i), 
+        // Get_NAME_LENGTH((R_INF), i)) == 0)
+        // {
+        //     fprintf(stderr, "\n3\n");
+        //     debug_print_overlap("m64016_190918_162737/179635219/ccs", &overlap_list, &R_INF, "third");
+        // }
+
+        // if(memcmp("m64016_190918_162737/179635219/ccs", Get_NAME((R_INF), i), 
+        // Get_NAME_LENGTH((R_INF), i)) == 0)
+        // {
+        //     fprintf(stderr, "\n3\n");
+        //     debug_print_overlap("m64016_190918_162737/130811282/ccs", &overlap_list, &R_INF, "third");
+        // }
+
+
         ///recover missing exact overlaps 
         reverse_complement(g_read.seq, g_read.length);
         for (j = 0; j < overlap_list.length; j++)
@@ -2991,18 +3171,60 @@ void* Final_overlap_calculate_heap_merge(void* arg)
         }
 
 
+        // if(memcmp("m64016_190918_162737/130811282/ccs", Get_NAME((R_INF), i), 
+        // Get_NAME_LENGTH((R_INF), i)) == 0)
+        // {
+        //     fprintf(stderr, "\n4\n");
+        //     debug_print_overlap("m64016_190918_162737/179635219/ccs", &overlap_list, &R_INF, "fourth");
+        // }
+
+        // if(memcmp("m64016_190918_162737/179635219/ccs", Get_NAME((R_INF), i), 
+        // Get_NAME_LENGTH((R_INF), i)) == 0)
+        // {
+        //     fprintf(stderr, "\n4\n");
+        //     debug_print_overlap("m64016_190918_162737/130811282/ccs", &overlap_list, &R_INF, "fourth");
+        // }
+
+
+
+
         if(R_INF.paf[i].is_fully_corrected)
         {
             potiental_matched_overlap_1++;
         }
 
-
-
-
         num_correct_base += 
-        push_final_overlaps(&(R_INF.paf[i]), R_INF.reverse_paf, &overlap_list, &g_read, &overlap_read);
+        push_final_overlaps(&(R_INF.paf[i]), R_INF.reverse_paf, 
+        &overlap_list, &g_read, &overlap_read, 1, 1);
+
+
+        push_final_overlaps(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, 
+        &overlap_list, &g_read, &overlap_read, 2, 0);
+
+
+        // if(memcmp("m64016_190918_162737/130811282/ccs", Get_NAME((R_INF), i), 
+        // Get_NAME_LENGTH((R_INF), i)) == 0)
+        // {
+        //     fprintf(stderr, "\n5\n");
+        //     debug_print_overlap("m64016_190918_162737/179635219/ccs", &overlap_list, &R_INF, "last");
+
+        //     debug_print_overlap(NULL, &overlap_list, &R_INF, "SET");
+        // }
+
+        // if(memcmp("m64016_190918_162737/179635219/ccs", Get_NAME((R_INF), i), 
+        // Get_NAME_LENGTH((R_INF), i)) == 0)
+        // {
+        //     fprintf(stderr, "\n5\n");
+        //     debug_print_overlap("m64016_190918_162737/130811282/ccs", &overlap_list, &R_INF, "last");
+
+        //     debug_print_overlap(NULL, &overlap_list, &R_INF, "SET");
+        // }
+
     }
 
+
+    
+
     finish_output_buffer();
 
     destory_buffer_sub_block(&current_sub_buffer);
@@ -3047,9 +3269,30 @@ void* Final_overlap_calculate_heap_merge(void* arg)
 
 
 
+        for (i = 0; i < R_INF.total_reads; i++)
+        {
+            cur_r_overlaps += R_INF.reverse_paf[i].length;
+            cur_overlaps += R_INF.paf[i].length;
+        }
+        fprintf(stderr, "pre_r_overlaps: %d, cur_r_overlaps: %d\n", 
+        pre_r_overlaps, cur_r_overlaps);
+
+        fprintf(stderr, "pre_overlaps: %d, cur_overlaps: %d\n", 
+        pre_overlaps, cur_overlaps);
+
+
+
+
+ 
+        // debug_info_of_specfic_read("m64016_190918_162737/130811282/ccs", 
+        // R_INF.paf, R_INF.reverse_paf, -1, "yyyy");
+
+        // debug_info_of_specfic_read("m64016_190918_162737/179635219/ccs", 
+        // R_INF.paf, R_INF.reverse_paf, -1, "yyyy");
+        
+              
     }
     pthread_mutex_unlock(&statistics);
-
     free(arg);
 }
 
@@ -3134,7 +3377,7 @@ void generate_overlaps(int last_round)
     fprintf(stdout, "Final overlaps have been calculated.\n");
     fprintf(stdout, "%-30s%18.2f\n\n", "Final overlaps calculation time:", Get_T() - start_time); 
 
-    ///Output_PAF();
+    Output_PAF();
     /**
     build_string_graph(MIN_OVERLAP_COVERAGE, R_INF.paf, R_INF.reverse_paf, R_INF.total_reads, R_INF.read_length, 
     MIN_OVERLAP_LEN, MAX_HANG_LEN, 3, 0.5, 0.7, 0.8, output_file_name, MAX_BUBBLE_DIST);

Correct.cpp

@@ -7,6 +7,7 @@
 #include "edlib.h"
 #include "Assembly.h"
 #include "CommandLines.h"
+#include "ksw2.h"
 
 long long T_total_match=0;
 long long T_total_unmatch=0;
@@ -14,6 +15,36 @@ long long T_total_mis=0;
 pthread_mutex_t debug_statistics ;
 
 
+
+
+void align(const char *tseq, const char *qseq, const int tl, const int ql, 
+const uint8_t *c, int sc_mch, int sc_mis, int gapo, int gape, int bandLen, int zdrop,
+int* max_q_pos, int* max_t, int* score)
+{
+	int i, a = sc_mch, b = sc_mis < 0? sc_mis : -sc_mis; // a>0 and b<0
+	int8_t mat[25] = { a,b,b,b,0, b,a,b,b,0, b,b,a,b,0, b,b,b,a,0, 0,0,0,0,0 };
+	uint8_t *ts, *qs;
+	ksw_extz_t ez;
+
+	memset(&ez, 0, sizeof(ksw_extz_t));
+	ts = (uint8_t*)malloc(tl);
+	qs = (uint8_t*)malloc(ql);
+	for (i = 0; i < tl; ++i) ts[i] = c[(uint8_t)tseq[i]]; // encode to 0/1/2/3
+	for (i = 0; i < ql; ++i) qs[i] = c[(uint8_t)qseq[i]];
+	///ksw_extz(0, ql, qs, tl, ts, 5, mat, gapo, gape, -1, -1, 0, &ez);
+    ksw_extz2_sse(0, ql, qs, tl, ts, 5, mat, gapo, gape, bandLen, zdrop, sc_mch, 0, &ez);
+    /**
+	for (i = 0; i < ez.n_cigar; ++i) // print CIGAR
+		printf("%d%c", ez.cigar[i]>>4, "MID"[ez.cigar[i]&0xf]);
+	putchar('\n');
+    **/
+	free(ez.cigar); free(ts); free(qs);
+    (*score) = ez.max;
+}
+
+
+
+
 void clear_Round2_alignment(Round2_alignment* h)
 {
     clear_Correct_dumy_pure(&(h->dumy));
@@ -950,6 +981,116 @@ inline double trim_error_rate(overlap_region_alloc* overlap_list, long long ID)
     return error_rate;
 }
 
+///error_rate should be 30%
+long long get_high_error(long long x_start, long long x_end, 
+long long y_start, long long y_end, long long y_id, long long y_strand, long long pre_threshold,
+long long n_steps, float error_rate, All_reads* R_INF, Correct_dumy* dumy,
+UC_Read* g_read)
+{
+    long long stepLen = (x_end - x_start + 1) / n_steps;
+    if((x_end - x_start + 1) % n_steps != 0)
+    {
+        stepLen++;
+    }
+    long long SubLen, SubWindowLen;
+    long long SubThreshold = THRESHOLD_MAX_SIZE;
+    int extra_begin, extra_end;
+    long long o_len;
+    long long T_error = 0;
+    y_start = y_start + pre_threshold;
+
+    while (x_start <= x_end)
+    {
+        SubLen = x_end - x_start + 1;
+        if(SubLen > stepLen)
+        {
+            SubLen = stepLen;
+        }
+
+
+        SubThreshold = SubLen * error_rate;
+        if(SubThreshold > THRESHOLD_MAX_SIZE)
+        {
+            SubThreshold = THRESHOLD_MAX_SIZE;
+        }
+
+        SubThreshold = Adjust_Threshold(SubThreshold, SubLen);
+
+
+        SubWindowLen = SubLen + (SubThreshold << 1);
+        if(determine_overlap_region(SubThreshold, y_start, y_id, SubWindowLen, R_INF,
+        &extra_begin, &extra_end, &y_start, &o_len) == 0)
+        {
+            T_error = T_error + (x_end - x_start + 1) * error_rate * 1.5;
+            break;
+        }
+
+        fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, R_INF, y_id, 
+        extra_begin, extra_end);
+
+        char* x_string = g_read->seq + x_start;
+        char* y_string = dumy->overlap_region;
+        int end_site;
+        unsigned int error;
+
+        end_site = Reserve_Banded_BPM(y_string, SubWindowLen, x_string, SubLen, SubThreshold, &error);
+                    
+        ///error等于-1说明没匹配
+        if (error!=(unsigned int)-1)
+        {
+            T_error = T_error + error;
+            y_start = y_start + end_site - extra_begin + 1;
+        }
+        else
+        {
+            T_error = T_error + SubLen * error_rate * 1.5;
+            y_start = y_start + SubThreshold - extra_begin + SubLen;
+        }
+        
+        x_start = x_start + SubLen;
+    }
+
+
+    return T_error;
+}
+
+inline double non_trim_error_rate(overlap_region_alloc* overlap_list, long long ID,
+All_reads* R_INF, Correct_dumy* dumy, UC_Read* g_read)
+{
+    long long tLen, tError,i, subWinLen, subWinNum;
+    
+    tLen = 0;
+    tError = 0;
+
+    subWinNum = overlap_list->list[ID].w_list_length;
+
+    
+    for (i = 0; i < subWinNum; i++)
+    {
+        subWinLen = overlap_list->list[ID].w_list[i].x_end - overlap_list->list[ID].w_list[i].x_start + 1;
+        tLen += subWinLen;
+
+        if(overlap_list->list[ID].w_list[i].y_end != -1)
+        {
+            tError += overlap_list->list[ID].w_list[i].error;
+        }
+        else
+        {
+            tError +=
+            get_high_error(overlap_list->list[ID].w_list[i].x_start, 
+            overlap_list->list[ID].w_list[i].x_end, overlap_list->list[ID].w_list[i].y_start,
+            overlap_list->list[ID].w_list[i].y_end, overlap_list->list[ID].y_id, 
+            overlap_list->list[ID].y_pos_strand, 
+            overlap_list->list[ID].w_list[i].error_threshold - overlap_list->list[ID].w_list[i].extra_begin,
+            3, 0.3, R_INF, dumy, g_read);
+        }
+    }
+
+    double error_rate = (double)(tError)/(double)(tLen);
+
+    return error_rate;
+}
+
 
 
 void mark_duplicate(overlap_region_alloc* overlap_list, All_reads* R_INF, 
@@ -3044,7 +3185,9 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN
                 }
             }
 
-            error_rate = trim_error_rate(overlap_list, j);
+            ///error_rate = trim_error_rate(overlap_list, j);
+            error_rate = non_trim_error_rate(overlap_list, j, R_INF, dumy, g_read);
+            
 
             ///if(error_rate <= 0.015)
             if(error_rate <= 0.03)
@@ -8639,6 +8782,19 @@ uint32_t* SNPs, long long SNPLen, haplotype_evdience_alloc* hap)
     {
         snpID = SNPs[i];
         
+        
+        ///if(overlap_list->list[0].x_id == 83735)
+        ///if(overlap_list->list[0].x_id == 83739)
+        // if(overlap_list->list[0].x_id == 1185538)
+        // {
+        //     fprintf(stderr, "SNPLen: %d, x_id: %d, hap->snp_stat[snpID].site: %d, occ_0: %d, occ_1: %d, occ_2: %d, overlap_num: %d\n", 
+        //     SNPLen, overlap_list->list[0].x_id, hap->snp_stat[snpID].site,
+        //     hap->snp_stat[snpID].occ_0, hap->snp_stat[snpID].occ_1, hap->snp_stat[snpID].occ_2,
+        //     hap->snp_stat[snpID].overlap_num);
+        // }
+        
+        
+
         ///check all overlaps
         for (j = 0; j < Get_SNP_Vector_Length((*hap)); j++)
         {
@@ -9578,6 +9734,22 @@ Correct_dumy* dumy)
 
 }
 
+int debug_print_snp_stat(char* name, haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, All_reads* R_INF)
+{
+    if(overlap_list->length > 0 &&
+       memcmp(name, Get_NAME((*R_INF), overlap_list->list[0].x_id), 
+       Get_NAME_LENGTH((*R_INF), overlap_list->list[0].x_id)) == 0)
+    {
+        fprintf(stderr, "\n%s, available_snp: %d\n", name, hap->available_snp);
+        int i;
+        for (i = 0; i < hap->available_snp; i++)
+        {
+            fprintf(stderr, "site: %d, occ_0: %d, occ_1: %d, occ_2: %d\n", 
+            hap->snp_stat[i].site, hap->snp_stat[i].occ_0,
+            hap->snp_stat[i].occ_1, hap->snp_stat[i].occ_2);
+        }
+    }  
+}
 
 int generate_haplotypes_DP(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, All_reads* R_INF, long long rLen, 
 int force_repeat)
@@ -9595,6 +9767,7 @@ int force_repeat)
         return 0;
     }
 
+    // debug_print_snp_stat("m64016_190918_162737/174131552/ccs", hap, overlap_list, R_INF);
 
     
     ///if hap->available_snp == 1, the following codes would have bugs
@@ -9817,6 +9990,8 @@ int force_repeat)
         }
     }
 
+    // debug_print_snp_stat("m64016_190918_162737/174131552/ccs", hap, overlap_list, R_INF);
+
 
     //if(hap->dp.max_snp_num > 0)
     if(hap->available_snp > 0)
@@ -10537,7 +10712,6 @@ void correct_overlap_back(overlap_region_alloc* overlap_list, All_reads* R_INF,
 }
 
 
-
 void print_overlap(char* name, long long readID, 
 overlap_region_alloc* overlap_list, All_reads* R_INF, int output_reads)
 {
@@ -10545,8 +10719,8 @@ overlap_region_alloc* overlap_list, All_reads* R_INF, int output_reads)
     Get_NAME_LENGTH((*R_INF),readID)) == 0)
     {
         long long i, j;
-        fprintf(stderr, "\n\n****************ref_read: %.*s****************\n", 
+        fprintf(stderr, "\n\n****************ref_read: %.*s, id: %d****************\n", 
-            Get_NAME_LENGTH((*R_INF),readID), Get_NAME((*R_INF),readID));
+            Get_NAME_LENGTH((*R_INF),readID), Get_NAME((*R_INF),readID), readID);
 
         fprintf(stderr, "\n###flag: 1\n");
 
@@ -10618,10 +10792,10 @@ overlap_region_alloc* overlap_list, All_reads* R_INF, int output_reads)
             fprintf(stderr, ">%.*s\n", Get_NAME_LENGTH((*R_INF),readID), Get_NAME((*R_INF),readID));
             fprintf(stderr, "%.*s\n", g_read.length, g_read.seq);
 
-            fprintf(stderr, "query_read:\n");
+            ///fprintf(stderr, "query_read:\n");
             for (i = 0; i < overlap_list->length; i++)
             {
-                fprintf(stderr, "i: %d\n", i);
+                ///fprintf(stderr, "i: %d\n", i);
                 recover_UC_Read(&g_read, R_INF, overlap_list->list[i].y_id);
                 fprintf(stderr, ">%.*s\n", 
                 Get_NAME_LENGTH((*R_INF),overlap_list->list[i].y_id),
@@ -10629,6 +10803,8 @@ overlap_region_alloc* overlap_list, All_reads* R_INF, int output_reads)
                 fprintf(stderr, "%.*s\n", g_read.length, g_read.seq);
             }
 
+            fprintf(stderr, "Has already output all related reads\n\n");
+
             destory_UC_Read(&g_read);
 
 
@@ -10684,13 +10860,14 @@ overlap_region_alloc* overlap_list, All_reads* R_INF, int output_reads)
 
 }
 
+
 void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, 
                         UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, Graph* g, Graph* DAGCon,
                         long long* matched_overlap_0, long long* matched_overlap_1, 
                         long long* potiental_matched_overlap_0, long long* potiental_matched_overlap_1,
                         Cigar_record* current_cigar, haplotype_evdience_alloc* hap, 
                         Round2_alignment* second_round, int force_repeat, int is_consensus,
-                        int* fully_cov, int* abnormal)
+                        int* fully_cov, int* abnormal, uint8_t* c2n)
 {
     reverse_complement(g_read->seq, g_read->length);
 
@@ -10735,19 +10912,13 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF,
     recalcate_window(overlap_list, R_INF, g_read, dumy, overlap_read);
 
 
-    // print_overlap("m64011_190329_072846/59507330/ccs", 
-    // overlap_list->list[0].x_id, overlap_list, R_INF);
-
 
     partition_overlaps(overlap_list, R_INF, g_read, dumy, hap, force_repeat);
 
 
-    // print_overlap("m64016_190918_162737/53545052/ccs", 
+    // print_overlap("m64016_190918_162737/174131552/ccs", 
     // overlap_list->list[0].x_id, overlap_list, R_INF, 1);
     
-
-
-    
     if(is_consensus)
     {
         generate_consensus(overlap_list, R_INF, g_read, dumy, g, DAGCon, current_cigar, second_round);

Correct.h

@@ -1230,7 +1230,7 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF,
                         long long* potiental_matched_overlap_0, long long* potiental_matched_overlap_1,
                         Cigar_record* current_cigar, haplotype_evdience_alloc* hap,
                         Round2_alignment* second_round, int force_repeat, int is_consensus,
-                        int* fully_cov, int* abnormal);
+                        int* fully_cov, int* abnormal, uint8_t* c2n);
 void init_Correct_dumy(Correct_dumy* list);
 void destory_Correct_dumy(Correct_dumy* list);
 void clear_Correct_dumy(Correct_dumy* list, overlap_region_alloc* overlap_list);

Hash_Table.cpp

@@ -643,11 +643,13 @@ int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_regi
         memset(list->list + (list->size/2), 0, sizeof(overlap_region)*(list->size/2));
     }
 
-    if (list->length!=0 && 
+    if (list->length!=0 && list->list[list->length - 1].y_id==tmp->y_id)
-    list->list[list->length - 1].y_id==tmp->y_id
-    )
     {    
-        if(list->list[list->length - 1].shared_seed >= tmp->shared_seed)
+        ///if(list->list[list->length - 1].shared_seed >= tmp->shared_seed)
+        if((list->list[list->length - 1].shared_seed > tmp->shared_seed)
+           ||
+           ((list->list[list->length - 1].shared_seed == tmp->shared_seed) && 
+           (list->list[list->length - 1].overlapLen <= tmp->overlapLen)))
         {
             return 0;
         }
@@ -1244,10 +1246,41 @@ void debug_chain(k_mer_hit* a, long long a_n, Chain_Data* dp)
     }
 }
 
+long long get_chainLen(long long x_beg, long long x_end, long long xLen, 
+long long y_beg, long long y_end, long long yLen)
+{
+    if(x_beg <= y_beg)
+    {
+        y_beg = y_beg - x_beg;
+        x_beg = 0;
+    }
+    else
+    {
+        x_beg = x_beg - y_beg;
+        y_beg = 0;
+    }
+
+    long long x_right_length = xLen - x_end - 1;
+    long long y_right_length = yLen - y_end - 1;
+
+
+    if(x_right_length <= y_right_length)
+    {
+        x_end = xLen - 1;
+        y_end = y_end + x_right_length;        
+    }
+    else
+    {
+        x_end = x_end + y_right_length;
+        y_end = yLen - 1;
+    }
+
+    return x_end - x_beg + 1;
+}
 
 ///double band_width_threshold = 0.05;
 void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* result, 
-double band_width_threshold)
+double band_width_threshold, int max_skip, int x_readLen, int y_readLen)
 {
     long long i, j;
     long long self_pos, pos, max_j, max_i, max_score, score, n_skip;
@@ -1294,31 +1327,36 @@ double band_width_threshold)
                 continue;
             }
 
+            ///min distance
             distance_min = distance_pos < distance_self_pos? distance_pos:distance_self_pos;
             score = distance_min < min_score? distance_min : min_score;
 
-            /**
-            log_distance_gap = distance_gap? ilog2_32(distance_gap) : 0;
-
-            score -= (long long)(distance_gap * 0.01 * min_score) + (log_distance_gap/2);
-            **/
             gap_rate = (double)((double)(total_indels)/(double)(total_self_length));
-            ///if the gap rate > 0.05, score will be negative
+            ///if the gap rate > 0.06, score will be negative
             score -= (long long)(gap_rate * score * band_width_penalty);
 
             score += dp->score[j];
 
+            ///find a new max score
             if(score > max_score)
             {
                 max_score = score;
                 max_j = j;
                 max_indels = total_indels;
                 max_self_length = total_self_length;
-                if (n_skip > 0)
+                /****************************may have bugs********************************/
+                n_skip = 0;
+                /****************************may have bugs********************************/
+            }/****************************may have bugs********************************/
+            else
             {
-                    n_skip--;
+                n_skip++;
+                if(n_skip > max_skip)
+                {
+                    break;
                 }
             }
+            /****************************may have bugs********************************/
         }
 
         dp->score[i] = max_score;
@@ -1332,15 +1370,33 @@ double band_width_threshold)
 
 
     
+
     max_score = -1;
     max_i = -1;
+    long long mini_xLen = x_readLen * 2 + 2, tmp_xLen;
     for (i = 0; i < a_n; ++i) 
     {
         if(dp->score[i] > max_score)
         {
             max_score = dp->score[i];
             max_i = i;
+            mini_xLen = get_chainLen(a[i].self_offset, a[i].self_offset, x_readLen, 
+            a[i].offset, a[i].offset, y_readLen);
         }
+        else if(dp->score[i] == max_score)
+        {
+            tmp_xLen = get_chainLen(a[i].self_offset, a[i].self_offset, x_readLen, 
+            a[i].offset, a[i].offset, y_readLen);
+
+            if(tmp_xLen < mini_xLen)
+            {
+                max_score = dp->score[i];
+                max_i = i;
+
+                mini_xLen = tmp_xLen;
+            }
+        }
+        
     }
 
 
@@ -1350,6 +1406,7 @@ double band_width_threshold)
     result->x_pos_e = a[i].self_offset;
     result->y_pos_e = a[i].offset;
     result->shared_seed = max_score;
+    result->overlapLen = mini_xLen;
 
     distance_self_pos = result->x_pos_e - a[i].self_offset;
     distance_pos = result->y_pos_e - a[i].offset;
@@ -1515,7 +1572,6 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_thresh
     long long tmp_pos_distance;
     long long tmp_self_pos_distance;
     long long constant_distance = 5;
-    double error_rate = 0.05;
 
 
     if (candidates->length == 0)
@@ -1563,8 +1619,79 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_thresh
         }
 
        
+        // if(
+        // ((memcmp("m64016_190918_162737/130811282/ccs", 
+        // Get_NAME((*R_INF), tmp_region.x_id), Get_NAME_LENGTH((*R_INF), tmp_region.x_id)) == 0) 
+        // && 
+        // (memcmp("m64016_190918_162737/179635219/ccs", 
+        // Get_NAME((*R_INF), tmp_region.y_id), Get_NAME_LENGTH((*R_INF), tmp_region.y_id)) == 0))
+        // ||
+        // ((memcmp("m64016_190918_162737/179635219/ccs", 
+        // Get_NAME((*R_INF), tmp_region.x_id), Get_NAME_LENGTH((*R_INF), tmp_region.x_id)) == 0) 
+        // && 
+        // (memcmp("m64016_190918_162737/130811282/ccs", 
+        // Get_NAME((*R_INF), tmp_region.y_id), Get_NAME_LENGTH((*R_INF), tmp_region.y_id)) == 0)))
+        // {
+        //     fprintf(stderr, "****************x_name: %.*s****************\n", 
+        //     Get_NAME_LENGTH((*R_INF), tmp_region.x_id), Get_NAME((*R_INF), tmp_region.x_id));
+        //     fprintf(stderr, "****************y_name: %.*s****************\n", 
+        //     Get_NAME_LENGTH((*R_INF), tmp_region.y_id), Get_NAME((*R_INF), tmp_region.y_id));
+            
+
+        //     k_mer_hit* k_list = candidates->list + sub_region_beg;
+        //     long long k_listLen = sub_region_end - sub_region_beg + 1;
+        //     long long k = 0;
+
+        //     fprintf(stderr, "k_listLen: %d\n", k_listLen);
+
+        //     for (k = 0; k < k_listLen; k++)
+        //     {
+        //         fprintf(stderr, "k: %d, readID: %d, strand: %d, offset: %d, self_offset: %d\n", 
+        //         k, k_list[k].readID, k_list[k].strand, k_list[k].offset, k_list[k].self_offset);
+        //     }
+        // }
+
+        
         chain_DP(candidates->list + sub_region_beg, 
-        sub_region_end - sub_region_beg + 1, &(candidates->chainDP), &tmp_region, band_width_threshold);
+        sub_region_end - sub_region_beg + 1, &(candidates->chainDP), &tmp_region, band_width_threshold,
+        50, Get_READ_LENGTH((*R_INF), tmp_region.x_id), Get_READ_LENGTH((*R_INF), tmp_region.y_id));
+
+        
+        // if(
+        // ((memcmp("m64016_190918_162737/130811282/ccs", 
+        // Get_NAME((*R_INF), tmp_region.x_id), Get_NAME_LENGTH((*R_INF), tmp_region.x_id)) == 0) 
+        // && 
+        // (memcmp("m64016_190918_162737/179635219/ccs", 
+        // Get_NAME((*R_INF), tmp_region.y_id), Get_NAME_LENGTH((*R_INF), tmp_region.y_id)) == 0))
+        // ||
+        // ((memcmp("m64016_190918_162737/179635219/ccs", 
+        // Get_NAME((*R_INF), tmp_region.x_id), Get_NAME_LENGTH((*R_INF), tmp_region.x_id)) == 0) 
+        // && 
+        // (memcmp("m64016_190918_162737/130811282/ccs", 
+        // Get_NAME((*R_INF), tmp_region.y_id), Get_NAME_LENGTH((*R_INF), tmp_region.y_id)) == 0)))
+        // {
+        //     fprintf(stderr, "****************x_name: %.*s****************\n", 
+        //     Get_NAME_LENGTH((*R_INF), tmp_region.x_id), Get_NAME((*R_INF), tmp_region.x_id));
+        //     fprintf(stderr, "****************y_name: %.*s****************\n", 
+        //     Get_NAME_LENGTH((*R_INF), tmp_region.y_id), Get_NAME((*R_INF), tmp_region.y_id));
+
+
+        //     k_mer_hit* k_list = candidates->list + sub_region_beg;
+        //     long long k_listLen = sub_region_end - sub_region_beg + 1;
+        //     long long k = 0;
+
+        //     fprintf(stderr, "k_listLen: %d\n", k_listLen);
+
+        //     for (k = 0; k < k_listLen; k++)
+        //     {
+        //         fprintf(stderr, "k: %d, score: %d, pre: %d, indels: %d, self_length: %d\n", 
+        //         k, candidates->chainDP.score[k], candidates->chainDP.pre[k], 
+        //         candidates->chainDP.indels[k], candidates->chainDP.self_length[k]);                
+        //     }
+        // }
+    
+
+
 
         ///自己和自己重叠的要排除
         ///if (tmp_region.x_id != tmp_region.y_id && tmp_region.shared_seed > 1)
@@ -2467,7 +2594,6 @@ int load_Total_Pos_Table(Total_Pos_Table* TCB, char* read_file_name)
     fread(&TCB->useful_k_mer, sizeof(TCB->useful_k_mer), 1, fp);
     fread(&TCB->total_occ, sizeof(TCB->total_occ), 1, fp);
 
-
     if (TCB->useful_k_mer+1)
     {
         TCB->k_mer_index = (uint64_t*)malloc(sizeof(uint64_t)*(TCB->useful_k_mer+1));

Hash_Table.h

@@ -120,6 +120,7 @@ typedef struct
     uint64_t y_pos_e;
     uint64_t y_pos_strand;
 
+    uint64_t overlapLen;
     uint64_t shared_seed;
     uint64_t align_length;
     ///uint64_t total_errors;

Levenshtein_distance.h

@@ -184,6 +184,9 @@ inline int Reserve_Banded_BPM
 	int i_last = i;
 	i = 0;
 
+	/****************************may have bugs********************************/
+	unsigned int ungap_error = (unsigned int)-1;
+	/****************************may have bugs********************************/
 
 	while (i < available_i)
 	{
@@ -196,7 +199,21 @@ inline int Reserve_Banded_BPM
 			*return_err = err;
 			return_site = site + i;
 		}
+
+		/****************************may have bugs********************************/
+		if(i == errthold)
+		{
+			ungap_error = err;
 		}
+		/****************************may have bugs********************************/
+	}
+
+	/****************************may have bugs********************************/
+	if((ungap_error<=errthold) && (ungap_error == (*return_err)))
+	{
+		return_site = site + errthold;
+	}
+	/****************************may have bugs********************************/
 
 	return return_site;
 
@@ -436,6 +453,11 @@ inline int Reserve_Banded_BPM_PATH
 	///int site = p_length - last_high - 1;
 	int site = t_length - 1;
 	int return_site = -1;
+
+	/****************************may have bugs********************************/
+	unsigned int ungap_error = (unsigned int)-1;
+	/****************************may have bugs********************************/
+
 	///p_length大部分情况下应该是t_length + 2 * errthold，这是i要小于last_high = 2 * errthold
 	///也就是p_length - t_length
 	///那么当p_length < t_length + 2 * errthold, available_i也应该是这个值
@@ -459,7 +481,16 @@ inline int Reserve_Banded_BPM_PATH
 			*return_err = err;
 			return_site = site + i;
 		}
+
+		/****************************may have bugs********************************/
+		if(i == errthold)
+		{
+			ungap_error = err;
 		}
+		/****************************may have bugs********************************/
+	}
+
+	
 
 
 	if ((*return_err) == (unsigned int)-1)
@@ -467,6 +498,14 @@ inline int Reserve_Banded_BPM_PATH
 		return return_site;
 	}
 
+
+	/****************************may have bugs********************************/
+	if((ungap_error<=errthold) && (ungap_error == (*return_err)))
+	{
+		return_site = site + errthold;
+	}
+	/****************************may have bugs********************************/
+	
     ////注意，这里p_length要矫正啊啊
 	///不矫正会出错
 	///因为p_length有可能不够
@@ -846,6 +885,13 @@ inline int Reserve_Banded_BPM_4_SSE_only(char *pattern1, char *pattern2, char *p
 
 	i = 0;
 
+	/****************************may have bugs********************************/
+	unsigned int ungap_error1 = (unsigned int)-1;
+	unsigned int ungap_error2 = (unsigned int)-1;
+	unsigned int ungap_error3 = (unsigned int)-1;
+	unsigned int ungap_error4 = (unsigned int)-1;
+	/****************************may have bugs********************************/
+
 
 	///p_length大部分情况下应该是t_length + 2 * errthold，这是i要小于last_high = 2 * errthold
 	///也就是p_length - t_length
@@ -891,7 +937,39 @@ inline int Reserve_Banded_BPM_4_SSE_only(char *pattern1, char *pattern2, char *p
 			return_sites[3] = site + i;
 			return_sites_error[3] = err4;
 		}
+
+		/****************************may have bugs********************************/
+		if(i == errthold)
+		{
+			ungap_error1 = err1;
+			ungap_error2 = err2;
+			ungap_error3 = err3;
+			ungap_error4 = err4;
 		}
+		/****************************may have bugs********************************/
+	}
+
+	/****************************may have bugs********************************/
+	if((ungap_error1<=errthold) && (ungap_error1 == return_sites_error[0]))
+	{
+		return_sites[0] = site + errthold;
+	}
+
+	if((ungap_error2<=errthold) && (ungap_error2 == return_sites_error[1]))
+	{
+		return_sites[1] = site + errthold;
+	}
+
+	if((ungap_error3<=errthold) && (ungap_error3 == return_sites_error[2]))
+	{
+		return_sites[2] = site + errthold;
+	}
+
+	if((ungap_error4<=errthold) && (ungap_error4 == return_sites_error[3]))
+	{
+		return_sites[3] = site + errthold;
+	}
+	/****************************may have bugs********************************/
 
 	return 1;
 }

Makefile

@@ -3,7 +3,7 @@ CXXFLAGS=	-g -O3 -msse4.2 -mpopcnt -fomit-frame-pointer -Winline -Wall
 CPPFLAGS=
 INCLUDES=
 OBJS=		Output.o CommandLines.o Process_Read.o Assembly.o kmer.o Hash_Table.o \
-			POA.o Correct.o Levenshtein_distance.o edlib.o Overlaps.o
+			POA.o Correct.o Levenshtein_distance.o edlib.o Overlaps.o ksw2_extz2_sse.o
 EXE=		ccs_assembly
 LIBS=		-lz -lpthread -lm
 
@@ -37,7 +37,7 @@ Assembly.o: Levenshtein_distance.h Output.h
 CommandLines.o: CommandLines.h ketopt.h
 Correct.o: Correct.h Hash_Table.h khash.h kmer.h Process_Read.h kseq.h
 Correct.o: Overlaps.h kvec.h kdq.h CommandLines.h Levenshtein_distance.h
-Correct.o: POA.h edlib.h Assembly.h
+Correct.o: POA.h edlib.h Assembly.h ksw2.h
 Hash_Table.o: Hash_Table.h khash.h kmer.h Process_Read.h kseq.h Overlaps.h
 Hash_Table.o: kvec.h kdq.h CommandLines.h Correct.h Levenshtein_distance.h
 Hash_Table.o: POA.h ksort.h
@@ -52,3 +52,4 @@ edlib.o: edlib.h
 kmer.o: kmer.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h CommandLines.h
 main.o: CommandLines.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h
 main.o: Assembly.h Levenshtein_distance.h edlib.h
+ksw2_extz2_sse.o: ksw2.h

Overlaps.cpp

@@ -884,6 +884,106 @@ long long mini_overlap_length, ma_sub_t** coverage_cut)
 
 
 
+void ma_hit_chimeric(int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
+long long n_read, uint64_t* readLen, ma_sub_t* coverage_cut)
+{
+    double startTime = Get_T();
+
+	int i, j, k, n_remove = 0;
+	kvec_t(uint32_t) b = {0,0,0};
+
+	for (i = 0; i < n_read; ++i) 
+    {
+
+        kv_resize(uint32_t, b, readLen[i]);
+        memset(b.a, 0, sizeof(uint32_t)*readLen[i]);
+
+        
+        
+        for (j = 0; j < sources[i].length; j++)
+        {
+            uint32_t qs, qe;
+            qs = Get_qs(sources[i].buffer[j]);
+            qe = Get_qe(sources[i].buffer[j]);
+
+
+            ///if(qe - qs < 1000) continue;
+            for (k = qs; k < qe; k++)
+            {
+                b.a[k]++;
+            }
+        }
+
+    
+        // for (j = 0; j < reverse_sources[i].length; j++)
+        // {
+        //     uint32_t qs, qe;
+        //     qs = Get_qs(reverse_sources[i].buffer[j]);
+        //     qe = Get_qe(reverse_sources[i].buffer[j]);
+
+        //     for (k = qs; k < qe; k++)
+        //     {
+        //         b.a[k]++;
+        //     }
+        // }
+    
+        
+        int left, right;
+        left = -1;
+        right = readLen[i];
+        for (k = 0; k < readLen[i]; k++)
+        {
+            if(b.a[k] < min_dp)
+            {
+                left = k;
+                break;
+            }
+        }
+        
+        for (k = readLen[i] - 1; k >= 0; k--)
+        {
+            if(b.a[k] < min_dp)
+            {
+                right = k;
+                break;
+            }
+        }
+
+        // if(i == 4616942 || i == 4024299 || i == 6135193)
+        // {
+        //     fprintf(stderr, "i: %d, left: %d, right: %d, rLen: %d, min_dp: %d\n", 
+        //     i, left, right, readLen[i], min_dp);
+        //     for (k = 0; k < readLen[i]; k++)
+        //     {
+        //         fprintf(stderr, "a[%d]: %d\n", k, b.a[k]);
+        //     }
+        // }
+    
+        if( (left <= right) && (left > 0) && (right < readLen[i] - 1))
+        {
+            coverage_cut[i].c = 1;
+            n_remove++;
+            /****************************may have bugs********************************/
+            coverage_cut[i].del = 1;
+            sources[i].length = 0;
+            /****************************may have bugs********************************/
+        }
+        else
+        {
+            coverage_cut[i].c = 0;
+        }
+	}
+
+	free(b.a);
+
+    fprintf(stderr, "[M::%s] takes %0.2f s, n_remove: %d\n\n", __func__, Get_T()-startTime, n_remove);
+}
+
+
+
+
+
+
 void ma_hit_cut(int min_dp, ma_hit_t_alloc* sources, long long n_read, uint64_t* readLen, 
 long long mini_overlap_length, ma_sub_t** coverage_cut)
 {
@@ -1488,6 +1588,7 @@ int max_hang, int min_ovlp)
     {
         ///if a read has been deleted, should we still add them?
 		asg_seq_set(g, i, coverage_cut[i].e - coverage_cut[i].s, coverage_cut[i].del);
+        g->seq[i].c = coverage_cut[i].c;
 	}
 
     g->seq_vis = (uint8_t*)calloc(g->n_seq*2, sizeof(uint8_t));
@@ -3889,15 +3990,6 @@ void* asg_arc_identify_simple_bubbles_pthread(void* arg)
 
         b.b.n = 0;
 
-        ///if(v%10000 == 0)
-        //if(v >= 18020000)
-        // if(v == 18021291)
-        // {
-        //     fprintf(stderr, "0 v: %d, thr_ID: %d\n", v, thr_ID);
-        //     fflush(stderr);
-        // }
-        
-
         if(g->seq_vis[v] != 1)
         {
             ///if(detect_bubble_with_bubbles(g, v, &w, &l, &b, (uint32_t)-1))
@@ -6294,7 +6386,6 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen)
 
 
         /****************************may have bugs********************************/
-        /**
         if(to_del_l != kv)
         {
             b_f.n = 0;
@@ -6384,7 +6475,317 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen)
         }
 
         terminal:
-        **/
+        /****************************may have bugs********************************/
+
+
+        if(to_del_l != kv) continue;
+
+
+
+
+
+        uint32_t convex1;
+        long long l1;
+
+
+
+
+        ////forward bubble
+        to_del_l = 0;
+        for (i = 0; i < b_f.n; i++)
+        {
+            if(b_f.a[i] == b_f.a[0])
+            {
+                to_del_l = 1;
+            } 
+            else
+            {
+                to_del_l = 0;
+                break;
+            }
+        }
+        //check the length
+        if(to_del_l == 0 && asg_bub_end_finder_with_del_advance(g, 
+        b_f.a, b_f.n, max_dist, &bub, 0, (u_int32_t)-1, &sink)==1)
+        {
+            to_del_l = 1;
+        }
+        if(to_del_l == 0 && detect_mul_bubble_end_with_bubbles(g, b_f.a, b_f.n, &convex1, &l1, NULL))
+        {
+            to_del_l = 1;
+        }
+
+        ///if(v>>1 == 4581428) fprintf(stderr, "to_del_l: %d, b_f.n: %d\n", to_del_l, b_f.n);
+
+        ///if(to_del_l == 0) continue;
+
+
+
+
+
+
+
+        ////backward bubble
+        to_del_r = 0;
+        for (i = 0; i < b_r.n; i++)
+        {
+            if(b_r.a[i] == b_r.a[0])
+            {
+                to_del_r = 1;
+            } 
+            else
+            {
+                to_del_r = 0;
+                break;
+            }
+        }
+        if(to_del_r == 0 && asg_bub_end_finder_with_del_advance
+        (g, b_r.a, b_r.n, max_dist, &bub, 1, v^1, &sink)==1)
+        {
+            to_del_r = 1;
+        }
+        if(to_del_r == 0 && detect_mul_bubble_end_with_bubbles(g, b_r.a, b_r.n, &convex1, &l1, NULL))
+        {
+            to_del_r = 1;
+        }
+
+        
+        // if(v>>1 == 4581428)
+        // {
+        //     fprintf(stderr, "to_del_l: %d, b_f.n: %d\n", to_del_l, b_f.n);
+        //     asg_bub_end_finder_with_del_advance_debug(g, b_r.a, b_r.n, max_dist, &bub, 1, v^1);
+        // }
+        
+        // if(v>>1 == 7318036)
+        // {
+        //     fprintf(stderr, "to_del_l: %d, to_del_r: %d, b_f.n: %d\n", to_del_l, to_del_r, b_f.n);
+        // } 
+
+
+
+		if (to_del_l && to_del_r)
+        {
+            for (i = 0; i < nv; ++i) 
+            {
+			    if (av[i].del) continue;
+                ///fprintf(stderr, "%.*s\n", Get_NAME_LENGTH((R_INF), v>>1), Get_NAME((R_INF), v>>1));
+                ++n_cut;
+                av[i].del = 1;
+                asg_arc_del(g, av[i].v^1, av[i].ul>>32^1, 1);
+		    }
+
+        }
+    }
+    
+    free(b.a); free(b_f.a); free(b_r.a);
+    free(bub.a); free(bub.S.a); free(bub.T.a); free(bub.b.a); free(bub.e.a);
+	if (n_cut) 
+    {
+		asg_cleanup(g);
+		asg_symm(g);
+	}
+	fprintf(stderr, "[M::%s] removed %d false overlaps\n", __func__, n_cut);
+    fprintf(stderr, "[M::%s] takes %0.2f s\n\n", __func__, Get_T()-startTime);
+	return n_cut;
+}
+
+
+int asg_arc_del_short_false_link_advance(asg_t *g, float drop_ratio, float o_drop_ratio, int max_dist, 
+ma_hit_t_alloc* reverse_sources, long long miniedgeLen)
+{
+    double startTime = Get_T();
+    
+    kvec_t(uint64_t) b;
+    memset(&b, 0, sizeof(b));
+
+
+    kvec_t(uint32_t) b_f;
+    memset(&b_f, 0, sizeof(b_f));
+
+    kvec_t(uint32_t) b_r;
+    memset(&b_r, 0, sizeof(b_r));
+
+
+	uint32_t v, w, n_vtx = g->n_seq * 2, n_cut = 0;
+    uint32_t sink;
+
+    buf_t bub;
+    if (!g->is_symm) asg_symm(g);
+    memset(&bub, 0, sizeof(buf_t));
+    bub.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t));
+    
+    for (v = 0; v < n_vtx; ++v) 
+    {
+        if(g->seq_vis[v] == 0)
+        {
+            asg_arc_t *av = asg_arc_a(g, v);
+            uint32_t nv = asg_arc_n(g, v);
+            if(nv == 1 && asg_arc_n(g, v^1) == 1) continue;
+
+            uint64_t t_ol = 0;
+            long long i;
+            for (i = 0; i < nv; ++i)
+            {
+                t_ol += av[i].ol;
+            }
+            kv_push(uint64_t, b, (uint64_t)(t_ol << 32 | v));  
+        }
+	}
+
+    fprintf(stderr, "[M::%s] %lld unsorted pending overlaps\n", __func__, b.n);
+
+    radix_sort_arch64(b.a, b.a + b.n);
+
+    fprintf(stderr, "[M::%s] %lld sorted pending overlaps\n", __func__, b.n);
+
+    
+
+    uint32_t min_edge;
+    
+    
+    long long k, t;
+    for (k = 0; k < b.n; k++)
+    {
+        ///v is the node
+        v = (uint32_t)b.a[k];
+        if (g->seq[v>>1].del) continue;
+        uint32_t nv = asg_arc_n(g, v), nw, to_del_l, to_del_r;
+        if (nv < 2) continue;
+        uint32_t kv = get_real_length(g, v, NULL), kw;
+        if (kv < 2) continue;
+        uint32_t i;
+        asg_arc_t *av = asg_arc_a(g, v), *aw;
+
+        b_f.n = 0;
+        b_r.n = 0;
+        to_del_l = 0;
+        for (i = 0; i < nv; i++)
+        {
+            if (av[i].del) continue;
+            
+            w = av[i].v^1;
+            nw = asg_arc_n(g, w);
+            if(nw < 2) break;
+            kw = get_real_length(g, w, NULL);
+            if(kw < 2) break;
+
+            kv_push(uint32_t, b_f, av[i].v);
+            kv_push(uint32_t, b_r, w);
+
+            aw = asg_arc_a(g, w);
+            min_edge = (u_int32_t)-1;
+            for (t = 0; t < nw; t++)
+            {
+                if(aw[t].del) continue;
+                if((aw[t].v>>1) == (v>>1)) continue;
+                if(aw[t].ol < min_edge) min_edge = aw[t].ol;
+                ///kv_push(uint32_t, b_r, aw[t].v);
+            }
+
+            if(av[i].ol < min_edge * drop_ratio) to_del_l++;
+        }
+
+
+
+
+
+
+
+
+
+
+
+
+        /****************************may have bugs********************************/
+        if(to_del_l != kv)
+        {
+            b_f.n = 0;
+            b_r.n = 0;
+            to_del_l = 0;
+
+            for (i = 0; i < nv; i++)
+            {
+                if (av[i].del) continue;
+
+                w = av[i].v^1;
+                nw = asg_arc_n(g, w);
+                if(nw < 2) break;
+                kw = get_real_length(g, w, NULL);
+                if(kw < 2) break;
+
+                kv_push(uint32_t, b_f, av[i].v);
+                kv_push(uint32_t, b_r, w);
+
+                aw = asg_arc_a(g, w);
+                min_edge = (u_int32_t)-1;
+                for (t = 0; t < nw; t++)
+                {
+                    if(aw[t].del) continue;
+                    if((aw[t].v>>1) == (v>>1)) continue;
+                    if(aw[t].ol < min_edge) min_edge = aw[t].ol;
+                }
+
+                if(av[i].ol < min_edge * o_drop_ratio) to_del_l++;
+            }
+
+            if(to_del_l == kv)
+            {
+                ///forward
+                to_del_l = 1;
+                for (i = 1; i < b_f.n; i++)
+                {
+                    if(check_if_diploid(b_f.a[0], b_f.a[i], g, reverse_sources, miniedgeLen) == 1)
+                    {
+                        to_del_l++;
+                    }
+                }
+
+                ///backward
+                if(to_del_l != kv && b_r.n >= 2)
+                {
+                    to_del_l = 0;
+                    uint32_t w0, w1;
+
+
+                    w = b_r.a[0];
+                    kw = get_real_length(g, w, NULL);
+                    if(kw != 2) goto terminal;
+                    aw = asg_arc_a(g, w);
+                    nw = asg_arc_n(g, w);
+                    for (t = 0; t < nw; t++)
+                    {
+                        if(aw[t].del) continue;
+                        if((aw[t].v>>1) == (v>>1)) continue;
+                        w0 = aw[t].v;
+                    }
+                    to_del_l = 1;
+
+
+                    for (i = 1; i < b_r.n; i++)
+                    {
+                        w = b_r.a[i];
+                        kw = get_real_length(g, w, NULL);
+                        if(kw != 2) goto terminal;
+                        aw = asg_arc_a(g, w);
+                        nw = asg_arc_n(g, w);
+                        for (t = 0; t < nw; t++)
+                        {
+                            if(aw[t].del) continue;
+                            if((aw[t].v>>1) == (v>>1)) continue;
+                            w1 = aw[t].v;
+                        }
+
+                        if(check_if_diploid(w0, w1, g, reverse_sources, miniedgeLen) == 1)
+                        {
+                            to_del_l++;
+                        }
+                    }
+
+                }
+            }
+        }
+
+        terminal:
         /****************************may have bugs********************************/
 
 
@@ -6472,10 +6873,11 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen)
 
 		if (to_del_l && to_del_r)
         {
+            fprintf(stderr, "%.*s\n", Get_NAME_LENGTH((R_INF), v>>1), Get_NAME((R_INF), v>>1));
             for (i = 0; i < nv; ++i) 
             {
 			    if (av[i].del) continue;
-                ///fprintf(stderr, "%.*s\n", Get_NAME_LENGTH((R_INF), v>>1), Get_NAME((R_INF), v>>1));
+                
                 ++n_cut;
                 av[i].del = 1;
                 asg_arc_del(g, av[i].v^1, av[i].ul>>32^1, 1);
@@ -7856,7 +8258,7 @@ ma_hit_t_alloc* reverse_sources, int id, char* command)
         }
     }
 
-
+    fflush(stderr);
 
     
 }
@@ -8871,7 +9273,50 @@ int asg_arc_del_orthology_multiple_way(asg_t *g, ma_hit_t_alloc* reverse_sources
 
 
 
+int asg_arc_del_chimeric_read(asg_t *g, long long miniedgeLen)
+{
+    double startTime = Get_T();
+    ///the reason is that each read has two direction (query->target, target->query)
+    uint32_t v, v_max, v_maxLen, w, n_vtx = g->n_seq * 2, n_reduced = 0, n_arc;
 
+    for (v = 0; v < n_vtx; ++v) 
+    {
+        ///if (g->seq_vis[v] != 0) continue;
+        if (g->seq[v>>1].del) continue;
+        if (g->seq[v>>1].c == 0) continue;
+        ///fprintf(stderr, "v>>1: %d\n", v>>1);
+        if((get_real_length(g, v, NULL) == 0) || (get_real_length(g, v^1, NULL) == 0))
+        {
+            continue;
+        }
+
+
+        uint32_t convex1, convex2, flag1, flag2;
+        long long l1, l2, ll;
+        flag1 = detect_single_path_with_dels(g, v, &convex1, &l1, NULL);
+        if(flag1 == END_TIPS || flag1 == LONG_TIPS) continue;
+        flag2 = detect_single_path_with_dels(g, v^1, &convex2, &l2, NULL);
+        if(flag2 == END_TIPS || flag2 == LONG_TIPS) continue;
+        ll = l1 + l2 - 1;
+        if(ll <= miniedgeLen)
+        {
+            ///fprintf(stderr, "***v>>1: %d\n", v>>1);
+            asg_seq_del(g, v>>1);
+            n_reduced++;
+        }
+    }
+
+    if (n_reduced) {
+        asg_cleanup(g);
+        asg_symm(g);
+    }
+
+    fprintf(stderr, "[M::%s] removed %d chimeric reads\n", 
+    __func__, n_reduced);
+    fprintf(stderr, "[M::%s] takes %0.2f s\n\n", __func__, Get_T()-startTime);
+
+    return n_reduced;
+}
 
 
 uint32_t detect_single_path_with_dels_by_length
@@ -9016,10 +9461,38 @@ long long asg_arc_del_self_circle_untig(asg_t *g, long long circleLen)
 }
 
 
+void output_unitig_graph_without_small_bubbles(asg_t *sg, ma_sub_t* coverage_cut, 
+char* output_file_name, long long n_read, long long bubble_dist, long long tipsLen)
+{
+    asg_cut_tip(sg, tipsLen);
+    asg_pop_bubble(sg, bubble_dist);
+    asg_cut_tip(sg, tipsLen);
+    
+    ma_ug_t *ug = NULL;
+    ug = ma_ug_gen(sg);
+    ma_ug_seq(ug, &R_INF, coverage_cut, n_read);
+
+    fprintf(stdout, "Writing unitig GFA to disk ...... \n");
+    char* gfa_name = (char*)malloc(strlen(output_file_name)+35);
+    sprintf(gfa_name, "%s.no_s_bub.gfa", output_file_name);
+    FILE* output_file = fopen(gfa_name, "w");
+    ma_ug_print(ug, &R_INF, coverage_cut, output_file);
+    fclose(output_file);
+    
+    
+    sprintf(gfa_name, "%s.simple.no_s_bub.gfa", output_file_name);
+    output_file = fopen(gfa_name, "w");
+    ma_ug_print_simple(ug, &R_INF, coverage_cut, output_file);
+    fclose(output_file);
+
+    free(gfa_name);
+    ma_ug_destroy(ug);
+}
+
+
 void output_contig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, long long n_read, long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long circleLen,
 ma_hit_t_alloc* reverse_sources, long long miniedgeLen)
 {
-
     asg_cut_tip(sg, tipsLen);
     // asg_pop_bubble(sg, bubble_dist);
     // asg_arc_del_self_circle_untig(sg, circleLen);
@@ -9038,12 +9511,9 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen)
 
         cur_cons = sg->n_seq + sg->n_arc;
     }
-    /**
+
-    fprintf(stderr, "pop bubbles: %d, sg->n_seq + sg->n_arc: %d\n", 
+    asg_arc_identify_simple_bubbles_multi(sg, 1);
-    asg_pop_bubble(sg, bubble_dist), sg->n_seq + sg->n_arc);
+    asg_arc_del_short_false_link(sg, 0.6, 0.85, bubble_dist, reverse_sources, MAX_SHORT_TIPS);
-    fprintf(stderr, "pop bubbles: %d, sg->n_seq + sg->n_arc: %d\n", 
-    asg_pop_bubble(sg, bubble_dist), sg->n_seq + sg->n_arc);
-    **/
 
     ///asg_arc_del_self_circle_untig(sg, circleLen);
     
@@ -9071,6 +9541,23 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen)
     ma_ug_destroy(ug);
 }
 
+
+int output_tips(asg_t *g, const All_reads *RNF)
+{
+	uint32_t v, w, n_vtx = g->n_seq * 2, n_reduced = 0;
+    for (v = 0; v < n_vtx; ++v)
+    {
+        if (g->seq[v>>1].del) continue;
+
+        if(asg_arc_n(g, v) == 0)
+        {
+            fprintf(stderr, "%.*s\n",
+             Get_NAME_LENGTH((*RNF), v>>1), 
+             Get_NAME((*RNF), v>>1));
+        }
+    }
+}
+
 void build_string_graph_without_clean(
 int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
 long long n_read, uint64_t* readLen, 
@@ -9086,24 +9573,21 @@ char* output_file_name, long long bubble_dist, int read_graph, int write)
     }
 
 
-    // debug_info_of_specfic_read("m64016_190918_162737/92668450/ccs", 
+    // debug_info_of_specfic_read("m64016_190918_162737/179635219/ccs", 
     // sources, reverse_sources, -1, "init");
 
-    // debug_info_of_specfic_read("m64016_190918_162737/53545052/ccs", 
+    // debug_info_of_specfic_read("m64016_190918_162737/130811282/ccs", 
     // sources, reverse_sources, -1, "init");
     
 
     
+    
     ma_sub_t* coverage_cut;
     ///normalize_ma_hit_t(sources, n_read);
     normalize_ma_hit_t_single_side(sources, n_read);
 
 
-    // debug_info_of_specfic_read("m64016_190918_162737/92668450/ccs", 
-    // sources, reverse_sources, -1, "normalize");
     
-    // debug_info_of_specfic_read("m64016_190918_162737/53545052/ccs", 
-    // sources, reverse_sources, -1, "normalize");
     
 
 
@@ -9119,7 +9603,13 @@ char* output_file_name, long long bubble_dist, int read_graph, int write)
     
     
     
-    // debug_info_of_specfic_read("m64016_190918_162737/49678749/ccs", 
+    // debug_info_of_specfic_read("m64013_190322_203854/82051959/ccs", 
+    // sources, reverse_sources, -1, "clean");
+
+    // debug_info_of_specfic_read("m64013_190322_203854/74385680/ccs", 
+    // sources, reverse_sources, -1, "clean");
+
+    // debug_info_of_specfic_read("m64011_190329_072846/80545633/ccs", 
     // sources, reverse_sources, -1, "clean");
 
 
@@ -9127,6 +9617,9 @@ char* output_file_name, long long bubble_dist, int read_graph, int write)
     
 
     ma_hit_sub(min_dp, sources, n_read, readLen, mini_overlap_length, &coverage_cut);
+
+    ma_hit_chimeric(1, sources, reverse_sources, n_read, readLen, coverage_cut);
+
     ma_hit_cut(min_dp, sources, n_read, readLen, mini_overlap_length, &coverage_cut);
     ///it seems we do not need ma_hit_flt
     ma_hit_flt(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length);
@@ -9136,7 +9629,11 @@ char* output_file_name, long long bubble_dist, int read_graph, int write)
 
     
 
+    // debug_info_of_specfic_read("m64016_190918_162737/179635219/ccs", 
+    // sources, reverse_sources, -1, "contain");
 
+    // debug_info_of_specfic_read("m64016_190918_162737/130811282/ccs", 
+    // sources, reverse_sources, -1, "contain");
     
     // debug_info_of_specfic_read("m64016_190918_162737/72220752/ccs", 
     // sources, reverse_sources, -1, "contain");
@@ -9311,16 +9808,35 @@ char* output_file_name, long long bubble_dist, int read_graph, int write)
     asg_arc_del_triangular_directly(sg, MAX_SHORT_TIPS, reverse_sources);
 
 
+
+
+
+    ///asg_arc_identify_simple_bubbles_multi(sg, 0);
+    // asg_arc_del_chimeric_read(sg, MAX_SHORT_TIPS*2);
+    // asg_cut_tip(sg, MAX_SHORT_TIPS);
+    
+
+
+
+    
     asg_arc_identify_simple_bubbles_multi(sg, 0);
     asg_arc_del_orthology_multiple_way(sg, reverse_sources, 0.4, MAX_SHORT_TIPS);
     asg_cut_tip(sg, MAX_SHORT_TIPS);
 
     
+
+
+
     asg_arc_identify_simple_bubbles_multi(sg, 0);
     asg_arc_del_too_short_overlaps(sg, 2000, min_ovlp_drop_ratio, reverse_sources, MAX_SHORT_TIPS);
     asg_cut_tip(sg, MAX_SHORT_TIPS);
 
 
+    /**
+    asg_arc_identify_simple_bubbles_multi(sg, 1);
+    asg_arc_del_short_false_link_advance(sg, 0.6, 0.85, bubble_dist, reverse_sources, MAX_SHORT_TIPS);
+    **/
+    
     
 
 
@@ -9374,9 +9890,17 @@ char* output_file_name, long long bubble_dist, int read_graph, int write)
     // debug_info_of_specfic_node("m64016_190918_162737/141297762/ccs", sg);
     
     out:
+    ///output_tips(sg, &R_INF);
+
+
     output_unitig_graph(sg, coverage_cut, output_file_name, n_read);
     output_read_graph(sg, coverage_cut, output_file_name, n_read);
 
+    /****************************may have bugs********************************/
+    output_unitig_graph_without_small_bubbles(sg, coverage_cut, output_file_name, n_read, 
+    100000, MAX_SHORT_TIPS);
+    /****************************may have bugs********************************/
+
     output_contig_graph(sg, coverage_cut, output_file_name, n_read, 10000000, MAX_SHORT_TIPS, 0.1, 20,
     reverse_sources, MAX_SHORT_TIPS);
     ///output_contig_graph(sg, coverage_cut, output_file_name, n_read, 10000000);

Overlaps.h

@@ -9,7 +9,8 @@
 ///#define MIN_OVERLAP_LEN 500
 ///#define MIN_OVERLAP_LEN 50
 #define MIN_OVERLAP_LEN 50
-#define MIN_OVERLAP_COVERAGE 1
+///#define MIN_OVERLAP_COVERAGE 1
+#define MIN_OVERLAP_COVERAGE 0
 #define MAX_HANG_LEN 1000
 #define MAX_HANG_PRE 0.8
 #define GAP_FUZZ 1000
@@ -81,6 +82,7 @@ void debug_normalize_ma_hit_t(ma_hit_t_alloc* sources, long long num_sources);
 
 typedef struct {
 	uint32_t s:31, del:1, e;
+	uint8_t c;
 } ma_sub_t;
 
 void ma_hit_sub(int min_dp, ma_hit_t_alloc* sources, long long n_read, uint64_t* readLen, 
@@ -104,6 +106,7 @@ typedef struct {
 
 typedef struct {
 	uint32_t len:31, del:1;
+	uint8_t c;
 } asg_seq_t;
 
 typedef struct {
@@ -324,4 +327,7 @@ long long clean_round, float min_ovlp_drop_ratio, float max_ovlp_drop_ratio,
 float corase_ovlp_drop_ratio, char* output_file_name, long long bubble_dist, int read_graph,
 int write);
 
+void debug_info_of_specfic_read(char* name, ma_hit_t_alloc* sources, 
+ma_hit_t_alloc* reverse_sources, int id, char* command);
+
 #endif

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.cpp

@@ -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__