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

lack final alignment

Browse Source
This commit is contained in:
chhylp123
2022-10-01 19:53:02 -04:00
parent 20a7dac82a
commit bbd3394b4a
10 changed files with 1397 additions and 165 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Assembly.cpp
+2
View File
@@ -468,6 +468,7 @@ ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov, int is_ug)
kv_init(b->k_flag.a);
kv_init(b->sp);
init_bit_extz_t(&(b->exz), 31);
init_bit_extz_t(&(b->exz1), 31);
if(!is_ug) b->ab = ha_abuf_init();
else b->abl = ha_abufl_init();
if (!b->is_final) {
@@ -497,6 +498,7 @@ void ha_ovec_destroy(ha_ovec_buf_t *b)
kv_destroy(b->k_flag.a);
kv_destroy(b->sp);
destroy_bit_extz_t(&(b->exz));
destroy_bit_extz_t(&(b->exz1));
if (!b->is_final) {
destory_Cigar_record(&b->cigar1);
destory_Graph(&b->POA_Graph);
Assembly.h
+1 -1
View File
@@ -40,7 +40,7 @@ typedef struct {
overlap_region tmp_region;
ma_utg_v *ua;
st_mt_t sp;
bit_extz_t exz;
bit_extz_t exz, exz1;
} ha_ovec_buf_t;
int ha_assemble(void);
Correct.cpp
+888 -130
View File
File diff suppressed because it is too large Load Diff
Correct.h
+2 -2
View File
@@ -20,7 +20,7 @@
#define ERROR_RATE 1.25
#define UL_TOPN 50
#define SGAP 16
#define MAX_LGAP(ql) ((((ql)*0.33)>2000)?((ql)*0.33):2000)
#define MAX_LGAP(ql) ((((ql)*0.2)<256)?((ql)*0.2):256)
#define WINDOW_MAX_SIZE (WINDOW + (int)(1.0 / HA_MIN_OV_DIFF) + 3) // TODO: why 1/max_ov_diff?
@@ -1136,7 +1136,7 @@ void correct_ul_overlap(overlap_region_alloc* overlap_list, const ul_idx_t *uref
double max_ov_diff_ec, long long winLen, void *km);
void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *uref, char *qstr,
uint64_t ql, UC_Read* qu, UC_Read* tu, Correct_dumy* dumy, bit_extz_t *exz,
haplotype_evdience_alloc* hap, kvec_t_u64_warp* v_idx,
bit_extz_t *exz1, haplotype_evdience_alloc* hap, kvec_t_u64_warp* v_idx,
double e_rate, int64_t wl, kv_ul_ov_t *aln, int64_t sid, uint64_t hpc_k, void *km);
void ul_lalign_old_ed(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *uref, char *qstr,
Hash_Table.cpp
+283
View File
@@ -185,6 +185,60 @@ int ovlp_chain_gen(overlap_region_alloc* ol, overlap_region* t, int64_t xl, int6
return 1;
}
int ovlp_chain_qgen(overlap_region_alloc* ol, overlap_region* t, int64_t xl, int64_t yl, int64_t apend_be, k_mer_hit* hit, int64_t n_hit)
{
if (ol->length + 1 > ol->size) {
uint64_t sl = ol->size;
ol->size = ol->length + 1;
kroundup64(ol->size);
REALLOC(ol->list, ol->size);
/// need to set new space to be 0
memset(ol->list + sl, 0, sizeof(overlap_region)*(ol->size - sl));
}
if ((ol->length!=0) && (ol->list[ol->length-1].y_id==t->y_id)) {
if((ol->list[ol->length-1].shared_seed > t->shared_seed) ||
((ol->list[ol->length-1].shared_seed == t->shared_seed) &&
(ol->list[ol->length-1].overlapLen <= t->overlapLen))) {
return 0;
} else {
ol->length--;
}
}
int64_t xr, yr;
if(t->x_pos_s <= t->y_pos_s) {
t->y_pos_s -= t->x_pos_s; t->x_pos_s = 0;
} else {
t->x_pos_s -= t->y_pos_s; t->y_pos_s = 0;
}
xr = xl-t->x_pos_e-1; yr = yl-t->y_pos_e-1;
if(xr <= yr) {
t->x_pos_e = xl-1; t->y_pos_e += xr;
} else {
t->y_pos_e = yl-1; t->x_pos_e += yr;
}
overlap_region *o = &(ol->list[ol->length++]);
o->shared_seed = t->shared_seed;
o->align_length = 0;
o->is_match = 0;
o->non_homopolymer_errors = 0;
o->strong = 0;
o->x_id = t->x_id;
o->y_id = t->y_id;
o->x_pos_strand = 0;///always 0
o->y_pos_strand = t->x_pos_strand;
o->x_pos_e = t->x_pos_e; o->x_pos_s = t->x_pos_s;
o->y_pos_e = t->y_pos_e; o->y_pos_s = t->y_pos_s;
///debug
// debug_cigar(&(t->f_cigar), o, apend_be, hit, n_hit);
return 1;
}
int ovlp_chain_gen_fcigar(overlap_region_alloc* ol, overlap_region* t, int64_t xl, int64_t yl, int64_t apend_be, k_mer_hit* hit, int64_t n_hit)
{
if (ol->length + 1 > ol->size) {
@@ -1359,6 +1413,58 @@ int32_t lchain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, double bw_thres)
return n_a;
}
int32_t lchain_qcheck(k_mer_hit *a, int32_t n_a, Chain_Data *dp, double bw_thres)
{
int32_t i, tot_g = 0, sc, dg, dq, dr, dd, span;
double bw_pen;
if (n_a == 0) return -1;
if (n_a > 1) {
if ((a[0].self_offset >= a[n_a-1].self_offset)||(a[0].offset >= a[n_a-1].offset)) return -1;
dq = (int32_t)a[n_a-1].self_offset - (int32_t)a[0].self_offset;
dr = (int32_t)a[n_a-1].offset - (int32_t)a[0].offset;
dd = ((dq>=dr)? (dq-dr): (dr-dq));//gap
dg = ((dq>=dr)? (dr): (dq));///len
if (dg == 0 || dd > (dg*bw_thres)) return -1;
}
for (i = 1; i < n_a; ++i) {///a[] is sorted by self_offset
if(a[i-1].self_offset >= a[i].self_offset) break;
if(a[i-1].offset >= a[i].offset) break;
}
if (i < n_a) return -1;
bw_pen = 1.0 / bw_thres;
dp->score[0] = normal_w((a[0].cnt&(0xffu)), (a[0].cnt>>8)); dp->pre[0] = -1;
for (i = 1; i < n_a; ++i) {
dq = (int32_t)a[i].self_offset - (int32_t)a[i-1].self_offset;
dr = (int32_t)a[i].offset - (int32_t)a[i-1].offset;
dd = ((dq>=dr)? (dq-dr): (dr-dq));//gap
dg = ((dq>=dr)? (dr): (dq));///len
if(dg == 0) break;
tot_g += dd;
if (dd > THRESHOLD_MAX_SIZE && dd > (dg*bw_thres)) break;
span = a[i].cnt&(0xffu);
sc = dg < span? dg : span;
sc = normal_w(sc, ((int32_t)(a[i].cnt>>8)));
sc -= (int32_t)((((double)dd)/((double)dg))*bw_pen*((double)sc));///bw_pen is 20 for HiFi
dp->score[i] = dp->score[i-1] + sc;
dp->pre[i] = i - 1;
}
if (i < n_a) return -1;
if(n_a > 1) {
dq = (int32_t)a[n_a-1].self_offset - (int32_t)a[0].self_offset;
dr = (int32_t)a[n_a-1].offset - (int32_t)a[0].offset;
dg = ((dq>=dr)? (dr): (dq));///len
dd = tot_g;///gap
if (dd > (dg*bw_thres)) return -1;
}
return n_a;
}
inline int32_t cal_bw(const k_mer_hit *ai, const k_mer_hit *aj, double bw_rate, int64_t sf_l, int64_t ot_l)
{
///ai is the suffix of aj
@@ -1624,6 +1730,180 @@ uint64_t lchain_dp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, ov
return cL;
}
uint64_t lchain_qdp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, overlap_region* res,
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 *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);
n_skip = 0; max_j = end_j = -1;
if ((i-st) > max_iter) st = i-max_iter;
for (j = i - 1; j >= st; --j) {
sc = comput_sc_ch(&a[i], &a[j], bw_rate, chn_pen_gap, chn_pen_skip, xl, yl);
if (sc == INT32_MIN) continue;
sc += f[j];
if (sc > max_f) {
max_f = sc, max_j = j;
if (n_skip > 0) --n_skip;
} else if (t[j] == (int32_t)i) {
if (++n_skip > max_skip)
break;
}
if (p[j] >= 0) t[p[j]] = i;
}
end_j = j;
if (max_ii < 0 || ((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 (max < f[j]) {
max = f[j], max_ii = j;
}
}
}
if (max_ii >= 0 && max_ii < end_j) {///just have a try with a[i]<->a[max_ii]
tmp = comput_sc_ch(&a[i], &a[max_ii], bw_rate, chn_pen_gap, chn_pen_skip, xl, yl);
if (tmp != INT32_MIN && max_f < tmp + f[max_ii])
max_f = tmp + f[max_ii], max_j = max_ii;
}
f[i] = max_f; p[i] = max_j;
if ((max_ii < 0) || (((((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:
///a[] has been sorted by self_offset
i = msc_i;
res->x_pos_s = res->x_pos_e = a[i].self_offset;
res->y_pos_s = res->y_pos_e = a[i].offset;
res->shared_seed = msc;
cL = 0;
while (i >= 0) {
t[cL++] = i; msc_i = i; i = p[i];
}
res->x_pos_s = a[t[cL-1]].self_offset;
res->y_pos_s = a[t[cL-1]].offset;
res->overlapLen = get_chainLen(res->x_pos_s, res->x_pos_e, xl, res->y_pos_s, res->y_pos_e, yl);
for (i = 0; i < cL; i++) des[i] = a[t[cL-i-1]];
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)
{
if(a_n <= 0) return 0;
int64_t *p, *t, max_f, n_skip, st, max_j, sc, msc, msc_i, dq, dr, dd;
int32_t *f; int64_t i, j, cL = 0;
resize_Chain_Data(dp, a_n, NULL);
t = dp->tmp; f = dp->score; p = dp->pre; msc = msc_i = -1;
for (i = 1, f[0] = 0, p[0] = -1, msc_i = a_n - 1; i < a_n; i++) {
j = i-1;
dq = (int64_t)(a[i].self_offset) - (int64_t)(a[j].self_offset);
dr = (int64_t)(a[i].offset) - (int64_t)(a[j].offset);
dd = dr > dq? dr - dq : dq - dr;//gap
if(dd <= long_gap || dq > max_dis) {
p[i] = i - 1; f[i] = i;
} else {
break;
}
}
if(i >= a_n) goto ss_kip;
memset(t, 0, (a_n*sizeof((*t))));
f[0] = 0; p[0] = -1;
for (i = 1, st = 0; i < a_n; ++i) {
max_f = INT32_MIN; n_skip = 0; max_j = -1;
if ((i-st) > max_iter) st = i-max_iter;
///i-1
j = i - 1;
dq = (int64_t)(a[i].self_offset) - (int64_t)(a[j].self_offset);
dr = (int64_t)(a[i].offset) - (int64_t)(a[j].offset);
dd = dr > dq? dr - dq : dq - dr;//gap
if(dd <= long_gap) dd = 0;
sc = f[j] - dd;
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;
///[st, i-2]
for (--j; (j >= st) && (a[i].self_offset <= (max_dis + a[j].self_offset)); --j) {
dq = (int64_t)(a[i].self_offset) - (int64_t)(a[j].self_offset);
dr = (int64_t)(a[i].offset) - (int64_t)(a[j].offset);
dd = dr > dq? dr - dq : dq - dr;//gap
if(dd <= long_gap) dd = 0;
sc = f[j] - dd;
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;
}
f[i] = max_f; p[i] = max_j;
}
i = a_n-1; msc = f[i]; msc_i = i;
for (j = i-1; (j >= 0) && (a[i].self_offset <= (max_dis + a[j].self_offset)); --j) {
if(msc < f[j] && p[j] >= 0) {///hold at least two hits in th final chain
msc = f[j]; msc_i = j;
}
}
ss_kip:
///a[] has been sorted by self_offset
i = msc_i;
cL = 0;
while (i >= 0) {
t[cL++] = 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(des) {
for (i = 0; i < cL; i++) des[i] = a[t[i]];
}
return cL;
}
inline int64_t hit_long_gap(k_mer_hit *a, k_mer_hit *b, int64_t max_lgap, double small_bw_rate, int64_t min_small_bw)
{
@@ -1680,10 +1960,13 @@ int64_t filter_bad_seed_dp(k_mer_hit *sk, k_mer_hit *ek, k_mer_hit* a, int64_t a
uint64_t lchain_dp_trace(k_mer_hit* a, int64_t a_n, int64_t max_lgap, double sgap_rate, int64_t sgap)
{
// fprintf(stderr, "[M::%s::] a_n::%ld\n", __func__, a_n);
if(a_n <= 0) return 0;
int64_t i, st, occ = 0;
for (i = 1, st = 0; i <= a_n; ++i) {
// fprintf(stderr, "[M::%s::i->%ld] q::%u, t::%u, cnt::%u, readID::%u\n", __func__,
// i-1, a[i-1].self_offset, a[i-1].offset, a[i-1].cnt, a[i-1].readID);
if((i == a_n) || (is_alnw(a[i]))) {///[st, i)
if(i > st) {
occ += filter_bad_seed_dp((st>0)?&(a[st-1]):NULL, (i<a_n)?&(a[i]):NULL, a, a_n, max_lgap, sgap_rate, sgap);
Hash_Table.h
+6 -1
View File
@@ -208,5 +208,10 @@ int append_utg_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_
#define is_pri_aln(a) ((((a).readID) == ((uint32_t)(0x7fffffff)))||((a).cnt >= (a).readID))
uint64_t lchain_dp_trace(k_mer_hit* a, int64_t a_n, int64_t max_lgap, double sgap_rate, int64_t sgap);
uint64_t lchain_qdp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, overlap_region* res,
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);
int ovlp_chain_qgen(overlap_region_alloc* ol, overlap_region* t, int64_t xl, int64_t yl, int64_t apend_be, k_mer_hit* hit, int64_t n_hit);
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);
#endif
Levenshtein_distance.h
+22 -15
View File
@@ -542,8 +542,14 @@ inline uint32_t pop_trace(asg16_v *res, uint32_t i, uint16_t *c, uint32_t *len)
// #define MAX_L 2500
///511 -> 32 64-bits
#define MAX_E 1023
#define MAX_L 5120
#define MAX_CNS_E 1023
#define MAX_CNS_L 3072
#define FORCE_CNS_L 256
#define MAX_SIN_E 2047
#define MAX_SIN_L 10000
#define FORCE_SIN_L 512
typedef uint64_t w_sig;
#define bitw (6)
@@ -1226,7 +1232,7 @@ inline void ed_band_cal_extension_##sf##_0_w(char *pstr, int32_t pn, char *tstr,
{\
poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/\
if(k >= 0) {\
if(k == 0) {\
if(tmp_e == INT32_MAX) {\
tmp_e = err;\
for ((k) = 0; (poff) < (pe); (poff)++, (k)++) {\
bd = (k>>bitw); k_bd = (k&bitz);\
@@ -1316,7 +1322,7 @@ inline void ed_band_cal_extension_##sf##_1_w(char *pstr, int32_t pn, char *tstr,
{\
poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/\
if(k >= 0) {\
if(k == 0) {\
if(tmp_e == INT32_MAX) {\
tmp_e = err;\
for ((k) = 0; (poff) < (pe); (poff)++, (k)++) {\
bd = (k>>bitw); k_bd = (k&bitz);\
@@ -1572,7 +1578,7 @@ inline void ed_band_cal_extension_##sf##_0_w_trace(char *pstr, int32_t pn, char
if(!done) {\
poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/\
if(k >= 0) {\
if(k == 0) {\
if(tmp_e == INT32_MAX) {\
tmp_e = err;\
for ((k) = 0; (poff) < (pe); (poff)++, (k)++) {\
bd = (k>>bitw); k_bd = (k&bitz);\
@@ -1690,7 +1696,7 @@ inline void ed_band_cal_extension_##sf##_1_w_trace(char *pstr, int32_t pn, char
if(!done) {\
poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/\
if(k >= 0) {\
if(k == 0) {\
if(tmp_e == INT32_MAX) {\
tmp_e = err;\
for ((k) = 0; (poff) < (pe); (poff)++, (k)++) {\
bd = (k>>bitw); k_bd = (k&bitz);\
@@ -2157,7 +2163,7 @@ inline void ed_band_cal_extension_infi_0_w(char *pstr, int32_t pn, char *tstr, i
{
poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/
if(k >= 0) {
if(k == 0) {
if(tmp_e == INT32_MAX) {
tmp_e = err;
for ((k) = 0; (poff) < (pe); (poff)++, (k)++) {
bd = (k>>bitw); k_bd = (k&bitz);
@@ -2251,7 +2257,7 @@ inline void ed_band_cal_extension_infi_1_w(char *pstr, int32_t pn, char *tstr, i
{
poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/
if(k >= 0) {
if(k == 0) {
if(tmp_e == INT32_MAX) {
tmp_e = err;
for ((k) = 0; (poff) < (pe); (poff)++, (k)++) {
bd = (k>>bitw); k_bd = (k&bitz);
@@ -2540,7 +2546,7 @@ inline void ed_band_cal_extension_infi_0_w_trace(char *pstr, int32_t pn, char *t
if(!done) {//diff
poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/
if(k >= 0) {
if(k == 0) {
if(tmp_e == INT32_MAX) {
tmp_e = err;
for ((k) = 0; (poff) < (pe); (poff)++, (k)++) {
bd = (k>>bitw); k_bd = (k&bitz);
@@ -2664,7 +2670,7 @@ inline void ed_band_cal_extension_infi_1_w_trace(char *pstr, int32_t pn, char *t
if(!done) {//diff
poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/
if(k >= 0) {
if(k == 0) {
if(tmp_e == INT32_MAX) {
tmp_e = err;
for ((k) = 0; (poff) < (pe); (poff)++, (k)++) {
bd = (k>>bitw); k_bd = (k&bitz);
@@ -3038,7 +3044,7 @@ inline void ed_band_cal_extension_64_0_w(char *pstr, int32_t pn, char *tstr, int
{
poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/
if(k >= 0) {
if(k == 0) {
if(tmp_e == INT32_MAX) {
tmp_e = err;
for (k = 0; poff < pe; poff++, (k)++) {
tmp_e += ((VP>>k)&(1ULL)); tmp_e -= ((VN>>k)&(1ULL));
@@ -3111,7 +3117,7 @@ inline void ed_band_cal_extension_64_1_w(char *pstr, int32_t pn, char *tstr, int
{
poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/
if(k >= 0) {
if(k == 0) {
if(tmp_e == INT32_MAX) {
tmp_e = err;
for (k = 0; poff < pe; poff++, (k)++) {
tmp_e += ((VP>>k)&(1ULL)); tmp_e -= ((VN>>k)&(1ULL));
@@ -3337,7 +3343,7 @@ inline void ed_band_cal_extension_64_0_w_trace(char *pstr, int32_t pn, char *tst
if(!done) {//diff
poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/
if(k >= 0) {
if(k == 0) {
if(tmp_e == INT32_MAX) {
tmp_e = err;
for (k = 0; poff < pe; poff++, (k)++) {
tmp_e += ((VP>>k)&(1ULL)); tmp_e -= ((VN>>k)&(1ULL));
@@ -3351,9 +3357,10 @@ inline void ed_band_cal_extension_64_0_w_trace(char *pstr, int32_t pn, char *tst
if(tmp_e <= (*ez).thre && tmp_e < (*ez).err) {
(*ez).err = tmp_e; (*ez).pe = pe; (*ez).te = i;
}
// fprintf(stderr, "i::%d, (*ez).err::%d, (*ez).pe::%d, (*ez).te::%d, thre::%d, err::%d, tmp_e::%d, poff::%d, k::%d, pe::%d\n",
// i, (*ez).err, (*ez).pe, (*ez).te, thre, err, tmp_e, i-thre, i+thre-pe, pe);
}
}
// fprintf(stderr, "i::%d, tmp_e::%d, err::%d\n", i, tmp_e, err);
Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1;
++i; ++i_bd;
@@ -3446,7 +3453,7 @@ inline void ed_band_cal_extension_64_1_w_trace(char *pstr, int32_t pn, char *tst
if(!done) {//diff
poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/
if(k >= 0) {
if(k == 0) {
if(tmp_e == INT32_MAX) {
tmp_e = err;
for (k = 0; poff < pe; poff++, (k)++) {
tmp_e += ((VP>>k)&(1ULL)); tmp_e -= ((VN>>k)&(1ULL));
anchor.cpp
+178 -2
View File
@@ -17,8 +17,11 @@ typedef struct { // this struct is not strictly necessary; we can use k_mer_pos
#define an_key1(a) ((a).srt)
#define an_key2(a) ((a).self_off)
#define an_key3(a) ((a).other_off)
KRADIX_SORT_INIT(ha_an1, anchor1_t, an_key1, 8)
KRADIX_SORT_INIT(ha_an2, anchor1_t, an_key2, 4)
KRADIX_SORT_INIT(ha_an3, anchor1_t, an_key3, 4)
#define oreg_xs_lt(a, b) (((uint64_t)(a).x_pos_s<<32|(a).x_pos_e) < ((uint64_t)(b).x_pos_s<<32|(b).x_pos_e))
KSORT_INIT(or_xs, overlap_region, oreg_xs_lt)
@@ -872,6 +875,102 @@ void *ha_flt_tab, ha_pt_t *ha_idx, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_
cl->length = ab->n_a;
}
void minimizers_qgen(ha_abufl_t *ab, char* rs, int64_t rl, uint64_t mz_w, uint64_t mz_k, Candidates_list *cl, kvec_t_u8_warp* k_flag,
void *ha_flt_tab, ha_pt_t *ha_idx, All_reads* rdb, const ul_idx_t *udb, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ,
uint32_t *low_occ)
{
// fprintf(stderr, "+[M::%s]\n", __func__);
uint64_t i, k, l, max_cnt = UINT32_MAX, min_cnt = 0; int n, j; ha_mzl_t *z; seedl_t *s;
if(high_occ) {
max_cnt = (*high_occ);
if(max_cnt < 2) max_cnt = 2;
}
if(low_occ) {
min_cnt = (*low_occ);
if(min_cnt < 2) min_cnt = 2;
}
clear_Candidates_list(cl); ab->mz.n = 0, ab->n_a = 0;
// get the list of anchors
mz2_ha_sketch(rs, rl, mz_w, mz_k, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0, NULL);
// minimizer of queried read
if (ab->mz.m > ab->old_mz_m) {
ab->old_mz_m = ab->mz.m;
REALLOC(ab->seed, ab->old_mz_m);
}
for (i = 0, ab->n_a = 0; i < ab->mz.n; ++i) {
ab->seed[i].a = ha_ptl_get(ha_idx, ab->mz.a[i].x, &n);
ab->seed[i].n = n;
ab->n_a += n;
}
if (ab->n_a > ab->m_a) {
ab->m_a = ab->n_a;
REALLOC(ab->a, ab->m_a);
}
for (i = 0, k = 0; i < ab->mz.n; ++i) {
///z is one of the minimizer
z = &ab->mz.a[i]; s = &ab->seed[i];
for (j = 0; j < s->n; ++j) {
const ha_idxposl_t *y = &s->a[j];
anchor1_t *an = &ab->a[k++];
uint8_t rev = z->rev == y->rev? 0 : 1;
an->other_off = rev?((uint32_t)-1)-1-(y->pos+1-y->span):y->pos;
an->self_off = z->pos;
///an->cnt: cnt<<8|span
an->cnt = s->n; if(an->cnt > ((uint32_t)(0xffffffu))) an->cnt = 0xffffffu;
an->cnt <<= 8; an->cnt |= ((z->span <= ((uint32_t)(0xffu)))?z->span:((uint32_t)(0xffu)));
an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->self_off;
}
}
// copy over to _cl_
if (ab->m_a >= (uint64_t)cl->size) {
cl->size = ab->m_a;
REALLOC(cl->list, cl->size);
}
k_mer_hit *p; uint64_t tid = (uint64_t)-1, tl = (uint64_t)-1;
radix_sort_ha_an1(ab->a, ab->a + ab->n_a);
for (k = 1, l = 0; k <= ab->n_a; ++k) {
if (k == ab->n_a || ab->a[k].srt != ab->a[l].srt) {
if (k-l>1) radix_sort_ha_an3(ab->a+l, ab->a+k);
if((ab->a[l].srt>>33)!=tid) {
tid = ab->a[l].srt>>33;
tl = rdb?Get_READ_LENGTH((*rdb), tid):udb->ug->u.a[tid].len;
}
for (i = l; i < k; i++) {
p = &cl->list[i];
p->readID = ab->a[i].srt>>33;
p->strand = (ab->a[i].srt>>32)&1;
if(!(p->strand)) {
p->offset = ab->a[i].other_off;
} else {
p->offset = ((uint32_t)-1)-ab->a[i].other_off;
p->offset = tl-p->offset;
}
p->self_offset = ab->a[i].self_off;
if(((ab->a[i].cnt>>8) < max_cnt) && ((ab->a[i].cnt>>8) > min_cnt)){
p->cnt = 1;
} else if((ab->a[i].cnt>>8) <= min_cnt) {
p->cnt = 2;
} else{
p->cnt = 1 + (((ab->a[i].cnt>>8) + (max_cnt<<1) - 1)/(max_cnt<<1));
p->cnt = pow(p->cnt, 1.1);
}
if(p->cnt > ((uint32_t)(0xffffffu))) p->cnt = 0xffffffu;
p->cnt <<= 8; p->cnt |= (((uint32_t)(0xffu))&(ab->a[i].cnt));
}
l = k;
}
}
cl->length = ab->n_a;
}
void inline reverse_k_mer_hit(k_mer_hit *a, uint64_t a_n, uint64_t xl, uint64_t yl)
{
uint64_t z, han = a_n>>1; k_mer_hit *ai, *aj, ka;
@@ -1009,6 +1108,81 @@ void lchain_gen(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, uin
ks_introsort_or_xs(ol->length, ol->list);
}
void lchain_qgen(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, uint64_t rl, All_reads* rdb,
const ul_idx_t *udb, uint32_t apend_be, overlap_region* tf, uint64_t max_n_chain,
int64_t max_skip, int64_t max_iter, int64_t max_dis, double chn_pen_gap, double chn_pen_skip, double bw_rate, int64_t quick_check, uint32_t gen_off)
{
// fprintf(stderr, "+[M::%s]\n", __func__);
uint64_t i, k, l, m, sm, cn = cl->length; overlap_region *r; ///srt = 0
clear_overlap_region_alloc(ol);
clear_fake_cigar(&(tf->f_cigar));
for (l = 0, k = 1, m = 0; k <= cn; k++) {
if((k == cn) || (cl->list[k].readID != cl->list[l].readID)
|| (cl->list[k].strand != cl->list[l].strand)) {
if(cl->list[l].readID != rid) {
tf->x_id = rid;
tf->x_pos_strand = cl->list[l].strand;
tf->y_id = cl->list[l].readID;
tf->y_pos_strand = 0;///always 0
// fprintf(stderr, "+[M::%s] l::%lu, k::%lu\n", __func__, l, k);
sm = lchain_qdp(cl->list+l, k-l, cl->list+m, &(cl->chainDP), tf, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_rate,
rl, rdb?Get_READ_LENGTH((*rdb), (*tf).y_id):udb->ug->u.a[(*tf).y_id].len, quick_check);
// assert(sm > 0);
if(ovlp_chain_qgen(ol, tf, rl, rdb?Get_READ_LENGTH((*rdb), (*tf).y_id):udb->ug->u.a[(*tf).y_id].len, apend_be, cl->list+m, sm)) {
r = &(ol->list[ol->length-1]); r->non_homopolymer_errors = m;
// if(tf->y_id == 66 || tf->y_id == 66) {
// fprintf(stderr, "\n[M::%s::] utg%.6dl(%c), i::%lu\n",
// __func__, (int32_t)tf->y_id+1, "+-"[tf->x_pos_strand], m);
// }
// reset_k_mer_hit(cl->list+m, sm, rl, rdb?Get_READ_LENGTH((*rdb), r->y_id):udb->ug->u.a[r->y_id].len, r->y_pos_strand, &(ol->length));
for (i = 0; i < sm; i++) {
cl->list[m+i].readID = ol->length;
// if(tf->y_id == 126) fprintf(stderr, "[M::%s::qoff->%u::toff->%u]\n", __func__, cl->list[m+i].self_offset, cl->list[m+i].offset);
}
if(gen_off) gen_fake_cigar(&(r->f_cigar), r, apend_be, cl->list+m, sm);
m += sm;
}
}
l = k;
}
}
cl->length = m;
k = ol->length;
if (ol->length > max_n_chain) {
int32_t w, n[4], s[4]; overlap_region t;
n[0] = n[1] = n[2] = n[3] = 0, s[0] = s[1] = s[2] = s[3] = 0;
ks_introsort_or_ss(ol->length, ol->list); ///srt = 1;
for (i = 0; i < ol->length; ++i) {
r = &(ol->list[i]);
w = ha_ov_type(r, rl);
++n[w];
if (((uint64_t)n[w]) == max_n_chain) s[w] = r->shared_seed;
}
if (s[0] > 0 || s[1] > 0 || s[2] > 0 || s[3] > 0) {
// n[0] = n[1] = n[2] = n[3] = 0;
for (i = 0, k = 0; i < ol->length; ++i) {
r = &(ol->list[i]);
w = ha_ov_type(r, rl);
// ++n[w];
// if (((int)n[w] <= max_n_chain) || (r->shared_seed >= s[w] && s[w] >= (asm_opt.k_mer_length<<1))) {
if (r->shared_seed >= s[w]) {
if (k != i) {
t = ol->list[k];
ol->list[k] = ol->list[i];
ol->list[i] = t;
}
++k;
}
}
ol->length = k;
}
}
ks_introsort_or_xs(ol->length, ol->list);
}
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)
{
double div, pen_gap, pen_skip, tmp;
@@ -1031,7 +1205,9 @@ void ul_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t
extern ha_pt_t *ha_idx;
int64_t max_skip, max_iter, max_dis, quick_check; double chn_pen_gap, chn_pen_skip;
set_lchain_dp_op(is_accurate, mz_k, &max_skip, &max_iter, &max_dis, &chn_pen_gap, &chn_pen_skip, &quick_check);
minimizers_gen(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, dbg_ct, sp, high_occ, low_occ);
lchain_gen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off);
// minimizers_gen(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, dbg_ct, sp, high_occ, low_occ);
minimizers_qgen(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, NULL, uref, dbg_ct, sp, high_occ, low_occ);
// lchain_gen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off);
lchain_qgen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off);
///no need to sort here, overlap_list has been sorted at lchain_gen
}
inter.cpp
+3 -3
View File
@@ -6050,7 +6050,7 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call
// if(s->id+i!=41927 && s->id+i!=47072 && s->id+i!=67641 && s->id+i!=90305 && s->id+i!=698342 && s->id+i!=329421) {
// return;
// }
if((s->id+i!=1) /**&& (s->id+i!=44) && (s->id+i!=948)**/) return;
// if((s->id+i!=291) /**&& (s->id+i!=44) && (s->id+i!=948)**/) return;
// fprintf(stderr, "\n[M::%s] rid::%ld, len::%lu, name::%.*s\n", __func__, s->id+i, s->len[i],
// (int32_t)UL_INF.nid.a[s->id+i].n, UL_INF.nid.a[s->id+i].a);
@@ -6073,7 +6073,7 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call
// memset(&b->self_read, 0, sizeof(b->self_read));
ul_lalign(&b->olist, &b->clist, s->uu, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read,
&b->correct, &b->exz, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, NULL, s->id+i, s->opt->k, NULL);
&b->correct, &b->exz, &b->exz1, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, NULL, s->id+i, s->opt->k, NULL);
// ul_lalign_old_ed(&b->olist, &b->clist, s->uu, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read,
// &b->correct, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, 1, NULL);
ton = b->olist.length;//all alignments pass similary check
@@ -6092,7 +6092,7 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call
copy_asg_arr(b->hap.snp_srt, b0); copy_asg_arr(s->sps[tid], b1); copy_asg_arr(b->r_buf.a, b2);
ul_lalign(&b->olist, &b->clist, s->uu, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read,
&b->correct, &b->exz, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, &(bl->lo), s->id+i, s->opt->k, NULL);
&b->correct, &b->exz, &b->exz1, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, &(bl->lo), s->id+i, s->opt->k, NULL);
// ul_lalign_old_ed(&b->olist, &b->clist, s->uu, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read,
// &b->correct, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, 0, NULL);
}
main.cpp
+12 -11
View File
@@ -12,16 +12,18 @@ int main(int argc, char *argv[])
yak_reset_realtime();
init_opt(&asm_opt);
if (!CommandLine_process(argc, argv, &asm_opt)) return 0;
/**
bit_extz_t exz, exz64; init_bit_extz_t(&exz, 2); init_bit_extz_t(&exz64, 2);
char *pstr = "AAATATT", *tsrt = "AATATTTT"; int32_t thre = 1;
ed_band_cal_semi_64_w_absent_diag((char*)pstr, strlen(pstr), (char*)tsrt, strlen(tsrt), thre, 0, &exz);
fprintf(stderr, "\n[M::%s::] exz.err::%d, exz.ps::%d, exz.pe::%d, exz.ts::%d, exz.te::%d\n", __func__,
exz.err, exz.ps, exz.pe, exz.ts, exz.te);
ed_band_cal_extension_64_0_w((char*)pstr, strlen(pstr), (char*)tsrt, strlen(tsrt), thre, &exz);
fprintf(stderr, "\n[M::%s::] exz.err::%d, exz.ps::%d, exz.pe::%d, exz.ts::%d, exz.te::%d\n", __func__,
exz.err, exz.ps, exz.pe, exz.ts, exz.te);
// bit_extz_t exz, exz64; init_bit_extz_t(&exz, 2); init_bit_extz_t(&exz64, 2);
// char *pstr = "GACCCAG", *tsrt = "GTTGTTAATTCCAT"; int32_t thre = 14; clear_align(exz); clear_align(exz64);
// ed_band_cal_extension_64_0_w_trace((char*)pstr, strlen(pstr), (char*)tsrt, strlen(tsrt), thre, &exz);
// // // ed_band_cal_semi_64_w_absent_diag((char*)pstr, strlen(pstr), (char*)tsrt, strlen(tsrt), thre, 0, &exz);
// fprintf(stderr, "\n[M::%s::] exz.err::%d, exz.ps::%d, exz.pe::%d, exz.ts::%d, exz.te::%d\n", __func__,
// exz.err, exz.ps, exz.pe, exz.ts, exz.te);
// cigar_check((char*)pstr, (char*)tsrt, &(exz));
// ed_band_cal_extension_64_0_w((char*)pstr, strlen(pstr), (char*)tsrt, strlen(tsrt), thre, &exz);
// fprintf(stderr, "\n[M::%s::] exz.err::%d, exz.ps::%d, exz.pe::%d, exz.ts::%d, exz.te::%d\n", __func__,
// exz.err, exz.ps, exz.pe, exz.ts, exz.te);
// ed_band_cal_semi_infi_w((char*)pstr, strlen(pstr), (char*)tsrt, strlen(tsrt), thre, NULL, &exz);
// ed_band_cal_semi_64_w((char*)pstr, strlen(pstr), (char*)tsrt, strlen(tsrt), thre, &exz64);
// fprintf(stderr, "\n[M::%s::] exz.err::%d, exz64.err::%d, exz.ps::%d, exz64.ps::%d, exz.pe::%d, exz64.pe::%d, exz.ts::%d, exz64.ts::%d, exz.te::%d, exz64.te::%d\n", __func__,
@@ -47,8 +49,7 @@ int main(int argc, char *argv[])
// ed_band_cal_extension_128bit((char *)"AAGTTTA", 7, (char *)"CCTTTTTT", 8, 4, &exz);
// ed_band_cal_extension_128bit((char *)"AA", 2, (char *)"ACTTTTTT", 8, 1, &exz);
// fprintf(stderr, "ed_extension::%d, pe::%d, te::%d\n", exz.err, exz.pe, exz.te);
exit(1);
**/
// exit(1);