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

0.20.0-r631

Browse Source
This commit is contained in:
chhylp123
2024-10-09 17:09:10 -04:00
parent d5f8a8a6c0
commit e8b18560a5
14 changed files with 7710 additions and 356 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Assembly.cpp
+57 -62
View File
@@ -970,86 +970,57 @@ void prt_dbg_rs(FILE *fp, Debug_reads* x, uint64_t round)
destory_UC_Read(&g_read);
}
void ha_ec(int64_t round)
void ha_ec(int64_t round, int num_pround, int des_idx, uint64_t *tot_b, uint64_t *tot_e)
{
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;
int hom_cov, het_cov, r_out = 0;
ha_flt_tab_hp = ha_idx_hp = NULL; (*tot_b) = (*tot_e) = 0;
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
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
asm_opt.hom_cov = hom_cov; asm_opt.het_cov = het_cov;
}
///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;
cal_ec_r(asm_opt.thread_num, round, num_pround, R_INF.total_reads, (round == (asm_opt.number_of_round-1))?1:0, tot_b, tot_e);
// exit(1);
// if (r_out) write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name);
if(des_idx) {
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);
// 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);
// sl_ec_r(asm_opt.thread_num, R_INF.total_reads);
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();
// Output_corrected_reads();
// exit(1);
}
@@ -1927,6 +1898,25 @@ void ha_overlap_final(void)
asm_opt.het_cov = het_cov;
}
void ha_ec_ff(int renew_idx)
{
int hom_cov, het_cov;
ha_flt_tab_hp = ha_idx_hp = NULL;
if(ha_idx && renew_idx) {
ha_pt_destroy(ha_idx); ha_idx = NULL;
}
if(!ha_idx) {
ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, 0, &R_INF, &hom_cov, &het_cov); // build the index
asm_opt.hom_cov = hom_cov; asm_opt.het_cov = het_cov;
}
cal_ov_r(asm_opt.thread_num, R_INF.total_reads, renew_idx);
ha_pt_destroy(ha_idx); ha_idx = NULL;
}
static void worker_ov_utg(void *data, long i, int tid)
{
ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
@@ -2016,7 +2006,7 @@ int ha_assemble(void)
// debug_mc_gg_t(MC_NAME, 0, 0);
// quick_debug_phasing(MC_NAME);
extern void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o);
int r, hom_cov = -1, ovlp_loaded = 0;
int r, 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)) {
ovlp_loaded = 1;
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> loaded corrected reads and overlaps from disk\n", __func__, yak_realtime(), yak_cpu_usage());
@@ -2042,24 +2032,29 @@ 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
tot_b = tot_e = 0;
// ha_overlap_and_correct(r);
ha_ec(r);
ha_ec(r, asm_opt.number_of_pround, (r<asm_opt.number_of_round-1)?1:0, &tot_b, &tot_e);
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__,
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);
fprintf(stderr, "[M::%s] # bases: %lu; # corrected bases: %lu\n", __func__, tot_b, tot_e);
// 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.flag & HA_F_WRITE_EC) Output_corrected_reads();
// overlap between corrected reads
ha_opt_reset_to_round(&asm_opt, asm_opt.number_of_round);
ha_overlap_final();
fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> found overlaps for the final round\n", __func__, yak_realtime(),
yak_cpu_usage(), yak_peakrss_in_gb());
ha_print_ovlp_stat(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads);
// ha_overlap_final();
ha_ec_ff(1/**0**/);
fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> found overlaps for the final round\n", __func__, yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb());
// fprintf(stderr, "\n[M::%s::%.3f*%.2f@%.3fGB] ==> found overlaps for the final round\n", __func__, yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb());
// ha_print_ovlp_stat(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads);
ha_ft_destroy(ha_flt_tab);
if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
ha_triobin(&asm_opt);
// exit(1);
}
if(ovlp_loaded == 2) ovlp_loaded = 0;
ha_opt_update_cov_min(&asm_opt, asm_opt.hom_cov, MIN_N_CHAIN);
CommandLines.cpp
+3 -2
View File
@@ -18,8 +18,8 @@ static ko_longopt_t long_options[] = {
{ "write-paf", ko_no_argument, 302 },
{ "write-ec", ko_no_argument, 303 },
{ "skip-triobin", ko_no_argument, 304 },
{ "max-od-ec", ko_no_argument, 305 },
{ "max-od-final", ko_no_argument, 306 },
{ "max-od-ec", ko_required_argument, 305 },
{ "max-od-final", ko_required_argument, 306 },
{ "ex-list", ko_required_argument, 307 },
{ "ex-iter", ko_required_argument, 308 },
{ "hom-cov", ko_required_argument, 309 },
@@ -249,6 +249,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->load_index_from_disk = 1;
asm_opt->write_index_to_disk = 1;
asm_opt->number_of_round = 3;
asm_opt->number_of_pround = 0/**3**/;
asm_opt->adapterLen = 0;
asm_opt->clean_round = 4;
///asm_opt->small_pop_bubble_size = 100000;
CommandLines.h
+2 -1
View File
@@ -5,7 +5,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.19.9-r616"
#define HA_VERSION "0.20.0-r631"
#define VERBOSE 0
@@ -74,6 +74,7 @@ typedef struct {
int load_index_from_disk;
int write_index_to_disk;
int number_of_round;
int number_of_pround;
int adapterLen;
int clean_round;
int roundID;
Correct.cpp
+1256 -69
View File
File diff suppressed because it is too large Load Diff
Correct.h
+7 -5
View File
@@ -1389,7 +1389,13 @@ const ul_idx_t *uref, char* qstr, UC_Read *tu, overlap_region_alloc *ol, overlap
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);
void 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);
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);
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);
#define ovlp_id(x) ((x).tn)
#define ovlp_min_wid(x) ((x).ts)
@@ -1400,9 +1406,5 @@ void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_all
#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
+47 -4
View File
@@ -1512,6 +1512,34 @@ inline int32_t comput_sc_ch(const k_mer_hit *ai, const k_mer_hit *aj, double bw_
return sc;
}
inline int32_t comput_sc_ch_ec(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
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);
///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*(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
@@ -1999,6 +2027,9 @@ int64_t *p, int64_t *t, int32_t *f, int32_t *ii, int64_t *plus, int64_t *msc, in
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(a_n && a[0].readID == 3125488) {
// fprintf(stderr, "[M::%s::] ii::[%ld,%ld)(%c), is_srt::%ld, chn_pen_gap::%f, chn_pen_skip::%f, bw_rate::%f\n", __func__, l, k, "+-"[a[l].strand], is_srt, chn_pen_gap, chn_pen_skip, bw_rate);
// }
if(is_srt) {
plus0 = 0; msc0 = msc_i0 = INT32_MIN; movl0 = INT32_MAX; ddt = 0;
@@ -2016,6 +2047,9 @@ int64_t *p, int64_t *t, int32_t *f, int32_t *ii, int64_t *plus, int64_t *msc, in
dr = (int64_t)(ai->offset) - (int64_t)(aj->offset);
if(dr <= 0) break;
dd = dr > dq? dr - dq : dq - dr;//gap
// if(a_n && a[0].readID == 3125488) {
// fprintf(stderr, "%ld,", dd);
// }
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);
@@ -2024,7 +2058,10 @@ int64_t *p, int64_t *t, int32_t *f, int32_t *ii, int64_t *plus, int64_t *msc, in
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;
///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*(double)dg);
sc -= (int32_t)lin_pen;
}
@@ -2036,6 +2073,10 @@ int64_t *p, int64_t *t, int32_t *f, int32_t *ii, int64_t *plus, int64_t *msc, in
if(f[z] < plus0) plus0 = f[z];
}
// if(a_n && a[0].readID == 3125488) {
// fprintf(stderr, "\n");
// fprintf(stderr, "[M::%s::] msc0::%ld, msc_i0::%ld, (%c)\n", __func__, msc0, msc_i0, "+-"[a[l].strand]);
// }
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)) {
@@ -2087,7 +2128,9 @@ uint64_t lchain_qdp_mcopy_fast(Candidates_list *cl, int64_t a_idx, int64_t a_n,
msc = msc_i = INT32_MIN; movl = INT32_MAX; plus = 0; si = 0; ei = a_n;
memset(t, 0, (a_n*sizeof((*t))));
}
// if(a_n && a[0].readID == 3125488) {
// fprintf(stderr, "[M::%s::] si::%ld, ei::%ld, a_n::%ld\n", __func__, si, ei, a_n);
// }
for (i = st = si, max_ii = -1; i < ei; ++i) {
max_f = a[i].cnt&(0xffu);
n_skip = 0; max_j = end_j = -1;
@@ -2095,7 +2138,7 @@ uint64_t lchain_qdp_mcopy_fast(Candidates_list *cl, int64_t a_idx, int64_t a_n,
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);
sc = comput_sc_ch_ec(&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) {
@@ -2119,7 +2162,7 @@ uint64_t lchain_qdp_mcopy_fast(Candidates_list *cl, int64_t a_idx, int64_t a_n,
}
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);
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;
}
Hash_Table.h
+1
View File
@@ -9,6 +9,7 @@
#define WINDOW 375
#define WINDOW_BOUNDARY 375
#define WINDOW_HC 775
#define WINDOW_HC_FAST 512
///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
Levenshtein_distance.h
+185 -8
View File
@@ -548,24 +548,201 @@ inline int32_t pop_trace_back(asg16_v *res, int32_t i, uint16_t *c, uint32_t *le
return i;
}
///compact functions
#define pop_trac_bpc(in, rc, rb, rl) do { \
(rc) = ((in)>>14);\
if((rc) == 1 || (rc) == 2) {(rb) = (((in)>>12)&3); (rl) = ((in)&(0xfff));}\
else {(rl) = ((in)&(0x3fff));}\
} while (0)
inline void push_trace_bp(asg16_v *res, uint16_t c, uint16_t b, uint32_t len, uint32_t is_append)
{
uint16_t p;
uint16_t p, c0, b0, len0, mm;
if((is_append) && (res->n)) {
b0 = b;
pop_trac_bpc(res->a[res->n-1], c0, b0, len0);
if((c == c0) && (b == b0)) {
res->n--; len += len0;
}
}
mm = (0x3fff); c0 = c; c <<= 14;
if(c0 == 1 || c0 == 2) {
mm = (0xfff); c += ((b&3) << 12);
}
while (len >= mm) {
p = (c + mm); kv_push(uint16_t, *res, p); len -= mm;
}
c <<= 14;
while (len >= (0x3fff)) {
p = (c + (0x3fff)); kv_push(uint16_t, *res, p); len -= (0x3fff);
}
// fprintf(stderr, "[M::%s] c::%u, len::%u\n", __func__, c, len);
if(len) {
p = (c + len); kv_push(uint16_t, *res, p);
}
}
inline uint32_t pop_trace_bp(asg16_v *res, uint32_t i, uint16_t *c, uint16_t *b, uint32_t *len)
{
(*c) = (res->a[i]>>14);
if((*c) == 1 || (*c) == 2) {
(*b) = ((res->a[i]>>12)&3);
(*len) = (res->a[i]&(0xfff));
} else {
(*b) = (uint16_t)-1;
(*len) = (res->a[i]&(0x3fff));
}
uint32_t sl; uint16_t sb;
for (i++; (i < res->n) && ((*c) == (res->a[i]>>14)); i++) {
if((*c) == 1 || (*c) == 2) {
sb = ((res->a[i]>>12)&3); sl = (res->a[i]&(0xfff));
} else {
sb = (uint16_t)-1; sl = (res->a[i]&(0x3fff));
}
if((*b) != sb) break;
(*len) += sl;
}
return i;
}
inline int64_t pop_trace_bp_rev(asg16_v *res, int64_t i, uint16_t *c, uint16_t *b, uint32_t *len)
{
(*c) = (res->a[i]>>14);
if((*c) == 1 || (*c) == 2) {
(*b) = ((res->a[i]>>12)&3);
(*len) = (res->a[i]&(0xfff));
} else {
(*b) = (uint16_t)-1;
(*len) = (res->a[i]&(0x3fff));
}
uint32_t sl; uint16_t sb;
for (i--; (i >= 0) && ((*c) == (res->a[i]>>14)); i--) {
if((*c) == 1 || (*c) == 2) {
sb = ((res->a[i]>>12)&3); sl = (res->a[i]&(0xfff));
} else {
sb = (uint16_t)-1; sl = (res->a[i]&(0x3fff));
}
if((*b) != sb) break;
(*len) += sl;
}
return i;
}
///full functions
#define pop_trac_bpc_f(in, rc, rbq, rbt, rl) do { \
(rc) = ((in)>>14);\
if((rc) == 2 || (rc) == 3) {(rbt) = (((in)>>12)&3); (rl) = ((in)&(0xfff));}\
else if((rc) == 1) {(rbt) = (((in)>>12)&3); (rbq) = (((in)>>10)&3); (rl) = ((in)&(0x3ff));}\
else {(rl) = ((in)&(0x3fff));}\
} while (0)
inline void push_trace_bp_f(asg16_v *res, uint16_t c, uint16_t bq, uint16_t bt, uint32_t len, uint32_t is_append)
{
uint16_t p, c0 = c, bq0, bt0, len0, mm;
if(c == 3) {
bt = bq; bq = (uint16_t)-1;
}
if((is_append) && (res->n)) {
bq0 = bq; bt0 = bt;
pop_trac_bpc_f(res->a[res->n-1], c0, bq0, bt0, len0);
if((c == c0) && (bq == bq0) && (bt == bt0)) {
res->n--; len += len0;
}
}
c0 = c; c <<= 14;
if(c0 == 2 || c0 == 3) {
mm = (0xfff); c += ((bt&3) << 12);
} else if(c0 == 1) {
mm = (0x3ff); c += ((bt&3) << 12); c += ((bq&3) << 10);
} else {
mm = (0x3fff);
}
while (len >= mm) {
p = (c + mm); kv_push(uint16_t, *res, p); len -= mm;
}
// fprintf(stderr, "[M::%s] c::%u, len::%u\n", __func__, c, len);
if(len) {
p = (c + len); kv_push(uint16_t, *res, p);
}
}
inline uint32_t pop_trace_bp_f(asg16_v *res, uint32_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; uint16_t sbq, sbt;
for (i++; (i < res->n) && ((*c) == (res->a[i]>>14)); i++) {
sbq = sbt = (uint16_t)-1;
if((*c) == 2 || (*c) == 3) {
sbt = ((res->a[i]>>12)&3);
sl = (res->a[i]&(0xfff));
} else if((*c) == 1) {
sbt = ((res->a[i]>>12)&3);
sbq = ((res->a[i]>>10)&3);
sl = (res->a[i]&(0x3ff));
} else {
sl = (res->a[i]&(0x3fff));
}
if((*bq) != sbq || (*bt) != sbt) break;
(*len) += sl;
}
if((*c) == 3) {
(*bq) = (*bt); (*bt) = (uint16_t)-1;
}
return i;
}
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; uint16_t sbq, sbt;
for (i--; (i >= 0) && ((*c) == (res->a[i]>>14)); i--) {
sbq = sbt = (uint16_t)-1;
if((*c) == 2 || (*c) == 3) {
sbt = ((res->a[i]>>12)&3);
sl = (res->a[i]&(0xfff));
} else if((*c) == 1) {
sbt = ((res->a[i]>>12)&3);
sbq = ((res->a[i]>>10)&3);
sl = (res->a[i]&(0x3ff));
} else {
sl = (res->a[i]&(0x3fff));
}
if((*bq) != sbq || (*bt) != sbt) break;
(*len) += sl;
}
if((*c) == 3) {
(*bq) = (*bt); (*bt) = (uint16_t)-1;
}
return i;
}
///511 -> 16 64-bits
// #define MAX_E 511
// #define MAX_L 2500
Overlaps.h
+4
View File
@@ -1242,4 +1242,8 @@ uint64_t trans_sec_cut0(kv_u_trans_t *ta, asg64_v *srt, uint32_t id, double sec_
void clean_u_trans_t_idx_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* src, ug_rid_cov_t *in);
void gen_ug_rid_cov_t_by_ovlp(kv_u_trans_t *ta, ug_rid_cov_t *cc);
#define UC_Read_resize(v, s) do {\
if ((v).size<(s)) {REALLOC((v).seq,(s));(v).size=(s);}\
} while (0)
#endif
Process_Read.h
+2
View File
@@ -107,6 +107,8 @@ typedef struct
#define CHAIN_MATCH 1
#define CHAIN_UNMATCH 0.334
#define NEC 1
typedef struct
{
uint64_t** N_site;
anchor.cpp
+1505 -45
View File
File diff suppressed because it is too large Load Diff
ecovlp.cpp
+4625 -144
View File
File diff suppressed because it is too large Load Diff
ecovlp.h
+15 -16
View File
@@ -5,60 +5,59 @@
#include <stdint.h>
#include "Hash_Table.h"
#include "Process_Read.h"
#include "kdq.h"
KDQ_INIT(uint32_t)
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 {size_t n, m, nou; 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_arc_v arc;
}cns_t;
typedef struct {
size_t n, m;
cns_t *a;
uint32_t si, ei;
uint32_t si, ei, off, bn, bb0, bb1, cns_g_wl;
kdq_t(uint32_t) *q;
}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;
// int64_t num_read_base, num_correct_base, num_recorrect_base;
uint64_t cnt[6], rr;
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;
cns_gfa cns;
} ec_ovec_buf_t0;
typedef struct {
ec_ovec_buf_t0 *a;
uint32_t n;
uint32_t n, rev;
} ec_ovec_buf_t;
ec_ovec_buf_t* gen_ec_ovec_buf_t(uint32_t n, uint32_t is_final, uint32_t save_ov);
ec_ovec_buf_t* gen_ec_ovec_buf_t(uint32_t n);
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);
void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, uint64_t is_sv, uint64_t *tot_b, uint64_t *tot_e);
void sl_ec_r(uint64_t n_thre, uint64_t n_a);
void cal_ov_r(uint64_t n_thre, uint64_t n_a, uint64_t new_idx);
#endif
htab.h
+1
View File
@@ -125,6 +125,7 @@ int adj_m_peak_hom(int m_peak_hom, int max_i, int max2_i, int max3_i, int *peak_
void print_hist_lines(int n_cnt, int start_cnt, const int64_t *cnt);
void debug_adapter(const hifiasm_opt_t *asm_opt, All_reads *rs);
inline int mz_low_b(int peak_hom, int peak_het)
{
int low_freq = 2;