xxay/hifiasm
1
0
Fork 0
mirror of https://github.com/chhylp123/hifiasm.git synced 2026-10-03 03:28:11 +08:00
Code Issues Packages Projects Releases Wiki Activity

new ec model

Browse Source
This commit is contained in:
chhylp123
2024-06-30 19:53:51 -04:00
parent 70fd9a0b1f
commit d5f8a8a6c0
12 changed files with 2746 additions and 52 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Assembly.cpp
+94 -4
View File
@@ -12,6 +12,7 @@
#include "kthread.h"
#include "rcut.h"
#include "kalloc.h"
#include "ecovlp.h"
void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres,
int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp);
@@ -595,9 +596,9 @@ static void worker_ovec(void *data, long i, int tid)
{
ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
int fully_cov, abnormal;
// if(i != 33) return;
// if(i != 12578) return;
// fprintf(stderr, "[M::%s-beg] rid->%ld\n", __func__, i);
// if (memcmp("m64012_190920_173625/88015004/ccs", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
// if (memcmp("7897e875-76e5-42c8-bc37-94b370c4cc8d", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
// fprintf(stderr, "[M::%s-beg] rid->%ld\n", __func__, i);
// } else {
// return;
@@ -606,6 +607,9 @@ static void worker_ovec(void *data, long i, int tid)
ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist,
0.02, asm_opt.max_n_chain, 1, NULL/**&(b->k_flag)**/, &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL, &(b->sp));
// prt_chain(&b->olist);
// return;
clear_Cigar_record(&b->cigar1);
clear_Round2_alignment(&b->round2);
@@ -967,6 +971,88 @@ void prt_dbg_rs(FILE *fp, Debug_reads* x, uint64_t round)
}
void ha_ec(int64_t round)
{
int i, hom_cov, het_cov, r_out = 0;
ec_ovec_buf_t *b = NULL;
ha_ecsave_buf_t *e = NULL;
ha_flt_tab_hp = ha_idx_hp = NULL;
if((ha_idx == NULL)&&(asm_opt.flag & HA_F_VERBOSE_GFA)&&(round == asm_opt.number_of_round - 1)) r_out = 1;
if(asm_opt.required_read_name) init_Debug_reads(&R_INF_FLAG, asm_opt.required_read_name); // for debugging only
// overlap and correct reads
b = gen_ec_ovec_buf_t(asm_opt.thread_num, 0, (round == asm_opt.number_of_round - 1));
// CALLOC(b, asm_opt.thread_num);
// for (i = 0; i < asm_opt.thread_num; ++i)
// b[i] = ha_ovec_init(0, (round == asm_opt.number_of_round - 1),0);
if(ha_idx) hom_cov = asm_opt.hom_cov;
if(ha_idx == NULL) ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, 0, &R_INF, &hom_cov, &het_cov); // build the index
///debug_adapter(&asm_opt, &R_INF);
if (round == 0 && ha_flt_tab == 0) // then asm_opt.hom_cov hasn't been updated
ha_opt_update_cov(&asm_opt, hom_cov);
het_cnt = NULL;
if(round == asm_opt.number_of_round-1 && asm_opt.is_dbg_het_cnt) CALLOC(het_cnt, R_INF.total_reads);
/**
// fprintf(stderr, "[M::%s-start]\n", __func__);
if (asm_opt.required_read_name)
kt_for(asm_opt.thread_num, worker_ovec_related_reads, b, R_INF.total_reads);
else
kt_for(asm_opt.thread_num, worker_ovec, b, R_INF.total_reads);///debug_for_fix
// fprintf(stderr, "[M::%s-end]\n", __func__);
**/
cal_ec_multiple(b, asm_opt.thread_num, R_INF.total_reads);
// kt_for(asm_opt.thread_num, worker_hap_ec, b, R_INF.total_reads);///debug_for_fix
if (r_out) write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name);
ha_pt_destroy(ha_idx);
ha_idx = NULL;
if(het_cnt) {
print_het_cnt_log(het_cnt); free(het_cnt); het_cnt = NULL;
}
/**
// collect statistics
for (i = 0; i < asm_opt.thread_num; ++i) {
asm_opt.num_bases += b[i]->num_read_base;
asm_opt.num_corrected_bases += b[i]->num_correct_base;
asm_opt.num_recorrected_bases += b[i]->num_recorrect_base;
asm_opt.mem_buf += ha_ovec_mem(b[i], NULL);
ha_ovec_destroy(b[i]);
}
free(b);
**/
destroy_ec_ovec_buf_t(b);
exit(1);
if (asm_opt.required_read_name) prt_dbg_rs(R_INF_FLAG.fp_r0, &R_INF_FLAG, 0); // for debugging only
// save corrected reads to R_INF
CALLOC(e, asm_opt.thread_num);
for (i = 0; i < asm_opt.thread_num; ++i) {
init_UC_Read(&e[i].g_read);
e[i].first_round_read_size = e[i].second_round_read_size = 50000;
CALLOC(e[i].first_round_read, e[i].first_round_read_size);
CALLOC(e[i].second_round_read, e[i].second_round_read_size);
}
kt_for(asm_opt.thread_num, worker_ec_save, e, R_INF.total_reads);
for (i = 0; i < asm_opt.thread_num; ++i) {
destory_UC_Read(&e[i].g_read);
free(e[i].first_round_read);
free(e[i].second_round_read);
}
free(e);
if (asm_opt.required_read_name) prt_dbg_rs(R_INF_FLAG.fp_r1, &R_INF_FLAG, 1); // for debugging only
if (asm_opt.required_read_name) destory_Debug_reads(&R_INF_FLAG), exit(0); // for debugging only
///debug_print_pob_regions();
}
void ha_overlap_and_correct(int round)
{
int i, hom_cov, het_cov, r_out = 0;
@@ -992,12 +1078,15 @@ void ha_overlap_and_correct(int round)
het_cnt = NULL;
if(round == asm_opt.number_of_round-1 && asm_opt.is_dbg_het_cnt) CALLOC(het_cnt, R_INF.total_reads);
// fprintf(stderr, "[M::%s-start]\n", __func__);
// double tt0 = yak_realtime_0();
if (asm_opt.required_read_name)
kt_for(asm_opt.thread_num, worker_ovec_related_reads, b, R_INF.total_reads);
else
kt_for(asm_opt.thread_num, worker_ovec, b, R_INF.total_reads);///debug_for_fix
// fprintf(stderr, "[M::%s-end]\n", __func__);
// fprintf(stderr, "[M::%s::%.3f] ==> chaining\n", __func__, yak_realtime_0()-tt0);
// exit(1);
if (r_out) write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name);
ha_pt_destroy(ha_idx);
ha_idx = NULL;
@@ -1953,7 +2042,8 @@ int ha_assemble(void)
assert(asm_opt.number_of_round > 0);
for (r = ha_idx?asm_opt.number_of_round-1:0; r < asm_opt.number_of_round; ++r) {
ha_opt_reset_to_round(&asm_opt, r); // this update asm_opt.roundID and a few other fields
ha_overlap_and_correct(r);
// ha_overlap_and_correct(r);
ha_ec(r);
fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> corrected reads for round %d\n", __func__, yak_realtime(),
yak_cpu_usage(), yak_peakrss_in_gb(), r + 1);
fprintf(stderr, "[M::%s] # bases: %lld; # corrected bases: %lld; # recorrected bases: %lld\n", __func__,
Correct.cpp
+1166 -44
View File
File diff suppressed because it is too large Load Diff
Correct.h
+17
View File
@@ -1388,4 +1388,21 @@ int64_t get_rid_backward_cigar_err(rtrace_iter *it, ul_ov_t *aln, kv_rtrace_t *t
const ul_idx_t *uref, char* qstr, UC_Read *tu, overlap_region_alloc *ol, overlap_region *o,
bit_extz_t *exz, double e_rate, int64_t qs);
void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf);
void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, UC_Read* tu, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, int64_t bd, int64_t wl, int64_t ql, uint8_t occ_thres);
#define ovlp_id(x) ((x).tn)
#define ovlp_min_wid(x) ((x).ts)
#define ovlp_max_wid(x) ((x).te)
#define ovlp_cur_wid(x) ((x).qn)
#define ovlp_cur_xoff(x) ((x).qs)
#define ovlp_cur_yoff(x) ((x).ts)
#define ovlp_cur_ylen(x) ((x).te)
#define ovlp_cur_coff(x) ((x).qe)
#define ovlp_bd(x) ((x).sec)
#define UC_Read_resize(v, s) do {\
if ((v).size<(s)) {REALLOC((v).seq,(s));(v).size=(s);}\
} while (0)
#endif
Hash_Table.cpp
+246
View File
@@ -1986,6 +1986,252 @@ uint64_t lchain_qdp_mcopy(Candidates_list *cl, int64_t a_idx, int64_t a_n, int64
return cL;
}
void quick_ck_lchain(k_mer_hit* a, int64_t a_n, int64_t xl, int64_t yl, double chn_pen_gap, double chn_pen_skip, double bw_rate,
int64_t *p, int64_t *t, int32_t *f, int32_t *ii, int64_t *plus, int64_t *msc, int64_t *msc_i, int64_t *movl, int64_t *si, int64_t *ei)
{
if(a_n <= 0) return;
int64_t l, k, is_srt = 1, z; k_mer_hit *ai, *aj;
int64_t dq, dr, dd, dg, q_span, sc, csc, ddt;
int64_t plus0, msc0, msc_i0, movl0; double lin_pen, a_pen;
*plus = 0; *msc = *msc_i = INT32_MIN; *movl = INT32_MAX; *si = 0; *ei = a_n;
for (k = 1, l = 0; k <= a_n; k++) {
if(k == a_n || a[k].strand != a[l].strand) {
t[k-1] = 0; ii[k-1] = 0;
if(is_srt) {
plus0 = 0; msc0 = msc_i0 = INT32_MIN; movl0 = INT32_MAX; ddt = 0;
p[l] = -1; f[l] = a[l].cnt&(0xffu);
if(f[l] >= msc0) {msc0 = f[l]; msc_i0 = l;}///difference
if(f[l] < plus0) plus0 = f[l];
for (z = l + 1; z < k; z++) {
///roughly same to comput_sc_ch(&a[z], &a[z-1])
ai = &a[z]; aj = &a[z-1];
dq = (int64_t)(ai->self_offset) - (int64_t)(aj->self_offset);
if(dq <= 0) break;
dr = (int64_t)(ai->offset) - (int64_t)(aj->offset);
if(dr <= 0) break;
dd = dr > dq? dr - dq : dq - dr;//gap
if((dd > 16) && (dd > cal_bw(&(a[z]), &(a[z-1]), bw_rate, xl, yl))) break;
dg = dr < dq? dr : dq;//len
q_span = ai->cnt&(0xffu);
sc = q_span < dg? q_span : dg;
sc = normal_w(sc, ((int32_t)(ai->cnt>>8)));
if (dd || (dg > q_span && dg > 0)) {
lin_pen = (chn_pen_gap*(double)dd);
a_pen = ((double)(sc))*((((double)dd)/((double)dg))/bw_rate);
if(lin_pen > a_pen) lin_pen = a_pen;
lin_pen += (chn_pen_skip*(double)dg);
sc -= (int32_t)lin_pen;
}
sc += f[z-1]; csc = a[z].cnt&(0xffu); if(sc < csc) break;
p[z] = z - 1; f[z] = sc; ddt += dd;
if(f[z] >= msc0) {msc0 = f[z]; msc_i0 = z;}///difference
if(f[z] < plus0) plus0 = f[z];
}
if((z >= k) && (msc_i0 == (k - 1))) {
if((k - l >= 2) && (ddt > 16) && (ddt > cal_bw(&(a[k-1]), &(a[l]), bw_rate, xl, yl))) msc_i0 = INT32_MIN;
if(msc_i0 == (k - 1)) {
if(msc0 >= (*msc)) {
movl0 = get_chainLen(a[msc_i0].self_offset, a[msc_i0].self_offset, xl, a[msc_i0].offset, a[msc_i0].offset, yl);
if(msc0 > (*msc) || movl0 < (*movl)) {
*msc = msc0; *msc_i = msc_i0; *movl = movl0;
}
}
if(plus0 < (*plus)) *plus = plus0;
if((*ei) > k) {
(*si) = k;
} else {
(*ei) = l;
}
}
}
}
l = k; is_srt = 1;
} else {
if((a[k].self_offset <= a[k-1].self_offset) || (a[k].offset <= a[k-1].offset)) is_srt = 0;
t[k-1] = 0; ii[k-1] = 0;
}
}
}
uint64_t lchain_qdp_mcopy_fast(Candidates_list *cl, int64_t a_idx, int64_t a_n, int64_t des_idx,
Chain_Data* dp, overlap_region_alloc* res, int64_t max_skip, int64_t max_iter,
int64_t max_dis, double chn_pen_gap, double chn_pen_skip, double bw_rate,
uint32_t xid, int64_t xl, int64_t yl, int64_t quick_check, uint32_t apend_be,
int64_t gen_cigar, int64_t enable_mcopy, double mcopy_rate, int64_t mcopy_khit_cutoff,
int64_t khit_n)
{
if(a_n <= 0) return 0;
int64_t *p, *t, max_f, n_skip, st, max_j, end_j, sc, msc, msc_i, max_ii, ovl, movl, plus = 0, min_sc, ch_n, si, ei;
int32_t *f, max, tmp, *ii; int64_t i, k, j, cL = 0; k_mer_hit* a; k_mer_hit* des; k_mer_hit *swap; overlap_region *z;
resize_Chain_Data(dp, a_n, NULL); ch_n = 1; // int64_t bw; bw = ((xl < yl)?xl:yl); bw *= bw_rate;
t = dp->tmp; f = dp->score; p = dp->pre; ii = dp->occ;
a = cl->list + a_idx; des = cl->list + des_idx;
// if(a_n && a[0].readID == 0) {
// fprintf(stderr, "---[M::%s::utg%.6dl::%c]\n",
// __func__, (int32_t)a[0].readID+1, "+-"[a[0].strand]);
// }
if(quick_check) {
quick_ck_lchain(a, a_n, xl, yl, chn_pen_gap, chn_pen_skip, bw_rate, p, t, f, ii, &plus, &msc, &msc_i, &movl, &si, &ei);
} else {
msc = msc_i = INT32_MIN; movl = INT32_MAX; plus = 0; si = 0; ei = a_n;
memset(t, 0, (a_n*sizeof((*t))));
}
for (i = st = si, max_ii = -1; i < ei; ++i) {
max_f = a[i].cnt&(0xffu);
n_skip = 0; max_j = end_j = -1;
if ((i-st) > max_iter) st = i-max_iter;
while (a[i].strand != a[st].strand) ++st;
for (j = i - 1; j >= st; --j) {
sc = comput_sc_ch(&a[i], &a[j], bw_rate, chn_pen_gap, chn_pen_skip, xl, yl);
if (sc == INT32_MIN) continue;
sc += f[j];
if (sc > max_f) {
max_f = sc, max_j = j;
if (n_skip > 0) --n_skip;
} else if (t[j] == (int32_t)i) {
if (++n_skip > max_skip)
break;
}
if (p[j] >= 0) t[p[j]] = i;
}
end_j = j;
if ((max_ii<0) || (a[i].self_offset>a[max_ii].self_offset+max_dis) || (a[i].strand!=a[max_ii].strand)) {
max = INT32_MIN; max_ii = -1;
for (j=i-1; (j>=st) && (a[i].self_offset<=max_dis+a[j].self_offset)&&(a[i].strand==a[j].strand); --j) {
if (max < f[j]) {
max = f[j], max_ii = j;
}
}
}
if ((max_ii >= 0) && (max_ii < end_j) && (a[i].strand == a[max_ii].strand)) {///just have a try with a[i]<->a[max_ii]
tmp = comput_sc_ch(&a[i], &a[max_ii], bw_rate, chn_pen_gap, chn_pen_skip, xl, yl);
if (tmp != INT32_MIN && max_f < tmp + f[max_ii])
max_f = tmp + f[max_ii], max_j = max_ii;
}
f[i] = max_f; p[i] = max_j;
if ((max_ii < 0) || ((a[i].self_offset<=max_dis+a[max_ii].self_offset)&&(a[i].strand==a[max_ii].strand)&&(f[max_ii]<f[i]))) {
max_ii = i;
}
if(f[i] >= msc) {
ovl = get_chainLen(a[i].self_offset, a[i].self_offset, xl, a[i].offset, a[i].offset, yl);
if(f[i] > msc || ovl < movl) {
msc = f[i]; msc_i = i; movl = ovl;
}
}
if(f[i] < plus) plus = f[i];
ii[i] = 0;///for mcopy, not here
// if(a_n && a[0].readID == 0) {
// fprintf(stderr, "i::%ld[M::%s::utg%.6dl::%c] x::%u, y::%u, st::%ld, max_ii::%ld, f[i]::%d, p[i]::%ld, msc_i::%ld, msc::%ld, movl::%ld\n",
// i, __func__, (int32_t)a[i].readID+1, "+-"[a[i].strand],
// a[i].self_offset, a[i].offset, st, max_ii, f[i], p[i], msc_i, msc, movl);
// }
}
for (i = msc_i, cL = 0; i >= 0; i = p[i]) { ii[i] = 1; t[cL++] = i;}///label the best chain
if((movl < xl) && enable_mcopy/**(movl < yl)**/) {
if(cL >= mcopy_khit_cutoff) {///if there are too few k-mers, disable mcopy
msc -= plus; min_sc = msc*mcopy_rate/**0.2**/; ii[msc_i] = 0;
for (i = ch_n = 0; i < a_n; ++i) {///make all f[] positive
f[i] -= plus; if(i >= ch_n) t[i] = 0;
if((!(ii[i])) && (f[i] >= min_sc)) {
t[ch_n] = ((uint64_t)f[i])<<32; t[ch_n] += (i<<1); ch_n++;
}
}
if(ch_n > 1) {
int64_t n_v, n_v0, ni, n_u, n_u0 = res->length;
radix_sort_hc64i(t, t + ch_n);
for (k = ch_n-1, n_v = n_u = 0; k >= 0; --k) {
n_v0 = n_v;
for (i = ((uint32_t)t[k])>>1; i >= 0 && (t[i]&1) == 0; ) {
ii[n_v++] = i; t[i] |= 1; i = p[i];
}
if(n_v0 == n_v) continue;
sc = (i<0?(t[k]>>32):((t[k]>>32)-f[i]));
if(sc >= min_sc) {
kv_pushp_ol(overlap_region, (*res), &z);
push_ovlp_chain_qgen(z, xid, xl, yl, sc+plus, &(a[ii[n_v-1]]), &(a[ii[n_v0]]));
///mcopy_khit_cutoff <= 1: disable the mcopy_khit_cutoff filtering, for the realignment
if((mcopy_khit_cutoff <= 1) || ((z->x_pos_e+1-z->x_pos_s) <= (movl<<2))) {
z->align_length = n_v-n_v0; z->x_id = n_v0;
n_u++;
} else {///non-best is too long
res->length--; n_v = n_v0;
}
} else {
n_v = n_v0;
}
}
if(n_u > 1) ks_introsort_or_sss(n_u, res->list + n_u0);
res->length = n_u0 + filter_non_ovlp_xchains(res->list + n_u0, n_u, &n_v);
n_u = res->length;
if(n_u > n_u0 + 1) {
kv_resize_cl(k_mer_hit, (*cl), (n_v+cl->length));
a = cl->list + a_idx; des = cl->list + des_idx; swap = cl->list + cl->length;
for (k = n_u0, i = n_v0 = n_v = 0; k < n_u; k++) {
z = &(res->list[k]);
z->non_homopolymer_errors = des_idx + i;
n_v0 = z->x_id; ni = z->align_length;
for (j = 0; j < ni; j++, i++) {
///k0 + (ni - j - 1)
swap[i] = a[ii[n_v0 + (ni- j - 1)]];
swap[i].readID = k;
}
z->x_id = xid;
if(gen_cigar) gen_fake_cigar(&(z->f_cigar), z, apend_be, swap+i-ni, ni);
if(!khit_n) z->align_length = 0;
}
memcpy(des, swap, i*sizeof((*swap))); //assert(i == ch_n);
// fprintf(stderr, "[M::%s::msc->%ld] msc_k_hits::%u, cL::%ld, min_sc::%ld, best_sc::%ld, n_u0_sc::%d, mcopy_rate::%f, # chains::%ld\n",
// __func__, msc, res->list[n_u0].align_length, cL, min_sc, msc+plus, res->list[n_u0].shared_seed,
// mcopy_rate, n_u-n_u0);
} else if(n_u == n_u0 + 1) {
z = &(res->list[n_u0]); k = n_u0; i = 0;
z->non_homopolymer_errors = des_idx + i;
n_v0 = z->x_id; ni = z->align_length;
for (j = 0; j < ni; j++, i++) {
///k0 + (ni - j - 1)
des[i] = a[ii[n_v0 + (ni- j - 1)]];
des[i].readID = k;
}
z->x_id = xid;
if(gen_cigar) gen_fake_cigar(&(z->f_cigar), z, apend_be, des+i-ni, ni);
if(!khit_n) z->align_length = 0;
}
return i;
} else {
msc += plus; i = msc_i; cL = 0;
while (i >= 0) {t[cL++] = i; i = p[i];}
}
}
}
///a[] has been sorted by self_offset
// i = msc_i; cL = 0;
// while (i >= 0) {t[cL++] = i; i = p[i];}
kv_pushp_ol(overlap_region, (*res), &z);
push_ovlp_chain_qgen(z, xid, xl, yl, msc, &(a[t[cL-1]]), &(a[t[0]]));
for (i = 0; i < cL; i++) {des[i] = a[t[cL-i-1]]; des[i].readID = res->length-1;}
z->non_homopolymer_errors = des_idx;
if(gen_cigar) gen_fake_cigar(&(z->f_cigar), z, apend_be, des, cL);
if(khit_n) z->align_length = cL;
return cL;
}
#define rev_khit(an, xl, yl) do { \
Hash_Table.h
+8
View File
@@ -8,6 +8,7 @@
#define WINDOW 375
#define WINDOW_BOUNDARY 375
#define WINDOW_HC 775
///for one side, the first or last WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY bases should not be corrected
#define WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY 25
#define THRESHOLD 15
@@ -236,4 +237,11 @@ uint64_t lchain_qdp_mcopy(Candidates_list *cl, int64_t a_idx, int64_t a_n, int64
int64_t gen_cigar, int64_t enable_mcopy, double mcopy_rate, int64_t mcopy_khit_cutoff,
int64_t khit_n);
uint64_t lchain_qdp_mcopy_fast(Candidates_list *cl, int64_t a_idx, int64_t a_n, int64_t des_idx,
Chain_Data* dp, overlap_region_alloc* res, int64_t max_skip, int64_t max_iter,
int64_t max_dis, double chn_pen_gap, double chn_pen_skip, double bw_rate,
uint32_t xid, int64_t xl, int64_t yl, int64_t quick_check, uint32_t apend_be,
int64_t gen_cigar, int64_t enable_mcopy, double mcopy_rate, int64_t mcopy_khit_cutoff,
int64_t khit_n);
#endif
Levenshtein_distance.h
+20 -2
View File
@@ -548,6 +548,24 @@ inline int32_t pop_trace_back(asg16_v *res, int32_t i, uint16_t *c, uint32_t *le
return i;
}
inline void push_trace_bp(asg16_v *res, uint16_t c, uint16_t b, uint32_t len, uint32_t is_append)
{
uint16_t p;
if((is_append) && (res->n)) {
}
c <<= 14;
while (len >= (0x3fff)) {
p = (c + (0x3fff)); kv_push(uint16_t, *res, p); len -= (0x3fff);
}
if(len) {
p = (c + len); kv_push(uint16_t, *res, p);
}
}
///511 -> 16 64-bits
// #define MAX_E 511
// #define MAX_L 2500
@@ -601,7 +619,7 @@ inline uint32_t cigar_check(char *pstr, char *tstr, bit_extz_t *ez)
if(c == 0) {
for (k=0;(k<cl)&&(pstr[pi]==tstr[ti]);k++,pi++,ti++);
if(k!=cl) {
fprintf(stderr, "ERROR-d-0\n");
fprintf(stderr, "ERROR-d-0, pi::%d, ti::%d\n", pi, ti);
return 0;
}
} else {
@@ -609,7 +627,7 @@ inline uint32_t cigar_check(char *pstr, char *tstr, bit_extz_t *ez)
if(c == 1) {
for (k=0;(k<cl)&&(pstr[pi]!=tstr[ti]);k++,pi++,ti++);
if(k!=cl) {
fprintf(stderr, "ERROR-d-1\n");
fprintf(stderr, "ERROR-d-1, pi::%d, ti::%d\n", pi, ti);
return 0;
}
} else if(c == 2) {///more p
Makefile
+3 -2
View File
@@ -6,7 +6,7 @@ CPPFLAGS=
INCLUDES=
OBJS= CommandLines.o Process_Read.o Assembly.o Hash_Table.o \
POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o Purge_Dups.o \
htab.o hist.o sketch.o anchor.o extract.o sys.o hic.o rcut.o horder.o \
htab.o hist.o sketch.o anchor.o extract.o sys.o hic.o rcut.o horder.o ecovlp.o\
tovlp.o inter.o kalloc.o gfa_ut.o gchain_map.o
EXE= hifiasm
LIBS= -lz -lpthread -lm
@@ -40,7 +40,7 @@ depend:
Assembly.o: Assembly.h CommandLines.h Process_Read.h Overlaps.h kvec.h kdq.h
Assembly.o: Hash_Table.h htab.h POA.h Correct.h Levenshtein_distance.h
Assembly.o: kthread.h
Assembly.o: kthread.h ecovlp.h
CommandLines.o: CommandLines.h ketopt.h
Correct.o: Correct.h Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h
Correct.o: kdq.h CommandLines.h Levenshtein_distance.h POA.h Assembly.h
@@ -58,6 +58,7 @@ Process_Read.o: Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
Purge_Dups.o: ksort.h Purge_Dups.h kvec.h kdq.h Overlaps.h Hash_Table.h
Purge_Dups.o: htab.h Process_Read.h CommandLines.h Correct.h
Purge_Dups.o: Levenshtein_distance.h POA.h kthread.h
ecovlp.o: Hash_Table.h Process_Read.h Overlaps.h kthread.h
Trio.o: khashl.h kthread.h kseq.h Process_Read.h Overlaps.h kvec.h kdq.h
Trio.o: CommandLines.h htab.h
anchor.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
anchor.cpp
+261
View File
@@ -979,6 +979,103 @@ uint32_t *low_occ)
cl->length = ab->n_a;
}
void minimizers_qgen0(ha_abuf_t *ab, char* rs, int64_t rl, uint64_t mz_w, uint64_t mz_k, Candidates_list *cl, kvec_t_u8_warp* k_flag,
void *ha_flt_tab, ha_pt_t *ha_idx, All_reads* rdb, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ)
{
// fprintf(stderr, "+[M::%s]\n", __func__);
uint64_t i, k, l, max_cnt = UINT32_MAX, min_cnt = 0; int n, j; ha_mz1_t *z; seed1_t *s;
if(high_occ) {
max_cnt = (*high_occ);
if(max_cnt < 2) max_cnt = 2;
}
if(low_occ) {
min_cnt = (*low_occ);
if(min_cnt < 2) min_cnt = 2;
}
clear_Candidates_list(cl); ab->mz.n = 0, ab->n_a = 0;
// get the list of anchors
mz1_ha_sketch(rs, rl, mz_w, mz_k, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0, NULL);
// minimizer of queried read
if (ab->mz.m > ab->old_mz_m) {
ab->old_mz_m = ab->mz.m;
REALLOC(ab->seed, ab->old_mz_m);
}
for (i = 0, ab->n_a = 0; i < ab->mz.n; ++i) {
ab->seed[i].a = ha_pt_get(ha_idx, ab->mz.a[i].x, &n);
ab->seed[i].n = n;
ab->n_a += n;
}
if (ab->n_a > ab->m_a) {
ab->m_a = ab->n_a;
REALLOC(ab->a, ab->m_a);
}
for (i = 0, k = 0; i < ab->mz.n; ++i) {
///z is one of the minimizer
z = &ab->mz.a[i]; s = &ab->seed[i];
for (j = 0; j < s->n; ++j) {
const ha_idxpos_t *y = &s->a[j];
anchor1_t *an = &ab->a[k++];
uint8_t rev = z->rev == y->rev? 0 : 1;
an->other_off = rev?((uint32_t)-1)-1-(y->pos+1-y->span):y->pos;
an->self_off = z->pos;
///an->cnt: cnt<<8|span
an->cnt = s->n; if(an->cnt > ((uint32_t)(0xffffffu))) an->cnt = 0xffffffu;
an->cnt <<= 8; an->cnt |= ((z->span <= ((uint32_t)(0xffu)))?z->span:((uint32_t)(0xffu)));
an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->self_off;
}
}
// copy over to _cl_
if (ab->m_a >= (uint64_t)cl->size) {
cl->size = ab->m_a;
REALLOC(cl->list, cl->size);
}
k_mer_hit *p; uint64_t tid = (uint64_t)-1, tl = (uint64_t)-1;
radix_sort_ha_an1(ab->a, ab->a + ab->n_a);
for (k = 1, l = 0; k <= ab->n_a; ++k) {
if (k == ab->n_a || ab->a[k].srt != ab->a[l].srt) {
if (k-l>1) radix_sort_ha_an3(ab->a+l, ab->a+k);
if((ab->a[l].srt>>33)!=tid) {
tid = ab->a[l].srt>>33;
tl = Get_READ_LENGTH((*rdb), tid);
// tl = rdb?Get_READ_LENGTH((*rdb), tid):udb->ug->u.a[tid].len;
}
for (i = l; i < k; i++) {
p = &cl->list[i];
p->readID = ab->a[i].srt>>33;
p->strand = (ab->a[i].srt>>32)&1;
if(!(p->strand)) {
p->offset = ab->a[i].other_off;
} else {
p->offset = ((uint32_t)-1)-ab->a[i].other_off;
p->offset = tl-p->offset;
}
p->self_offset = ab->a[i].self_off;
if(((ab->a[i].cnt>>8) < max_cnt) && ((ab->a[i].cnt>>8) > min_cnt)){
p->cnt = 1;
} else if((ab->a[i].cnt>>8) <= min_cnt) {
p->cnt = 2;
} else{
p->cnt = 1 + (((ab->a[i].cnt>>8) + (max_cnt<<1) - 1)/(max_cnt<<1));
p->cnt = pow(p->cnt, 1.1);
}
if(p->cnt > ((uint32_t)(0xffffffu))) p->cnt = 0xffffffu;
p->cnt <<= 8; p->cnt |= (((uint32_t)(0xffu))&(ab->a[i].cnt));
}
l = k;
}
}
cl->length = ab->n_a;
}
void gen_pair_chain(ha_abufl_t *ab, uint64_t rid, st_mt_t *tid, uint64_t tid_n, ha_mzl_t *in, uint64_t in_n, ha_mzl_t *idx, int64_t idx_n, uint64_t mzl_cutoff)
{
if(!tid_n) return;
@@ -1615,6 +1712,155 @@ void lchain_qgen_mcopy(Candidates_list* cl, overlap_region_alloc* ol, uint32_t r
for (i = 0; i < ol->length; ++i) ol->list[i].align_length = 0;
}
void lchain_qgen_mcopy_fast(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, uint64_t rl, All_reads* rdb,
uint32_t apend_be, uint64_t max_n_chain, int64_t max_skip, int64_t max_iter,
int64_t max_dis, double chn_pen_gap, double chn_pen_skip, double bw_rate, int64_t quick_check,
uint32_t gen_off, int64_t enable_mcopy, double mcopy_rate, uint32_t chain_cutoff, uint32_t mcopy_khit_cut, st_mt_t *sp)
{
// fprintf(stderr, "+[M::%s]\n", __func__);
uint64_t i, k, l, m, cn = cl->length, yid, ol0, lch; overlap_region *r, t; ///srt = 0
clear_overlap_region_alloc(ol);
for (l = 0, k = 1, m = 0, lch = 0; k <= cn; k++) {
if((k == cn) || (cl->list[k].readID != cl->list[l].readID)) {
if(cl->list[l].readID != rid) {
yid = cl->list[l].readID; ol0 = ol->length;
m += lchain_qdp_mcopy_fast(cl, l, k-l, m, &(cl->chainDP), ol, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_rate,
rid, rl, Get_READ_LENGTH((*rdb), yid), quick_check, apend_be, gen_off, enable_mcopy, mcopy_rate, mcopy_khit_cut, 1);
if((chain_cutoff >= 2) && (!lch)) {
for (i = ol0; (i<ol->length) && (!lch); i++) {
if(ol->list[i].align_length < chain_cutoff) lch = 1;
}
}
}
l = k;
}
}
cl->length = m;
// for (k = 0; k < ol->length; k++) {
// fprintf(stderr, "---[M::%s::utg%.6dl] q[%d, %d), t[%d, %d), khit_off::%u\n", __func__,
// (int32_t)ol->list[k].y_id+1, ol->list[k].x_pos_s, ol->list[k].x_pos_e+1,
// ol->list[k].y_pos_s, ol->list[k].y_pos_e+1, ol->list[k].non_homopolymer_errors);
// }
k = ol->length;
if (ol->length > max_n_chain) {
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(ol->length, ol->list);
for (i = 0; i < ol->length; ++i) {
r = &(ol->list[i]);
w = ha_ov_type(r, rl);
++n[w];
if (((uint64_t)n[w]) == max_n_chain) s[w] = r->shared_seed;
}
if (s[0] > 0 || s[1] > 0 || s[2] > 0 || s[3] > 0) {
// n[0] = n[1] = n[2] = n[3] = 0;
for (i = 0, k = 0, lch = 0; i < ol->length; ++i) {
r = &(ol->list[i]);
w = ha_ov_type(r, rl);
// ++n[w];
// if (((int)n[w] <= max_n_chain) || (r->shared_seed >= s[w] && s[w] >= (asm_opt.k_mer_length<<1))) {
if (r->shared_seed >= s[w]) {
if (k != i) {
t = ol->list[k];
ol->list[k] = ol->list[i];
ol->list[i] = t;
}
if(ol->list[k].align_length < chain_cutoff) lch = 1;
++k;
}
}
ol->length = k;
}
}
ks_introsort_or_xs(ol->length, ol->list);
if(lch) {
//@brief r485
uint64_t zs, ze, rs, re, ob, os, oe, ocn, pp, kn, ms, me; int64_t osc;
for (i = l = 0, cn = cl->length; i < ol->length; ++i) {
if(ol->list[i].align_length < chain_cutoff) {
zs = ol->list[i].x_pos_s; ze = ol->list[i].x_pos_e + 1;
ob = (ze - zs)*OFL; if(ob < 16) ob = 16;
osc = ol->list[i].shared_seed*CH_SC;
ocn = ol->list[i].align_length<<CH_OCC;
for (k = 0; (k < ol->length) && (ze > ol->list[k].x_pos_s); k++) {
if(ol->list[k].align_length < chain_cutoff) continue;
if(ol->list[k].align_length < ocn) continue;
if(ol->list[k].shared_seed < osc) continue;
rs = ol->list[k].x_pos_s; re = ol->list[k].x_pos_e + 1;
os = ((rs>=zs)?rs:zs); oe = ((re<=ze)?re:ze);
if((oe > os) && (oe - os) >= ob) {
m = ol->list[k].non_homopolymer_errors;
pp = cl->list[m].readID; kn = 0;
for (; (m < cn) && (cl->list[m].readID == pp) && (kn < ocn); m++) {
me = cl->list[m].self_offset; ms = me - (cl->list[m].cnt&(0xffu));
if((ms >= os) && (me <= oe)) kn++;
}
if(kn >= ocn) break;
}
}
if((k < ol->length) && (ze > ol->list[k].x_pos_s)) continue;
}
if (l != i) {
t = ol->list[l];
ol->list[l] = ol->list[i];
ol->list[i] = t;
}
l++;
}
// fprintf(stderr, "+[M::%s] rid::%u, ol->length0::%lu, ol->length1::%lu\n", __func__, rid, ol->length, l);
ol->length = l;
/**
//@brief r484
for (i = sp->n = 0; i < ol->length; ++i) {
if(ol->list[i].align_length < chain_cutoff) continue;
os = ol->list[i].x_pos_s; oe = ol->list[i].x_pos_e + 1;
if((sp->n) && (((uint32_t)sp->a[sp->n-1]) >= os)) {
if(oe > ((uint32_t)sp->a[sp->n-1])) {
oe = oe - ((uint32_t)sp->a[sp->n-1]);
sp->a[sp->n-1] += oe;
}
} else {
os = (os<<32)|oe; kv_push(uint64_t, *sp, os);
}
}
for (i = k = 0; i < ol->length; ++i) {
if(ol->list[i].align_length < chain_cutoff) {///ol has been sorted by x_pos_s
r = &(ol->list[i]); rs = r->x_pos_s; re = r->x_pos_e + 1;
rl = re - rs; ovl = 0;
for (m = 0; (m < sp->n) && (re > (sp->a[m]>>32)); m++) {
os = ((rs>=(sp->a[m]>>32))? rs:(sp->a[m]>>32));
oe = ((re<=((uint32_t)sp->a[m]))? re:((uint32_t)sp->a[m]));
if(oe > os) {
ovl += (oe - os); if(ovl >= (rl*0.95)) break;
}
}
if(ovl >= (rl*0.95)) continue;
}
if (k != i) {
t = ol->list[k];
ol->list[k] = ol->list[i];
ol->list[i] = t;
}
ol->list[k++].align_length = 0;
}
// fprintf(stderr, "+[M::%s] ol->length0::%lu, ol->length1::%lu\n", __func__, ol->length, k);
ol->length = k;
**/
}
/**else {
for (i = 0; i < ol->length; ++i) ol->list[i].align_length = 0;
}
**/
for (i = 0; i < ol->length; ++i) ol->list[i].align_length = 0;
}
inline uint64_t special_lchain(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, uint64_t rl, All_reads* rdb,
const ul_idx_t *udb, uint32_t apend_be, int64_t max_skip, int64_t max_iter, int64_t max_dis, double chn_pen_gap,
double chn_pen_skip, double bw_rate, int64_t quick_check, uint32_t gen_off, double mcopy_rate, uint32_t mcopy_khit_cut,
@@ -1802,6 +2048,21 @@ void ul_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t
lchain_qgen_mcopy(cl, overlap_list, rid, rl, NULL, uref, apend_be, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off, mcopy_rate, chain_cutoff, mcopy_khit_cut, sp);
}
void h_ec_lchain(ha_abuf_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, All_reads *rref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres,
int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate, uint32_t gen_off, int64_t enable_mcopy, double mcopy_rate, uint32_t chain_cutoff, uint32_t mcopy_khit_cut)
{
extern void *ha_flt_tab;
extern ha_pt_t *ha_idx;
int64_t max_skip, max_iter, max_dis, quick_check; double chn_pen_gap, chn_pen_skip;
set_lchain_dp_op(is_accurate, mz_k, &max_skip, &max_iter, &max_dis, &chn_pen_gap, &chn_pen_skip, &quick_check);
// minimizers_gen(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, dbg_ct, sp, high_occ, low_occ);
minimizers_qgen0(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, rref, dbg_ct, sp, high_occ, low_occ);
// lchain_gen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off);
// lchain_qgen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off);
///no need to sort here, overlap_list has been sorted at lchain_gen
lchain_qgen_mcopy_fast(cl, overlap_list, rid, rl, rref, apend_be, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off, enable_mcopy, mcopy_rate, chain_cutoff, mcopy_khit_cut, sp);
}
int64_t ug_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, const ul_idx_t *uref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, double bw_thres_sec,
int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate,
uint32_t gen_off, double mcopy_rate, uint32_t mcopy_khit_cut, uint32_t is_hpc, ha_mzl_t *res, uint64_t res_n, ha_mzl_t *idx, uint64_t idx_n, uint64_t mzl_cutoff, uint64_t chain_cutoff, kv_u_trans_t *kov)
ecovlp.cpp
+860
View File
@@ -0,0 +1,860 @@
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include "Correct.h"
#include "Process_Read.h"
#include "ecovlp.h"
#include "kthread.h"
#include "htab.h"
#define HA_KMER_GOOD_RATIO 0.333
#define E_KHIT 31
#define generic_key(x) (x)
KRADIX_SORT_INIT(ec16, uint16_t, generic_key, 2)
KRADIX_SORT_INIT(ec32, uint32_t, generic_key, 4)
KRADIX_SORT_INIT(ec64, uint64_t, generic_key, 8)
void h_ec_lchain(ha_abuf_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, All_reads *rref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres,
int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate, uint32_t gen_off, int64_t enable_mcopy, double mcopy_rate, uint32_t chain_cutoff, uint32_t mcopy_khit_cut);
ec_ovec_buf_t* gen_ec_ovec_buf_t(uint32_t n, uint32_t is_final, uint32_t save_ov)
{
uint32_t k; ec_ovec_buf_t0 *z = NULL;
ec_ovec_buf_t *p = NULL; CALLOC(p, 1);
p->n = n; CALLOC(p->a, p->n);
for (k = 0; k < p->n; k++) {
z = &(p->a[k]);
z->is_final = !!is_final; z->save_ov = !!save_ov;
init_UC_Read(&z->self_read);
init_UC_Read(&z->ovlp_read);
init_Candidates_list(&z->clist);
init_overlap_region_alloc(&z->olist);
init_fake_cigar(&(z->tmp.f_cigar));
memset(&(z->tmp.w_list), 0, sizeof(z->tmp.w_list));
CALLOC(z->tmp.w_list.a, 1); z->tmp.w_list.n = z->tmp.w_list.m = 1;
// kv_init(z->b_buf.a);
kv_init(z->r_buf.a);
kv_init(z->k_flag.a);
kv_init(z->sp);
kv_init(z->pidx);
kv_init(z->v64);
kv_init(z->v32);
kv_init(z->v16);
init_bit_extz_t(&(z->exz), 31);
z->ab = ha_abuf_init();
if (!z->is_final) {
init_Cigar_record(&z->cigar);
// init_Graph(&b->POA_Graph);
// init_Graph(&b->DAGCon);
init_Correct_dumy(&z->correct);
InitHaplotypeEvdience(&z->hap);
}
}
return p;
}
void destroy_cns_gfa(cns_gfa *p)
{
size_t k;
for (k = 0; k < p->m; k++) {
kv_destroy(p->a[k].in);
kv_destroy(p->a[k].ou);
}
free(p->a);
}
void destroy_ec_ovec_buf_t(ec_ovec_buf_t *p)
{
uint32_t k; ec_ovec_buf_t0 *z = NULL;
for (k = 0; k < p->n; k++) {
z = &(p->a[k]);
destory_UC_Read(&z->self_read);
destory_UC_Read(&z->ovlp_read);
destory_Candidates_list(&z->clist);
destory_overlap_region_alloc(&z->olist);
destory_fake_cigar(&(z->tmp.f_cigar));
free(z->tmp.w_list.a); free(z->tmp.w_list.c.a);
kv_destroy(z->r_buf.a);
kv_destroy(z->k_flag.a);
kv_destroy(z->sp);
kv_destroy(z->pidx);
kv_destroy(z->v64);
kv_destroy(z->v32);
kv_destroy(z->v16);
destroy_bit_extz_t(&(z->exz));
ha_abuf_destroy(z->ab);
if (!z->is_final) {
destory_Cigar_record(&z->cigar);
// destory_Graph(&b->POA_Graph);
// destory_Graph(&b->DAGCon);
destory_Correct_dumy(&z->correct);
destoryHaplotypeEvdience(&z->hap);
}
destroy_cns_gfa(&(z->cns));
asm_opt.num_bases += z->num_read_base;
asm_opt.num_corrected_bases += z->num_correct_base;
asm_opt.num_recorrected_bases += z->num_recorrect_base;
// asm_opt.mem_buf += ha_ovec_mem(b[i], NULL);
}
free(p->a); free(p);
fprintf(stderr, "[M::%s-chains] #->%lld\n", __func__, asm_opt.num_bases);
fprintf(stderr, "[M::%s-passed-chains-0] #->%lld\n", __func__, asm_opt.num_corrected_bases);
fprintf(stderr, "[M::%s-cis-chains-1] #->%lld\n", __func__, asm_opt.num_recorrected_bases);
}
void prt_chain(overlap_region_alloc *o)
{
uint64_t k;
for (k = 0; k < o->length; k++) {
fprintf(stderr, "[M::%s]\t#%u\tlen::%lu\t%u\t%u\t%c\t#%u\tlen::%lu\t%u\t%u\tsc::%d\taln::%u\terr::%u\n", __func__, o->list[k].x_id, Get_READ_LENGTH(R_INF, o->list[k].x_id), o->list[k].x_pos_s, o->list[k].x_pos_e+1, "+-"[o->list[k].y_pos_strand],
o->list[k].y_id, Get_READ_LENGTH(R_INF, o->list[k].y_id), o->list[k].y_pos_s, o->list[k].y_pos_e+1, o->list[k].shared_seed, o->list[k].align_length, o->list[k].non_homopolymer_errors);
}
}
overlap_region *fetch_aux_ovlp(overlap_region_alloc* ol) /// exactly same to gen_aux_ovlp
{
if (ol->length + 1 > ol->size) {
uint64_t sl = ol->size;
ol->size = ol->length + 1;
kroundup64(ol->size);
REALLOC(ol->list, ol->size);
/// need to set new space to be 0
memset(ol->list + sl, 0, sizeof(overlap_region)*(ol->size - sl));
}
return &(ol->list[ol->length+1]);
}
typedef struct {
ul_ov_t *c_idx;
asg64_v *idx;
int64_t i, i0, srt_n, rr;
uint64_t mms, mme;
} cc_idx_t;
///[s, e)
int64_t extract_sub_cigar_mm(overlap_region *z, int64_t s, int64_t e, ul_ov_t *p, uint64_t *ct)
{
int64_t wk = ovlp_cur_wid(*p), xk = ovlp_cur_xoff(*p), yk = ovlp_cur_yoff(*p), ck = ovlp_cur_coff(*p), os, oe, t;
bit_extz_t ez; int64_t bd = ovlp_bd(*p), s0, e0;
s0 = ((int64_t)(z->w_list.a[wk].x_start)) + bd;
e0 = ((int64_t)(z->w_list.a[wk].x_end)) + 1 - bd;
if(s < s0) s = s0; if(e > e0) e = e0;///exclude boundary
if(s >= e) return -1;
os = MAX(s, s0); oe = MIN(e, e0);
if(oe <= os) return -1;
set_bit_extz_t(ez, (*z), wk);
if(!ez.cigar.n) return -1;
int64_t cn = ez.cigar.n, op; int64_t ws, we, ovlp;
if((ck < 0) || (ck > cn)) {//(*ck) == cn is allowed
ck = 0; xk = ez.ts; yk = ez.ps;
}
while (ck > 0 && xk >= s) {///x -> t; y -> p; first insertion and then match/mismatch
--ck;
op = ez.cigar.a[ck]>>14;
if(op!=2) xk -= (ez.cigar.a[ck]&(0x3fff));
if(op!=3) yk -= (ez.cigar.a[ck]&(0x3fff));
}
//some cigar will span s or e
while (ck < cn && xk < e) {//[s, e)
ws = xk;
op = ez.cigar.a[ck]>>14;
///op == 3: -> x; op == 2: -> y;
if(op!=2) xk += (ez.cigar.a[ck]&(0x3fff));
if(op!=3) yk += (ez.cigar.a[ck]&(0x3fff));
ck++; we = xk;
os = MAX(s, ws); oe = MIN(e, we);
ovlp = ((oe>os)? (oe-os):0);
if(op != 2) {
if(!ovlp) continue;
} else {///ws == we
if(ws < s || ws >= e) continue;
}
if(op == 0) {
for (t = os + 1; t < oe; t++) {
ct[(t-s)<<1]++; ct[(t-s)<<1] += ((uint64_t)(0x100000000));
ct[((t-s)<<1)+1]++; ct[((t-s)<<1)+1] += ((uint64_t)(0x100000000));
}
t = os;
if(t < oe) {
ct[(t-s)<<1]++; ct[(t-s)<<1] += ((uint64_t)(0x100000000));
if(os > ws) {
ct[((t-s)<<1)+1]++; ct[((t-s)<<1)+1] += ((uint64_t)(0x100000000));
}
}
} else if(op!=2) {
for (t = os + 1; t < oe; t++) {
ct[(t-s)<<1]++;
ct[((t-s)<<1)+1]++;
}
t = os;
if(t < oe) {
ct[(t-s)<<1]++;
if(os > ws) {
ct[((t-s)<<1)+1]++;
}
}
} else {
ct[((ws-s)<<1)+1]++; ///ct[((ws-s)<<1)+1] += ((uint64_t)(0x100000000));
}
}
ovlp_cur_xoff(*p) = xk; ovlp_cur_yoff(*p) = yk; ovlp_cur_coff(*p) = ck; ovlp_cur_ylen(*p) = 0;
return 1;
}
#define simp_vote_len 6
///[s, e)
uint32_t extract_sub_cigar_ii(overlap_region *z, All_reads *rref, int64_t s, int64_t e, UC_Read* tu, ul_ov_t *p)
{
int64_t wk = ovlp_cur_wid(*p), xk = ovlp_cur_xoff(*p), yk = ovlp_cur_yoff(*p), ck = ovlp_cur_coff(*p), os, oe, ol;
bit_extz_t ez; int64_t bd = ovlp_bd(*p), s0, e0, ii[2], it[2]; uint32_t res = (uint32_t)-1;
s0 = ((int64_t)(z->w_list.a[wk].x_start)) + bd;
e0 = ((int64_t)(z->w_list.a[wk].x_end)) + 1 - bd;
if(s < s0) s = s0; if(e > e0) e = e0;///exclude boundary
if(s > e) return -1;///it is possible s == e
os = MAX(s, s0); oe = MIN(e, e0);
if(oe < os) return -1;///it is possible os == oe
set_bit_extz_t(ez, (*z), wk);
if(!ez.cigar.n) return -1;
int64_t cn = ez.cigar.n; uint16_t op; int64_t ws, we, wts, wte, ovlp, cc = 0, cci;
if((ck < 0) || (ck > cn)) {//(*ck) == cn is allowed
ck = 0; xk = ez.ts; yk = ez.ps;
}
while (ck > 0 && xk >= s) {///x -> t; y -> p; first insertion and then match/mismatch
--ck;
op = ez.cigar.a[ck]>>14;
if(op!=2) xk -= (ez.cigar.a[ck]&(0x3fff));
if(op!=3) yk -= (ez.cigar.a[ck]&(0x3fff));
}
char cm[4]; cm[0] = 'M'; cm[1] = 'S'; cm[2] = 'I'; cm[3] = 'D';
//some cigar will span s or e
ii[0] = ii[1] = it[0] = it[1] = -1; res = cc = 0;
while (ck < cn && xk < e) {//[s, e)
ws = xk; wts = yk;
op = ez.cigar.a[ck]>>14; ol = (ez.cigar.a[ck]&(0x3fff));
///op == 3: -> x; op == 2: -> y;
if(op!=2) xk += ol;
if(op!=3) yk += ol;
ck++; we = xk; wte = yk;
// if(s == 10480) {
// fprintf(stderr, "[%ld, %ld)\t%c\n", ws, we, cm[op]);
// }
os = MAX(s, ws); oe = MIN(e, we);
ovlp = ((oe>os)? (oe-os):0);
if(op != 2) {
if(!ovlp) continue;
} else {///ws == we
if(ws < s || ws >= e) continue;
}
if(ii[0] == -1) {
ii[0] = os;
if(op < 2) {
it[0] = os - ws + wts;
} else {///op == 2: more y; p == 3: more x
it[0] = wts;
}
}
ii[1] = oe;
if(op < 2) {
it[1] = oe - ws + wts;
} else {///op == 2: more y; p == 3: more x
it[1] = wte;
}
if(op != 2) ol = oe-os;
cc += ol;
fprintf(stderr, "%ld%c", ol, cm[op]);
if(cc <= simp_vote_len) {
for (cci = 0; cci < ol; cci++) {
res <<= 2; res |= op;
}
}
}
while (ck < cn && xk <= e) {//[s, e)
ws = xk; wts = yk;
op = ez.cigar.a[ck]>>14; ol = (ez.cigar.a[ck]&(0x3fff));
if(op != 2) break;
yk += (ez.cigar.a[ck]&(0x3fff));
ck++; we = xk; wte = yk;
if(ws >= s && ws <= e) {
if(ii[0] == -1) {
ii[0] = os; it[0] = wts;
}
ii[1] = oe; it[1] = wte;
cc += ol;
fprintf(stderr, "%ld%c", ol, cm[op]);
if(cc <= simp_vote_len) {
for (cci = 0; cci < ol; cci++) {
res <<= 2; res |= op;
}
}
}
}
fprintf(stderr, "\tx::[%ld, %ld)\ty::[%ld, %ld)\tcc::%ld\n", ii[0], ii[1], it[0], it[1], cc);
if((cc <= simp_vote_len)
&& (ii[1] >= ii[0]) && (ii[1] - ii[0] <= simp_vote_len)
&& (it[1] >= it[0]) && (it[1] - it[0] <= simp_vote_len)) {
// ii[0] = ii[0] - s; ii[1] = e - ii[1];
if((ii[0] == s) && (ii[1] == e)) {
op = cc; op <<= 12; res |= op;
char *ystr = NULL; res <<= 16; cc = it[1] - it[0]; op = 0;
if(cc > 0) {
UC_Read_resize(*tu, (it[1] - it[0])); ystr = tu->seq;
recover_UC_Read_sub_region(ystr, it[0], (it[1] - it[0]), z->y_pos_strand, rref, z->y_id);
for (cci = 0; cci < cc; cci++) {
op <<= 2; op |= seq_nt6_table[(uint32_t)(ystr[cci])];
}
}
res |= op;
op = it[1] - it[0]; op <<= 12; res |= op;
} else {
res = (uint32_t)-1;
}
} else {
res = (uint32_t)-1;
}
ovlp_cur_xoff(*p) = xk; ovlp_cur_yoff(*p) = yk; ovlp_cur_coff(*p) = ck; ovlp_cur_ylen(*p) = 0;
return res;
}
uint64_t iter_cc_idx_t(overlap_region* ol, cc_idx_t *z, int64_t s, int64_t e, uint64_t is_reduce, uint64_t is_insert, uint64_t **ra)
{
int64_t rm_n, q[2], os, oe; ul_ov_t *cp; uint64_t m; *ra = NULL;
// if(s == 15816 && e == 15819) {
// fprintf(stderr, "[M::%s] is_reduce::%lu\n", __func__, is_reduce);
// }
if(is_reduce) {
for (m = rm_n = z->srt_n; m < z->idx->n; m++) {
cp = &(z->c_idx[z->idx->a[m]]);
// if(s == 15816 && e == 15819) {
// fprintf(stderr, "-0-[M::%s] ii::%lu, ii0::%ld\n", __func__, z->idx->a[m], z->i0);
// }
q[0] = ol[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_start+ovlp_bd(*cp);
q[1] = ol[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_end+1-ovlp_bd(*cp);
os = MAX(q[0], s); oe = MIN(q[1], e);
if((oe > os) || ((is_insert) && (s == e) && (s >= q[0]) && (s <= q[1]))) {
z->idx->a[rm_n++] = z->idx->a[m];
}
}
z->idx->n = rm_n;
}
for (; z->i < z->srt_n; ++z->i) {
cp = &(z->c_idx[(uint32_t)z->idx->a[z->i]]);
// if(s == 15816 && e == 15819) {
// fprintf(stderr, "-1-[M::%s] ii::%u, ii0::%ld\n", __func__, (uint32_t)z->idx->a[z->i], z->i0);
// }
q[0] = ol[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_start+ovlp_bd(*cp);
q[1] = ol[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_end+1-ovlp_bd(*cp);
if(q[0] > e) break;
if((!is_insert) && (q[0] >= e)) break;
os = MAX(q[0], s); oe = MIN(q[1], e);
if((oe > os) || ((is_insert) && (s == e) && (s >= q[0]) && (s <= q[1]))) {
kv_push(uint64_t, *(z->idx), ((uint32_t)z->idx->a[z->i]));
}
}
(*ra) = z->idx->a + z->srt_n;
return z->idx->n - z->srt_n;
}
void debug_inter0(overlap_region* ol, ul_ov_t *c_idx, uint64_t *idx, int64_t idx_n, uint64_t *res, int64_t res_n, int64_t s, int64_t e, uint64_t is_insert, uint64_t is_hard_check, const char *cmd)
{
ul_ov_t *cp; int64_t q[2], a_n = 0, i, k = 0, os, oe;
for (i = 0; i < idx_n; i++) {
cp = &(c_idx[(uint32_t)idx[i]]);
q[0] = ol[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_start+ovlp_bd(*cp);
q[1] = ol[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_end+1-ovlp_bd(*cp);
// fprintf(stderr, "%s[M::%s] tid::%u\t%.*s\twid::%u\tq::[%u, %u)\terr::%d\toerr::%u\n", cmd, __func__, ol[ovlp_id(*cp)].y_id, (int)Get_NAME_LENGTH(R_INF, ol[ovlp_id(*cp)].y_id), Get_NAME(R_INF, ol[ovlp_id(*cp)].y_id),
// ovlp_cur_wid(*cp), ol[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_start, ol[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_end+1, ol[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].error, ol[ovlp_id(*cp)].non_homopolymer_errors);
os = MAX(q[0], s); oe = MIN(q[1], e);
if((oe > os) || ((is_insert) && (s == e) && (s >= q[0]) && (s <= q[1]))) {
a_n++;
// if(!(((uint32_t)idx[i]) == res[k])) {
// fprintf(stderr, "[M::%s] a_n::%ld\tres_n::%ld\ts::%ld\te::%ld\ti::%ld\tk::%ld\n", __func__, a_n, res_n, s, e, i, k);
// }
if(is_hard_check) {
assert(((uint32_t)idx[i]) == res[k++]);
} else {
for (; (k < res_n) && (((uint32_t)idx[i]) != res[k]); k++);
assert(k < res_n);
}
}
}
// if(a_n != res_n) {
// fprintf(stderr, "[M::%s] a_n::%ld\tres_n::%ld\ts::%ld\te::%ld\tidx_n::%ld\n", __func__, a_n, res_n, s, e, idx_n);
// }
if(is_hard_check) {
assert(a_n == res_n);
} else {
assert(a_n <= res_n);
}
}
void prt_cigar0(uint64_t in, int64_t len)
{
int64_t k; uint64_t mp;
char cm[4]; cm[0] = 'M'; cm[1] = 'S'; cm[2] = 'I'; cm[3] = 'D';
for (k = 0; k < len; k++) {
mp = len - 1 - k; mp <<= 1;
fprintf(stderr, "%c", cm[(in >> mp)&3]);
}
fprintf(stderr, "\n");
}
void prt_bp0(uint64_t in, int64_t len)
{
int64_t k; uint64_t mp;
char cm[4]; cm[0] = 'A'; cm[1] = 'C'; cm[2] = 'G'; cm[3] = 'T';
for (k = 0; k < len; k++) {
mp = len - 1 - k; mp <<= 1;
fprintf(stderr, "%c", cm[(in >> mp)&3]);
}
fprintf(stderr, "\n");
}
uint64_t cns_gen0(overlap_region* ol, All_reads *rref, uint64_t s, uint64_t e, UC_Read* tu, cc_idx_t *idx, uint64_t occ_tot, double occ_max, asg32_v* b32, uint32_t *rc)
{
if(e > s + simp_vote_len) return 0;///too long
uint64_t *id_a = NULL, id_n, an = 0, oc[2]; b32->n = 0; uint32_t m, *a = NULL;
id_n = iter_cc_idx_t(ol, idx, s, e, idx->rr, ((s==e)?1:0), &id_a);
// debug_inter0(ol, idx->c_idx, idx->idx->a + idx->i0, idx->srt_n - idx->i0, id_a, id_n, s, e, ((s==e)?1:0), 0, "-1-");
uint64_t k, l, q[2], os, oe; ul_ov_t *p; overlap_region *z; idx->rr = 0;
fprintf(stderr, "[M::%s] [%lu, %lu) id_n::%lu\n", __func__, s, e, id_n);
for (k = 0; k < id_n; k++) {
p = &(idx->c_idx[id_a[k]]); z = &(ol[ovlp_id(*p)]);
q[0] = z->w_list.a[ovlp_cur_wid(*p)].x_start+ovlp_bd(*p);
q[1] = z->w_list.a[ovlp_cur_wid(*p)].x_end+1-ovlp_bd(*p);
// fprintf(stderr, "[M::%s] tid::%u\t%.*s\twid::%u\tq::[%u, %u)\terr::%d\toerr::%u\n", __func__, ol[ovlp_id(*p)].y_id, (int)Get_NAME_LENGTH(R_INF, ol[ovlp_id(*p)].y_id), Get_NAME(R_INF, ol[ovlp_id(*p)].y_id),
// ovlp_cur_wid(*p), ol[ovlp_id(*p)].w_list.a[ovlp_cur_wid(*p)].x_start, ol[ovlp_id(*p)].w_list.a[ovlp_cur_wid(*p)].x_end+1, ol[ovlp_id(*p)].w_list.a[ovlp_cur_wid(*p)].error, ol[ovlp_id(*p)].non_homopolymer_errors);
if(q[1] <= e) idx->rr = 1;
os = MAX(q[0], s); oe = MIN(q[1], e);
if((oe > os) || ((s == e) && (s >= q[0]) && (s <= q[1]))) {
// if(oe >= os) {
///[-4-][-12-][-4-][-12-]
///[cigar_len][cigar][base_len][base]
m = extract_sub_cigar_ii(z, rref, os, oe, tu, p); an++;
if(m != ((uint32_t)-1)) {///no gap in both sides
kv_push(uint32_t, *b32, m);
}
}
}
oc[0] = b32->n; oc[1] = an + 1; //+1 for the reference read
fprintf(stderr, "-0-[M::%s] oc[0]::%lu, oc[1]::%lu\n", __func__, oc[0], oc[1]);
if(((oc[0] > (oc[1]*occ_max)) && (oc[0] > (oc[1]-oc[0])) && (oc[1] >= occ_tot) && (oc[0] > 1))) {
radix_sort_ec32(b32->a, b32->a+b32->n); an = 0;
for (k = 1, l = 0; k <= b32->n; ++k) {
if (k == b32->n || b32->a[k] != b32->a[l]) {
if(k - l > an) {
an = k - l; a = b32->a + l;
}
l = k;
}
}
oc[0] = an;
fprintf(stderr, "-1-[M::%s] oc[0]::%lu, oc[1]::%lu\n", __func__, oc[0], oc[1]);
if(((oc[0] > (oc[1]*occ_max)) && (oc[0] > (oc[1]-oc[0])) && (oc[1] >= occ_tot) && (oc[0] > 1))) {
(*rc) = a[0];
// prt_cigar0((a[0]<<4)>>20, a[0]>>28);
// prt_bp0((a[0]<<20)>>20, (a[0]<<16)>>28);
return 1;
}
}
return 0;
}
void push_correct0(window_list *idx, window_list_alloc *res, uint32_t len0, uint32_t rc)
{
if(len0 != ((uint32_t)-1)) {
;
} else if(rc != ((uint32_t)-1)) {
uint32_t cc = (rc<<4)>>20, cn = rc>>28, ck = 0, cs, cp;
uint32_t bc = (rc<<20)>>20, bn = (rc<<16)>>28, bk = 0, bs, bp;
for (ck = 0; ck < cn; ck++) {
cs = (cn-1-ck)<<1; cp = (cc>>cs)&3;
bp = (uint32_t)-1;
if(cp != 3) {
bs = (bn-1-bk)<<1; bp = (bc>>bs)&3;
bk++;
}
push_trace_bp(((asg16_v *)(&(res->c))), cp, bp, 1, ((idx->clen>0)?1:0));
fprintf(stderr, "%c", cm[(in >> mp)&3]);
}
}
}
void push_cns_anchor(overlap_region* ol, All_reads *rref, uint64_t s, uint64_t e, UC_Read* tu, cc_idx_t *idx, overlap_region *aux_o, uint64_t is_tail, uint64_t occ_tot, double occ_max, asg32_v* b32)
{
if((!is_tail) && (s >= e)) return;//if s >= e && is_tail = 1, gen the cns of the last a few bases -> s = e = ql
fprintf(stderr, "\n[M::%s] [%lu, %lu)\n", __func__, s, e);
window_list *p = NULL; uint64_t e0 = 0; uint32_t rc;
if(aux_o->w_list.n > 0) {
p = &(aux_o->w_list.a[aux_o->w_list.n-1]);
e0 = p->x_end+1;
///make sure e > s
}
assert(s >= e0);
if((((!is_tail) && (s > 0)) || ((is_tail) && (s > e0)))
&& (cns_gen0(ol, rref, e0, s, tu, idx, occ_tot, occ_max, b32, &rc))) {///CNS in between
push_correct0(p, &(aux_o->w_list), (uint32_t)-1, rc);
} else {
}
kv_pushp(window_list, aux_o->w_list, &p);
p->x_start = s; p->x_end = e-1;
// p->y_start = exz->ps; p->y_end = exz->pe;
}
uint64_t wcns_vote(overlap_region* ol, All_reads *rref, char* qstr, UC_Read* tu, uint64_t *id_a, uint64_t id_n, uint64_t s, uint64_t e, ul_ov_t *c_idx, cc_idx_t *occ, uint64_t occ_tot, double occ_exact, overlap_region *aux_o, asg32_v* b32)
{
uint64_t k, q[2], rr = 0, os, oe, wl, oc[2], fI; ul_ov_t *p, *gp; overlap_region *z;
// uint64_t *ct = occ->idx->a;///occ->idx->a[0, wl<<1)
for (k = 0; k < id_n; k++) {
p = &(c_idx[id_a[k]]); z = &(ol[ovlp_id(*p)]);
q[0] = z->w_list.a[ovlp_cur_wid(*p)].x_start+ovlp_bd(*p);
q[1] = z->w_list.a[ovlp_cur_wid(*p)].x_end+1-ovlp_bd(*p);
if(q[1] <= e) rr = 1;
os = MAX(q[0], s); oe = MIN(q[1], e);
if(oe > os) {
///prepare for CNS
gp = &(occ->c_idx[id_a[k]]);
ovlp_cur_xoff(*gp) = ovlp_cur_xoff(*p); ovlp_cur_yoff(*gp) = ovlp_cur_yoff(*p); ovlp_cur_coff(*gp) = ovlp_cur_coff(*p); ovlp_cur_ylen(*gp) = ovlp_cur_ylen(*p);
// assert(ovlp_cur_wid(*p) == ovlp_cur_wid(*gp));
// assert(ovlp_id(*p) == ovlp_id(*gp));
// fprintf(stderr, "[M::%s] tid::%u\t%.*s\twid::%u\tq::[%u, %u)\tos::%lu\toe::%lu\n", __func__, ol[ovlp_id(*p)].y_id, (int)Get_NAME_LENGTH(R_INF, ol[ovlp_id(*p)].y_id), Get_NAME(R_INF, ol[ovlp_id(*p)].y_id),
// ovlp_cur_wid(*p), ol[ovlp_id(*p)].w_list.a[ovlp_cur_wid(*p)].x_start, ol[ovlp_id(*p)].w_list.a[ovlp_cur_wid(*p)].x_end+1, os, oe);
extract_sub_cigar_mm(z, os, oe, p, occ->idx->a + os - s);
}
}
wl = e - s;
os = occ->mms; oe = occ->mme;
// fprintf(stderr, "[M::%s] s::%lu\te::%lu\n", __func__, s, e);
for (k = 0; k < wl; k++) {
//+1 for the reference read
oc[0] = (occ->idx->a[(k<<1)]>>32) + 1;
oc[1] = ((uint32_t)occ->idx->a[(k<<1)]) + 1;
// fprintf(stderr, "-0-p::%lu\toc[0]::%lu\toc[1]::%lu\tgoc[0]::%lu\tgoc[1]::%u\n", s + k, oc[0], oc[1], (ct[(k<<1)+1]>>32) + 1, ((uint32_t)ct[(k<<1)+1]) + 1);
// if(oc[1] < occ_tot || oc[0] <= 1) {
// ct[(k<<1)] = ct[(k<<1)+1] = 0;
// continue;
// }
// fprintf(stderr, "-1-p::%lu\toc[0]::%lu\toc[1]::%lu\n", s + k, oc[0], oc[1]);
if((oc[0] > (oc[1]*occ_exact)) && (oc[0] > (oc[1]-oc[0])) && (oc[1] >= occ_tot) && (oc[0] > 1)) {
///note: there might be insertions at q[k-1, k], insead if q[k, k+1]
fI = 1;
///make sure there is no insertion
//+1 for the reference read
oc[0] = (occ->idx->a[(k<<1)+1]>>32) + 1;
oc[1] = ((uint32_t)occ->idx->a[(k<<1)+1]) + 1;
if(((oc[0] > (oc[1]*occ_exact)) && (oc[0] > (oc[1]-oc[0])) && (oc[1] >= occ_tot) && (oc[0] > 1))) fI = 0;
if(fI) {
// fprintf(stderr, "-1-p::%lu\toc[0]::%lu\toc[1]::%u\tgoc[0]::%lu\tgoc[1]::%u\n", s + k, (occ->idx->a[(k<<1)]>>32) + 1, ((uint32_t)occ->idx->a[(k<<1)]) + 1, (occ->idx->a[(k<<1)+1]>>32) + 1, ((uint32_t)occ->idx->a[(k<<1)+1]) + 1);
if(oe > os && os != ((uint64_t)-1)) {///push previous intervals
push_cns_anchor(ol, rref, os, oe, tu, occ, aux_o, 0, occ_tot, occ_exact, b32);
}
os = oe = (uint64_t)-1;
}
//+1 for the reference read
oc[0] = (occ->idx->a[(k<<1)]>>32) + 1;
oc[1] = ((uint32_t)occ->idx->a[(k<<1)]) + 1;
if((s+k) == oe) {
oe++;
} else {
if(oe > os && os != ((uint64_t)-1)) {///push previous intervals
push_cns_anchor(ol, rref, os, oe, tu, occ, aux_o, 0, occ_tot, occ_exact, b32);
}
os = s+k; oe = s+k+1;
}
} else {
// fprintf(stderr, "-2-p::%lu\toc[0]::%lu\toc[1]::%u\tgoc[0]::%lu\tgoc[1]::%u\n", s + k, (occ->idx->a[(k<<1)]>>32) + 1, ((uint32_t)occ->idx->a[(k<<1)]) + 1, (occ->idx->a[(k<<1)+1]>>32) + 1, ((uint32_t)occ->idx->a[(k<<1)+1]) + 1);
if(oe > os && os != ((uint64_t)-1)) {///push previous intervals
push_cns_anchor(ol, rref, os, oe, tu, occ, aux_o, 0, occ_tot, occ_exact, b32);
}
os = oe = (uint64_t)-1;
}
occ->idx->a[(k<<1)] = occ->idx->a[(k<<1)+1] = 0;
}
occ->mms = occ->mme = (uint64_t)-1;
if(oe > os && os != ((uint64_t)-1)) {
occ->mms = os; occ->mme = oe;
}
return rr;
}
void print_debug_ovlp_cigar(overlap_region_alloc* ol, asg64_v* idx, kv_ul_ov_t *c_idx)
{
uint64_t k, ci; uint32_t cl; ul_ov_t *cp; bit_extz_t ez; uint16_t c; char cm[4];
cm[0] = 'M'; cm[1] = 'S'; cm[2] = 'I'; cm[3] = 'D';
for (k = 0; k < idx->n; k++) {
cp = &(c_idx->a[(uint32_t)idx->a[k]]);
fprintf(stderr, "**********[M::%s] tid::%u\t%.*s\twid::%u\tq::[%u, %u)\terr::%d\toerr::%u**********\n", __func__, ol->list[ovlp_id(*cp)].y_id, (int)Get_NAME_LENGTH(R_INF, ol->list[ovlp_id(*cp)].y_id), Get_NAME(R_INF, ol->list[ovlp_id(*cp)].y_id),
ovlp_cur_wid(*cp), ol->list[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_start, ol->list[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_end+1, ol->list[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].error, ol->list[ovlp_id(*cp)].non_homopolymer_errors);
set_bit_extz_t(ez, ol->list[ovlp_id(*cp)], ovlp_cur_wid(*cp)); ci = 0;
while (ci < ez.cigar.n) {
ci = pop_trace(&(ez.cigar), ci, &c, &cl);
fprintf(stderr, "%u%c", cl, cm[c]);
}
fprintf(stderr, "\n");
}
}
void wcns_gen(overlap_region_alloc* ol, All_reads *rref, UC_Read* qu, UC_Read* tu, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, int64_t bd, uint64_t wl, int64_t ql, uint64_t occ_tot, double occ_exact, overlap_region *aux_o, asg32_v* b32)
{
int64_t on = ol->length, k, i, zwn, q[2];
uint64_t m, *ra, rn; overlap_region *z; ul_ov_t *cp;
for (k = idx->n = c_idx->n = 0; k < on; k++) {
z = &(ol->list[k]); zwn = z->w_list.n;
if((!zwn) || (z->is_match != 1)) continue;
for (i = 0; i < zwn; i++) {
if(is_ualn_win(z->w_list.a[i])) continue;
q[0] = z->w_list.a[i].x_start; q[1] = z->w_list.a[i].x_end;
q[0] += bd; q[1] -= bd;
if(q[1] >= q[0]) {
m = ((uint64_t)q[0]); m <<= 32;
m += c_idx->n; kv_push(uint64_t, *idx, m);
kv_pushp(ul_ov_t, *c_idx, &cp);
ovlp_id(*cp) = k; ///ovlp id
// ovlp_min_wid(*cp) = i; ///beg id of windows
// ovlp_max_wid(*cp) = i; ///end id of windows
ovlp_cur_wid(*cp) = i; ///cur id of windows
ovlp_cur_xoff(*cp) = z->w_list.a[i].x_start; ///cur xpos
ovlp_cur_yoff(*cp) = z->w_list.a[i].y_start; ///cur xpos
ovlp_cur_ylen(*cp) = 0;
ovlp_cur_coff(*cp) = 0; ///cur cigar off in cur window
ovlp_bd(*cp) = bd;
}
}
}
int64_t srt_n = idx->n, s, e, t, rr; i = 0;
radix_sort_ec64(idx->a, idx->a+idx->n);
for (k = 1, i = 0; k < srt_n; k++) {
if (k == srt_n || (idx->a[k]>>32) != (idx->a[i]>>32)) {
if(k - i > 1) {
for (t = i; t < k; t++) {
cp = &(c_idx->a[(uint32_t)idx->a[t]]);
// s = ol->list[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_start+ovlp_bd(*cp);
// assert(s == (int64_t)(idx->a[i]>>32));
m = ol->list[ovlp_id(*cp)].w_list.a[ovlp_cur_wid(*cp)].x_end+1-ovlp_bd(*cp);
m <<= 32; m += ((uint32_t)idx->a[t]); idx->a[t] = m;
// fprintf(stderr, "[M::%s] s::%ld\tsi::%lu\n", __func__, s, (idx->a[i]>>32));
}
radix_sort_ec64(idx->a + i, idx->a + k);
}
i = k;
}
}
print_debug_ovlp_cigar(ol, idx, c_idx);
///second index
kv_resize(ul_ov_t, *c_idx, (c_idx->n<<1));
ul_ov_t *idx_a = NULL, *idx_b = NULL;
idx_a = c_idx->a; idx_b = c_idx->a;
memcpy(idx_b, idx_a, c_idx->n * (sizeof((*(idx_a)))));
kv_resize(uint64_t, *buf, ((wl<<1) + idx->n)); buf->n = ((wl<<1) + idx->n);
memcpy(buf->a + (wl<<1), idx->a, idx->n * (sizeof((*(idx->a)))));
memset(buf->a, 0, (wl<<1)*(sizeof((*(idx->a)))));
cc_idx_t ii_a, ii_b; memset(&ii_a, 0, sizeof(ii_a)); memset(&ii_b, 0, sizeof(ii_b));
ii_a.c_idx = idx_a; ii_a.idx = idx; ii_a.i = ii_a.i0 = 0; ii_a.srt_n = ii_a.idx->n; ii_a.mms = ii_a.mme = (uint64_t)-1;
ii_b.c_idx = idx_b; ii_b.idx = buf; ii_b.i = ii_b.i0 = (wl<<1); ii_b.srt_n = ii_b.idx->n; ii_b.mms = ii_b.mme = (uint64_t)-1;
s = 0; e = wl; e = ((e<=ql)?e:ql); rr = 0;
aux_o->w_list.n = aux_o->w_list.c.n = 0; ///for cigar
for (; s < ql; ) {
rn = iter_cc_idx_t(ol->list, &ii_a, s, e, rr, 0, &ra);
// debug_inter0(ol->list, ii_a.c_idx, ii_a.idx->a + ii_a.i0, ii_a.srt_n - ii_a.i0, ra, rn, s, e, 0, 1, "-0-");
rr = wcns_vote(ol->list, rref, qu->seq, tu, ra, rn, s, e, ii_a.c_idx, &ii_b, occ_tot, occ_exact, aux_o, b32);
s += wl; e += wl; e = ((e<=ql)?e:ql);
}
if(ii_b.mme > ii_b.mms && ii_b.mms != (uint64_t)-1) {
push_cns_anchor(ol->list, rref, ii_b.mms, ii_b.mme, tu, &ii_b, aux_o, 0, occ_tot, occ_exact, b32);
}
push_cns_anchor(ol->list, rref, ql, ql, tu, &ii_b, aux_o, 1, occ_tot, occ_exact, b32);
}
static void worker_hap_ec(void *data, long i, int tid)
{
ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]);
uint32_t high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO);
uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO;
overlap_region *aux_o = NULL; asg64_v buf0;
// if (memcmp("m64012_190920_173625/7210046/ccs", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
// fprintf(stderr, "[M::%s-beg] rid->%ld\n", __func__, i);
// } else {
// return;
// }
if(i != 596/**1024**/) return;
recover_UC_Read(&b->self_read, &R_INF, i);
h_ec_lchain(b->ab, i, b->self_read.seq, b->self_read.length, asm_opt.mz_win, asm_opt.k_mer_length, &R_INF, &b->olist, &b->clist, 0.02, asm_opt.max_n_chain, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 0, 2, 2, UINT32_MAX);
b->num_read_base += b->olist.length;
aux_o = fetch_aux_ovlp(&b->olist);///must be here
gen_hc_r_alin(&b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_ov_diff_ec, WINDOW_HC, i, E_KHIT/**asm_opt.k_mer_length**/, 1, &b->v16);
// fprintf(stderr, "\n[M::%s] rid::%ld\t%.*s\tlen::%lld\tocc::%lu\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i),
// Get_NAME(R_INF, i), b->self_read.length, b->olist.length);
b->num_correct_base += b->olist.length;
copy_asg_arr(buf0, b->sp);
rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1);
copy_asg_arr(b->sp, buf0);
copy_asg_arr(buf0, b->sp);
wcns_gen(&b->olist, &R_INF, &b->self_read, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, 512, b->self_read.length, 3, 0.500001, aux_o, &b->v32);
copy_asg_arr(b->sp, buf0);
uint32_t k;
for (k = 0; k < b->olist.length; k++) {
if(b->olist.list[k].is_match == 1) b->num_recorrect_base++;
}
// exit(1);
// prt_chain(&b->olist);
// ul_map_lchain(b->abl, (uint32_t)-1, s->seq[i], s->len[i], s->opt->w, s->opt->k, s->uu, &b->olist, &b->clist, s->opt->bw_thres,
// s->opt->max_n_chain, 1, NULL, &(b->tmp_region), NULL, &(b->sp), &high_occ, NULL, 0, 1, 0.2/**0.75**/, 2, 3);
/**
int fully_cov, abnormal;
// if(i != 12578) return;
// fprintf(stderr, "[M::%s-beg] rid->%ld\n", __func__, i);
// if (memcmp("7897e875-76e5-42c8-bc37-94b370c4cc8d", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
// fprintf(stderr, "[M::%s-beg] rid->%ld\n", __func__, i);
// } else {
// return;
// }
ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist,
0.02, asm_opt.max_n_chain, 1, NULL, &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL, &(b->sp));
clear_Cigar_record(&b->cigar1);
clear_Round2_alignment(&b->round2);
correct_overlap(&b->olist, &R_INF, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon,
&b->cigar1, &b->hap, &b->round2, &b->r_buf, &(b->tmp_region.w_list), 0, 1, &fully_cov, &abnormal);
b->num_read_base += b->self_read.length;
b->num_correct_base += b->correct.corrected_base;
b->num_recorrect_base += b->round2.dumy.corrected_base;
push_cigar(R_INF.cigars, i, &b->cigar1);
push_cigar(R_INF.second_round_cigar, i, &b->round2.cigar);
R_INF.paf[i].is_fully_corrected = 0;
if (fully_cov) {
if (get_cigar_errors(&b->cigar1) == 0 && get_cigar_errors(&b->round2.cigar) == 0)
R_INF.paf[i].is_fully_corrected = 1;
}
R_INF.paf[i].is_abnormal = abnormal;
R_INF.trio_flag[i] = AMBIGU;
///need to be fixed in r305
// if(ha_idx_hp == NULL)
// {
// R_INF.trio_flag[i] += collect_hp_regions(&b->olist, &R_INF, &(b->k_flag), RESEED_HP_RATE, Get_READ_LENGTH(R_INF, i), NULL);
// }
if (R_INF.trio_flag[i] != AMBIGU || b->save_ov) {
int is_rev = (asm_opt.number_of_round % 2 == 0);
push_overlaps(&(R_INF.paf[i]), &b->olist, 1, &R_INF, is_rev);
push_overlaps(&(R_INF.reverse_paf[i]), &b->olist, 2, &R_INF, is_rev);
}
if(het_cnt) het_cnt[i] = get_het_cnt(&b->hap);
// fprintf(stderr, "[M::%s-end] rid->%ld\n", __func__, i);
**/
}
void cal_ec_multiple(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a)
{
double tt0 = yak_realtime_0();
kt_for(n_thre, worker_hap_ec, b, n_a);///debug_for_fix
fprintf(stderr, "[M::%s-reads] #->%lu\n", __func__, n_a);
fprintf(stderr, "[M::%s::%.3f] ==> chaining\n", __func__, yak_realtime_0()-tt0);
}
ecovlp.h
+64
View File
@@ -0,0 +1,64 @@
#ifndef __ECOVLP_PARSER__
#define __ECOVLP_PARSER__
#define __STDC_LIMIT_MACROS
#include <stdint.h>
#include "Hash_Table.h"
#include "Process_Read.h"
typedef struct {
uint32_t v:31, f:1;
uint32_t sc;
} cns_arc;
typedef struct {size_t n, m; cns_arc *a;} cns_arc_v;
typedef struct {
// uint16_t c:2, t:2, f:1, sc:3;
uint32_t c:2, f:1, sc:29;
cns_arc_v in, ou;
}cns_t;
typedef struct {
size_t n, m;
cns_t *a;
uint32_t si, ei;
}cns_gfa;
typedef struct {
int is_final, save_ov;
// chaining and overlapping related buffers
UC_Read self_read, ovlp_read;
Candidates_list clist;
overlap_region_alloc olist;
overlap_region tmp;
ha_abuf_t *ab;
// error correction related buffers
int64_t num_read_base, num_correct_base, num_recorrect_base;
Cigar_record cigar;
Correct_dumy correct;
haplotype_evdience_alloc hap;
bit_extz_t exz;
// asg32_v v32;
kv_ul_ov_t pidx;
asg64_v v64;
asg32_v v32;
asg16_v v16;
kvec_t_u64_warp r_buf;
kvec_t_u8_warp k_flag;
st_mt_t sp;
cns_gfa cns;
} ec_ovec_buf_t0;
typedef struct {
ec_ovec_buf_t0 *a;
uint32_t n;
} ec_ovec_buf_t;
ec_ovec_buf_t* gen_ec_ovec_buf_t(uint32_t n, uint32_t is_final, uint32_t save_ov);
void destroy_ec_ovec_buf_t(ec_ovec_buf_t *p);
void cal_ec_multiple(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a);
void prt_chain(overlap_region_alloc *o);
#endif
htab.h
+1
View File
@@ -108,6 +108,7 @@ uint64_t ha_abufl_mem(const ha_abufl_t *ab);
double yak_cputime(void);
void yak_reset_realtime(void);
double yak_realtime_0(void);
double yak_realtime(void);
long yak_peakrss(void);
double yak_peakrss_in_gb(void);
sys.cpp
+6
View File
@@ -31,6 +31,12 @@ double yak_realtime(void)
return yak_realtime_core() - yak_realtime0;
}
double yak_realtime_0(void)
{
return yak_realtime_core();
}
long yak_peakrss(void)
{
struct rusage r;