diff --git a/Assembly.cpp b/Assembly.cpp index b71ef4e..f707edd 100644 --- a/Assembly.cpp +++ b/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,25 +1288,22 @@ overlap_region_alloc* overlap_list, UC_Read* x_read, UC_Read* y_read) - - - if(overlap_list->list[i].y_pos_strand == 0) + if(test_exact == 1) { - recover_UC_Read(y_read, &R_INF, overlap_list->list[i].y_id); + if(overlap_list->list[i].y_pos_strand == 0) + { + recover_UC_Read(y_read, &R_INF, overlap_list->list[i].y_id); + } + else + { + recover_UC_Read_RC(y_read, &R_INF, overlap_list->list[i].y_id); + } + + 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); } - else - { - recover_UC_Read_RC(y_read, &R_INF, overlap_list->list[i].y_id); - } - - 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(¤t_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, - ¤t_cigar, &hap, &second_round, 0, 1, &fully_cov, &abnormal); + ¤t_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,9 +3009,25 @@ 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"); + // } + + + overlap_list.mapped_overlaps_length = 0; inner_j = 0; @@ -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(¤t_sub_buffer); @@ -3045,11 +3267,32 @@ void* Final_overlap_calculate_heap_merge(void* arg) fprintf(stderr, "fully corrected reads: %llu\n", total_potiental_matched_overlap_1); + + + 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); diff --git a/Correct.cpp b/Correct.cpp index 36a7998..746ad3d 100644 --- a/Correct.cpp +++ b/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) @@ -8638,6 +8781,19 @@ uint32_t* SNPs, long long SNPLen, haplotype_evdience_alloc* hap) for (i = 0; i < SNPLen; i++) { 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,8 +9767,9 @@ 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 ///filter snps that are highly likly false if(hap->available_snp > 1) @@ -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", - Get_NAME_LENGTH((*R_INF),readID), Get_NAME((*R_INF),readID)); + fprintf(stderr, "\n\n****************ref_read: %.*s, id: %d****************\n", + 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,18 +10912,12 @@ 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) { diff --git a/Correct.h b/Correct.h index 0b6172c..3ef973a 100644 --- a/Correct.h +++ b/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); diff --git a/Hash_Table.cpp b/Hash_Table.cpp index 1b2f7e1..d1d6036 100644 --- a/Hash_Table.cpp +++ b/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 && - list->list[list->length - 1].y_id==tmp->y_id - ) + if (list->length!=0 && 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++; + if(n_skip > max_skip) { - n_skip--; + break; } } + /****************************may have bugs********************************/ } dp->score[i] = max_score; @@ -1331,16 +1369,34 @@ double band_width_threshold) ///debug_chain(a, a_n, dp); + 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) @@ -1562,9 +1618,80 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_thresh continue; } + + // 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)); diff --git a/Hash_Table.h b/Hash_Table.h index 18e592e..8b48caf 100644 --- a/Hash_Table.h +++ b/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; diff --git a/Levenshtein_distance.h b/Levenshtein_distance.h index 81ffbe0..ea0bd0e 100644 --- a/Levenshtein_distance.h +++ b/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,8 +199,22 @@ 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,13 +481,30 @@ 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) { return return_site; } + + + /****************************may have bugs********************************/ + if((ungap_error<=errthold) && (ungap_error == (*return_err))) + { + return_site = site + errthold; + } + /****************************may have bugs********************************/ ////注意，这里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,8 +937,40 @@ 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; } diff --git a/Makefile b/Makefile index b1b7770..c83c10b 100644 --- a/Makefile +++ b/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 \ No newline at end of file diff --git a/Overlaps.cpp b/Overlaps.cpp index 9ffce29..2fd624e 100644 --- a/Overlaps.cpp +++ b/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_pop_bubble(sg, bubble_dist), sg->n_seq + sg->n_arc); - 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_identify_simple_bubbles_multi(sg, 1); + asg_arc_del_short_false_link(sg, 0.6, 0.85, bubble_dist, reverse_sources, MAX_SHORT_TIPS); ///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"); @@ -9310,16 +9807,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); diff --git a/Overlaps.h b/Overlaps.h index 8054576..9538ee8 100644 --- a/Overlaps.h +++ b/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 \ No newline at end of file diff --git a/ksw2.h b/ksw2.h new file mode 100644 index 0000000..213c27f --- /dev/null +++ b/ksw2.h @@ -0,0 +1,177 @@ +#ifndef KSW2_H_ +#define KSW2_H_ + +#include + +#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 +#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 diff --git a/ksw2_extz2_sse.cpp b/ksw2_extz2_sse.cpp new file mode 100644 index 0000000..02bb4c2 --- /dev/null +++ b/ksw2_extz2_sse.cpp @@ -0,0 +1,305 @@ +#include +#include +#include "ksw2.h" + +#ifdef __SSE2__ +#include + +#ifdef KSW_SSE2_ONLY +#undef __SSE4_1__ +#endif + +#ifdef __SSE4_1__ +#include +#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__