From dbc72b8c2de0a82044dde4046f36b5a7af4d002f Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Tue, 28 Dec 2021 21:48:41 -0500 Subject: [PATCH] has bug multiple threads --- Overlaps.cpp | 32 +-- Process_Read.cpp | 506 +++++++++++++++++++++++++++++++++++++---------- Process_Read.h | 10 +- htab.cpp | 2 +- inter.cpp | 57 +++--- 5 files changed, 447 insertions(+), 160 deletions(-) diff --git a/Overlaps.cpp b/Overlaps.cpp index 5980665..20f0bd2 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -31102,9 +31102,14 @@ char *get_outfile_name(char* output_file_name) return buf; } -void create_ul_info() +void create_ul_info(ma_hit_t_alloc* sources, int max_hang, int min_ovlp, long long gap_fuzz) { ug_opt_t opt; memset(&opt, 0, sizeof(opt)); + opt.sources = sources; + opt.max_hang = max_hang; + opt.min_ovlp = min_ovlp; + opt.gap_fuzz = gap_fuzz; + ul_load(&opt); exit(1); } @@ -31117,30 +31122,6 @@ float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name, long long bubble_dist, int read_graph, R_to_U* ruIndex, asg_t **sg_ptr, ma_sub_t **coverage_cut_ptr, int debug_g) { - if(asm_opt.ar) create_ul_info(); - - - - - - - - - - - - - - - - - - - - - - - char *o_file = get_outfile_name(output_file_name); ma_sub_t *coverage_cut = *coverage_cut_ptr; asg_t *sg = *sg_ptr; @@ -31166,6 +31147,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g) { memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads*sizeof(uint8_t)); } + if(asm_opt.ar) create_ul_info(sources, max_hang_length, mini_overlap_length, gap_fuzz); ///print_binned_reads(sources, n_read, coverage_cut); clean_weak_ma_hit_t(sources, reverse_sources, n_read); diff --git a/Process_Read.cpp b/Process_Read.cpp index d9dcde3..38fb881 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -28,7 +28,7 @@ uint8_t seq_nt6_table[256] = { char bit_t_seq_table[256][4] = {{0}}; char bit_t_seq_table_rc[256][4] = {{0}}; char s_H[5] = {'A', 'C', 'G', 'T', 'N'}; -char rc_Table[5] = {'T', 'G', 'C', 'A', 'N'}; +char rc_Table[6] = {'T', 'G', 'C', 'A', 'N', 'N'}; void init_All_reads(All_reads* r) { @@ -341,87 +341,92 @@ void init_UC_Read(UC_Read* r) } -void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID) +void recover_UC_Read_sub_region(char* r, int64_t start_pos, int64_t length, uint8_t strand, All_reads* R_INF, int64_t ID) { - long long readLen = Get_READ_LENGTH((*R_INF), ID); + int64_t readLen = Get_READ_LENGTH((*R_INF), ID); + int64_t end_pos = start_pos + length - 1, begLen, tailLen, offset, mn, src_i, des_i, i; uint8_t* src = Get_READ((*R_INF), ID); - long long i; - long long copyLen; - long long end_pos = start_pos + length - 1; + if(strand == 0) { + offset = start_pos&3; + begLen = 4-offset; + if(begLen > length) begLen = length; + tailLen = (length-begLen)&3; + mn = (length - begLen - tailLen)>>2; + src_i = start_pos; des_i = 0; i = 0; - if (strand == 0) - { - i = start_pos; - copyLen = 0; - - long long initLen = start_pos % 4; - - if (initLen != 0) - { - memcpy(r, bit_t_seq_table[src[i>>2]] + initLen, 4 - initLen); - copyLen = copyLen + 4 - initLen; - i = i + copyLen; + if(begLen > 0) { + memcpy(r+des_i, bit_t_seq_table[src[src_i>>2]]+offset, begLen); + des_i += begLen; src_i += begLen; } - while (copyLen < length) - { - memcpy(r+copyLen, bit_t_seq_table[src[i>>2]], 4); - copyLen = copyLen + 4; - i = i + 4; + for (i = 0; i < mn; i++) { + memcpy(r+des_i, bit_t_seq_table[src[src_i>>2]], 4); + des_i += 4; src_i += 4; } - if (R_INF->N_site[ID]) - { - for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) - { - if ((long long)R_INF->N_site[ID][i] >= start_pos && (long long)R_INF->N_site[ID][i] <= end_pos) - { + if(tailLen > 0) { + memcpy(r+des_i, bit_t_seq_table[src[src_i>>2]], tailLen); + des_i += tailLen; src_i += tailLen; + } + + if (R_INF->N_site[ID]) { + for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) { + if ((long long)R_INF->N_site[ID][i] >= start_pos && (long long)R_INF->N_site[ID][i] <= end_pos) { r[R_INF->N_site[ID][i] - start_pos] = 'N'; } - else if((long long)R_INF->N_site[ID][i] > end_pos) - { + else if((long long)R_INF->N_site[ID][i] > end_pos) { break; } } } } - else - { + else { start_pos = readLen - start_pos - 1; end_pos = readLen - end_pos - 1; - ///start_pos > end_pos - i = start_pos; - copyLen = 0; - long long initLen = (start_pos + 1) % 4; + begLen = (start_pos+1)&3; + offset = 4 - begLen; + if(begLen > length) begLen = length; + tailLen = (length-begLen)&3; + mn = (length - begLen - tailLen)>>2; + src_i = start_pos; des_i = 0; i = 0; - if (initLen != 0) - { - memcpy(r, bit_t_seq_table_rc[src[i>>2]] + 4 - initLen, initLen); - copyLen = copyLen + initLen; - i = i - initLen; + if(begLen > 0) { + memcpy(r+des_i, bit_t_seq_table_rc[src[src_i>>2]]+offset, begLen); + des_i += begLen; src_i -= begLen; } - while (copyLen < length) - { - memcpy(r+copyLen, bit_t_seq_table_rc[src[i>>2]], 4); - copyLen = copyLen + 4; - i = i - 4; + for (i = 0; i < mn; i++) { + memcpy(r+des_i, bit_t_seq_table_rc[src[src_i>>2]], 4); + des_i += 4; src_i -= 4; } - if (R_INF->N_site[ID]) - { - long long offset = readLen - start_pos - 1; + if(tailLen > 0) { + memcpy(r+des_i, bit_t_seq_table_rc[src[src_i>>2]], tailLen); + des_i += tailLen; src_i -= tailLen; + } - for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) - { - if ((long long)R_INF->N_site[ID][i] >= end_pos && (long long)R_INF->N_site[ID][i] <= start_pos) - { + /** + if (R_INF->N_site[ID]) { + offset = readLen - start_pos - 1; + for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) { + if ((long long)R_INF->N_site[ID][i] >= end_pos && (long long)R_INF->N_site[ID][i] <= start_pos) { r[readLen - R_INF->N_site[ID][i] - 1 - offset] = 'N'; } - else if((long long)R_INF->N_site[ID][i] > start_pos) - { + else if((long long)R_INF->N_site[ID][i] > start_pos) { + break; + } + } + } + **/ + if (R_INF->N_site[ID]) { + start_pos = readLen - start_pos - 1; end_pos = readLen - end_pos - 1; + for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) { + offset = readLen - R_INF->N_site[ID][i] - 1; + if(offset >= start_pos && offset <= end_pos) { + r[offset - start_pos] = 'N'; + } else if(offset < start_pos) { break; } } @@ -676,6 +681,7 @@ void reverse_complement(char* pattern, uint64_t length) for (i = 0; i < end; i++) { + index = length - i - 1; k = pattern[index]; pattern[index] = RC_CHAR(pattern[i]); @@ -743,8 +749,9 @@ void destory_Debug_reads(Debug_reads* x) void init_all_ul_t(all_ul_t *x, All_reads *hR) { memset(x, 0, sizeof(*x)); x->hR = hR; x->mm = 0x7fffffff; + init_aux_table(); } -void destory_all_ul_t(all_ul_t *x, All_reads *hR) { +void destory_all_ul_t(all_ul_t *x) { uint64_t i; for (i = 0; i < x->n; i++) { free(x->a[i].n_n); free(x->a[i].N_site.a); @@ -753,7 +760,7 @@ void destory_all_ul_t(all_ul_t *x, All_reads *hR) { free(x->a); } -void ha_encode_base(uint8_t* dest, char* src, uint64_t src_l, N_t *nn) +void ha_encode_base(uint8_t* dest, char* src, uint64_t src_l, N_t *nn, uint64_t nn_offset) { uint64_t i = 0; uint64_t dest_i = 0; @@ -765,28 +772,28 @@ void ha_encode_base(uint8_t* dest, char* src, uint64_t src_l, N_t *nn) c = seq_nt6_table[(uint8_t)src[i]]; if (c >= 4) { - c = 0; kv_push(uint32_t, *nn, i); + c = 0; kv_push(uint32_t, *nn, i+nn_offset); } i++; tmp = tmp | (c<<6); c = seq_nt6_table[(uint8_t)src[i]]; if (c >= 4) { - c = 0; kv_push(uint32_t, *nn, i); + c = 0; kv_push(uint32_t, *nn, i+nn_offset); } i++; tmp = tmp | (c<<4); c = seq_nt6_table[(uint8_t)src[i]]; if (c >= 4) { - c = 0; kv_push(uint32_t, *nn, i); + c = 0; kv_push(uint32_t, *nn, i+nn_offset); } i++; tmp = tmp | (c<<2); c = seq_nt6_table[(uint8_t)src[i]]; if (c >= 4) { - c = 0; kv_push(uint32_t, *nn, i); + c = 0; kv_push(uint32_t, *nn, i+nn_offset); } i++; tmp = tmp | c; @@ -803,7 +810,7 @@ void ha_encode_base(uint8_t* dest, char* src, uint64_t src_l, N_t *nn) while (i < src_l) { c = seq_nt6_table[(uint8_t)src[i]]; if (c >= 4) { - c = 0; kv_push(uint32_t, *nn, i); + c = 0; kv_push(uint32_t, *nn, i+nn_offset); } i++; tmp = tmp | (c << shift); @@ -816,14 +823,33 @@ void ha_encode_base(uint8_t* dest, char* src, uint64_t src_l, N_t *nn) } #define B4L(x) (((x)>>2)+(((x)&3)?1:0)) +void push_subblock_original_bases(char* str, all_ul_t *x, ul_vec_t *p, uint32_t s, uint32_t e, uint32_t subLen)//for debug +{ + uc_block_t *b = NULL; + uint32_t qs = s, qe = e; + while (qs < e) { + qe = qs + subLen; if(qe > e) qe = e; + kv_pushp(uc_block_t, p->bb, &b); + b->hid = x->mm; b->rev = 0; + b->qs = qs; b->qe = qe; + b->ts = p->r_base.n; b->te = b->ts + B4L(b->qe - b->qs); + kv_resize(uint8_t, p->r_base, b->te); p->r_base.n = b->te; + ha_encode_base(p->r_base.a+b->ts, str+b->qs, b->qe-b->qs, &(p->N_site), b->qs); + + qs = qe; + } +} + void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str, int64_t str_l, ul_ov_t *o, int64_t on) { int64_t i, mine, maxs, ovlp, end; ul_vec_t *p = NULL; ul_ov_t *z = NULL, *zp = NULL; uc_block_t *b = NULL; + if(rid) fprintf(stderr, "rid:%lu\n", *rid); if(rid == NULL) { kv_pushp(ul_vec_t, *x, &p); memset(p, 0, sizeof(*p)); + fprintf(stderr, "x->n:%u\n", x->n); } else { p = &(x->a[*rid]); @@ -854,7 +880,7 @@ void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str, b->qs = end; b->qe = b->qs - ovlp; b->ts = p->r_base.n; b->te = b->ts + B4L(-ovlp); kv_resize(uint8_t, p->r_base, b->te); p->r_base.n = b->te; - ha_encode_base(p->r_base.a+b->ts, str+b->qs, b->qe-b->qs, &(p->N_site)); + ha_encode_base(p->r_base.a+b->ts, str+b->qs, b->qe-b->qs, &(p->N_site), b->qs); } ///push ovlp bases @@ -876,72 +902,348 @@ void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str, b->qs = end; b->qe = str_l; b->ts = p->r_base.n; b->te = b->ts + B4L(b->qe - b->qs); kv_resize(uint8_t, p->r_base, b->te); p->r_base.n = b->te; - ha_encode_base(p->r_base.a+b->ts, str+b->qs, b->qe-b->qs, &(p->N_site)); + ha_encode_base(p->r_base.a+b->ts, str+b->qs, b->qe-b->qs, &(p->N_site), b->qs); + // push_subblock_original_bases(str, x, p, end, str_l, 321);//for debug } + // char *sst = NULL; CALLOC(sst, str_l);//for debug + // retrieve_ul_t(NULL, sst, x, rid?*rid:x->n-1, 0, 0, -1); + // if(memcmp(sst, str, str_l)) { + // fprintf(stderr, "ap-Wrong read, id: %ld, [%d, %ld)\n", (int64_t)(rid?*rid:x->n-1), 0, str_l); + // for (i = 0; i < str_l; i++) { + // if(sst[i] != str[i]) fprintf(stderr, "[%ld] input:%c, decompress:%c\n", i, str[i], sst[i]); + // } + // } + // free(sst); } } -void retrieve_ul_t(UC_Read* r, all_ul_t *ref, uint64_t ID, uint8_t strand) { - ul_vec_t *p = &(ref->a[ID]); +void retrieve_ul_t(UC_Read* i_r, char *i_s, all_ul_t *ref, uint64_t ID, uint8_t strand, int64_t s, int64_t l) { + ul_vec_t *p = &(ref->a[ID]); + if(l < 0) l = p->rlen; uc_block_t *b = NULL; - char *a = NULL; + char *r = NULL; uint8_t *src = NULL; - int64_t k, i, a_n, l_chr, idx; - r->length = p->rlen; r->RID = ID; + int64_t k, i, a_n, e = s + l, ssp, sep, sl, rts, rte; + int64_t offset, begLen, tailLen, src_i, des_i; + if(i_r) { + i_r->length = l; i_r->RID = ID; + if(i_r->length > i_r->size) { + i_r->size = i_r->length; + i_r->seq = (char*)realloc(i_r->seq,sizeof(char)*(i_r->size)); + } + r = i_r->seq; + } + if(i_s) r = i_s; + + - if(r->length + 4 > r->size) { - r->size = r->length + 4; - r->seq = (char*)realloc(r->seq,sizeof(char)*(r->size)); - } - - a = NULL; if(strand == 0) { - for (k = 0; k < (int64_t)p->bb.n; k++) { - b = &(p->bb.a[k]); a = r->seq + b->qs; a_n = b->qe - b->qs; src = p->r_base.a + b->ts; + for (k = 0, des_i = 0; k < (int64_t)p->bb.n; k++) { + b = &(p->bb.a[k]); + if(b->qe <= s) continue; + if(b->qs >= e) break; + src = p->r_base.a + b->ts; + ssp = MAX(s, b->qs) - b->qs; + sep = MIN(e, b->qe) - b->qs; + sl = sep - ssp; + if(b->hid == ref->mm){///original bases - i = 0; - while (i < a_n) { - memcpy(a+i, bit_t_seq_table[src[i>>2]], 4); - i += 4; + offset = ssp&3; + begLen = 4-offset; + if(begLen > sl) begLen = sl; + tailLen = (sl-begLen)&3; + a_n = (sl - begLen - tailLen)>>2; + + src_i = ssp; i = 0; + if(begLen > 0) { + memcpy(r+des_i, bit_t_seq_table[src[src_i>>2]]+offset, begLen); + des_i += begLen; src_i += begLen; + } + + for (i = 0; i < a_n; i++) { + memcpy(r+des_i, bit_t_seq_table[src[src_i>>2]], 4); + des_i += 4; src_i += 4; + } + + if(tailLen > 0) { + memcpy(r+des_i, bit_t_seq_table[src[src_i>>2]], tailLen); + des_i += tailLen; src_i += tailLen; } } else {///ovlps - recover_UC_Read_sub_region(a, b->ts, b->te - b->ts, b->rev, ref->hR, b->hid); + if(b->rev == 0) { + recover_UC_Read_sub_region(r+des_i, b->ts + ssp, sep - ssp, b->rev, ref->hR, b->hid); + } else{ + rts = b->ts + (b->qe - sep); rte = rts + sep - ssp; + recover_UC_Read_sub_region(r+des_i, Get_READ_LENGTH((*ref->hR), b->hid) - rte, + sep - ssp, b->rev, ref->hR, b->hid); + } + des_i += sep - ssp; + } + } + for (k = 0; k < (int64_t)p->N_site.n; k++) { + if(p->N_site.a[k] >= s && p->N_site.a[k] < e){ + r[p->N_site.a[k]-s] = 'N'; + } + else if(p->N_site.a[k] >= e) { + break; } } - for (k = 0; k < (int64_t)p->N_site.n; k++) r->seq[p->N_site.a[k]] = 'N'; } else { - for (k = 0; k < (int64_t)p->bb.n; k++) { - b = &(p->bb.a[k]); a = r->seq + r->length - b->qe; a_n = b->qe - b->qs; src = p->r_base.a + b->ts; + sep = p->rlen - s; + ssp = p->rlen - e; + s = ssp; e = sep; + ///[s, e) + for (k = ((int64_t)p->bb.n)-1, des_i = 0; k >= 0; k--) { + b = &(p->bb.a[k]); + if(b->qe <= s) break; + if(b->qs >= e) continue; + src = p->r_base.a + b->ts; + ssp = MAX(s, b->qs) - b->qs; + sep = MIN(e, b->qe) - b->qs; + sl = sep - ssp; + ///[ssp, sep) + if(b->hid == ref->mm){///original bases - l_chr = (b->qe - b->qs)&((uint32_t)3); - i = ((b->qe - b->qs)>>2)+(l_chr!= 0)-1; - idx = 0; - if(l_chr > 0) { - memcpy(a + idx, bit_t_seq_table_rc[src[i]]+4-l_chr, l_chr); - idx = l_chr; i--; + begLen = sep&3; + offset = 4 - begLen; + if(begLen > sl) begLen = sl; + tailLen = (sl-begLen)&3; + a_n = (sl - begLen - tailLen)>>2; + src_i = sep-1; i = 0; + + if(begLen > 0) { + memcpy(r+des_i, bit_t_seq_table_rc[src[src_i>>2]]+offset, begLen); + des_i += begLen; src_i -= begLen; } - while (i >= 0) { - memcpy(a + idx, bit_t_seq_table_rc[src[i]], 4); - i--; idx += 4; + + for (i = 0; i < a_n; i++) { + memcpy(r+des_i, bit_t_seq_table_rc[src[src_i>>2]], 4); + des_i += 4; src_i -= 4; + } + + if(tailLen > 0) { + memcpy(r+des_i, bit_t_seq_table_rc[src[src_i>>2]], tailLen); + des_i += tailLen; src_i -= tailLen; } } else {///ovlps - recover_UC_Read_sub_region(a, b->ts, b->te - b->ts, b->rev==strand?0:1, ref->hR, b->hid); + if(b->rev == 0) {///b->rev != strand + rts = b->ts + ssp; rte = rts + sep - ssp; + recover_UC_Read_sub_region(r+des_i, Get_READ_LENGTH((*ref->hR), b->hid) - rte, sep - ssp, 1, ref->hR, b->hid); + } else {///b->rev == strand + rts = b->ts + (b->qe - sep); rte = rts + sep - ssp; + recover_UC_Read_sub_region(r+des_i, rts, sep - ssp, 0, ref->hR, b->hid); + } + des_i += sep - ssp; } } - for (k = 0; k < (int64_t)p->N_site.n; k++) r->seq[r->length - p->N_site.a[k] - 1] = 'N'; + + sep = p->rlen - s; + ssp = p->rlen - e; + s = ssp; e = sep; + for (k = 0; k < (int64_t)p->N_site.n; k++) { + sl = p->rlen - p->N_site.a[k] - 1; + if(sl >= s && sl < e) r[sl-s] = 'N'; + else if(sl < s) { + break; + } + } } } -void retrieve_u_seq(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand) + +void retrieve_u_seq(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand, int64_t s, int64_t l) { + if(u->m == 0 || u->n == 0) return; + if(l < 0) l = u->len; + char *r = NULL, *a = NULL; + int64_t e = s + l, ssp, sep, rs, re, des_i; + uint64_t k, rId, ori, r_l; if(i_r) { + i_r->length = l; i_r->RID = 0; + if(i_r->length > i_r->size) { + i_r->size = i_r->length; + i_r->seq = (char*)realloc(i_r->seq,sizeof(char)*(i_r->size)); + } + r = i_r->seq; + } + if(i_s) r = i_s; + + if(strand == 1) { + sep = u->len - s; + ssp = u->len - e; + s = ssp; e = sep; + } + for (k = l = des_i = 0; k < u->n; k++) { + rId = u->a[k]>>33; + ori = u->a[k]>>32&1; + r_l = (uint32_t)u->a[k]; + if(r_l == 0) continue; + ssp = l; sep = l + r_l; + l += r_l; + if(sep <= s) continue; + if(ssp >= e) break; + rs = MAX(ssp, s); re = MIN(sep, e); + a = r + des_i; des_i += re - rs; + recover_UC_Read_sub_region(a, rs-ssp, re-rs, ori, &R_INF, rId); + } + if(strand == 1) { + char t; + re = (e - s); + l = re>>1; + for (k = 0; k < (uint64_t)l; k++) { + des_i = re - k - 1; + t = r[des_i]; + r[des_i] = RC_CHAR(r[k]); + r[k] = RC_CHAR(t); + } + if(re&1) r[l] = RC_CHAR(r[l]); + } +} + +void produce_u_seq(char* r, ma_utg_t *u, UC_Read *buf) +{ + if(u->m == 0 || u->n == 0) return; + uint32_t j, k, l = 0; + uint32_t rId, ori, start, eLen, readLen; + char *readS = NULL; + memset(r, 'N', u->len); + for (j = 0; j < u->n; ++j) { + rId = u->a[j]>>33; + ///uId = i; + ori = u->a[j]>>32&1; + start = l; + eLen = (uint32_t)u->a[j]; + l += eLen; + + if(eLen == 0) continue; + recover_UC_Read(buf, &R_INF, rId); + + readS = buf->seq; + readLen = Get_READ_LENGTH(R_INF, rId); + if (!ori) // forward strand + { + for (k = 0; k < eLen; k++) + { + r[start + k] = readS[k]; + } + } + else + { + for (k = 0; k < eLen; k++) + { + uint8_t c = (uint8_t)readS[readLen - 1 - k]; + r[start + k] = c >= 128? 'N' : RC_CHAR(c); + } + } + } +} + +void debug_retrieve_rc_sub(all_ul_t *ref, const All_reads *R_INF, ma_utg_v *u, uint32_t n_step) +{ + uint64_t i, step, s, e, occ; + UC_Read f, r; + init_UC_Read(&f); init_UC_Read(&r); + kvec_t(char) ss; kv_init(ss); + if(ref) { + for (i = 0, occ = 0; i < ref->n; i++) { + retrieve_ul_t(&f, NULL, ref, i, 0, 0, -1); + retrieve_ul_t(&r, NULL, ref, i, 1, 0, -1); + + kv_resize(char, ss, ref->a[i].rlen); + ss.n = ref->a[i].rlen; + memcpy(ss.a, r.seq, ss.n); + reverse_complement(ss.a, ss.n); + if(memcmp(ss.a, f.seq, ss.n)) { + fprintf(stderr, "1-Wrong whole reverse-read, id: %lu\n", i); + // for (s = 0; s < ss.n; s++) { + // if(ss.a[s] != f.seq[s]) { + // fprintf(stderr,"s:%lu, ss.a[s]:%c, f.seq[s]:%c, r.seq[rs]:%c\n", s, ss.a[s], f.seq[s], r.seq[ref->a[i].rlen - s - 1]); + // } + // } + } + + step = ss.n/n_step; + if(step <= 0) step = 1; + + for (s = 0; s < ss.n; s += step) { + e = MIN(s+step, ss.n); + retrieve_ul_t(NULL, ss.a, ref, i, 0, s, e-s); + if(memcmp(ss.a, f.seq + s, e - s)) { + fprintf(stderr, "1-Wrong sub forward-read, id: %lu, [%lu, %lu)\n", i, s, e); + } + + retrieve_ul_t(NULL, ss.a, ref, i, 1, s, e-s); + if(memcmp(ss.a, r.seq + s, e - s)) { + fprintf(stderr, "1-Wrong sub reverse-read, id: %lu, [%lu, %lu)\n", i, s, e); + // uint64_t dk; char *tf = ss.a; char *rf = r.seq + s; + // for (dk = 0; dk < e - s; dk++) { + // if(tf[dk] != rf[dk]) { + // fprintf(stderr,"dk:%lu, tf[dk]:%c, rf[dk]:%c\n", + // dk, tf[dk], rf[dk]); + // } + // } + } + occ++; + } + } + fprintf(stderr, "[M::%s::# checking: %lu] ==> all_ul_t\n", __func__, occ); } - if (r->length + 4 > r->size) - { - r->size = r->length + 4; - r->seq = (char*)realloc(r->seq,sizeof(char)*(r->size)); + if(R_INF) { + for (i = 0, occ = 0; i < R_INF->total_reads; i++) { + recover_UC_Read(&f, R_INF, i); + recover_UC_Read_RC(&r, (All_reads*)R_INF, i); + kv_resize(char, ss, Get_READ_LENGTH((*R_INF), i)); + ss.n = Get_READ_LENGTH((*R_INF), i); + memcpy(ss.a, r.seq, ss.n); + reverse_complement(ss.a, ss.n); + if(memcmp(ss.a, f.seq, ss.n)) fprintf(stderr, "2-Wrong whole reverse-read, id: %lu\n", i); + + step = ss.n/n_step; + if(step <= 0) step = 1; + for (s = 0; s < ss.n; s += step) { + e = MIN(s+step, ss.n); + recover_UC_Read_sub_region(ss.a, s, e-s, 0, (All_reads*)R_INF, i); + if(memcmp(ss.a, f.seq + s, e - s)) fprintf(stderr, "2-Wrong sub forward-read, id: %lu, [%lu, %lu)\n", i, s, e); + + recover_UC_Read_sub_region(ss.a, s, e-s, 1, (All_reads*)R_INF, i); + if(memcmp(ss.a, r.seq + s, e - s)) fprintf(stderr, "2-Wrong sub reverse-read, id: %lu, [%lu, %lu)\n", i, s, e); + occ++; + } + } + fprintf(stderr, "[M::%s::# checking: %lu] ==> All_reads\n", __func__, occ); } + + if(u) { + for (i = 0, occ = 0; i < u->n; i++) { + kv_resize(char, ss, u->a[i].len); ss.n = u->a[i].len; + retrieve_u_seq(&f, NULL, &(u->a[i]), 0, 0, -1); + produce_u_seq(ss.a, &(u->a[i]), &r); + if(memcmp(ss.a, f.seq, ss.n)) fprintf(stderr, "4-Wrong whole reverse-read, id: %lu\n", i); + retrieve_u_seq(&r, NULL, &(u->a[i]), 1, 0, -1); + + memcpy(ss.a, r.seq, ss.n); + reverse_complement(ss.a, ss.n); + if(memcmp(ss.a, f.seq, ss.n)) fprintf(stderr, "3-Wrong whole reverse-read, id: %lu\n", i); + + step = ss.n/n_step; + if(step <= 0) step = 1; + + for (s = 0; s < ss.n; s += step) { + e = MIN(s+step, ss.n); + retrieve_u_seq(NULL, ss.a, &(u->a[i]), 0, s, e-s); + if(memcmp(ss.a, f.seq + s, e - s)) fprintf(stderr, "3-Wrong sub forward-read, id: %lu, [%lu, %lu)\n", i, s, e); + + retrieve_u_seq(NULL, ss.a, &(u->a[i]), 1, s, e-s); + if(memcmp(ss.a, r.seq + s, e - s)) fprintf(stderr, "3-Wrong sub reverse-read, id: %lu, [%lu, %lu)\n", i, s, e); + occ++; + } + } + fprintf(stderr, "[M::%s::# checking: %lu] ==> ma_utg_v\n", __func__, occ); + } + + destory_UC_Read(&f); destory_UC_Read(&r); + kv_destroy(ss); } \ No newline at end of file diff --git a/Process_Read.h b/Process_Read.h index 49a7986..d77017d 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -29,7 +29,7 @@ extern uint8_t seq_nt6_table[256]; extern char bit_t_seq_table[256][4]; extern char bit_t_seq_table_rc[256][4]; extern char s_H[5]; -extern char rc_Table[5]; +extern char rc_Table[6]; #define RC_CHAR(x) rc_Table[seq_nt6_table[(uint8_t)x]] @@ -198,7 +198,7 @@ void ha_compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_sit void init_UC_Read(UC_Read* r); void recover_UC_Read(UC_Read* r, const All_reads *R_INF, uint64_t ID); void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID); -void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID); +void recover_UC_Read_sub_region(char* r, int64_t start_pos, int64_t length, uint8_t strand, All_reads* R_INF, int64_t ID); void destory_UC_Read(UC_Read* r); void reverse_complement(char* pattern, uint64_t length); void write_All_reads(All_reads* r, char* read_file_name); @@ -210,8 +210,10 @@ void destory_Debug_reads(Debug_reads* x); void recover_UC_sub_Read(UC_Read* i_r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID); void init_all_ul_t(all_ul_t *x, All_reads *hR); -void destory_all_ul_t(all_ul_t *x, All_reads *hR); +void destory_all_ul_t(all_ul_t *x); void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str, int64_t str_l, ul_ov_t *o, int64_t on); -void retrieve_ul_t(UC_Read* r, all_ul_t *ref, uint64_t ID, uint8_t strand); +void retrieve_ul_t(UC_Read* i_r, char *i_s, all_ul_t *ref, uint64_t ID, uint8_t strand, int64_t s, int64_t l); +void retrieve_u_seq(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand, int64_t s, int64_t l); +void debug_retrieve_rc_sub(all_ul_t *ref, const All_reads *R_INF, ma_utg_v *u, uint32_t n_step); #endif diff --git a/htab.cpp b/htab.cpp index 97b3d18..6d36cf5 100644 --- a/htab.cpp +++ b/htab.cpp @@ -733,7 +733,7 @@ static void *sf##_worker_count(void *data, int step, void *in) /** callback for }\ MALLOC(s->seq[s->n_seq], u->len);\ if(u->s) memcpy(s->seq[s->n_seq], u->s, u->len);\ - else memcpy(s->seq[s->n_seq], u->s, u->len);\ + else retrieve_u_seq(NULL, s->seq[s->n_seq], u, 0, 0, -1);\ s->len[s->n_seq++] = u->len;\ ++p->n_seq;\ s->sum_len += u->len;\ diff --git a/inter.cpp b/inter.cpp index e3ab5ea..02616f1 100644 --- a/inter.cpp +++ b/inter.cpp @@ -288,10 +288,10 @@ typedef struct { // data structure for each step in kt_pipeline() int n, m, sum_len; uint64_t *len, id; char **seq; - ha_mzl_v *mzs; - st_mt_t *sps; - mg_gchains_t **gcs; - mg_tbuf_t **buf; + ha_mzl_v *mzs;///useless + st_mt_t *sps;///useless + mg_gchains_t **gcs;///useless + mg_tbuf_t **buf;///useless ha_ovec_buf_t **hab; } utepdat_t; @@ -2320,6 +2320,7 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac } **/ for (i = 0; i < p->n_thread; ++i) ha_ovec_destroy(s->hab[i]); + free(s->hab); // free(s->buf); free(s->mzs); free(s->sps); return s; } @@ -2337,13 +2338,14 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac **/ rid = s->id + i; append_ul_t(&UL_INF, &rid, NULL, 0, s->seq[i], s->len[i], NULL, 0); + // fprintf(stderr, "%.*s\n", (int)s->len[i], s->seq[i]); free(s->seq[i]); p->total_base += s->len[i]; } ///debug /** free(s->gcs); **/ - free(s); + free(s->len); free(s->seq); free(s); } return 0; } @@ -2380,7 +2382,7 @@ int print_ul_rs(all_ul_t *U_INF) init_UC_Read(&ur); for (i = 0; i < U_INF->n; i++) { p = &(U_INF->a[i]); - retrieve_ul_t(&ur, U_INF, i, 0); + retrieve_ul_t(&ur, NULL, U_INF, i, 0, 0, -1); fprintf(stderr, ">%s\n", p->n_n); fprintf(stderr, "%.*s\n", (int)ur.length, ur.seq); } @@ -2870,19 +2872,19 @@ void ul_resolve(ma_ug_t *ug, const asg_t *rg, const ug_opt_t *uopt, int hap_n) uidx_destory(); } -int ul_v_call(mg_idxopt_t *opt, const ug_opt_t *uopt, const asg_t *rg, const enzyme *fn, void *ha_flt_tab, ha_pt_t *ha_idx, ma_ug_t *ug) +int ul_v_call(mg_idxopt_t *opt, const ug_opt_t *uopt, const enzyme *fn, void *ha_flt_tab, ha_pt_t *ha_idx, ma_ug_t *ug) { uldat_t sl; memset(&sl, 0, sizeof(sl)); sl.ha_flt_tab = ha_flt_tab; sl.ha_idx = ha_idx; sl.opt = opt; - sl.chunk_size = 500000000; + sl.chunk_size = 100000000; sl.n_thread = asm_opt.thread_num; sl.ug = ug; - sl.rg = rg; sl.uopt = uopt; scall_ul_pipeline(&sl, fn); - print_ul_rs(&UL_INF); + // print_ul_rs(&UL_INF); + debug_retrieve_rc_sub(&UL_INF, &R_INF, &(ug->u), 100); // if(!load_ul_hits(&sl.hits, &sl.nn, asm_opt.output_file_name)) { // scall_ul_pipeline(&sl, fn); // write_ul_hits(&sl.hits, &sl.nn, asm_opt.output_file_name); @@ -2891,13 +2893,13 @@ int ul_v_call(mg_idxopt_t *opt, const ug_opt_t *uopt, const asg_t *rg, const enz return 1; } -ma_ug_t *dedup_HiFis(ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, uint64_t n_read, int64_t readLen) +ma_ug_t *dedup_HiFis(ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, int64_t gap_fuzz) { - uint64_t i, k, qn, tn; + uint64_t i, k, qn, tn, n_read = R_INF.total_reads; int32_t r; asg_arc_t t, *p = NULL; asg_t *rg = asg_init(); - rg->m_seq = n_read; MALLOC(rg->seq, rg->m_seq); + rg->m_seq = rg->n_seq = n_read; MALLOC(rg->seq, rg->m_seq); for (i = 0; i < n_read; ++i) rg->seq[i].len = Get_READ_LENGTH(R_INF, i); for (i = 0; i < n_read; i++) { @@ -2906,6 +2908,8 @@ ma_ug_t *dedup_HiFis(ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, ui if(!src[i].buffer[k].el) continue; qn = Get_qn(src[i].buffer[k]); tn = Get_tn(src[i].buffer[k]); if(rg->seq[qn].del || rg->seq[tn].del) continue; + if((Get_qe(src[i].buffer[k]) - Get_qs(src[i].buffer[k])) < min_ovlp) continue; + if((Get_te(src[i].buffer[k]) - Get_ts(src[i].buffer[k])) < min_ovlp) continue; r = ma_hit2arc(&(src[i].buffer[k]), rg->seq[qn].len, rg->seq[tn].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); if (r == MA_HT_QCONT) { rg->seq[qn].del = 1; @@ -2922,6 +2926,8 @@ ma_ug_t *dedup_HiFis(ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, ui if(!src[i].buffer[k].el) continue; qn = Get_qn(src[i].buffer[k]); tn = Get_tn(src[i].buffer[k]); if(rg->seq[qn].del || rg->seq[tn].del) continue; + if((Get_qe(src[i].buffer[k]) - Get_qs(src[i].buffer[k])) < min_ovlp) continue; + if((Get_te(src[i].buffer[k]) - Get_ts(src[i].buffer[k])) < min_ovlp) continue; r = ma_hit2arc(&(src[i].buffer[k]), rg->seq[qn].len, rg->seq[tn].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); if (r >= 0) { p = asg_arc_pushp(rg); @@ -2931,34 +2937,29 @@ ma_ug_t *dedup_HiFis(ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, ui } asg_cleanup(rg); asg_symm(rg); - - - - + asg_arc_del_trans(rg, gap_fuzz); ma_ug_t *ug = NULL; ug = ma_ug_gen(rg); + asg_destroy(rg); - - - - - - asg_destroy(rg); ma_ug_destroy(ug); + for (i = k = 0; i < ug->u.n; i++) k += ug->u.a[i].len; + fprintf(stderr, "[M::%s::] # unitigs: %lu, # bases: %lu\n", __func__, (uint64_t)ug->u.n, k); + return ug; } void ul_load(const ug_opt_t *uopt) { fprintf(stderr, "[M::%s::] ==> UL\n", __func__); mg_idxopt_t opt; - ma_ug_t *ug = NULL; + ma_ug_t *ug = dedup_HiFis(uopt->sources, uopt->min_ovlp, uopt->max_hang, uopt->gap_fuzz); int cutoff = asm_opt.hom_cov * asm_opt.high_factor; init_aux_table(); init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff); + /** int exist = (asm_opt.load_index_from_disk? uidx_load(&ha_flt_tab, &ha_idx, asm_opt.output_file_name) : 0); if(exist == 0) uidx_l_build(ug, &opt, cutoff); if(exist == 0) uidx_write(ha_flt_tab, ha_idx, asm_opt.output_file_name); - - - - ul_v_call(&opt, uopt, /**rg**/NULL, asm_opt.ar, ha_flt_tab, ha_idx, /**ug**/NULL); + **/ + ul_v_call(&opt, uopt, asm_opt.ar, ha_flt_tab, ha_idx, ug); + ma_ug_destroy(ug); destory_all_ul_t(&UL_INF); } \ No newline at end of file