diff --git a/Correct.cpp b/Correct.cpp index 34c4450..5bcc376 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -290,7 +290,7 @@ int* r_extra_begin, int* r_extra_end, long long* r_y_start, long long* r_y_lengt ///since Window_Len == x_len + (threshold << 1) if(y_start < 0 || currentIDLen <= y_start || - currentIDLen - y_start + 2 * threshold + THRESHOLD_MAX_SIZE < Window_Len) + currentIDLen - y_start + 2 * threshold + THRESHOLD_MAX_SIZE < Window_Len)///if ylen is too small { return 0; } @@ -1770,9 +1770,9 @@ inline int move_gap_greedy(char* path, int path_i, int path_length, char* x, int /** * - GGCG-TGTGCCTGT + GGCG-TGTGCCTGT [Y] * - GGCAATGTGCCTGT + GGCAATGTGCCTGT [X] * 00013000000000 **/ @@ -1782,7 +1782,7 @@ inline int move_gap_greedy(char* path, int path_i, int path_length, char* x, int char oper = path[path_i]; - if(oper == 3)///there are more x + if(oper == 3)///there are more x; path[] is in reverse order { path_i++; y_i--; @@ -3197,7 +3197,7 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea // z->w_list.a[i].x_start, z->w_list.a[i].x_end, w_e); // } assert(z->w_list.a[i].x_end == w_e); - total_y_start = z->w_list.a[i].y_end + 1 - z->w_list.a[i].extra_begin; + total_y_start = z->w_list.a[i].y_end + 1 - z->w_list.a[i].extra_begin;///since y_end has extra_begin for (k = w_id + 1; k < nw; k++) { if(w_idx[k] != (uint64_t)-1) break; w_s = w_e + 1; @@ -9287,7 +9287,7 @@ int insert_snp_ee(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a uint32_t is_homopolymer = (uint32_t)-1; if(occ_0 == 0 || diff <= 1) return 0; for (i = m = 0; i < 4; i++) { - if(occ_1[i] >= 2){ + if(occ_1[i] >= 2){///Improve: it would be better to double check the bases around this snp at different overlapped reads to make sure this snp is real, instead of by chance if(!km) kv_pushp(SnpStats, h->snp_stat, &p); else kv_pushp_km(km, SnpStats, h->snp_stat, &p); p->id = h->snp_stat.n-1; @@ -9674,7 +9674,7 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, int flag = 0; - while(get_Window(&w_inf, &window_start, &window_end) && flag != -2) + while(get_Window(&w_inf, &window_start, &window_end) && flag != -2)///Improve: discard alignments with too many diff early { dumy->length = 0; dumy->lengthNT = 0; @@ -9700,9 +9700,12 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, // recalcate_window(overlap_list, R_INF, g_read, dumy, overlap_read); // partition_overlaps(overlap_list, R_INF, g_read, dumy, hap, force_repeat); + ///Improve: we should also use the base quality to avoid additional base-pair dp or even for the chaining + //////Improve: we may also consider to do quick conesensus among candidate overlapped reads to make sure if inforamtive sites are real recalcate_window_advance(overlap_list, R_INF, NULL, g_read, dumy, overlap_read, v_idx, w_inf.window_length, asm_opt.max_ov_diff_ec, asm_opt.max_ov_diff_final); // fprintf(stderr, "[M::%s-beg] occ[0]->%lu, occ[1]->%lu, occ[2]->%lu, occ[3]->%lu\n", __func__, // ovlp_occ(overlap_list, 0), ovlp_occ(overlap_list, 1), ovlp_occ(overlap_list, 2), ovlp_occ(overlap_list, 3)); + ///Improve: only keeps specific number of trans overlaps to save memory partition_overlaps_advance(overlap_list, R_INF, g_read, overlap_read, dumy, hap, force_repeat); // fprintf(stderr, "[M::%s-after] occ[0]->%lu, occ[1]->%lu, occ[2]->%lu, occ[3]->%lu\n", __func__, // ovlp_occ(overlap_list, 0), ovlp_occ(overlap_list, 1), ovlp_occ(overlap_list, 2), ovlp_occ(overlap_list, 3)); diff --git a/Hash_Table.cpp b/Hash_Table.cpp index 0626305..1108834 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -1134,7 +1134,7 @@ void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_r { continue; } - + ///Improve: direct filter out too sparse chain before chaining chain_DP(candidates->list + sub_region_beg, sub_region_end - sub_region_beg + 1, &(candidates->chainDP), f_cigar, band_width_threshold, 25, /**Get_READ_LENGTH((*R_INF), (*f_cigar).x_id)**/readLength, diff --git a/anchor.cpp b/anchor.cpp index 24b553a..cd7895e 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -126,7 +126,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg ha_mz1_t *z = &ab->mz.a[i]; seed1_t *s = &ab->seed[i]; for (j = 0; j < s->n; ++j) { - const ha_idxpos_t *y = &s->a[j]; + const ha_idxpos_t *y = &s->a[j];///Improve: filter out y->rid == rid anchor1_t *an = &ab->a[k++]; uint8_t rev = z->rev == y->rev? 0 : 1; an->other_off = y->pos; @@ -152,7 +152,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg cl->size = ab->m_a; REALLOC(cl->list, cl->size); } - for (k = 0; k < ab->n_a; ++k) { + for (k = 0; k < ab->n_a; ++k) {///Improve: direct use ab->a instead of cl->list k_mer_hit *p = &cl->list[k]; p->readID = ab->a[k].srt >> 33; p->strand = ab->a[k].srt >> 32 & 1; @@ -184,7 +184,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg } #endif - if ((int)overlap_list->length > max_n_chain) { + if ((int)overlap_list->length > max_n_chain) {///Improve: use HEAP to directly avoid chaining with too small scores int32_t w, n[4], s[4]; n[0] = n[1] = n[2] = n[3] = 0, s[0] = s[1] = s[2] = s[3] = 0; ks_introsort_or_ss(overlap_list->length, overlap_list->list);