xxay/hifiasm
1
0
Fork 0
mirror of https://github.com/chhylp123/hifiasm.git synced 2026-09-30 14:38:12 +08:00
Code Issues Packages Projects Releases Wiki Activity

regen_scb

Browse Source
This commit is contained in:
chhylp123
2026-05-03 03:34:23 -04:00
parent f5078f7b23
commit 7884b5ad88
15 changed files with 2177 additions and 133 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Assembly.cpp
+12 -5
View File
@@ -1026,7 +1026,7 @@ void ha_ec(int64_t round, int num_pround, int des_idx, uint64_t *tot_b, uint64_t
write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name);
if((asm_opt.flag & HA_F_VERBOSE_GFA) && (asm_opt.bin_only == 1)) exit(1);///just for debug
}
if (w_tmp) tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name, round, asm_opt.number_of_round, 0);
if (w_tmp) tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &scb, &asm_opt, asm_opt.output_file_name, round, asm_opt.number_of_round, 0);
// Output_corrected_fastq();
@@ -2077,7 +2077,7 @@ int ha_assemble(void)
// quick_debug_phasing(MC_NAME);
extern void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o);
int r, r0 = -1, hom_cov = -1, ovlp_loaded = 0; uint64_t tot_b, tot_e;
if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) {
if ((asm_opt.load_index_from_disk) && (asm_opt.dbg_ec_rr < 0) && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) {
ovlp_loaded = 1;
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> loaded corrected reads and overlaps from disk\n", __func__, yak_realtime(), yak_cpu_usage());
if (asm_opt.extract_list) {
@@ -2096,12 +2096,18 @@ int ha_assemble(void)
if((asm_opt.flag & HA_F_VERBOSE_GFA)) load_pt_index(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name), load_ct_index(&ha_ct_table, asm_opt.output_file_name);
r = ha_idx?asm_opt.number_of_round-1:0;
if((!ha_idx) && (asm_opt.restart)) {
for (r = asm_opt.number_of_round - 1; r >= 0; --r) {
if(tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name, r, asm_opt.number_of_round, 1)) {
r = asm_opt.number_of_round - 1;
if(asm_opt.dbg_ec_rr >= 0) r = asm_opt.dbg_ec_rr;
for (; r >= 0; --r) {
if(tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &scb, &asm_opt, asm_opt.output_file_name, r, asm_opt.number_of_round, 1)) {
load_ct_index(&ha_ct_table, asm_opt.output_file_name); r0 = r;
break;
}
}
if((asm_opt.dbg_ec_rr >= 0) && (asm_opt.dbg_ec_rr != r)) {
fprintf(stderr, "[E::%s] no matching debug error-correction bins found\n", __func__);
exit(1);
}
if(r < 0) r = 0;
}
@@ -2123,6 +2129,7 @@ int ha_assemble(void)
// fprintf(stderr, "[M::%s] # bases: %lld; # corrected bases: %lld; # recorrected bases: %lld\n", __func__,
// asm_opt.num_bases, asm_opt.num_corrected_bases, asm_opt.num_recorrected_bases);
// fprintf(stderr, "[M::%s] size of buffer: %.3fGB\n", __func__, asm_opt.mem_buf / 1073741824.0);
if(asm_opt.dbg_ec_rr >= 0) exit(1);
}
if (asm_opt.flag & HA_F_WRITE_EC) {
if(asm_opt.is_sc) Output_corrected_fastq();
@@ -2149,7 +2156,7 @@ int ha_assemble(void)
build_string_graph_without_clean(asm_opt.min_overlap_coverage, R_INF.paf, R_INF.reverse_paf,
R_INF.total_reads, R_INF.read_length, asm_opt.min_overlap_Len, asm_opt.max_hang_Len, asm_opt.clean_round,
asm_opt.gap_fuzz, asm_opt.min_drop_rate, asm_opt.max_drop_rate, asm_opt.output_file_name, asm_opt.large_pop_bubble_size, 0, !ovlp_loaded);
destory_All_reads(&R_INF); if(asm_opt.dbg_bam) destroy_cc_v(&scb);
destory_All_reads(&R_INF); /**if(asm_opt.dbg_bam)**/ destroy_cc_v(&scb);
return 0;
}
CommandLines.cpp
+5
View File
@@ -94,6 +94,7 @@ static ko_longopt_t long_options[] = {
{ "hyb-syn", ko_required_argument, 376},
{ "simd-m", ko_required_argument, 377},
{ "del-hf", ko_no_argument, 378},
{ "dbg-rr", ko_required_argument, 379},
// { "path-round", ko_required_argument, 348},
{ 0, 0, 0 }
};
@@ -444,6 +445,8 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->simd_mm = -1;
asm_opt->del_hf = 0;
asm_opt->dbg_ec_rr = -1;
}
void destory_enzyme(enzyme* f)
@@ -1121,6 +1124,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
asm_opt->simd_mm = atoi(opt.arg);
} else if (c == 378) {
asm_opt->del_hf = 1;
} else if (c == 379) {
asm_opt->dbg_ec_rr = atoi(opt.arg);
} else if (c == 'l') { ///0: disable purge_dup; 1: purge containment; 2: purge overlap
asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg);
}
CommandLines.h
+3 -1
View File
@@ -5,7 +5,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.25.1-r910"
#define HA_VERSION "0.25.1-r920"
#define VERBOSE 0
@@ -204,6 +204,8 @@ typedef struct {
int8_t simd_mm;
int8_t del_hf;
int64_t dbg_ec_rr;
} hifiasm_opt_t;
extern hifiasm_opt_t asm_opt;
Correct.cpp
+997 -47
View File
File diff suppressed because it is too large Load Diff
Correct.h
+3 -2
View File
@@ -1456,8 +1456,9 @@ void gen_hc_r_alin_adv(gen_hc_aln_t *ez);
uint64_t gen_hc_r_alin_nec(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, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max, uint64_t sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, asg64_v *kp, uint8_t *hpf);
void gen_hc_r_alin_nec_adv(gen_hc_aln_t *ez);
uint64_t gen_hc_r_alin_re(overlap_region* z, Candidates_list *cl, char* qstr, uint64_t ql, char* tstr, uint64_t tl, 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/**, uint8_t is_dbg**/, uint64_t rid, uint64_t hpc_len, uint64_t std_bs, Chain_Data *dp, asg8_v *q8, asg8_v *t8, uint8_t lindel, uint64_t tcut, uint64_t site_sc, int64_t h0_w, asg32_v *b32,
int64_t hap_cov_match, int64_t hap_cov_unmatch, int64_t het_cov_a, int64_t hom_cov_a, int64_t n_hap, double hf_rate);
uint64_t gen_hc_r_alin_self(overlap_region* z, Candidates_list *cl, char* qstr, uint64_t ql, char* tstr, uint64_t tl, 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, asg16_v* scc);
void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, asg16_v *qc0, 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/**, uint8_t is_dbg**/, uint64_t rid, uint64_t hpc_len, uint64_t std_bs, Chain_Data *dp, asg8_v *q8, asg8_v *t8, uint8_t lindel, uint64_t tcut, uint64_t site_sc, int64_t h0_w, asg32_v *b32,
int64_t hap_cov_match, int64_t hap_cov_unmatch, int64_t het_cov_a, int64_t hom_cov_a, int64_t n_hap, double hf_rate, overlap_region *rchn, int64_t rref_len);
void set_exact_exz(bit_extz_t *exz, int64_t qs, int64_t qe, int64_t ts, int64_t te);
void push_alnw(overlap_region *aux_o, bit_extz_t *exz);
void cal_exz_global(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez);
Hash_Table.cpp
+481
View File
@@ -1540,6 +1540,63 @@ inline int32_t comput_sc_ch_ec(const k_mer_hit *ai, const k_mer_hit *aj, double
return sc;
}
inline int32_t comput_sc_ch_ec_global(const k_mer_hit *ai, const k_mer_hit *aj, double bw_rate, double chn_pen_gap, double chn_pen_skip, int64_t sl, int64_t ol, uint8_t no_adj)
{
///ai is the suffix of aj
int32_t dq, dr, dd, dg, q_span, sc; double dg_of;
dq = (int64_t)(ai->self_offset) - (int64_t)(aj->self_offset);
dr = (int64_t)(ai->offset) - (int64_t)(aj->offset);
if((dq < 0) || (dr < 0)) return INT32_MIN;
if((no_adj) && ((dq == 0) || (dr == 0))) return INT32_MIN;
dd = dr > dq? dr - dq : dq - dr;//gap
if((dd > 16) && (dd > cal_bw(ai, aj, bw_rate, sl, ol))) return INT32_MIN;
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)) {
double lin_pen, a_pen;
lin_pen = (chn_pen_gap*(double)dd);
dg_of = (dg>0)?((double)dg):(0.333333);
a_pen = ((double)(sc))*((((double)dd)/dg_of)/bw_rate);
///for long gap
// if(lin_pen > a_pen) lin_pen = a_pen;
if(dd < 4) lin_pen = ((lin_pen > a_pen)?(a_pen):(lin_pen));
else lin_pen = ((lin_pen < a_pen)?(a_pen):(lin_pen));
lin_pen += (chn_pen_skip*dg_of);
sc -= (int32_t)lin_pen;
}
return sc;
}
inline int32_t comput_sc_ff_adv(const k_mer_hit *ai, const k_mer_hit *aj, double bw_rate, double chn_pen_gap, double chn_pen_skip, int64_t sl, int64_t ol)
{
///ai is the suffix of aj
int32_t dq, dr, dd, dg, q_span, sc;
dq = (int64_t)(ai->self_offset) - (int64_t)(aj->self_offset);
if(dq < 0) return INT32_MIN;
dr = (int64_t)(ai->offset) - (int64_t)(aj->offset);
if(dr < 0) return INT32_MIN;
dd = dr > dq? dr - dq : dq - dr;//gap
// if((dd > 16) && (dd > cal_bw(ai, aj, bw_rate, sl, ol))) return INT32_MIN;
dg = dr < dq? dr : dq;//len
if(dg <= 0) return INT32_MIN;
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)) {
double lin_pen, a_pen;
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;
}
return sc;
}
inline int32_t comput_sc_ff(const k_mer_hit *ai, const k_mer_hit *aj, double bw_rate, double chn_pen_gap, double chn_pen_skip, int64_t sl, int64_t ol)
{
///ai is the suffix of aj
@@ -2283,6 +2340,315 @@ uint64_t lchain_qdp_mcopy_fast(Candidates_list *cl, int64_t a_idx, int64_t a_n,
}
/**
void quick_ck_lchain_global(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, int32_t *gf, int64_t *gp, int64_t *si, int64_t *ei)
{
*si = 0; *ei = a_n; gf[0] = gf[1] = gp[0] = gp[1] = INT32_MIN;
if((a_n <= 0) || (xl <= 0) || (yl <= 0)) return;
int64_t l, k, is_srt = 1, z; k_mer_hit *ai, *aj, ft, fz;
int64_t dq, dr, dd, dg, q_span, sc, csc, ddt; uint8_t ff; double lin_pen, a_pen, dg_of;
ft.cnt = fz.cnt = 0xFFFFFF00u; ft.readID = fz.readID = a[0].readID;
ft.offset = ft.self_offset = 0; fz.offset = yl-1; fz.self_offset = xl-1;
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) {
ddt = 0; ff = 0; p[l] = f[l] = INT32_MIN;
ft.strand = a[l].strand;
aj = &ft; z = l; ai = &a[z];
dq = (int64_t)(ai->self_offset) - (int64_t)(aj->self_offset);
dr = (int64_t)(ai->offset) - (int64_t)(aj->offset);
dd = dr > dq? dr - dq : dq - dr;//gap
if((dd > 16) && (dd > cal_bw(ai, aj, bw_rate, xl, yl))) ff = 1;
if(!ff) {
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);
dg_of = ((dg>0)?(dg):(0.333333));
a_pen = ((double)(sc))*((((double)dd)/dg_of)/bw_rate);
///for long gap
// if(lin_pen > a_pen) lin_pen = a_pen;
if(dd < 4) lin_pen = ((lin_pen > a_pen)?(a_pen):(lin_pen));
else lin_pen = ((lin_pen < a_pen)?(a_pen):(lin_pen));
lin_pen += (chn_pen_skip*dg_of);
sc -= (int32_t)lin_pen;
}
csc = a[z].cnt&(0xffu); if(sc < csc) ff = 1;
if(!ff) {
p[z] = -1; f[z] = sc; ddt += dd;
}
}
if(!ff) {
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);
dg_of = ((dg>0)?(dg):(0.333333));
a_pen = ((double)(sc))*((((double)dd)/dg_of)/bw_rate);
///for long gap
// if(lin_pen > a_pen) lin_pen = a_pen;
if(dd < 4) lin_pen = ((lin_pen > a_pen)?(a_pen):(lin_pen));
else lin_pen = ((lin_pen < a_pen)?(a_pen):(lin_pen));
lin_pen += (chn_pen_skip*dg_of);
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(z < k) ff = 1;
}
if(!ff) {
fz.strand = a[l].strand;
ai = &fz; aj = &(a[k-1]);
dq = (int64_t)(ai->self_offset) - (int64_t)(aj->self_offset);
dr = (int64_t)(ai->offset) - (int64_t)(aj->offset);
dd = dr > dq? dr - dq : dq - dr;//gap
if(((dd > 16) && (dd > cal_bw(ai, aj, bw_rate, xl, yl)))||
(((ddt + dd) > 16) && ((ddt + dd) > cal_bw(&fz, &ft, bw_rate, xl, yl)))) {
ff = 1;
} else {
dg = dr < dq? dr : dq;//len
sc = q_span = 0;
if (dd || (dg > q_span && dg > 0)) {
lin_pen = (chn_pen_gap*(double)dd);
dg_of = ((dg>0)?(dg):(0.333333));
a_pen = ((double)(sc))*((((double)dd)/dg_of)/bw_rate);
///for long gap
// if(lin_pen > a_pen) lin_pen = a_pen;
if(dd < 4) lin_pen = ((lin_pen > a_pen)?(a_pen):(lin_pen));
else lin_pen = ((lin_pen < a_pen)?(a_pen):(lin_pen));
lin_pen += (chn_pen_skip*dg_of);
sc -= (int32_t)lin_pen;
}
sc += f[k-1]; ///csc = a[z].cnt&(0xffu); if(sc < csc) break;
gp[fz.strand] = k-1; gf[fz.strand] = sc; ///ddt += dd;
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_global_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 khit_n)
{
if(a_n <= 0) return 0;
int64_t *p, *t, *gp, max_f, n_skip, st, max_j, end_j, sc, max_ii, ovl, min_sc, ch_n, si, ei;
int32_t *f, max, tmp, *ii, *gf; int64_t i, k, j, cL = 0; k_mer_hit* a; k_mer_hit* des; k_mer_hit *swap, ft; overlap_region *z;
resize_Chain_Data(dp, a_n + 2, NULL); ch_n = 1;
t = dp->tmp; f = dp->score; p = dp->pre; ii = dp->occ; gp = p + a_n; gf = f + a_n; gf[0] = gf[1] = gp[0] = gp[1] = INT32_MIN;
a = cl->list + a_idx; des = cl->list + des_idx;
if(quick_check) {
quick_ck_lchain_global(a, a_n, xl, yl, chn_pen_gap, chn_pen_skip, bw_rate, p, t, f, ii, gf, gp, &si, &ei);
} else {
si = 0; ei = a_n; memset(t, 0, (a_n*sizeof((*t))));
}
ft.cnt = 0xFFFFFF00u; ft.readID = a[0].readID; ft.strand = 0; ft.offset = ft.self_offset = 0;
for (i = st = si, max_ii = -1; i < ei; ++i) {
///max_f = a[i].cnt&(0xffu);
ft.strand = a[i].strand;
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_ec_global(&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_ec(&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 == 27105 || a[0].readID == 7603)) {///r833
// fprintf(stderr, "i::%ld[M::%s::rid->%u::%c] q::%u, t::%u, st::%ld, max_ii::%ld, f[i]::%d, p[i]::%ld, msc_i::%ld, msc::%ld, movl::%ld\n",
// i, __func__, a[i].readID, "+-"[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(mcopy_num > 1) {
// if(a[0].readID == 4412344) {
// fprintf(stderr, "[M::%s::] msc::%ld, cL::%ld\n", __func__, msc, cL);
// }
if(cL >= mcopy_khit_cutoff) {///if there are too few k-mers, disable mcopy
msc -= plus; min_sc = msc*mcopy_rate; 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)) {///!(ii[i]): skip the best chain
t[ch_n] = ((uint64_t)f[i])<<32; t[ch_n] += (i<<1); ch_n++;
}
}
// if(a[0].readID == 4412344) {
// fprintf(stderr, "[M::%s::] msc::%ld, min_sc::%ld, cL::%ld, ch_n::%ld, mcopy_num::%ld\n", __func__, msc, min_sc, cL, ch_n, mcopy_num);
// }
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 && n_u < mcopy_num; --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(a[0].readID == 4412344) {
// fprintf(stderr, "+[M::%s::] sc::%ld, n_a::%ld\n", __func__, sc, n_v-n_v0);
// }
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]]));
// if(a[0].readID == 4412344) {
// fprintf(stderr, "-[M::%s::] sc::%ld, n_a::%ld, q::[%u,%u), t::[%u,%u), %c\n", __func__, sc, n_v-n_v0, z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1, "+-"[z->y_pos_strand]);
// }
///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))) {
if((!n_u) || (n_v - n_v0 > 1)) {
z->align_length = n_v-n_v0; z->x_id = n_v0;
n_u++;
} else {///non-best is tiny
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 { \
(an).self_offset = (xl)-1-((an).self_offset+1-((an).cnt&((uint32_t)(0xffu)))); \
(an).offset = (yl)-1-((an).offset+1-((an).cnt&((uint32_t)(0xffu))));\
@@ -2453,6 +2819,121 @@ uint64_t lchain_qdp_fix(k_mer_hit* a, int64_t a_n, Chain_Data* dp, int64_t max_s
}
uint64_t lchain_qdp_fix_adv(k_mer_hit *a, int64_t a_n, Chain_Data* dp, 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 xl, int64_t yl, int64_t quick_check,
int64_t left_fix, int64_t right_fix, k_mer_hit *res)
{
if(a_n <= 0) return 0;
int64_t *p, *t, max_f, n_skip, st, max_j, end_j, sc, msc, msc_i, bw, max_ii, ovl, movl;
int32_t *f, max, tmp; int64_t i, j, ret, cL = 0;
resize_Chain_Data(dp, a_n, NULL);
t = dp->tmp; f = dp->score; p = dp->pre;
bw = ((xl < yl)?xl:yl); bw *= bw_rate;
msc = msc_i = -1; movl = INT32_MAX;
if(quick_check) {
ret = lchain_qcheck(a, a_n, dp, bw_rate);
if (ret > 0) {
a_n = ret; msc_i = a_n-1; msc = f[msc_i];
goto skip_ldp;
}
}
memset(t, 0, (a_n*sizeof((*t))));
for (i = st = 0, max_ii = -1; i < a_n; ++i) {
max_f = a[i].cnt&(0xffu); if(left_fix && i > 0) max_f = INT32_MIN;
n_skip = 0; max_j = end_j = -1;
if ((i-st) > max_iter) st = i-max_iter;
for (j = i - 1; j >= 0; --j) {
if(left_fix && f[j] == INT32_MIN)continue;
sc = comput_sc_ff_adv(&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) {
if((max_j != -1) || (left_fix == 0)) break;
}
}
if (p[j] >= 0) t[p[j]] = i;
///put it here will allow at least one prefix no matter max_dis
///this is special for gap filling, not for chaining
if (a[i].self_offset > (max_dis + a[j].self_offset)) {
if((max_j != -1)) break;
}
if (j < st) {
if((max_j != -1) || (left_fix == 0)) break;
}
}
end_j = j;
if (max_ii < 0 || ((int64_t)a[i].self_offset) - ((int64_t)a[max_ii].self_offset) > max_dis) {
max = INT32_MIN; max_ii = -1;
for (j = i - 1; (j >= st) && ((((int64_t)a[i].self_offset)-((int64_t)a[j].self_offset))<=max_dis); --j) {
if ((f[j] != INT32_MIN) && (max < f[j])) {
max = f[j], max_ii = j;
}
}
}
if ((max_ii >= 0) && (max_ii < end_j) && (f[max_ii] != INT32_MIN)) {///just have a try with a[i]<->a[max_ii]
tmp = comput_sc_ff_adv(&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;
}
p[i] = max_j; f[i] = max_f;
if ((max_ii < 0) || (((((int64_t)a[i].self_offset)-((int64_t)a[max_ii].self_offset))<=max_dis) && (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;
}
}
}
skip_ldp:
if(right_fix && f[a_n-1] == INT32_MIN) return 0;
if(right_fix) msc_i = a_n-1;
///a[] has been sorted by self_offset
i = msc_i; cL = 0;
while (i >= 0) {
t[cL++] = i; msc_i = i; i = p[i];
}
n_skip = cL>>1;
for (i = 0; i < n_skip; i++) {
msc_i = t[i]; t[i] = t[cL-i-1]; t[cL-i-1] = msc_i;
}
if((cL > 0) && (right_fix) && (t[cL-1] != (a_n-1))) {
cL = 0;
}
if((cL > 0) && (left_fix) && (t[0] != 0)) {
cL = 0;
}
if(cL > 0 && res) {
for (i = 0; i < cL; i++) {
res[i] = a[t[i]];
}
}
return cL;
}
uint64_t lchain_refine(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp,
int64_t max_skip, int64_t max_iter, int64_t max_dis, int64_t long_gap)
{
Hash_Table.h
+10
View File
@@ -234,6 +234,10 @@ uint64_t lchain_qdp_fix(k_mer_hit* a, int64_t a_n, Chain_Data* dp, int64_t max_s
int64_t max_iter, int64_t max_dis, double chn_pen_gap, double chn_pen_skip,
double bw_rate, int64_t xl, int64_t yl, int64_t quick_check,
int64_t left_fix, int64_t right_fix);
uint64_t lchain_qdp_fix_adv(k_mer_hit* a, int64_t a_n, Chain_Data* dp, 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 xl, int64_t yl, int64_t quick_check,
int64_t left_fix, int64_t right_fix, k_mer_hit *res);
uint64_t lchain_simple(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp,
int64_t max_skip, int64_t max_iter);
uint64_t lchain_simple0(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, int64_t max_skip, int64_t max_iter);
@@ -253,6 +257,12 @@ uint64_t lchain_qdp_mcopy_fast(Candidates_list *cl, int64_t a_idx, int64_t a_n,
int64_t gen_cigar, int64_t enable_mcopy, double mcopy_rate, int64_t mcopy_khit_cutoff,
int64_t khit_n);
uint64_t lchain_qdp_global_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 khit_n);
#define kv_pushp_ol(type, v, p) do { \
if ((v).length == (v).size) { \
(v).list = (type*)realloc((v).list, sizeof(type)*((v).size?((v).size<<1):(2))); \
Levenshtein_distance.h
+22 -11
View File
@@ -712,19 +712,30 @@ inline uint32_t pop_trace_bp_f(asg16_v *res, uint32_t i, uint16_t *c, uint16_t *
inline int64_t pop_trace_bp_rev_f(asg16_v *res, int64_t i, uint16_t *c, uint16_t *bq, uint16_t *bt, uint32_t *len)
{
(*c) = (res->a[i]>>14); (*bq) = (*bt) = (uint16_t)-1;
if((*c) == 2 || (*c) == 3) {
(*bt) = ((res->a[i]>>12)&3);
(*len) = (res->a[i]&(0xfff));
} else if((*c) == 1) {
(*bt) = ((res->a[i]>>12)&3);
(*bq) = ((res->a[i]>>10)&3);
(*len) = (res->a[i]&(0x3ff));
} else {
(*len) = (res->a[i]&(0x3fff));
uint32_t sl = 1;
(*c) = (*bq) = (*bt) = (uint16_t)-1; (*len) = 0;
if(i >= ((int64_t)res->n)) {
i = ((int64_t)res->n) - 1; sl = 0;
}
if(i >= 0) {
(*c) = (res->a[i]>>14);
(*bq) = (*bt) = (uint16_t)-1;
if((*c) == 2 || (*c) == 3) {
(*bt) = ((res->a[i]>>12)&3);
(*len) = (res->a[i]&(0xfff));
} else if((*c) == 1) {
(*bt) = ((res->a[i]>>12)&3);
(*bq) = ((res->a[i]>>10)&3);
(*len) = (res->a[i]&(0x3ff));
} else {
(*len) = (res->a[i]&(0x3fff));
}
}
if(sl == 0) {
(*len) = 0; i = res->n;
}
uint32_t sl; uint16_t sbq, sbt;
uint16_t sbq, sbt;
for (i--; (i >= 0) && ((*c) == (res->a[i]>>14)); i--) {
sbq = sbt = (uint16_t)-1;
if((*c) == 2 || (*c) == 3) {
Overlaps.cpp
+2 -2
View File
@@ -23747,7 +23747,7 @@ void write_all_data_to_disk(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sou
sprintf(gfa_name, "%s.ovlp.reverse", output_file_name);
write_ma_hit_ts(reverse_sources, RNF->total_reads, gfa_name);
if(asm_opt.dbg_bam) {
/**if(asm_opt.dbg_bam)**/ {
sprintf(gfa_name, "%s.rec", output_file_name);
write_cc_v(&scb, gfa_name);
}
@@ -23794,7 +23794,7 @@ int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_s
return 0;
}
if(asm_opt.dbg_bam) {
/**if(asm_opt.dbg_bam)**/ {
sprintf(gfa_name, "%s.rec", output_file_name);
load_cc_v(&scb, gfa_name); ///write_ec_reads("lec.raw.fa", &scb, 0);
}
Process_Read.cpp
+24 -7
View File
@@ -125,6 +125,11 @@ void write_All_reads(All_reads* r, char* read_file_name)
}
}
mm = 3;
fwrite(&mm, sizeof(mm), 1, fp);
fwrite(&(r->tr[0]), sizeof(r->tr[0]), 1, fp);
fwrite(&(r->tr[1]), sizeof(r->tr[1]), 1, fp);
free(index_name);
fflush(fp);
fclose(fp);
@@ -235,11 +240,18 @@ int load_All_reads(All_reads* r, char* read_file_name)
MALLOC(r->rsc[i], (r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0));
f_flag += fread(r->rsc[i], sizeof(uint8_t), (r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0), fp);
}
mm = 0;
if(!feof(fp)) fread(&mm, sizeof(mm), 1, fp);
}
if(mm == 3) {
fread(&(r->tr[0]), sizeof(r->tr[0]), 1, fp);
fread(&(r->tr[1]), sizeof(r->tr[1]), 1, fp);
}
}
free(index_name);
fclose(fp);
fprintf(stderr, "Reads has been loaded.\n");
@@ -283,12 +295,17 @@ uint8_t load_cc_v(cc_v* r, char* read_file_name)
///typedef struct {size_t n, m; asg16_v *a; uint8_t *f; uint16_t *er; uint64_t bid;} cc_v;
asg16_v *z; int f_flag = 0;
uint64_t k, rn; uint32_t zn; f_flag += fread(&rn, sizeof(rn), 1, fp);
r->n = r->m = rn; MALLOC(r->a, r->n);
for (k = 0; k < r->n; k++) {
z = &(r->a[k]);
f_flag += fread(&zn, sizeof(zn), 1, fp);
z->n = z->m = zn; MALLOC(z->a, z->n);
f_flag += fread(z->a, sizeof((*(z->a))), zn, fp);
r->n = r->m = rn;
if(rn == 0) {
free(r->a); r->a = NULL;
} else {
MALLOC(r->a, r->n);
for (k = 0; k < r->n; k++) {
z = &(r->a[k]);
f_flag += fread(&zn, sizeof(zn), 1, fp);
z->n = z->m = zn; MALLOC(z->a, z->n);
f_flag += fread(z->a, sizeof((*(z->a))), zn, fp);
}
}
fflush(fp); fclose(fp);
anchor.cpp
+150
View File
@@ -4014,6 +4014,156 @@ void get_pi_ec_chain(ha_abuf_t *ab, uint64_t rid, uint64_t rl, uint32_t tid, cha
}
uint64_t gen_srt_chain(ha_abuf_t *ab, ha_mz1_t *rfa, uint64_t *rfi, uint64_t rfn, ha_mz1_t *qra, uint64_t *qri, uint64_t qrn, uint32_t *high_occ, uint32_t *low_occ, Candidates_list *cl, uint32_t qrid)
{
uint64_t i, k, l, max_cnt = UINT32_MAX, min_cnt = 0, ri, qi, sn, x, rz, qz, n_a0, l0, k0, i0; anchor1_t *an; k_mer_hit *p;
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;
}
radix_sort_anc64(rfi, rfi + rfn);
radix_sort_anc64(qri, qri + qrn);
i = ab->n_a = 0; n_a0 = l0 = i0 = 0; k0 = 1;
for (k = 1, l = 0; k <= rfn; ++k) {
if (k == rfn || (rfi[k]>>32) != (rfi[l]>>32)) {
x = rfi[l]>>32;
for (; i < qrn && (qri[i]>>32) < x; i++);
if(ab->n_a < ab->m_a) {
n_a0 = ab->n_a; l0 = l; i0 = i; k0 = k;
}
if(i < qrn && (qri[i]>>32) == x) {
// sn = 0;
// if(ab->n_a < ab->m_a) sn = ab->seed[l].n;///ha_pt_cnt(ha_idx, ra[l].x);
for (qi = i; qi < qrn && (qri[qi]>>32) == x; qi++) {
qz = (uint32_t)qri[qi];
for (ri = l; ri < k; ri++) {
rz = (uint32_t)rfi[ri];
if(qra[qz].rev != rfa[rz].rev) continue;
if(ab->n_a < ab->m_a) {
an = &(ab->a[ab->n_a++]);
an->other_off = qra[qz].pos;
an->self_off = rfa[rz].pos;
///an->cnt: cnt<<8|span
an->cnt = ab->seed[rz].n; if(an->cnt > ((uint32_t)(0xffffffu))) an->cnt = 0xffffffu;
an->cnt <<= 8; an->cnt |= ((rfa[rz].span <= ((uint32_t)(0xffu)))?rfa[rz].span:((uint32_t)(0xffu)));
an->srt = (((uint64_t)(an->self_off))<<32)|((uint64_t)(an->other_off));
} else {
ab->n_a++;
}
}
}
i = qi;
}
l = k;
}
}
if (ab->n_a > ab->m_a) {
ab->m_a = ab->n_a; REALLOC(ab->a, ab->m_a);
k = k0; i = i0; l = l0; ab->n_a = n_a0;
for (; k <= rfn; ++k) {
if (k == rfn || (rfi[k]>>32) != (rfi[l]>>32)) {
x = rfi[l]>>32;
for (; i < qrn && (qri[i]>>32) < x; i++);
if(ab->n_a < ab->m_a) {
n_a0 = ab->n_a; l0 = l; i0 = i; k0 = k;
}
if(i < qrn && (qri[i]>>32) == x) {
// sn = 0;
// if(ab->n_a < ab->m_a) sn = ab->seed[l].n;///ha_pt_cnt(ha_idx, ra[l].x);
for (qi = i; qi < qrn && (qri[qi]>>32) == x; qi++) {
qz = (uint32_t)qri[qi];
for (ri = l; ri < k; ri++) {
rz = (uint32_t)rfi[ri];
if(qra[qz].rev != rfa[rz].rev) continue;
if(ab->n_a < ab->m_a) {
an = &(ab->a[ab->n_a++]);
an->other_off = qra[qz].pos;
an->self_off = rfa[rz].pos;
///an->cnt: cnt<<8|span
an->cnt = ab->seed[rz].n; if(an->cnt > ((uint32_t)(0xffffffu))) an->cnt = 0xffffffu;
an->cnt <<= 8; an->cnt |= ((rfa[rz].span <= ((uint32_t)(0xffu)))?rfa[rz].span:((uint32_t)(0xffu)));
an->srt = (((uint64_t)(an->self_off))<<32)|((uint64_t)(an->other_off));
} else {
ab->n_a++;
}
}
}
i = qi;
}
l = k;
}
}
}
// copy over to _cl_
sn = ab->n_a + cl->length;
if (sn > (uint64_t)cl->size) {
cl->size = sn;
REALLOC(cl->list, cl->size);
}
radix_sort_ha_an1(ab->a, ab->a + ab->n_a);
for (i = 0; i < ab->n_a; i++) {
p = &cl->list[cl->length++];
p->readID = qrid;
p->strand = 0;
p->offset = ab->a[i].other_off;
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));
}
// cl->length = ab->n_a;
return ab->n_a;
}
void gen_self_global_chain(ha_abuf_t *ab, Candidates_list *cl, uint32_t rid, uint64_t rl, uint32_t tid, char *ts, uint64_t tl, uint64_t mz_w, uint64_t mz_k, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, asg64_v *ix,
uint8_t is_accurate, double bw_thres, int apend_be, uint32_t gen_off, int64_t enable_mcopy, double mcopy_rate, uint32_t mcopy_khit_cut, overlap_region_alloc *ores)
{
extern void *ha_flt_tab;
uint64_t k, rn = ab->mz.n, m = cl->length;
mz1_ha_sketch(ts, tl, 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);
kv_resize(uint64_t, *ix, ab->mz.n);
for (k = 0; k < rn; k++) {
ix->a[k] = ab->mz.a[k].x;
ix->a[k] <<= 32; ix->a[k] |= k;
}
for (; k < ab->mz.n; k++) {
ix->a[k] = ab->mz.a[k].x;
ix->a[k] <<= 32; ix->a[k] |= (k-rn);
}
ix->n = ab->mz.n;
if(gen_srt_chain(ab, ab->mz.a, ix->a, rn, ab->mz.a + rn, ix->a + rn, ab->mz.n - rn, high_occ, low_occ, cl, tid)) {
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);
m += lchain_qdp_mcopy_fast(cl, m, cl->length - m, m, &(cl->chainDP), ores, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres,
rid, rl, tl, quick_check, apend_be, gen_off, enable_mcopy, mcopy_rate, mcopy_khit_cut, 1);
cl->length = m;
}
ab->mz.n = rn;
}
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
+445 -50
View File
@@ -13,6 +13,7 @@
#define del_cns_arc(z, arc_i) ((z).arc.a[(arc_i)].v == CNS_DEL_E)
#define CNS_DEL_V (0x1fffffffu)
#define del_cns_nn(z, nn_i) ((z).a[(nn_i)].sc == CNS_DEL_V)
#define is_cns_bb(z, nn_i) ((nn_i) >= (z).bb0 && (nn_i) < (z).bb1)
#define REFRESH_N 128
#define COV_W 3072
#define COV_W_AC 512
@@ -242,6 +243,8 @@ uint64_t get_mz1(const char *str, int len, int w, int k, uint32_t rid, int is_hp
void get_pi_ec_chain(ha_abuf_t *ab, uint64_t rid, uint64_t rl, uint32_t tid, char* ts, uint64_t tl, uint64_t mz_w, uint64_t mz_k, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres,
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 mcopy_khit_cut,
int64_t max_skip, int64_t max_iter, int64_t max_dis, int64_t quick_check, double chn_pen_gap, double chn_pen_skip);
void gen_self_global_chain(ha_abuf_t *ab, Candidates_list *cl, uint32_t rid, uint64_t rl, uint32_t tid, char *ts, uint64_t tl, uint64_t mz_w, uint64_t mz_k, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, asg64_v *ix,
uint8_t is_accurate, double bw_thres, int apend_be, uint32_t gen_off, int64_t enable_mcopy, double mcopy_rate, uint32_t mcopy_khit_cut, overlap_region_alloc *ores);
void set_lchain_dp_op(uint32_t is_accurate, uint32_t mz_k, int64_t *max_skip, int64_t *max_iter, int64_t *max_dis, double *chn_pen_gap, double *chn_pen_skip, int64_t *quick_check);
void h_ec_lchain_re_gen_srt(ha_abuf_t *ab, ha_pt_t *ha_idx, overlap_region_alloc *olst, Candidates_list *cl);
uint64_t h_ec_lchain_re_gen_qry(ha_abuf_t *ab, uint64_t *k, uint64_t *l, uint64_t *i, uint64_t *idx_a, uint64_t idx_n, uint64_t *tid, uint64_t *trev);
@@ -1430,7 +1433,174 @@ void del_cns_g_nn(cns_gfa *cns, uint32_t v)
cns->a[v].arc.n = cns->a[v].arc.nou = 0;
cns->a[v].c = cns->a[v].f = 0; cns->a[v].sc = CNS_DEL_V;
/**cns->a[v].c =**/ cns->a[v].f = 0; cns->a[v].sc = CNS_DEL_V;
}
void merge_cns_g_in_adv(cns_gfa *cns, uint32_t v0, asg32_v* b32)
{
cns_t *av, *aw; uint32_t v, bp, vk, wk, wka, w, wn, nn, mn, wh, mn_k[2];
b32->n = 0;
kv_push(uint32_t, *b32, v0);
while (b32->n) {
v = b32->a[--b32->n];
if(del_cns_nn((*cns), v)) continue;
av = &((*cns).a[v]);
for (bp = 0; bp < 4; bp++) {
//nn: number of node; wh: weight
nn = wh = 0; mn = mn_k[0] = mn_k[1] = wka = (uint32_t)-1;
for (vk = av->arc.nou; vk < av->arc.n; vk++) {///in-edge of v
if(del_cns_arc((*av), vk)) continue;
w = av->arc.a[vk].v; aw = &((*cns).a[w]);
if(aw->c != bp) continue;
if(w == cns->si || w == cns->ei) continue;
for (wk = wn = 0; wk < aw->arc.nou; wk++) {///out-edge of w
if(del_cns_arc((*aw), wk)) continue;
wn++; wka = wk; if(wn > 1) break;
}
if(wn != 1) continue;
assert(aw->arc.a[wka].v == v);
///deal with out-edge of w
if((nn == 0) || is_cns_bb((*cns), w)) {///this is still risky as there might be multiple backbone
mn = w; mn_k[0] = vk; mn_k[1] = wka;
///not sure if we should set these edges as visited
// av->arc.a[vk].f = 1; aw->arc.a[wka].f = 1;
}
wh += aw->arc.a[wka].sc;
nn++;
}
if(nn > 1) {
for (vk = av->arc.nou; vk < av->arc.n; vk++) {///in-edge of v
if(del_cns_arc((*av), vk)) continue;
w = av->arc.a[vk].v; aw = &((*cns).a[w]);
if(aw->c != bp) continue;
if(w == cns->si || w == cns->ei) continue;
for (wk = wn = 0; wk < aw->arc.nou; wk++) {///out-edge of w
if(del_cns_arc((*aw), wk)) continue;
wn++; wka = wk; if(wn > 1) break;
}
if(wn != 1) continue;
assert(aw->arc.a[wka].v == v);
///deal with in-edge of w
///all edges to w, should be move to mn
if(mn != w) {///not sure if we should set these edges as visited; affect when nn == 0
for (wk = aw->arc.nou; wk < aw->arc.n; wk++) {
if(del_cns_arc((*aw), wk)) continue;
///previously, aw->arc.a[wk].v -> w
///currently, aw->arc.a[wk].v -> mn
/// if(nn == 0), then mn = w
gen_mm_cns_arc(cns, aw->arc.a[wk].v, mn, aw->arc.a[wk].sc/**(nn?(aw->arc.a[wk].sc):(0))**/, aw->arc.a[wk].f);///not sure if we should set these edges as visited
}
del_cns_g_nn(cns, w);
}
}
}
if(nn) {
aw = &((*cns).a[mn]);
av->arc.a[mn_k[0]].sc = aw->arc.a[mn_k[1]].sc = wh;
// merge_cns_g_in(cns_gfa *cns, uint32_t v, asg32_v* b32)
kv_push(uint32_t, *b32, mn);
}
}
}
}
void merge_cns_g_ou_adv(cns_gfa *cns, uint32_t v0, asg32_v* b32)
{
cns_t *av, *aw; uint32_t v, bp, vk, wk, wka, w, wn, nn, mn, wh, mn_k[2];
b32->n = 0;
kv_push(uint32_t, *b32, v0);
while (b32->n) {
v = b32->a[--b32->n];
if(del_cns_nn((*cns), v)) continue;
av = &((*cns).a[v]);
for (bp = 0; bp < 4; bp++) {
//nn: number of node; wh: weight
nn = wh = 0; mn = mn_k[0] = mn_k[1] = wka = (uint32_t)-1;
for (vk = 0; vk < av->arc.nou; vk++) {///ou-edge of v
if(del_cns_arc((*av), vk)) continue;
w = av->arc.a[vk].v; aw = &((*cns).a[w]);
if(aw->c != bp) continue;
if(w == cns->si || w == cns->ei) continue;
for (wk = aw->arc.nou, wn = 0; wk < aw->arc.n; wk++) {///in-edge of w
if(del_cns_arc((*aw), wk)) continue;
wn++; wka = wk; if(wn > 1) break;
}
if(wn != 1) continue;
assert(aw->arc.a[wka].v == v);
///deal with out-edge of w
if((nn == 0) || is_cns_bb((*cns), w)) {///this is still risky as there might be multiple backbone
mn = w; mn_k[0] = vk; mn_k[1] = wka;
///not sure if we should set these edges as visited
// av->arc.a[vk].f = 1; aw->arc.a[wka].f = 1;
}
wh += aw->arc.a[wka].sc;
nn++;
}
if(nn > 1) {
for (vk = 0; vk < av->arc.nou; vk++) {///ou-edge of v
if(del_cns_arc((*av), vk)) continue;
w = av->arc.a[vk].v; aw = &((*cns).a[w]);
if(aw->c != bp) continue;
if(w == cns->si || w == cns->ei) continue;
for (wk = aw->arc.nou, wn = 0; wk < aw->arc.n; wk++) {///in-edge of w
if(del_cns_arc((*aw), wk)) continue;
wn++; wka = wk; if(wn > 1) break;
}
if(wn != 1) continue;
assert(aw->arc.a[wka].v == v);
///deal with in-edge of w
///all edges to w, should be move to mn
if(mn != w) {///not sure if we should set these edges as visited; affect when nn == 0
for (wk = 0; wk < aw->arc.nou; wk++) {
if(del_cns_arc((*aw), wk)) continue;
///previously, w -> aw->arc.a[wk].v
///currently, mn -> aw->arc.a[wk].v
/// if(nn == 0), then mn = w
gen_mm_cns_arc(cns, mn, aw->arc.a[wk].v, aw->arc.a[wk].sc/**(nn?(aw->arc.a[wk].sc):(0))**/, aw->arc.a[wk].f);///not sure if we should set these edges as visited
}
del_cns_g_nn(cns, w);
}
}
}
if(nn) {
aw = &((*cns).a[mn]);
// fprintf(stderr, "\n[M::%s] nn::%u, mn::%u\n", __func__, nn, mn);
// fprintf(stderr, "[M::%s] vi::%u, vn::%u\n", __func__, mn_k[0], (uint32_t)av->arc.n);
// fprintf(stderr, "[M::%s] wi::%u, wn::%u\n", __func__, mn_k[1], (uint32_t)aw->arc.n);
av->arc.a[mn_k[0]].sc = aw->arc.a[mn_k[1]].sc = wh;
// merge_cns_g_in(cns_gfa *cns, uint32_t v, asg32_v* b32)
kv_push(uint32_t, *b32, mn);
}
}
}
}
void merge_cns_g_in(cns_gfa *cns, uint32_t v0, asg32_v* b32)
@@ -1936,6 +2106,9 @@ uint64_t push_correct1_fhc_indel_exz(asg16_v *sc, int64_t sc0, window_list *idx,
uint64_t push_correct1_fhc(window_list *idx, window_list_alloc *res, cns_gfa *cns, char* qstr, UC_Read* tu, bit_extz_t *exz, asg32_v *rc, uint32_t bl, uint32_t rid)
{
// if(rid == 11206 && bl == 6) {
// fprintf(stderr, "[M::%s]\trc->n::%u\tbl::%u\trid::%u\n", __func__, (uint32_t)rc->n, bl, rid);
// }
// fprintf(stderr, "[M::%s]\trc->n::%u\tbl::%u\n", __func__, (uint32_t)rc->n, bl);
uint64_t nec = 0; uint32_t k, l, i, ff, sl, sk, bs = cns->off, be = bl + cns->off, bend = cns->off, is_i = 0, sc0 = res->c.n;///[bs, be)
if(rc->n) {///it is possible that rc->n == 0, which means there is a deletion
@@ -1983,13 +2156,17 @@ uint64_t push_correct1_fhc(window_list *idx, window_list_alloc *res, cns_gfa *cn
l = k;
}
}
// if(rid == 11206 && bl == 6) {
// fprintf(stderr, "[M::%s]\trc->n::%u\tbl::%u\trid::%u\tbend::%u\tbe::%u\n", __func__, (uint32_t)rc->n, bl, rid, bend, be);
// }
///push remaining deletion
if(be > bend) {
if(is_i && exz) {
nec += push_correct1_fhc_indel_exz(((asg16_v *)(&(res->c))), sc0, idx, cns, qstr, tu, exz, cns->off, 3, be - bend, bend-cns->off);
} else {
for (i = bend; i < be; i++) {
// fprintf(stderr, "%c(%u)\n", s_H[cns->a[i].c], i);
push_trace_bp_f(((asg16_v *)(&(res->c))), 3, cns->a[i].c, (uint16_t)-1, 1, ((idx->clen>0)?1:0));
idx->clen = res->c.n - idx->cidx; nec++;
}
@@ -2017,6 +2194,23 @@ uint64_t cns_gen_full0(overlap_region* ol, All_reads *rref, uint64_t s, uint64_t
}
init_cns_g(cns, qstr + s, e - s, hw, rid);
// uint64_t dk = 0;
// if(s <= 29788 && e > 29788) {
// fprintf(stderr, "[M::%s]\tqid::%u\tq::[%lu, %lu)\n", __func__, rid, s, e);
// fprintf(stderr, "-z-[M::%s]\toriginal::", __func__);
// for (dk = s; dk < e; dk++) {
// fprintf(stderr, "%c", qstr[dk]);
// }
// fprintf(stderr, "\n");
// fprintf(stderr, "-a-[M::%s]\tGraph::\t\t\t", __func__);
// for (dk = 2; dk < (*cns).n; dk++) {
// fprintf(stderr, "%c", s_H[(*cns).a[dk].c]);
// }
// fprintf(stderr, "\n");
// }
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, q[2], os, oe; ul_ov_t *p; overlap_region *z; idx->rr = 0;
@@ -2049,16 +2243,42 @@ uint64_t cns_gen_full0(overlap_region* ol, All_reads *rref, uint64_t s, uint64_t
}
}
// if(s <= 29788 && e > 29788) {
// fprintf(stderr, "-b-[M::%s]\tGraph::\t\t\t", __func__);
// for (dk = 2; dk < (*cns).n; dk++) {
// fprintf(stderr, "%c(del::%u)", s_H[(*cns).a[dk].c], del_cns_nn((*cns), dk));
// }
// fprintf(stderr, "\n");
// }
// fprintf(stderr, "-2-[M::%s] cns->n::%u\n", __func__, (uint32_t)cns->n);
// return;
refine_cns_g(cns, b32);
// if(s <= 29788 && e > 29788) {
// fprintf(stderr, "-c-[M::%s]\tGraph::\t\t\t", __func__);
// for (dk = 2; dk < (*cns).n; dk++) {
// fprintf(stderr, "%c(del::%u)", s_H[(*cns).a[dk].c], del_cns_nn((*cns), dk));
// }
// fprintf(stderr, "\n");
// }
// fprintf(stderr, "-3-[M::%s] cns->n::%u\n", __func__, (uint32_t)cns->n);
gseq_cns_g(cns, b32, e - s);
// if(s <= 29788 && e > 29788) {
// fprintf(stderr, "-d-[M::%s]\tGraph::\t\t\t", __func__);
// for (dk = 2; dk < (*cns).n; dk++) {
// fprintf(stderr, "%c(del::%u)", s_H[(*cns).a[dk].c], del_cns_nn((*cns), dk));
// }
// fprintf(stderr, "\n");
// }
// fprintf(stderr, "-4-[M::%s] cns->n::%u\n", __func__, (uint32_t)cns->n);
// nec += push_correct1(ridx, res, cns, b32, e - s);
@@ -2392,7 +2612,7 @@ uint64_t wcns_vote(overlap_region* ol, All_reads *rref, uint64_t q_hf, char* qst
}
///a) pass coverage check; b) no enough coverage
if((((oc[0] > (oc[1]*occ_exact)) && (oc[0] > (oc[1]-oc[0])) && (oc[1] >= occ_tot) && (oc[0] > 1) && ((!hf_idx) || ((ow[0] > (ow[1]*occ_exact)) && (ow[0] > (ow[1] - ow[0])))))) || (oc[1] < occ_tot)) {
fI = 0;
fI = 0;///keep q itself
}
if(fI) {
@@ -2434,8 +2654,6 @@ uint64_t wcns_vote(overlap_region* ol, All_reads *rref, uint64_t q_hf, char* qst
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];
@@ -2533,7 +2751,7 @@ uint64_t wcns_gen(overlap_region_alloc* ol, All_reads *rref, uint64_t qid, UC_Re
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++) {
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++) {
@@ -2555,10 +2773,11 @@ uint64_t wcns_gen(overlap_region_alloc* ol, All_reads *rref, uint64_t qid, UC_Re
///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;
idx_a = c_idx->a; idx_b = c_idx->a + c_idx->n;///update here?
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);
///buf->a[wl<<1, idx->n): this is for sorted index
memcpy(buf->a + (wl<<1), idx->a, idx->n * (sizeof((*(idx->a)))));
memset(buf->a, 0, (wl<<1)*(sizeof((*(idx->a)))));
if(hf_idx) {
@@ -4340,15 +4559,157 @@ static void worker_init_ec_step(void *data, long i, int tid)
void update_scb(All_reads *R_INF, asg16_v *scc, asg16_v *scb, asg16_v *scb_res, UC_Read *qu, UC_Read *tu, asg64_v *srt, bit_extz_t *exz, uint64_t rid);
uint64_t gen_ori_seq0(char *tstr, uint64_t tl, UC_Read *qu, asg16_v *sc, uint64_t rid);
uint8_t refresh_check_scc(overlap_region *z, int64_t sc_e, int64_t ql, int64_t tl, int64_t gap_bd, double gap_rate, double err_rate, int64_t err_diff_bd, double err_diff_rate)
{
int64_t k, wn = z->w_list.n, tot_g = 0, tot_e = 0, wq, wt;
if(wn <= 0) return 0;
tot_g += z->w_list.a[0].x_start; tot_g += ql - z->w_list.a[wn-1].x_end - 1;
tot_g += z->w_list.a[0].y_start; tot_g += tl - z->w_list.a[wn-1].y_end - 1;
for (k = 0; k < wn; k++) {
if(is_ualn_win(z->w_list.a[k])) {
if(z->w_list.a[k].extra_begin == -1 && z->w_list.a[k].extra_end == -1) {
return 0;
}
tot_e += (z->w_list.a[k].extra_begin*(-1)) - 1;
// fprintf(stderr, "+[M::%s]\tq::[%u,%u)\tt::[%u,%u)\ttot_e::%ld\n", __func__, z->w_list.a[k].x_start, z->w_list.a[k].x_end + 1, z->w_list.a[k].y_start, z->w_list.a[k].y_end + 1, tot_e);
} else {
tot_e += z->w_list.a[k].error;
// fprintf(stderr, "-[M::%s]\tq::[%u,%u)\tt::[%u,%u)\ttot_e::%ld\n", __func__, z->w_list.a[k].x_start, z->w_list.a[k].x_end + 1, z->w_list.a[k].y_start, z->w_list.a[k].y_end + 1, tot_e);
}
}
// fprintf(stderr, "[M::%s]\tq::[%u,%u)\tt::[%u,%u)\ttot_e::%ld\trtot_g::%ld\n", __func__, z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1, tot_e, tot_g);
if((tot_g > 0) && (tot_g > gap_bd)) {
if(tot_g > (ql*gap_rate)) return 0;
if(tot_g > (tl*gap_rate)) return 0;
}
wq = (z->w_list.a[wn-1].x_end+1) - z->w_list.a[0].x_start;
wt = (z->w_list.a[wn-1].y_end+1) - z->w_list.a[0].y_start;
if((tot_e > 0) && ((tot_e > (wq*err_rate)) || (tot_e > (wt*err_rate)))) return 0;
// tot_e += tot_g;
if((tot_e > 0) && ((tot_e > (ql*err_rate)) || (tot_e > (tl*err_rate)))) return 0;
if((tot_e > sc_e) && ((tot_e - sc_e) > err_diff_bd) && ((tot_e - sc_e) > (sc_e*err_diff_rate))) return 0;
return 1;
}
void dbg_prt_asg16_v_sc(uint64_t rid, asg16_v *sc)
{
uint64_t ck, qk, tk, tot_e = 0; uint32_t len; uint16_t c, bq, bt;
fprintf(stderr, "[M::%s]\trid::%lu\trlen::%lu\n", __func__, rid, Get_READ_LENGTH(R_INF, rid));
ck = qk = tk = 0;
while (ck < sc->n) {
ck = pop_trace_bp_f(sc, ck, &c, &bq, &bt, &len);
if(c != 2) qk += len;
if(c != 3) tk += len;
if(c!=0) tot_e += len;
// fprintf(stderr, "%u(%c)\t", len, "MSID"[c]);
}
// fprintf(stderr, "\n");
fprintf(stderr, "[M::%s]\trid::%lu\t#\talter::%lu\n", __func__, rid, tot_e);
}
uint8_t regen_scb(ha_abuf_t *ab, Candidates_list *cl, uint32_t rid, asg16_v *sc, UC_Read *ia, UC_Read *ob, asg64_v *srt,
uint64_t mz_w, uint64_t mz_k, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ,
uint8_t is_accurate, double bw_thres, int apend_be, uint32_t gen_off, overlap_region_alloc *ol, bit_extz_t *exz, double erate, uint64_t wl, asg16_v *buf,
overlap_region **r_aux_o, overlap_region **rchn, int64_t *rref_len)
{
// return 0;
// dbg_prt_asg16_v_sc(rid, sc);
uint64_t oln0 = ol->length, e0; (*rchn) = NULL; (*rref_len) = -1;
e0 = gen_ori_seq0(ia->seq, ia->length, ob, sc, rid);
// fprintf(stderr, "[M::%s]\tinitial\tlength::%ld\n", __func__, (int64_t)ob->length);
// fprintf(stderr, "\n[M::%s]\trid::%u\trlen::%ld\te0::%ld\n", __func__, rid, (int64_t)ia->length, e0);
if(e0 == 0) return 0;
cl->length = 0;
gen_self_global_chain(ab, cl, rid, ia->length, rid/**((uint32_t)-1)**/, ob->seq, ob->length, mz_w, mz_k, k_flag, dbg_ct, sp, high_occ, low_occ, srt, is_accurate, bw_thres, apend_be, gen_off, 1, -1, UINT32_MAX, ol);
(*r_aux_o) = fetch_aux_ovlp(ol, NULL);///refresh
if(ol->length > oln0) {
assert(ol->length == oln0 + 1);
if(gen_hc_r_alin_self(&(ol->list[oln0]), cl, ia->seq, ia->length, ob->seq, ob->length, exz, (*r_aux_o), erate, wl, rid, E_KHIT, 1, buf, sc)) {
if(refresh_check_scc(&(ol->list[oln0]), e0, ia->length, ob->length, 32, 0.012, 0.1, 32, 0.2)) {
(*rchn) = &(ol->list[oln0]); (*rref_len) = ob->length;
resize_UC_Read(ia, ia->length + ob->length); memcpy(ia->seq + ia->length, ob->seq, sizeof((*(ob->seq)))*ob->length);
ol->length = oln0;
return 1;
}
}
}
ol->length = oln0;
return 0;
}
void cmp_smp_ac(UC_Read *ref, asg16_v *scz, UC_Read *res, uint64_t rid)
{
uint64_t qn = ref->length; uint16_t c, bq, bt; uint32_t len, ck, qk, tk, tn, wq[2], wt[2], k/**, Nn = 0**/; char *qsr = NULL, *tsr = NULL;
ck = qk = tk = 0;
while (ck < scz->n) {
ck = pop_trace_bp_f(scz, ck, &c, &bq, &bt, &len);
if(c != 3) tk += len;
}
tn = tk; resize_UC_Read(ref, qn + tn);
qsr = ref->seq; tsr = ref->seq + qn;
ck = qk = tk = 0;
while (ck < scz->n) {
wq[0] = qk; wt[0] = tk;
ck = pop_trace_bp_f(scz, ck, &c, &bq, &bt, &len);
if(c != 2) qk += len;
if(c != 3) tk += len;
wq[1] = qk; wt[1] = tk;
// fprintf(stderr, "%u(%c)\tq::[%u,%u)\tbq::%u\tt::[%u,%u)\tbt::%u\n", len, "MSID"[c], wq[0], wq[1], bq, wt[0], wt[1], bt);
if(c == 0) {
for (; wq[0] < wq[1]; wq[0]++, wt[0]++) {
tsr[wt[0]] = qsr[wq[0]];
// if(p->a[wy[0]] == 'N') Nn++;
}
} else if(c == 1 || c == 2) {
for (k = wt[0]; k < wt[1]; k++) {
tsr[k] = s_H[bt];
// if(p->a[k] == 'N') Nn++;
}
}
}
gen_ori_seq0(tsr, tn, res, scz, rid);
assert(res->length == ((int64_t)qn));
if(memcmp(qsr, res->seq, qn) != 0) {
fprintf(stderr, "[M::%s]\trid::%lu(%.*s)\tres->length::%ld\tqn::%lu\n", __func__, rid, (int)Get_NAME_LENGTH(R_INF, rid), Get_NAME(R_INF, rid), (int64_t)res->length, qn);
for (k = 0; k < qn; k++) {
if(res->seq[k] != qsr[k]) {
fprintf(stderr, "[M::%s]\tk::%u\tqsr[%u]::%c\tres->seq[%u]::%c\n", __func__, k, k, qsr[k], k, res->seq[k]);
}
}
}
// assert(memcmp(qsr, res->seq, qn) == 0);
}
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; int64_t het_a, hom_a; ///gen_hc_aln_t ez;
overlap_region *aux_o = NULL/**, *rse_o = NULL**/; asg64_v buf0; uint32_t qlen = 0, qw = 0; uint64_t tot_b = 0; double tt0 = 0, tt1 = 0;
uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; int64_t het_a, hom_a, rl0 = -1; ///gen_hc_aln_t ez;
overlap_region *aux_o = NULL/**, *rse_o = NULL**/, *rcc = NULL; asg64_v buf0; uint32_t qlen = 0, qw = 0; uint64_t tot_b = 0; double tt0 = 0, tt1 = 0;
b->v8q.n = b->v8t.n = 0; set_ec_cov(asm_opt.het_cov, asm_opt.hom_cov, asm_opt.het_cov_set, asm_opt.polyploidy, het_a, hom_a);
// if((i != 733166) && (i != 858708) && (i != 858732) && (i != 859819) && (i != 859899) && (i != 863486) && (i != 872165) && (i != 899887) && (i != 902298) &&
// (i != 906808) && (i != 946173) && (i != 952685) && (i != 983977) && (i != 1000227) && (i != 1011228) && (i != 1042858) && (i != 1045860) && (i != 1118558) &&
// (i != 1143886) && (i != 1155956) && (i != 1159490) && (i != 1179151) && (i != 1180199) && (i != 1230524) && (i != 1232338) && (i != 1244031) && (i != 1268467) &&
@@ -4502,11 +4863,15 @@ static void worker_hap_ec(void *data, long i, int tid)
gen_reseed_re(&b->olist, &b->clist, aux_o, rse_o, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, asm_opt.mz_win, 19, i, asm_opt.max_ov_diff_ec, asm_opt.max_ov_diff_ec, &b->v16, R_INF.tqn, b->v8q.a);
copy_asg_arr(b->sp, buf0);
**/
if(scb.a[i].n) {
regen_scb(b->ab, &b->clist, i, &(scb.a[i]), &b->self_read, &b->ovlp_read, &b->v64, asm_opt.mz_win, asm_opt.k_mer_length, NULL, NULL, &(b->sp), &high_occ, &low_occ,
1, ((asm_opt.is_ont)?(0.05):(0.02)), 1, 1, &b->olist, &b->exz, ((asm_opt.max_ov_diff_ec>0.1)?(asm_opt.max_ov_diff_ec):(0.1)), (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), &b->v16, &aux_o, &rcc, &rl0);
}
copy_asg_arr(buf0, b->sp);
//site_sc: r765 -> r766: 1 -> 0
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/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), ((asm_opt.is_sc)?&(b->v8t):NULL)/**&(b->v8t)**/, (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32,
asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1.0);
rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &(scb.a[i]), &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), ((asm_opt.is_sc)?&(b->v8t):NULL)/**&(b->v8t)**/, (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32,
asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1.0, rcc, rl0);
copy_asg_arr(b->sp, buf0);
///for debug indel
// stderr_phase_ovlp(&b->olist);
@@ -4540,7 +4905,7 @@ static void worker_hap_ec(void *data, long i, int tid)
push_nec_re(aux_o, &(scc.a[i]));
// push_nec_re(aux_o, &(scb.a[i]));
if(asm_opt.dbg_bam) {
/**if(asm_opt.dbg_bam)**/ {
update_scb(&R_INF, &(scc.a[i]), &(scb.a[i]), &(b->v16), &b->self_read, &b->ovlp_read, &b->v64, &b->exz, i);
kv_resize(uint16_t, scb.a[i], b->v16.n); scb.a[i].n = b->v16.n; memcpy(scb.a[i].a, b->v16.a, b->v16.n*sizeof(*(scb.a[i].a)));
}
@@ -4630,7 +4995,6 @@ static void worker_hap_ec(void *data, long i, int tid)
//fprintf(stderr, "-[M::%s]\trid::%ld\t%.*s\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
}
void gen_ori_seq0(char *tstr, uint64_t tl, UC_Read *qu, asg16_v *sc, uint64_t rid);
static void worker_gfa_ec(void *data, long i, int tid)
{
@@ -4699,7 +5063,6 @@ static void worker_gfa_ec(void *data, long i, int tid)
refresh_gc_ovec_buf_t0(b, REFRESH_N);
}
void worker_hap_ec_back_dbg(void *data, long i, int tid)
{
ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]);
@@ -4817,8 +5180,8 @@ void worker_hap_ec_back_dbg(void *data, long i, int tid)
copy_asg_arr(buf0, b->sp);
//site_sc: r765 -> r766: 1 -> 0
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/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), /**((asm_opt.is_sc)?&(b->v8t):NULL)**/&(b->v8t), (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32,
asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1.0);
rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, NULL, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), /**((asm_opt.is_sc)?&(b->v8t):NULL)**/&(b->v8t), (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32,
asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1.0, NULL, -1);
copy_asg_arr(b->sp, buf0);
///for debug indel
// stderr_phase_ovlp(&b->olist);
@@ -4940,8 +5303,6 @@ void worker_hap_ec_back_dbg(void *data, long i, int tid)
//fprintf(stderr, "-[M::%s]\trid::%ld\t%.*s\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
}
static void worker_hap_ec_step(void *data, long i, int tid)
{
ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]); i += scc.bid;
@@ -5085,8 +5446,8 @@ static void worker_hap_ec_step(void *data, long i, int tid)
copy_asg_arr(buf0, b->sp);
//site_sc: r765 -> r766: 1 -> 0
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/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), /**((asm_opt.is_sc)?&(b->v8t):NULL)**/&(b->v8t), (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32,
asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1.0);
rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, NULL, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), /**((asm_opt.is_sc)?&(b->v8t):NULL)**/&(b->v8t), (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32,
asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1.0, NULL, -1);
copy_asg_arr(b->sp, buf0);
///for debug indel
// stderr_phase_ovlp(&b->olist);
@@ -5207,9 +5568,6 @@ static void worker_hap_ec_step(void *data, long i, int tid)
//fprintf(stderr, "-[M::%s]\trid::%ld\t%.*s\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
}
static void worker_hap_ec_ss(void *data, long i, int tid)
{
ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]);
@@ -5303,8 +5661,8 @@ static void worker_hap_ec_ss(void *data, long i, int tid)
copy_asg_arr(buf0, b->sp);
//site_sc: r765 -> r766: 1 -> 0
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/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), /**((asm_opt.is_sc)?&(b->v8t):NULL)**/&(b->v8t), (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32,
asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1.0);
rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, NULL, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), /**((asm_opt.is_sc)?&(b->v8t):NULL)**/&(b->v8t), (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32,
asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1.0, NULL, -1);
copy_asg_arr(b->sp, buf0);
///for debug indel
// stderr_phase_ovlp(&b->olist);
@@ -5335,13 +5693,12 @@ static void worker_hap_ec_ss(void *data, long i, int tid)
}
static void worker_hap_ec_hybrid(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); int64_t het_a, hom_a;
uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; double bw_h, bw_l, e_h, e_l;
gen_hc_aln_t ez; overlap_region *aux_o = NULL, *rse_o = NULL; asg64_v buf0, buf1; uint64_t qlen = 0, qw = 0, qid = i; //uint64_t sk[2], ek[2], fn, qid = i, nec;
uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; double bw_h, bw_l, e_h, e_l; int64_t rl0 = -1;
gen_hc_aln_t ez; overlap_region *aux_o = NULL, *rse_o = NULL, *rcc = NULL; asg64_v buf0, buf1; uint64_t qlen = 0, qw = 0, qid = i; //uint64_t sk[2], ek[2], fn, qid = i, nec;
if(qid < R_INF.tqn) {///ont
bw_h = 0.05; bw_l = 0.035; e_h = asm_opt.max_ov_diff_ec; e_l = (asm_opt.max_ov_diff_ec + asm_opt.max_ov_diff_ec_sec)/2;
} else { ///HiFi
@@ -5353,16 +5710,20 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid)
}
b->v8q.n = b->v8t.n = 0; set_ec_cov(asm_opt.het_cov, asm_opt.hom_cov, asm_opt.het_cov_set, asm_opt.polyploidy, het_a, hom_a);
// if(i != 10) return;
// if((i%16) != 0) return;
// if(i != 5966) return;
// if(i != 11206) return;
// fprintf(stderr, "-a-[M::%s] rid::%ld\n", __func__, i);
//id:i:21102
// if (memcmp("a59fab4a-892b-4ab7-bf4b-926bed57865b_1", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
//id:i:3504
// if (memcmp("485f7963-eeb4-4745-ab74-1be4d61460c3", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
// fprintf(stderr, "-a-[M::%s-beg] rid->%ld, rlen->%lu\n", __func__, i, Get_READ_LENGTH((R_INF),i));
// if (memcmp("c7ecbd6b-e09d-4042-93ac-2400839feaf6", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
// fprintf(stderr, "-a-[M::%s-beg] rid->%ld, rlen->%lu, scb.a[i].n::%u\n", __func__, i, Get_READ_LENGTH((R_INF),i), (uint32_t)scb.a[i].n);
// } else {
// return;
// }
@@ -5434,12 +5795,24 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid)
///for debug indel
// prt_ovlp_sam(&b->olist, &b->ovlp_read, b->self_read.seq, b->self_read.length);
if(scb.a[i].n) {
regen_scb(b->ab, &b->clist, i, &(scb.a[i]), &b->self_read, &b->ovlp_read, &b->v64, asm_opt.mz_win, asm_opt.k_mer_length, NULL, NULL, &(b->sp), &high_occ, &low_occ,
1, bw_h, 1, 1, &b->olist, &b->exz, ((e_h>0.1)?(e_h):(0.1)), (qid < R_INF.tqn)?(WINDOW_OHC):(WINDOW_HC), &b->v16, &aux_o, &rcc, &rl0);
// if(rcc) {
// fprintf(stderr, "[M::%s]\t%.*s(qid::%u)\tql::%ld\tq::[%u,\t%u)\t%c\t%.*s(tid::%u)\ttl::%ld\tt::[%u,\t%u)\terr::%u\n", __func__,
// (int32_t)Get_NAME_LENGTH(R_INF, rcc->x_id), Get_NAME(R_INF, rcc->x_id), rcc->x_id, (int64_t)b->self_read.length, rcc->x_pos_s, rcc->x_pos_e + 1, "+-"[rcc->y_pos_strand],
// (int32_t)Get_NAME_LENGTH(R_INF, rcc->y_id), Get_NAME(R_INF, rcc->y_id), rcc->y_id, rl0, rcc->y_pos_s, rcc->y_pos_e + 1, rcc->non_homopolymer_errors);
// } else {
// fprintf(stderr, "[M::%s]\tunalined\n", __func__);
// }
}
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/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), /**((asm_opt.is_sc)?&(b->v8q):NULL)**/&(b->v8q), ((asm_opt.is_sc)?&(b->v8t):NULL), (asm_opt.is_ont)?1:0, R_INF.tqn, 0/**1**/, HC0_W, &b->v32,
asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, (((double)R_INF.tr[1])/((double)(R_INF.tr[0] + R_INF.tr[1]))));
rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &(scb.a[i]), &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), /**((asm_opt.is_sc)?&(b->v8q):NULL)**/&(b->v8q), ((asm_opt.is_sc)?&(b->v8t):NULL), (asm_opt.is_ont)?1:0, R_INF.tqn, 0/**1**/, HC0_W, &b->v32,
asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, (((double)R_INF.tr[1])/((double)(R_INF.tr[0] + R_INF.tr[1]))), rcc, rl0);
copy_asg_arr(b->sp, buf0);
///for debug indel
// if(i == 23863) stderr_phase_ovlp(&b->olist);
// stderr_phase_ovlp(&b->olist);
// exit(1);
dedup_chains(&b->olist);
@@ -5449,11 +5822,17 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid)
R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, &buf1);
copy_asg_arr(b->sp, buf0); copy_asg_arr(b->hap.snp_srt, buf1);
push_nec_re(aux_o, &(scc.a[i]));
// if(DBG_TIME && dbg_a) {
// dbg_a[i].faln = b->cnt[1];
// }
push_nec_re(aux_o, &(scc.a[i]));
// cmp_smp_ac(&b->self_read, &(scc.a[i]), &b->ovlp_read, i);///for debug
// push_nec_re(aux_o, &(scb.a[i]));
if(asm_opt.dbg_bam) {
/**if(asm_opt.dbg_bam)**/ {
update_scb(&R_INF, &(scc.a[i]), &(scb.a[i]), &(b->v16), &b->self_read, &b->ovlp_read, &b->v64, &b->exz, i);
kv_resize(uint16_t, scb.a[i], b->v16.n); memcpy(scb.a[i].a, b->v16.a, b->v16.n);
kv_resize(uint16_t, scb.a[i], b->v16.n); memcpy(scb.a[i].a, b->v16.a, b->v16.n * sizeof((*(b->v16.a)))); scb.a[i].n = b->v16.n;
// fprintf(stderr, "-b-[M::%s-beg] rid->%ld, rlen->%lu, scb.a[i].n::%u\n", __func__, i, Get_READ_LENGTH((R_INF),i), (uint32_t)scb.a[i].n);
}
// if((asm_opt.is_ont) && is_chemical_r_qual(&b->olist, &b->v64, qlen, 1, 16, &(b->v8q), i)/**(is_uncorrected_read(&b->olist, &b->v64, qlen, 1600))**/) {
@@ -5607,8 +5986,8 @@ static void worker_hap_ec_hybrid_sync(void *data, long i, int tid)
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/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), /**((asm_opt.is_sc)?&(b->v8q):NULL)**/&(b->v8q), ((asm_opt.is_sc)?&(b->v8t):NULL), (asm_opt.is_ont)?1:0, R_INF.tqn, 0/**1**/, HC0_W, &b->v32,
asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, hf_rate);
rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, NULL, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), /**((asm_opt.is_sc)?&(b->v8q):NULL)**/&(b->v8q), ((asm_opt.is_sc)?&(b->v8t):NULL), (asm_opt.is_ont)?1:0, R_INF.tqn, 0/**1**/, HC0_W, &b->v32,
asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, hf_rate, NULL, -1);
copy_asg_arr(b->sp, buf0);
///for debug indel
// if(i == 23863) stderr_phase_ovlp(&b->olist);
@@ -8090,9 +8469,9 @@ overlap_region* h_ec_lchain_re3(ha_abuf_t *ab, uint32_t rid, UC_Read *qu, UC_Rea
}
void gen_ori_seq0(char *tstr, uint64_t tl, UC_Read *qu, asg16_v *sc, uint64_t rid)
uint64_t gen_ori_seq0(char *tstr, uint64_t tl, UC_Read *qu, asg16_v *sc, uint64_t rid)
{
uint64_t ck, qk, tk, k, wq[2], wt[2]; uint32_t len; uint16_t c, bq, bt; char *qstr = NULL;
uint64_t ck, qk, tk, k, wq[2], wt[2], tot_e = 0; uint32_t len; uint16_t c, bq, bt; char *qstr = NULL;
ck = qk = tk = 0;
while (ck < sc->n) {
@@ -8101,7 +8480,14 @@ void gen_ori_seq0(char *tstr, uint64_t tl, UC_Read *qu, asg16_v *sc, uint64_t ri
if(c != 2) qk += len;
if(c != 3) tk += len;
wq[1] = qk; wt[1] = tk;
if(c!=0) tot_e += len;
// if(rid == 24) {
// fprintf(stderr, "%u(%c)\t", len, "MSID"[c]);
// }
}
// if(rid == 24) {
// fprintf(stderr, "\n");
// }
// if(!(tk == tl)) {
// if(rid == 8) {
// fprintf(stderr, "[M::%s] rid::%lu, tk::%lu, tl::%lu\n", __func__, rid, tk, tl);
@@ -8116,6 +8502,9 @@ void gen_ori_seq0(char *tstr, uint64_t tl, UC_Read *qu, asg16_v *sc, uint64_t ri
// }
// }
// }
// if((!(tk == tl)) && (rid == 24)) {
// fprintf(stderr, "[M::%s]\trid::%lu\tqk::%lu\ttk::%lu\ttl::%lu\n", __func__, rid, qk, tk, tl);
// }
assert(tk == tl);
@@ -8135,6 +8524,8 @@ void gen_ori_seq0(char *tstr, uint64_t tl, UC_Read *qu, asg16_v *sc, uint64_t ri
}
// fprintf(stderr, "%u%c(%c)(x::[%lu,%ld))(y::[%lu,%ld))\n", len, cm[c], ((c==1)||(c==2))?(cc[bt]):('*'), wx[0], wx[1], wy[0], wy[1]); // s_H
}
return tot_e;
}
void gen_cc_fly(asg16_v *sc, char *qstr, uint64_t ql, char *tstr, uint64_t tl, bit_extz_t *exz, double e_rate, uint64_t maxn, uint64_t maxe)
@@ -8259,6 +8650,10 @@ void cal_updated_trace_len(asg16_v *sc, uint64_t *ql, uint64_t *tl)
*ql = qk; *tl = tk;
}
///qstr:: latest; tstr:: original; there is an intermidate string I between qstr and tstr
///tcc:: tstr -> I;
///qcc:: I -> qstr;
///tcc_res:: tstr -> qstr
void gen_updated_trace(asg16_v *qcc, asg16_v *tcc, asg16_v *tcc_res, char *qstr, uint64_t ql, char *tstr, uint64_t tl, asg64_v *srt, bit_extz_t *exz, uint64_t rid)
{
uint64_t k, ck, qk, tk, wq[2], wt[2], old_dp, dp, s, e, srt_n, si, ei, so, os, oe, *qd, *td, qs, qe, ts, te, q0, t0;
@@ -8466,8 +8861,8 @@ void update_scb(All_reads *R_INF, asg16_v *scc, asg16_v *scb, asg16_v *scb_res,
// if(i == 700) fprintf(stderr, "|%u%c(%c)(x::%u)(y::%u)", len, cm[c], ((c==1)||(c==2))?(cc[b]):('*'), wx[1], wy[1]); // s_H
}
qstr = tstr; ql = tl;
tstr = tu->seq; tl = tu->length;
qstr = tstr; ql = tl;///latest version
tstr = tu->seq; tl = tu->length;///orginal version
// fprintf(stderr, "\n[M::%s] ql::%lu, tl::%lu, rid::%lu\n", __func__, ql, tl, rid);
@@ -8575,8 +8970,8 @@ static void worker_hap_dc_ec0(void *data, long i, int tid)
b->cnt[0] += b->self_read.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/**, 1**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), /**((asm_opt.is_sc)?&(b->v8q):NULL)**/&(b->v8q), ((asm_opt.is_sc)?&(b->v8t):NULL), (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32,
asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1);
rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, NULL, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 1**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), /**((asm_opt.is_sc)?&(b->v8q):NULL)**/&(b->v8q), ((asm_opt.is_sc)?&(b->v8t):NULL), (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32,
asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1, NULL, -1);
copy_asg_arr(b->sp, buf0);
copy_asg_arr(buf0, b->sp);
@@ -9299,7 +9694,7 @@ uint64_t cal_ec_multiple_step(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a, u
// fprintf(stderr, "[M::%s] # corrected bases->%lu\n", __func__, num_correct);
// fprintf(stderr, "[M::%s::%.3f] running time\n", __func__, yak_realtime_0()-tt0);
fprintf(stderr, "[M::pec::%.3f] # bases: %lu; # corrected bases: %lu\n", yak_realtime_0()-tt0, num_base, num_correct);
exit(1);
// exit(1);
(*r_base) = num_base;
return num_correct;
@@ -9567,7 +9962,7 @@ void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, u
// prt_nel_ovlp(R_INF.paf, n_a);
// exit(1);
// dbg_write_ec_reads("ec12.fa", round, &scb, is_cr);
// dbg_write_ec_reads("ec12.fa", round, &scb, 0/**is_cr**/);
if((!is_sv) || (is_sv && is_cr)) {
kt_for(n_thre, worker_hap_post_rev, b, n_a);
@@ -9582,7 +9977,7 @@ void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, u
fprintf(stderr, "-4-[M::%s]\t# tqn::%lu, Ont base::%lu, # HiFi bases::%lu\n", __func__, R_INF.tqn, R_INF.tr[0], R_INF.tr[1]);
// dbg_write_ec_reads("ec16.fa", round, &scb, !is_cr);
// dbg_write_ec_reads("ec16.fa", round, &scb, 0/**!is_cr**/);
// exit(1);
// uint64_t z;
@@ -9753,12 +10148,12 @@ void cal_ov_r(uint64_t n_thre, uint64_t n_a, uint64_t new_idx)
b = gen_ec_ovec_buf_t(n_thre);
if(new_idx) {
// kt_for(n_thre, worker_hap_dc_ec, b, n_a);///update overlaps
destroy_cc_v(&scc); if(!asm_opt.dbg_bam) destroy_cc_v(&scb); destroy_cc_v(&sca);
destroy_cc_v(&scc); /**if(!asm_opt.dbg_bam) destroy_cc_v(&scb);**/ destroy_cc_v(&sca);
ha_print_ovlp_stat_0(b, n_thre, n_a);
} else {
ha_print_ovlp_stat_1(b, n_thre, n_a);
destroy_cc_v(&scc); if(!asm_opt.dbg_bam) destroy_cc_v(&scb); destroy_cc_v(&sca);
destroy_cc_v(&scc); /**if(!asm_opt.dbg_bam) destroy_cc_v(&scb);**/ destroy_cc_v(&sca);
}
destroy_ec_ovec_buf_t(b);
gfa_ut.cpp
+2 -1
View File
@@ -3053,13 +3053,14 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
// fprintf(stderr, "%.*s\tid::%u\tis_c::%u\n",
// (int)Get_NAME_LENGTH(R_INF, 10819), Get_NAME(R_INF, 10819), 10819, is_contain_r((*rI), 10819));
// debug_info_of_specfic_node("m64011_190830_220126/47516220/ccs", sg, rI, "beg-0");
// debug_info_of_specfic_node("bcb40bcc-d9cf-48e6-88ee-47ac3dde22ff", sg, rI, "beg-0");
// debug_info_of_specfic_node("c7ecbd6b-e09d-4042-93ac-2400839feaf6", sg, rI, "beg-0");
if(asm_opt.is_ont) asg_arc_cut_weak(sg, &bu, max_tip, 0.975, 0, is_ou, 0, 1, 16, UL_COV_THRES-1, 0, rev, NULL, NULL);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);///p_telo
// fprintf(stderr, "[M::%s] count_edges_v_w(sg, 49778, 49847)->%ld\n", __func__, count_edges_v_w(sg, 49778, 49847));
// if(is_ou) dedup_contain_g(uopt, sg);
// debug_info_of_specfic_node("c7ecbd6b-e09d-4042-93ac-2400839feaf6", sg, rI, "beg-1");
for (i = 0; i < clean_round; i++, drop += step) {
if(drop > max_ovlp_drop_ratio) drop = max_ovlp_drop_ratio;
if(is_ou) {
htab.cpp
+20 -6
View File
@@ -9,6 +9,7 @@
#include "ksort.h"
#include "htab.h"
#include "Process_Read.h"
#include "ecovlp.h"
#define YAK_COUNTER_BITS 12
#define YAK_N_COUNTS (1<<YAK_COUNTER_BITS)
@@ -1627,12 +1628,11 @@ void refresh_pt_idx(void **flt_tab, ha_pt_t **ha_idx, All_reads *r, hifiasm_opt_
}
uint64_t tmp_pt_pro(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, hifiasm_opt_t *opt, char *file_name, uint64_t rr, uint64_t tot_rr, uint64_t is_load)
uint64_t tmp_pt_pro(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, cc_v* rcc, hifiasm_opt_t *opt, char *file_name, uint64_t rr, uint64_t tot_rr, uint64_t is_load)
{
char* gfa_name = (char*)malloc(strlen(file_name)+64);
FILE *fp = NULL; int f_flag = 0; uint64_t rr0 = (uint64_t)-1, tot_rr0 = (uint64_t)-1;
if(is_load) {
sprintf(gfa_name, "%s.r%lu.ht.bin", file_name, rr);
fp = fopen(gfa_name, "r");
@@ -1648,19 +1648,33 @@ uint64_t tmp_pt_pro(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, hifiasm_
return 0;
}
sprintf(gfa_name, "%s.r%lu", file_name, rr);
if(!load_pt_index(r_flt_tab, r_ha_idx, r, opt, gfa_name)) {
destory_All_reads(r);
ha_pt_destroy(*r_ha_idx); (*r_ha_idx) = NULL;
ha_ft_destroy(*r_flt_tab); (*r_flt_tab) = NULL;
free(gfa_name);
return 0;
}
sprintf(gfa_name, "%s.r%lu.rec", file_name, rr);
if(!load_cc_v(rcc, gfa_name)) {
destory_All_reads(r);
ha_pt_destroy(*r_ha_idx); (*r_ha_idx) = NULL;
ha_ft_destroy(*r_flt_tab); (*r_flt_tab) = NULL;
destroy_cc_v(rcc);
free(gfa_name);
return 0;
}
fprintf(stderr, "[M::%s] restarting from round %lu\n", __func__, rr0);
} else {
sprintf(gfa_name, "%s.r%lu.rec", file_name, rr);
write_cc_v(rcc, gfa_name);
sprintf(gfa_name, "%s.r%lu", file_name, rr);
write_pt_index(*r_flt_tab, *r_ha_idx, r, opt, gfa_name);
sprintf(gfa_name, "%s.r%lu.ht.bin", file_name, rr);
fp = fopen(gfa_name, "w");
if (!fp) {
htab.h
+1 -1
View File
@@ -94,7 +94,7 @@ int uidx_load(void **r_flt_tab, ha_pt_t **r_ha_idx, char* file_name, ma_ug_t *ug
int write_ct_index(void *ct_idx, char* file_name);
int load_ct_index(void **ct_idx, char* file_name);
int query_ct_index(void* ct_idx, uint64_t hash);
uint64_t tmp_pt_pro(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, hifiasm_opt_t *opt, char *file_name, uint64_t rr, uint64_t tot_rr, uint64_t is_load);
uint64_t tmp_pt_pro(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, cc_v* rcc, hifiasm_opt_t *opt, char *file_name, uint64_t rr, uint64_t tot_rr, uint64_t is_load);
ha_abuf_t *ha_abuf_init_buf(void *km);
ha_abufl_t *ha_abufl_init_buf(void *km);