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

Merge pull request #380 from chhylp123/hifiasm_dev_debug

Browse Source 
Hifiasm dev debug
This commit is contained in:
chhylp123
2023-01-12 01:23:04 -05:00
committed by GitHub
parent 34948a093c 2bf5119bb5
commit 2969af129d
19 changed files with 3259 additions and 467 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.cpp
+16 -6
View File
@@ -51,6 +51,7 @@ static ko_longopt_t long_options[] = {
{ "dbg-het-cnt", ko_no_argument, 337},
{ "ul-tip", ko_required_argument, 338},
{ "low-het", ko_no_argument, 339},
{ "s-base", ko_required_argument, 340},
{ 0, 0, 0 }
};
@@ -120,7 +121,7 @@ void Print_H(hifiasm_opt_t* asm_opt)
fprintf(stderr, " Purge-dups:\n");
fprintf(stderr, " -l INT purge level. 0: no purging; 1: light; 2/3: aggressive [0 for trio; 3 for unzip]\n");
fprintf(stderr, " -s FLOAT similarity threshold for duplicate haplotigs [%g for -l1/-l2, %g for -l3]\n",
fprintf(stderr, " -s FLOAT similarity threshold for duplicate haplotigs in read-level [%g for -l1/-l2, %g for -l3]\n",
asm_opt->purge_simi_rate_l2, asm_opt->purge_simi_rate_l3);
fprintf(stderr, " -O INT min number of overlapped reads for duplicate haplotigs [%d]\n",
asm_opt->purge_overlap_len);
@@ -134,7 +135,9 @@ void Print_H(hifiasm_opt_t* asm_opt)
fprintf(stderr, " --h1 FILEs file names of Hi-C R1 [r1_1.fq,r1_2.fq,...]\n");
fprintf(stderr, " --h2 FILEs file names of Hi-C R2 [r2_1.fq,r2_2.fq,...]\n");
fprintf(stderr, " --seed INT RNG seed [%lu]\n", asm_opt->seed);
fprintf(stderr, " --s-base FLOAT\n");
fprintf(stderr, " similarity threshold for homology detection in base-level;\n");
fprintf(stderr, " -1 to disable [%.3g]; -s for read-level (see <Purge-dups>)\n", asm_opt->trans_base_rate_sec);
fprintf(stderr, " --n-weight INT\n");
fprintf(stderr, " rounds of reweighting Hi-C links [%d]\n", asm_opt->n_weight);
@@ -215,7 +218,8 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->purge_level_trio = 0;
asm_opt->purge_simi_rate_l2 = 0.75;
asm_opt->purge_simi_rate_l3 = 0.55;
asm_opt->trans_base_rate = 0.8;
asm_opt->trans_base_rate = 0.93;
asm_opt->trans_base_rate_sec = 0.5;
asm_opt->purge_overlap_len = 1;
///asm_opt->purge_overlap_len_hic = 50;
asm_opt->recover_atg_cov_min = -1024;
@@ -252,6 +256,10 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->ul_ec_round = 3;
asm_opt->is_dbg_het_cnt = 0;
asm_opt->is_low_het_ul = 0;
asm_opt->is_base_trans = 1;
asm_opt->is_read_trans = 1;
asm_opt->is_topo_trans = 1;
asm_opt->is_bub_trans = 1;
}
void destory_enzyme(enzyme* f)
@@ -771,7 +779,10 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
else if (c == 337) asm_opt->is_dbg_het_cnt = 1;
else if (c == 338) asm_opt->max_short_ul_tip = atol(opt.arg);
else if (c == 339) asm_opt->is_low_het_ul = 1;
else if (c == 'l')
else if (c == 340) {
asm_opt->trans_base_rate_sec = atof(opt.arg);
if(asm_opt->trans_base_rate_sec < 0) asm_opt->is_base_trans = 0;
} 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);
}
@@ -788,10 +799,9 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
return 1;
}
}
if(asm_opt->purge_level_primary > 2) asm_opt->purge_simi_thres = asm_opt->purge_simi_rate_l3;
else asm_opt->purge_simi_thres = asm_opt->purge_simi_rate_l2;
if (argc == opt.ind)
{
CommandLines.h
+6 -1
View File
@@ -4,7 +4,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.18.2-r467"
#define HA_VERSION "0.18.4-r496"
#define VERBOSE 0
@@ -101,6 +101,7 @@ typedef struct {
float purge_simi_rate_l3;
float purge_simi_thres;
float trans_base_rate;
float trans_base_rate_sec;
///float purge_simi_rate_hic;
@@ -130,6 +131,10 @@ typedef struct {
int32_t ul_ec_round;
uint8_t is_dbg_het_cnt;
uint8_t is_low_het_ul;
uint8_t is_base_trans;
uint8_t is_read_trans;
uint8_t is_topo_trans;
uint8_t is_bub_trans;
} hifiasm_opt_t;
extern hifiasm_opt_t asm_opt;
Correct.cpp
+172 -1
View File
@@ -3485,6 +3485,24 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea
///note!!! need notification
end_site = Reserve_Banded_BPM_PATH(y_string, Window_Len, x_string, x_len, threshold, &error, &real_y_start,
&(dumy->path_length), dumy->matrix_bit, dumy->path, p->error, p->y_end - y_start);
// if(!(error != (unsigned int)-1)) {
// fprintf(stderr, "[M::%s::]\tqid::%u\tqlen::%lu\tq::[%d,\t%d)\ttid::%u\ttlen::%lu\tt::[%d,\t%d)\te_beg::%d\te_end::%d\terr::%d\n", __func__,
// overlap_list->list[j].x_id, Get_READ_LENGTH((*rref), overlap_list->list[j].x_id),
// p->x_start, p->x_end+1,
// overlap_list->list[j].y_id, Get_READ_LENGTH((*rref), overlap_list->list[j].y_id),
// p->y_start, p->y_end+1,
// p->extra_begin, p->extra_end, p->error);
// fprintf(stderr, "[M::%s::]\tqcal_len::%lld\ttcal_len::%lld\tthres::%d\n", __func__,
// x_len, Window_Len, threshold);
// fprintf(stderr, "qid::%u\nqname::%.*s\n\t%.*s\n", overlap_list->list[j].x_id,
// (int32_t)Get_NAME_LENGTH((*rref), overlap_list->list[j].x_id),
// Get_NAME((*rref), overlap_list->list[j].x_id), (int32_t)x_len, x_string);
// fprintf(stderr, "tid::%u\ntname::%.*s\n\t%.*s\n", overlap_list->list[j].y_id,
// (int32_t)Get_NAME_LENGTH((*rref), overlap_list->list[j].y_id),
// Get_NAME((*rref), overlap_list->list[j].y_id), (int32_t)Window_Len, y_string);
// }
assert(error != (unsigned int)-1);
{
@@ -14918,6 +14936,13 @@ int64_t ql, int64_t tl, double e_rate, int64_t h_khit, int64_t mode, int64_t rid
ws = qs; if(ws < z->x_pos_s) ws = z->x_pos_s;
we = qe-1; if(we > z->x_pos_e) we = z->x_pos_e;
wsk = get_win_id_by_s(z, ((ws/wl)*wl), wl, NULL);
// if(rid == 7) {
// fprintf(stderr, "[M::%s::]\tutg%.6u%c\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\tq::[%ld,%ld)\tw::[%ld,%ld]\twn::%ld\twsk::%ld\n", __func__,
// z->x_id+1, "lc"[uref->ug->u.a[z->x_id].circ],
// z->x_pos_s, z->x_pos_e+1, "+-"[z->y_pos_strand],
// z->y_id+1, "lc"[uref->ug->u.a[z->y_id].circ],
// z->y_pos_s, z->y_pos_e+1, qs, qe, ws, we, wn, wsk);
// }
// assert((ws>=z->w_list.a[wsk].x_start) && (ws<=z->w_list.a[wsk].x_end));
// wek = get_win_id_by_e(z, ((we/wl)*wl), wl, NULL);
// assert((we>=z->w_list.a[wek].x_start) && (we<=z->w_list.a[wek].x_end));
@@ -18854,6 +18879,61 @@ kv_rtrace_t *trace, uint64_t ql, uint64_t rid, uint64_t oid, uint64_t khit, over
}
}
int64_t ul_raw_aln(overlap_region *z, Candidates_list *cl, const ul_idx_t *udb, char* qstr, UC_Read *tu,
bit_extz_t *exz, double e_rate, int64_t w_l, uint64_t ql, uint64_t rid, uint64_t khit, overlap_region *aux_o)
{
return_t_chain(z, cl);
int64_t ch_idx = z->shared_seed, ch_n;
int64_t i, tl, id = z->y_id, m; ul_ov_t ov;
k_mer_hit *ch_a = cl->list + ch_idx;
tl = udb->ug->u.a[id].len;
for (i = ch_idx; i < cl->length && cl->list[i].readID == cl->list[ch_idx].readID; i++); ch_n = i-ch_idx;
// on = fusion_chain_ovlp(z, ch_a, ch_n, ov, on, wl, ql, tl);
aux_o->w_list.n = aux_o->w_list.c.n = 0;
aux_o->y_id = z->y_id; aux_o->y_pos_strand = z->y_pos_strand;
aux_o->x_pos_s = z->x_pos_s; aux_o->x_pos_e = z->x_pos_e;
aux_o->y_pos_s = z->y_pos_s; aux_o->y_pos_e = z->y_pos_e;
memset(&ov, 0, sizeof(ov)); ov.sec = 3;
ov.ts = ov.te = (uint32_t)-1;
ov.qs = 0; ov.qn = (uint32_t)-1; //extension to left
ov.qe = ql; ov.tn = ch_n; //extension to right
ovlp_base_direct(z, ch_a, ch_n, &ov, w_l, udb, NULL, NULL, qstr, tu, exz, aux_o, e_rate, ql, tl, rid);
int64_t aux_n = aux_o->w_list.n;
for (i = 0; i < aux_n; i++) {
if(!(is_ualn_win(aux_o->w_list.a[i]))) continue;
//will overwrite ch_a; does not matter
rechain_aln(z, cl, aux_o, i, w_l, udb, NULL, NULL, qstr, tu, exz, e_rate, ql, tl, khit, rid);
}
if(((int64_t)aux_o->w_list.n) > aux_n) {
for (i = m = 0; i < ((int64_t)aux_o->w_list.n); i++) {
if((i < aux_n) && (is_ualn_win(aux_o->w_list.a[i]))) continue;
aux_o->w_list.a[m++] = aux_o->w_list.a[i];
}
aux_o->w_list.n = m;
radix_sort_window_list_xs_srt(aux_o->w_list.a, aux_o->w_list.a+aux_o->w_list.n);
}
///update z by aux_o
update_overlap_region(z, aux_o, ql, tl);
int64_t zwn = z->w_list.n, zerr = 0, zlen = z->x_pos_e+1-z->x_pos_s;
for (i = 0; i < zwn; i++) {
if(is_ualn_win(z->w_list.a[i])) zerr += z->w_list.a[i].x_end+1-z->w_list.a[i].x_start;
else zerr += z->w_list.a[i].error;
}
z->non_homopolymer_errors = zerr;
if(zerr >= zlen) return 0;
if(zerr <= 0 && zlen > 0) return 1;
if(zerr > (zlen*e_rate)) return 0;
return 1;
}
void dedup_ul_ov_t(kv_ul_ov_t *in)
{
ul_ov_t *a = in->a; int64_t k, l, m, z, r, a_n = in->n; uint64_t qo, to; double rr = 0.95;
@@ -18998,4 +19078,95 @@ void ul_raw_lalign_adv(overlap_region_alloc* ol, Candidates_list *cl, const ul_i
// copy_asg_arr(hap->snp_srt, iidx); copy_asg_arr(v_idx->a, buf); copy_asg_arr((*stb), buf1);
// }
}
**/
**/
uint64_t gen_nkhits(Candidates_list *cl, overlap_region *z)
{
uint64_t pid, kn; int64_t k;
pid = cl->list[z->shared_seed].readID;kn = 0;
for (k = z->shared_seed; k < cl->length && cl->list[k].readID == pid; k++) kn++;
return kn;
}
void ug_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *uref, const ug_opt_t *uopt,
char *qstr, uint64_t ql, UC_Read* qu, UC_Read* tu, Correct_dumy* dumy, bit_extz_t *exz,
overlap_region *aux_o, double e_rate, int64_t wl, int64_t sid, uint64_t khit,
uint64_t chain_cut, void *km)
{
uint64_t i, bs, k, ovl, pid, kl, kn; Window_Pool w; double err;
overlap_region t; overlap_region *z;
ol->mapped_overlaps_length = 0;
if(ol->length <= 0) return;
///base alignment
clear_Correct_dumy(dumy, ol, km); err = e_rate;
init_Window_Pool(&w, ql, wl, (int)(1.0/err));
bs = (w.window_length)+(THRESHOLD_MAX_SIZE<<1)+1;
resize_UC_Read(tu, bs<<1);
if(!aux_o) {
resize_UC_Read(qu, ql); qu->length = ql; memcpy(qu->seq, qstr, ql);
for (i = 0; i < ol->length; i++) {
z = &(ol->list[i]); ovl = z->x_pos_e+1-z->x_pos_s;
z->shared_seed = z->non_homopolymer_errors;///for index
pid = cl->list[z->shared_seed].readID; kl = cl->length; kn = 0;
for (k = z->shared_seed; k < kl && cl->list[k].readID == pid && kn < chain_cut; k++) kn++;
if(kn < chain_cut) continue;
if(!align_ul_ed_post_extz(z, uref, NULL, qu->seq, tu->seq, exz, err, w.window_length, -1, 0, km)) {
continue;
}
if(uref && simi_pass(ovl, z->align_length, uref?1:0, -1, NULL)) {
z->is_match = 3; ol->mapped_overlaps_length += z->align_length;
}
}
// if(uref && ol->mapped_overlaps_length > 0) {
// set_herror_win(ol, dumy, v_idx, err, ql, w.window_length);
// }
double e_max = err*1.5, rr; int64_t re;
for (i = k = 0; i < ol->length; i++) {
z = &(ol->list[i]); ovl = z->x_pos_e + 1 - z->x_pos_s;
rr = gen_extend_err_exz(z, uref, NULL, NULL, qu->seq, tu->seq, exz, NULL/**v_idx?v_idx->a.a:NULL**/, w.window_length, -1, err, (e_max+0.000001), &re);
z->is_match = 0;///must be here;
// if(sid == 7 && z->y_id == 135) {
// fprintf(stderr, "===utg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\talign_length::%u\terr::%f\trr::%f\twn::%u\n",
// z->x_id+1, "lc"[uref->ug->u.a[z->x_id].circ], uref->ug->u.a[z->x_id].len,
// z->x_pos_s, z->x_pos_e+1, "+-"[z->y_pos_strand],
// z->y_id+1, "lc"[uref->ug->u.a[z->y_id].circ], uref->ug->u.a[z->y_id].len,
// z->y_pos_s, z->y_pos_e+1, z->align_length, err, rr, (uint32_t)z->w_list.n/**gen_nkhits(cl, z)**/);
// }
if (rr <= err) {
if(k != i) {
t = ol->list[k];
ol->list[k] = ol->list[i];
ol->list[i] = t;
}
ol->list[k].is_match = 1; ol->list[k].non_homopolymer_errors = re;
k++;
}
}
ol->length = k;
if(ol->length <= 0) return;
} else {
if(ol->length <= 1) return;
for (i = k = 0; i < ol->length; i++) {
z = &(ol->list[i]);
if(!ul_raw_aln(z, cl, uref, qu->seq, tu, exz, err, wl, ql, sid, khit, aux_o)) {
z->is_match = 0;
continue;
}
if(k != i) {
t = ol->list[k];
ol->list[k] = ol->list[i];
ol->list[i] = t;
}
ol->list[k].is_match = 1;
k++;
}
ol->length = k;
}
}
Correct.h
+5
View File
@@ -1160,6 +1160,11 @@ void lchain_align(overlap_region_alloc* overlap_list, const ul_idx_t *uref,
kvec_t_u64_warp* v_idx, window_list_alloc* win_ciagr_buf,
int force_repeat, int is_consensus, int* fully_cov, int* abnormal,
double max_ov_diff_ec, long long winLen, void *km);
void ug_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *uref, const ug_opt_t *uopt,
char *qstr, uint64_t ql, UC_Read* qu, UC_Read* tu, Correct_dumy* dumy, bit_extz_t *exz,
overlap_region *aux_o, double e_rate, int64_t wl, int64_t sid, uint64_t khit, uint64_t chain_cut,
void *km);
/***
type:
0. match
Hash_Table.cpp
+81 -63
View File
@@ -1823,8 +1823,10 @@ uint64_t lchain_qdp_mcopy(Candidates_list *cl, int64_t a_idx, int64_t a_n, int64
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 gen_cigar, int64_t enable_mcopy, double mcopy_rate, int64_t mcopy_khit_cutoff,
int64_t khit_n)
{
if(a_n <= 0) return 0;
int64_t *p, *t, max_f, n_skip, st, max_j, end_j, sc, msc, msc_i, bw, max_ii, ovl, movl, plus = 0, min_sc, ch_n;
int32_t *f, max, tmp, *ii; int64_t i, k, j, cL = 0; k_mer_hit* a; k_mer_hit* des; k_mer_hit *swap; overlap_region *z;
resize_Chain_Data(dp, a_n, NULL);
@@ -1884,86 +1886,102 @@ uint64_t lchain_qdp_mcopy(Candidates_list *cl, int64_t a_idx, int64_t a_n, int64
}
}
if(f[i] < plus) plus = f[i];
ii[i] = 0;///for mcopy, not here
// if(a_n && a[0].readID == 0) {
// fprintf(stderr, "i::%ld[M::%s::utg%.6dl::%c] x::%u, y::%u, st::%ld, max_ii::%ld, f[i]::%d, p[i]::%ld, msc_i::%ld, msc::%ld, movl::%ld\n",
// i, __func__, (int32_t)a[i].readID+1, "+-"[a[i].strand],
// a[i].self_offset, a[i].offset, st, max_ii, f[i], p[i], msc_i, msc, movl);
// }
}
if((movl < xl) && (movl < yl)) {
msc -= plus; min_sc = msc*0.2;
for (i = ch_n = 0; i < a_n; ++i) {///make all f[] positive
f[i] -= plus; if(i >= ch_n) t[i] = 0;
if(f[i] >= min_sc) {
t[ch_n] = ((uint64_t)f[i])<<32; t[ch_n] += (i<<1); ch_n++;
}
}
int64_t n_v, n_v0, ni, n_u, n_u0 = res->length;
radix_sort_hc64i(t, t + ch_n);
for (k = ch_n-1, n_v = n_u = 0; k >= 0; --k) {
n_v0 = n_v;
for (i = ((uint32_t)t[k])>>1; i >= 0 && (t[i]&1) == 0; ) {
ii[n_v++] = i; t[i] |= 1; i = p[i];
for (i = msc_i, cL = 0; i >= 0; i = p[i]) { ii[i] = 1; t[cL++] = i;}///label the best chain
if((movl < xl) && enable_mcopy/**(movl < yl)**/) {
if(cL >= mcopy_khit_cutoff) {///if there are too few k-mers, disable mcopy
msc -= plus; min_sc = msc*mcopy_rate/**0.2**/; ii[msc_i] = 0;
for (i = ch_n = 0; i < a_n; ++i) {///make all f[] positive
f[i] -= plus; if(i >= ch_n) t[i] = 0;
if((!(ii[i])) && (f[i] >= min_sc)) {
t[ch_n] = ((uint64_t)f[i])<<32; t[ch_n] += (i<<1); ch_n++;
}
}
// if(a_n && a[0].readID == 0) {
// fprintf(stderr, "init_sc::%ld[M::%s::k->%ld] n_v0::%ld, n_v::%ld\n", t[k]>>32, __func__, k, n_v0, n_v);
// }
if(n_v0 == n_v) continue;
sc = (i<0?(t[k]>>32):((t[k]>>32)-f[i]));
if(sc >= min_sc) {
// if(a_n && a[0].readID == 0) {
// fprintf(stderr, "sc::%ld[M::%s::] n_v0::%ld, n_v::%ld, k::%ld, ch_n::%ld, msc::%ld, min_sc::%ld\n",
// sc, __func__, n_v0, n_v, k, ch_n, msc, 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]]));
z->align_length = n_v-n_v0; z->x_id = n_v0;
n_u++;
if(ch_n > 1) {
int64_t n_v, n_v0, ni, n_u, n_u0 = res->length;
radix_sort_hc64i(t, t + ch_n);
for (k = ch_n-1, n_v = n_u = 0; k >= 0; --k) {
n_v0 = n_v;
for (i = ((uint32_t)t[k])>>1; i >= 0 && (t[i]&1) == 0; ) {
ii[n_v++] = i; t[i] |= 1; i = p[i];
}
if(n_v0 == n_v) continue;
sc = (i<0?(t[k]>>32):((t[k]>>32)-f[i]));
if(sc >= min_sc) {
kv_pushp_ol(overlap_region, (*res), &z);
push_ovlp_chain_qgen(z, xid, xl, yl, sc+plus, &(a[ii[n_v-1]]), &(a[ii[n_v0]]));
///mcopy_khit_cutoff <= 1: disable the mcopy_khit_cutoff filtering, for the realignment
if((mcopy_khit_cutoff <= 1) || ((z->x_pos_e+1-z->x_pos_s) <= (movl<<2))) {
z->align_length = n_v-n_v0; z->x_id = n_v0;
n_u++;
} else {///non-best is too long
res->length--; n_v = n_v0;
}
} else {
n_v = n_v0;
}
}
if(n_u > 1) ks_introsort_or_sss(n_u, res->list + n_u0);
res->length = n_u0 + filter_non_ovlp_xchains(res->list + n_u0, n_u, &n_v);
n_u = res->length;
if(n_u > n_u0 + 1) {
kv_resize_cl(k_mer_hit, (*cl), (n_v+cl->length));
a = cl->list + a_idx; des = cl->list + des_idx; swap = cl->list + cl->length;
for (k = n_u0, i = n_v0 = n_v = 0; k < n_u; k++) {
z = &(res->list[k]);
z->non_homopolymer_errors = des_idx + i;
n_v0 = z->x_id; ni = z->align_length;
for (j = 0; j < ni; j++, i++) {
///k0 + (ni - j - 1)
swap[i] = a[ii[n_v0 + (ni- j - 1)]];
swap[i].readID = k;
}
z->x_id = xid;
if(gen_cigar) gen_fake_cigar(&(z->f_cigar), z, apend_be, swap+i-ni, ni);
if(!khit_n) z->align_length = 0;
}
memcpy(des, swap, i*sizeof((*swap))); //assert(i == ch_n);
// fprintf(stderr, "[M::%s::msc->%ld] msc_k_hits::%u, cL::%ld, min_sc::%ld, best_sc::%ld, n_u0_sc::%d, mcopy_rate::%f, # chains::%ld\n",
// __func__, msc, res->list[n_u0].align_length, cL, min_sc, msc+plus, res->list[n_u0].shared_seed,
// mcopy_rate, n_u-n_u0);
} else if(n_u == n_u0 + 1) {
z = &(res->list[n_u0]); k = n_u0; i = 0;
z->non_homopolymer_errors = des_idx + i;
n_v0 = z->x_id; ni = z->align_length;
for (j = 0; j < ni; j++, i++) {
///k0 + (ni - j - 1)
des[i] = a[ii[n_v0 + (ni- j - 1)]];
des[i].readID = k;
}
z->x_id = xid;
if(gen_cigar) gen_fake_cigar(&(z->f_cigar), z, apend_be, des+i-ni, ni);
if(!khit_n) z->align_length = 0;
}
return i;
} else {
n_v = n_v0;
msc += plus; i = msc_i; cL = 0;
while (i >= 0) {t[cL++] = i; i = p[i];}
}
}
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;
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;
// if(n_u > 1) {
// fprintf(stderr, "\nk::%ld[M::%s::utg%.6dl::%c] ni::%ld\n",
// i, __func__, (int32_t)des[i].readID+1, "+-"[des[i].strand], ni);
// }
for (j = 0; j < ni; j++, i++) {
///k0 + (ni - j - 1)
swap[i] = a[ii[n_v0 + (ni- j - 1)]];
swap[i].readID = k;
// if(a_n && a[0].readID == 0) {
// fprintf(stderr, "i::%ld[M::%s::utg%.6dl::%c] x::%u, y::%u, ni::%ld\n",
// i, __func__, (int32_t)swap[i].readID+1, "+-"[swap[i].strand],
// swap[i].self_offset, swap[i].offset, ni);
// }
}
z->x_id = xid;
if(gen_cigar) gen_fake_cigar(&(z->f_cigar), z, apend_be, swap+i-ni, ni);
z->align_length = 0;
}
memcpy(des, swap, i*sizeof((*swap))); //assert(i == ch_n);
return i;
}
///a[] has been sorted by self_offset
i = msc_i; cL = 0;
while (i >= 0) {t[cL++] = i; i = p[i];}
// 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;
}
Hash_Table.h
+2 -1
View File
@@ -233,6 +233,7 @@ uint64_t lchain_qdp_mcopy(Candidates_list *cl, int64_t a_idx, int64_t a_n, int64
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 gen_cigar, int64_t enable_mcopy, double mcopy_rate, int64_t mcopy_khit_cutoff,
int64_t khit_n);
#endif
Overlaps.cpp
+906 -205
View File
File diff suppressed because it is too large Load Diff
Overlaps.h
+3
View File
@@ -33,6 +33,7 @@
// #define ALTER_LABLE 2
// #define HAP_LABLE 4
#define HA_RE_UL_ID "re"
#define ug_ext_len 75000
#define Get_qn(RECORD) ((uint32_t)((RECORD).qns>>32))
#define Get_qs(RECORD) ((uint32_t)((RECORD).qns))
@@ -1135,6 +1136,8 @@ asg_t *i_read_sg, trans_chain* i_t_ch, uint32_t i_cBeg, uint32_t i_cEnd);
void extract_sub_overlaps(uint32_t i_tScur, uint32_t i_tEcur, uint32_t i_tSpre, uint32_t i_tEpre,
uint32_t tn, kv_u_trans_hit_t* ktb, uint32_t bn);
void clean_u_trans_t_idx(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g);
void clean_u_trans_t_idx_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g);
void clean_u_trans_t_idx_filter_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g);
uint32_t test_dbug(ma_ug_t* ug, FILE* fp);
void write_dbug(ma_ug_t* ug, FILE* fp);
int asg_arc_identify_simple_bubbles_multi(asg_t *g, bub_label_t* x, int check_cross);
anchor.cpp
+528 -13
View File
@@ -5,8 +5,12 @@
#include "ksort.h"
#include "Hash_Table.h"
#include "kalloc.h"
#include "Overlaps.h"
#define HA_KMER_GOOD_RATIO 0.333
#define OFL 0.95
#define CH_OCC 4
#define CH_SC 16
typedef struct { // this struct is not strictly necessary; we can use k_mer_pos instead, with modifications
uint64_t srt;
@@ -21,7 +25,8 @@ typedef struct { // this struct is not strictly necessary; we can use k_mer_pos
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 generic_key(x) (x)
KRADIX_SORT_INIT(anc64, uint64_t, generic_key, 8)
#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)
@@ -29,6 +34,9 @@ KSORT_INIT(or_xs, overlap_region, oreg_xs_lt)
#define oreg_ss_lt(a, b) ((a).shared_seed > (b).shared_seed) // in the decending order
KSORT_INIT(or_ss, overlap_region, oreg_ss_lt)
#define oreg_occ_lt(a, b) ((a).align_length > (b).align_length) // in the decending order
KSORT_INIT(or_occ, overlap_region, oreg_occ_lt)
#define oreg_id_lt(a, b) ((a).y_id < (b).y_id)
KSORT_INIT(or_id, overlap_region, oreg_id_lt)
@@ -971,6 +979,275 @@ uint32_t *low_occ)
cl->length = ab->n_a;
}
void gen_pair_chain(ha_abufl_t *ab, uint64_t rid, st_mt_t *tid, uint64_t tid_n, ha_mzl_t *in, uint64_t in_n, ha_mzl_t *idx, int64_t idx_n, uint64_t mzl_cutoff)
{
if(!tid_n) return;
uint64_t i, l, k, n, m, rev, x, tn, kn; ha_mzl_t *z, *p, *y; int64_t zi, ns, ne; anchor1_t *an;
for (i = tn = 0; i < in_n; ++i) {
///z is one of the minimizer
z = &in[i]; p = &(idx[z->x]); n = 0;
assert(z->pos == p->pos && z->rid == p->rid && z->span == p->span && z->rev == p->rev);
for (zi = z->x+1; zi < idx_n && idx[zi].x == p->x && n < mzl_cutoff; zi++) {
if(idx[zi].rid == rid) continue;
x = idx[zi].rid; x <<= 1; x |= ((uint64_t)((z->rev == idx[zi].rev)?0:1));
for (m = 0; m < tid_n; m++) {
if(tid->a[m] == x) {
n++; break;
}
}
}
for (zi = ((int64_t)z->x)-1; zi >= 0 && idx[zi].x == p->x && n < mzl_cutoff; zi--) {
if(idx[zi].rid == rid) continue;
x = idx[zi].rid; x <<= 1; x |= ((uint64_t)((z->rev == idx[zi].rev)?0:1));
for (m = 0; m < tid_n; m++) {
if(tid->a[m] == x) {
n++; break;
}
}
}
if((!n) || (n >= mzl_cutoff)) continue;
tn += n;
}
if(!tn) return;
ab->n_a += tn;
if (ab->n_a > ab->m_a) {
ab->m_a = ab->n_a;
REALLOC(ab->a, ab->m_a);
}
for (i = 0, k = ab->n_a - tn; i < in_n; ++i) {
///z is one of the minimizer
z = &in[i]; p = &(idx[z->x]); n = 0; ns = 0; ne = -1;
for (zi = z->x+1; zi < idx_n && idx[zi].x == p->x && n < mzl_cutoff; zi++) {
if(idx[zi].rid == rid) continue;
x = idx[zi].rid; x <<= 1; x |= ((uint64_t)((z->rev == idx[zi].rev)?0:1));
for (m = 0; m < tid_n; m++) {
if(tid->a[m] == x) {
n++; break;
}
}
}
ne = zi;
for (zi = ((int64_t)z->x)-1; zi >= 0 && idx[zi].x == p->x && n < mzl_cutoff; zi--) {
if(idx[zi].rid == rid) continue;
x = idx[zi].rid; x <<= 1; x |= ((uint64_t)((z->rev == idx[zi].rev)?0:1));
for (m = 0; m < tid_n; m++) {
if(tid->a[m] == x) {
n++; break;
}
}
}
ns = zi + 1;
if((!n) || (n >= mzl_cutoff)) continue;
n = ne - ns; if(n > ((uint32_t)(0xffffffu))) n = 0xffffffu; n<<=8;
for (zi = z->x+1; zi < idx_n && idx[zi].x == p->x; zi++) {
if(idx[zi].rid == rid) continue;
x = idx[zi].rid; x <<= 1; x |= ((uint64_t)((z->rev == idx[zi].rev)?0:1));
for (m = 0; m < tid_n; m++) {
if(tid->a[m] == x) break;
}
if(m >= tid_n) continue;
y = &(idx[zi]); an = &ab->a[k++];
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 = ((z->span <= ((uint32_t)(0xffu)))?z->span:((uint32_t)(0xffu))); an->cnt |= n;
an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->self_off;
}
for (zi = ((int64_t)z->x)-1; zi >= 0 && idx[zi].x == p->x; zi--) {
if(idx[zi].rid == rid) continue;
x = idx[zi].rid; x <<= 1; x |= ((uint64_t)((z->rev == idx[zi].rev)?0:1));
for (m = 0; m < tid_n; m++) {
if(tid->a[m] == x) break;
}
if(m >= tid_n) continue;
y = &(idx[zi]); an = &ab->a[k++];
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 = ((z->span <= ((uint32_t)(0xffu)))?z->span:((uint32_t)(0xffu))); an->cnt |= n;
an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->self_off;
}
}
assert(k == ab->n_a);
radix_sort_ha_an1(ab->a + ab->n_a - tn, ab->a + ab->n_a); kn = 0;
for (k = ab->n_a - tn + 1, l = ab->n_a - tn; k <= ab->n_a; ++k) {
if (k == ab->n_a || (ab->a[k].srt>>32) != (ab->a[l].srt>>32)) {
for (; kn < tid_n && tid->a[kn] < (ab->a[l].srt>>32); kn++);
if(kn < tid_n && tid->a[kn] == (ab->a[l].srt>>32)) kv_push(uint64_t, *tid, l);
for (; kn < tid_n && tid->a[kn] == (ab->a[l].srt>>32); kn++);
l = k;
}
}
}
void minimizers_qgen_input(ha_abufl_t *ab, uint64_t rid, 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, ha_mzl_t *in, uint64_t in_n, ha_mzl_t *idx, int64_t idx_n, uint64_t mzl_cutoff, uint64_t chain_cutoff, kv_u_trans_t *kov)
{
// fprintf(stderr, "+[M::%s]\n", __func__);
uint64_t i, k, l, max_cnt = UINT32_MAX, min_cnt = 0, n; ha_mzl_t *z, *p, *y; anchor1_t *an;
int64_t zi, ns, ne; uint8_t rev; k_mer_hit *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);
for (i = 0, ab->n_a = 0; i < in_n; ++i) {
///z is one of the minimizer
z = &in[i]; p = &(idx[z->x]); n = 0;
assert(z->pos == p->pos && z->rid == p->rid && z->span == p->span && z->rev == p->rev);
for (zi = z->x+1; zi < idx_n && idx[zi].x == p->x && n < mzl_cutoff; zi++) {
if(idx[zi].rid == rid) continue; n++;
}
for (zi = ((int64_t)z->x)-1; zi >= 0 && idx[zi].x == p->x && n < mzl_cutoff; zi--) {
if(idx[zi].rid == rid) continue; n++;
}
if((!n) || (n >= mzl_cutoff)) continue;
ab->n_a += n;
}
// if(rid == 0) {
// fprintf(stderr, "-0-[M::%s] ab->n_a::%lu, in_n::%lu\n", __func__, ab->n_a, in_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 < in_n; ++i) {
///z is one of the minimizer
z = &in[i]; p = &(idx[z->x]); n = 0; ns = 0; ne = -1;
for (zi = z->x+1; zi < idx_n && idx[zi].x == p->x && n < mzl_cutoff; zi++) {
if(idx[zi].rid == rid) continue; n++;
}
ne = zi;
for (zi = ((int64_t)z->x)-1; zi >= 0 && idx[zi].x == p->x && n < mzl_cutoff; zi--) {
if(idx[zi].rid == rid) continue; n++;
}
ns = zi + 1;
if((!n) || (n >= mzl_cutoff)) continue;
n = ne - ns; if(n > ((uint32_t)(0xffffffu))) n = 0xffffffu; n<<=8;
for (zi = z->x+1; zi < idx_n && idx[zi].x == p->x; zi++) {
if(idx[zi].rid == rid) continue;
y = &(idx[zi]); an = &ab->a[k++];
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 = ((z->span <= ((uint32_t)(0xffu)))?z->span:((uint32_t)(0xffu))); an->cnt |= n;
an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->self_off;
}
for (zi = ((int64_t)z->x)-1; zi >= 0 && idx[zi].x == p->x; zi--) {
if(idx[zi].rid == rid) continue;
y = &(idx[zi]); an = &ab->a[k++];
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 = ((z->span <= ((uint32_t)(0xffu)))?z->span:((uint32_t)(0xffu))); an->cnt |= n;
an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->self_off;
}
}
assert(k == ab->n_a);
uint64_t kn, spn; u_trans_t *ka;
kn = spn = 0; ka = NULL; sp->n = 0;
if(kov) {
kn = u_trans_n(*kov, rid); ka = u_trans_a(*kov, rid);
}
if(kn > 0) {
kv_resize(uint64_t, *sp, kn);
for (k = 0; k < kn; ++k) {
if(ka[k].del) continue;
l = ka[k].tn; l <<= 1; l |= ka[k].rev;
kv_push(uint64_t, *sp, l);
}
if(sp->n > 1) radix_sort_anc64(sp->a, sp->a + sp->n);
}
radix_sort_ha_an1(ab->a, ab->a + ab->n_a); kn = 0;
for (k = 1, l = n = 0, spn = sp->n; k <= ab->n_a; ++k) {
if (k == ab->n_a || (ab->a[k].srt>>32) != (ab->a[l].srt>>32)) {
if(k - l >= chain_cutoff) {
for (; kn < spn && sp->a[kn] < (ab->a[l].srt>>32); kn++);
if(kn < spn && sp->a[kn] == (ab->a[l].srt>>32)) kv_push(uint64_t, *sp, n);
for (; kn < spn && sp->a[kn] == (ab->a[l].srt>>32); kn++) sp->a[kn] = (uint64_t)-1;
for (i = l; i < k; i++) ab->a[n++] = ab->a[i];
}
l = k;
}
}
ab->n_a = n;
for (k = kn = 0; k < spn; k++) {
if(sp->a[k] == (uint64_t)-1) continue;
sp->a[kn++] = sp->a[k];
}
// sp->n = kn;
if(kn > 0) gen_pair_chain(ab, rid, sp, kn, in, in_n, idx, idx_n, mzl_cutoff);
if(spn > 0) {
for (k = spn, kn = 0; k < sp->n; k++) sp->a[kn++] = sp->a[k];
sp->n = kn;
}
// copy over to _cl_
if (ab->n_a >= (uint64_t)cl->size) {
cl->size = ab->n_a;
REALLOC(cl->list, cl->size);
}
uint64_t tid = (uint64_t)-1, tl = (uint64_t)-1;
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++) {
s = &cl->list[i];
s->readID = ab->a[i].srt>>33;
s->strand = (ab->a[i].srt>>32)&1;
if(!(s->strand)) {
s->offset = ab->a[i].other_off;
} else {
s->offset = ((uint32_t)-1)-ab->a[i].other_off;
s->offset = tl-s->offset;
}
s->self_offset = ab->a[i].self_off;
if(((ab->a[i].cnt>>8) < max_cnt) && ((ab->a[i].cnt>>8) > min_cnt)){
s->cnt = 1;
} else if((ab->a[i].cnt>>8) <= min_cnt) {
s->cnt = 2;
} else{
s->cnt = 1 + (((ab->a[i].cnt>>8) + (max_cnt<<1) - 1)/(max_cnt<<1));
s->cnt = pow(s->cnt, 1.1);
}
if(s->cnt > ((uint32_t)(0xffffffu))) s->cnt = 0xffffffu;
s->cnt <<= 8; s->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;
@@ -1192,18 +1469,250 @@ void lchain_qgen(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, ui
void lchain_qgen_mcopy(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, 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)
uint32_t gen_off, double mcopy_rate, uint32_t chain_cutoff, uint32_t mcopy_khit_cut, st_mt_t *sp)
{
// fprintf(stderr, "+[M::%s]\n", __func__);
uint64_t i, k, l, m, cn = cl->length, yid; overlap_region *r; ///srt = 0
uint64_t i, k, l, m, cn = cl->length, yid, ol0, lch; overlap_region *r, t; ///srt = 0
clear_overlap_region_alloc(ol);
for (l = 0, k = 1, m = 0; k <= cn; k++) {
for (l = 0, k = 1, m = 0, lch = 0; k <= cn; k++) {
if((k == cn) || (cl->list[k].readID != cl->list[l].readID)) {
if(cl->list[l].readID != rid) {
yid = cl->list[l].readID; ol0 = ol->length;
m += lchain_qdp_mcopy(cl, l, k-l, m, &(cl->chainDP), ol, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_rate,
rid, rl, rdb?Get_READ_LENGTH((*rdb), yid):udb->ug->u.a[yid].len, quick_check, apend_be, gen_off, 1, mcopy_rate, mcopy_khit_cut, 1);
if((chain_cutoff >= 2) && (!lch)) {
for (i = ol0; (i<ol->length) && (!lch); i++) {
if(ol->list[i].align_length < chain_cutoff) lch = 1;
}
}
}
l = k;
}
}
cl->length = m;
// for (k = 0; k < ol->length; k++) {
// fprintf(stderr, "---[M::%s::utg%.6dl] q[%d, %d), t[%d, %d), khit_off::%u\n", __func__,
// (int32_t)ol->list[k].y_id+1, ol->list[k].x_pos_s, ol->list[k].x_pos_e+1,
// ol->list[k].y_pos_s, ol->list[k].y_pos_e+1, ol->list[k].non_homopolymer_errors);
// }
k = ol->length;
if (ol->length > max_n_chain) {
int32_t w, n[4], s[4];
n[0] = n[1] = n[2] = n[3] = 0, s[0] = s[1] = s[2] = s[3] = 0;
ks_introsort_or_ss(ol->length, ol->list);
for (i = 0; i < ol->length; ++i) {
r = &(ol->list[i]);
w = ha_ov_type(r, rl);
++n[w];
if (((uint64_t)n[w]) == max_n_chain) s[w] = r->shared_seed;
}
if (s[0] > 0 || s[1] > 0 || s[2] > 0 || s[3] > 0) {
// n[0] = n[1] = n[2] = n[3] = 0;
for (i = 0, k = 0, lch = 0; i < ol->length; ++i) {
r = &(ol->list[i]);
w = ha_ov_type(r, rl);
// ++n[w];
// if (((int)n[w] <= max_n_chain) || (r->shared_seed >= s[w] && s[w] >= (asm_opt.k_mer_length<<1))) {
if (r->shared_seed >= s[w]) {
if (k != i) {
t = ol->list[k];
ol->list[k] = ol->list[i];
ol->list[i] = t;
}
if(ol->list[k].align_length < chain_cutoff) lch = 1;
++k;
}
}
ol->length = k;
}
}
ks_introsort_or_xs(ol->length, ol->list);
if(lch) {
//@brief r485
uint64_t zs, ze, rs, re, ob, os, oe, ocn, pp, kn, ms, me; int64_t osc;
for (i = l = 0, cn = cl->length; i < ol->length; ++i) {
if(ol->list[i].align_length < chain_cutoff) {
zs = ol->list[i].x_pos_s; ze = ol->list[i].x_pos_e + 1;
ob = (ze - zs)*OFL; if(ob < 16) ob = 16;
osc = ol->list[i].shared_seed*CH_SC;
ocn = ol->list[i].align_length<<CH_OCC;
for (k = 0; (k < ol->length) && (ze > ol->list[k].x_pos_s); k++) {
if(ol->list[k].align_length < chain_cutoff) continue;
if(ol->list[k].align_length < ocn) continue;
if(ol->list[k].shared_seed < osc) continue;
rs = ol->list[k].x_pos_s; re = ol->list[k].x_pos_e + 1;
os = ((rs>=zs)?rs:zs); oe = ((re<=ze)?re:ze);
if((oe > os) && (oe - os) >= ob) {
m = ol->list[k].non_homopolymer_errors;
pp = cl->list[m].readID; kn = 0;
for (; (m < cn) && (cl->list[m].readID == pp) && (kn < ocn); m++) {
me = cl->list[m].self_offset; ms = me - (cl->list[m].cnt&(0xffu));
if((ms >= os) && (me <= oe)) kn++;
}
if(kn >= ocn) break;
}
}
if((k < ol->length) && (ze > ol->list[k].x_pos_s)) continue;
}
if (l != i) {
t = ol->list[l];
ol->list[l] = ol->list[i];
ol->list[i] = t;
}
l++;
}
// fprintf(stderr, "+[M::%s] rid::%u, ol->length0::%lu, ol->length1::%lu\n", __func__, rid, ol->length, l);
ol->length = l;
/**
//@brief r484
for (i = sp->n = 0; i < ol->length; ++i) {
if(ol->list[i].align_length < chain_cutoff) continue;
os = ol->list[i].x_pos_s; oe = ol->list[i].x_pos_e + 1;
if((sp->n) && (((uint32_t)sp->a[sp->n-1]) >= os)) {
if(oe > ((uint32_t)sp->a[sp->n-1])) {
oe = oe - ((uint32_t)sp->a[sp->n-1]);
sp->a[sp->n-1] += oe;
}
} else {
os = (os<<32)|oe; kv_push(uint64_t, *sp, os);
}
}
for (i = k = 0; i < ol->length; ++i) {
if(ol->list[i].align_length < chain_cutoff) {///ol has been sorted by x_pos_s
r = &(ol->list[i]); rs = r->x_pos_s; re = r->x_pos_e + 1;
rl = re - rs; ovl = 0;
for (m = 0; (m < sp->n) && (re > (sp->a[m]>>32)); m++) {
os = ((rs>=(sp->a[m]>>32))? rs:(sp->a[m]>>32));
oe = ((re<=((uint32_t)sp->a[m]))? re:((uint32_t)sp->a[m]));
if(oe > os) {
ovl += (oe - os); if(ovl >= (rl*0.95)) break;
}
}
if(ovl >= (rl*0.95)) continue;
}
if (k != i) {
t = ol->list[k];
ol->list[k] = ol->list[i];
ol->list[i] = t;
}
ol->list[k++].align_length = 0;
}
// fprintf(stderr, "+[M::%s] ol->length0::%lu, ol->length1::%lu\n", __func__, ol->length, k);
ol->length = k;
**/
}
/**else {
for (i = 0; i < ol->length; ++i) ol->list[i].align_length = 0;
}
**/
for (i = 0; i < ol->length; ++i) ol->list[i].align_length = 0;
}
inline uint64_t special_lchain(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, uint64_t rl, All_reads* rdb,
const ul_idx_t *udb, uint32_t apend_be, int64_t max_skip, int64_t max_iter, int64_t max_dis, double chn_pen_gap,
double chn_pen_skip, double bw_rate, int64_t quick_check, uint32_t gen_off, double mcopy_rate, uint32_t mcopy_khit_cut,
st_mt_t *sp, uint64_t *si, uint64_t m, uint64_t l, uint64_t k)
{
uint64_t z, s, e, yid, ol0 = ol->length, ol1;
if(l >= k) return m;
yid = cl->list[l].readID;
for (z = l; z < k && cl->list[z].strand == cl->list[l].strand; z++);
if(z > l) {
s = l; e = z; ol1 = ol->length; // m0 = m;
m += lchain_qdp_mcopy(cl, s, e-s, m, &(cl->chainDP), ol, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_rate,
rid, rl, rdb?Get_READ_LENGTH((*rdb), yid):udb->ug->u.a[yid].len, quick_check, apend_be, gen_off, 0, mcopy_rate, mcopy_khit_cut, 0);
// if(rid == 7) {
// fprintf(stderr, "+[M::%s::utg%.6ul] inn::[%lu, %lu), (*si)::%lu, ol1::%lu, ol->length::%lu\n",
// __func__, rid+1, s, e, (*si), ol1, ol->length);
// }
if(((*si) < sp->n) && (sp->a[(*si)] == s) && (ol->length > ol1)) {
ol->list[ol->length-1].x_pos_strand = 1;
// fprintf(stderr, "+[M::%s] utg%.6u%c -> utg%.6u%c, n_khits0::%lu, n_khits::%lu\n", __func__, rid+1, "lc"[udb->ug->u.a[rid].circ],
// ol->list[ol->length-1].y_id+1, "lc"[udb->ug->u.a[ol->list[ol->length-1].y_id].circ], e-s, m-m0);
}
for (; (*si) < sp->n && sp->a[(*si)] == s; (*si)++);
}
if(z < k) {
s = z; e = k; ol1 = ol->length; // m0 = m;
m += lchain_qdp_mcopy(cl, s, e-s, m, &(cl->chainDP), ol, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_rate,
rid, rl, rdb?Get_READ_LENGTH((*rdb), yid):udb->ug->u.a[yid].len, quick_check, apend_be, gen_off, 0, mcopy_rate, mcopy_khit_cut, 0);
// if(rid == 7) {
// fprintf(stderr, "-[M::%s::utg%.6ul] inn::[%lu, %lu), (*si)::%lu, ol1::%lu, ol->length::%lu\n",
// __func__, rid+1, s, e, (*si), ol1, ol->length);
// }
if(((*si) < sp->n) && (sp->a[(*si)] == s) && (ol->length > ol1)) {
ol->list[ol->length-1].x_pos_strand = 1;
// fprintf(stderr, "-[M::%s] utg%.6u%c -> utg%.6u%c, n_khits0::%lu, n_khits::%lu\n", __func__, rid+1, "lc"[udb->ug->u.a[rid].circ],
// ol->list[ol->length-1].y_id+1, "lc"[udb->ug->u.a[ol->list[ol->length-1].y_id].circ], e-s, m-m0);
}
for (; (*si) < sp->n && sp->a[(*si)] == s; (*si)++);
}
if(ol->length > ol0) {
overlap_region *p = &(ol->list[ol0]); uint64_t pi = ol0; overlap_region t;
for (z = ol0; z < ol->length; z++) {
if((p->shared_seed < ol->list[z].shared_seed) ||
((p->shared_seed == ol->list[z].shared_seed) && (ol->list[z].x_pos_strand == 1))) {
p = &(ol->list[z]); pi = z;
}
}
// if(rid == 7) {
// fprintf(stderr, ">[M::%s::] utg%.6ul\t->\tutg%.6ul\tflag::%u\n",
// __func__, p->x_id+1, p->y_id+1, p->x_pos_strand);
// }
for (z = s = ol0; z < ol->length; z++) {
if((ol->list[z].x_pos_strand == 0) && (z != pi)) continue;
if(s != z) {
t = ol->list[s];
ol->list[s] = ol->list[z];
ol->list[z] = t;
}
s++;
}
ol->length = s;
}
return m;
}
void lchain_qgen_mcopy_input(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, 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, double bw_thres_sec,
int64_t quick_check, uint32_t gen_off, double mcopy_rate, uint32_t mcopy_khit_cut, st_mt_t *sp)
{
// fprintf(stderr, "+[M::%s]\n", __func__);
uint64_t i, k, l, m, cn = cl->length, yid, si; overlap_region *r; ///srt = 0
clear_overlap_region_alloc(ol);
// if(rid == 160) {
// fprintf(stderr, "[M::%s::utg%.6ul] sp->n::%u\n", __func__, rid+1, (uint32_t)sp->n);
// }
for (l = 0, k = 1, m = 0, si = 0; k <= cn; k++) {
if((k == cn) || (cl->list[k].readID != cl->list[l].readID)) {
if(cl->list[l].readID != rid) {
yid = cl->list[l].readID;
m += lchain_qdp_mcopy(cl, l, k-l, m, &(cl->chainDP), ol, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_rate,
rid, rl, rdb?Get_READ_LENGTH((*rdb), yid):udb->ug->u.a[yid].len, quick_check, apend_be, gen_off);
// if(rid == 160) {
// fprintf(stderr, "+[M::%s::utg%.6lul] [%lu, %lu), m::%lu, ol->length::%lu\n", __func__, yid+1, l, k, m, ol->length);
// }
for (; si < sp->n && sp->a[si] < l; si++);
if(si >= sp->n || sp->a[si] >= k) {///might be unmatched
m += lchain_qdp_mcopy(cl, l, k-l, m, &(cl->chainDP), ol, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_rate,
rid, rl, rdb?Get_READ_LENGTH((*rdb), yid):udb->ug->u.a[yid].len, quick_check, apend_be, gen_off, 0, mcopy_rate, mcopy_khit_cut, 0);
} else {
m = special_lchain(cl, ol, rid, rl, rdb, udb, apend_be, max_skip, max_iter, max_dis, chn_pen_gap,
chn_pen_skip, bw_thres_sec, quick_check, gen_off, mcopy_rate, mcopy_khit_cut, sp, &si, m, l, k);
}
// if(rid == 160) {
// fprintf(stderr, "-[M::%s::utg%.6lul] [%lu, %lu), m::%lu, ol->length::%lu\n", __func__, yid+1, l, k, m, ol->length);
// }
}
l = k;
}
@@ -1234,7 +1743,7 @@ void lchain_qgen_mcopy(Candidates_list* cl, overlap_region_alloc* ol, uint32_t r
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 ((r->shared_seed >= s[w]) || (r->x_pos_strand == 1)) {
if (k != i) {
t = ol->list[k];
ol->list[k] = ol->list[i];
@@ -1265,7 +1774,7 @@ void set_lchain_dp_op(uint32_t is_accurate, uint32_t mz_k, int64_t *max_skip, in
}
void ul_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,
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)
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 chain_cutoff, uint32_t mcopy_khit_cut)
{
extern void *ha_flt_tab;
extern ha_pt_t *ha_idx;
@@ -1276,18 +1785,24 @@ void ul_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t
// lchain_gen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off);
// lchain_qgen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off);
///no need to sort here, overlap_list has been sorted at lchain_gen
lchain_qgen_mcopy(cl, overlap_list, rid, rl, NULL, uref, apend_be, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off);
lchain_qgen_mcopy(cl, overlap_list, rid, rl, NULL, uref, apend_be, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off, mcopy_rate, chain_cutoff, mcopy_khit_cut, sp);
}
void 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,
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, uint32_t is_hpc)
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)
{
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);
if((!overlap_list) || (!cl)) {
ab->mz.n = 0;
mz2_ha_sketch(rs, rl, mz_w, mz_k, 0, is_hpc, &ab->mz, NULL, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0, NULL);
mz2_ha_sketch(rs, rl, mz_w, mz_k, rid, is_hpc, &ab->mz, NULL, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0, NULL);
if(res) memcpy(res, ab->mz.a, ab->mz.n * (sizeof((*(ab->mz.a)))));
return ab->mz.n;
} else {
sp->n = 0;
minimizers_qgen_input(ab, rid, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, NULL, uref, dbg_ct, sp, high_occ, low_occ, res, res_n, idx, idx_n, mzl_cutoff, chain_cutoff, kov);
lchain_qgen_mcopy_input(cl, overlap_list, rid, rl, NULL, uref, apend_be, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, bw_thres_sec, quick_check, gen_off, mcopy_rate, mcopy_khit_cut, sp);
return 0;
}
}
gchain_map.cpp
+2 -2
View File
@@ -18,8 +18,8 @@
#include "Assembly.h"
#include "gchain_map.h"
KSEQ_INIT(gzFile, gzread)
void ul_map_lchain(ha_abufl_t *ab, int64_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, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres,
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);
void ul_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,
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 chain_cutoff, uint32_t mcopy_khit_cut);
typedef struct { // global data structure for kt_pipeline()
const void *ha_flt_tab;
gfa_ut.cpp
+187 -26
View File
@@ -2477,8 +2477,8 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
// (int)Get_NAME_LENGTH(R_INF, 10809), Get_NAME(R_INF, 10809), 10809, is_contain_r((*rI), 10809));
// 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("m64012_190921_234837/111673711/ccs", sg, rI, "beg");
// debug_info_of_specfic_node("m64011_190830_220126/95028102/ccs", sg, rI, "beg");
// debug_info_of_specfic_node("m64011_190830_220126/47516220/ccs", sg, rI, "beg-0");
// debug_info_of_specfic_node("m64012_190920_173625/163644465/ccs", sg, rI, "beg-0");
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL);
// fprintf(stderr, "[M::%s] count_edges_v_w(sg, 49778, 49847)->%ld\n", __func__, count_edges_v_w(sg, 49778, 49847));
@@ -2491,12 +2491,11 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
// print_vw_edge(sg, 34156, 34090, "0");
// stats_chimeric(sg, src, &bu);
if(!is_ou) asg_iterative_semi_circ(sg, src, &bu, max_tip, 1);
asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 1);
asg_arc_cut_chimeric(sg, src, &bu, is_ou?ou_thres:(uint32_t)-1);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL);
// prt_specfic_sge(sg, 10531, 10519, "--1--");
asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0);
asg_arc_cut_inexact(sg, src, &bu, max_tip, is_ou, is_trio, ou_drop_rate/**, NULL**//**&dbg**/);
// debug_edges(&dbg, d, 2);
@@ -2519,7 +2518,6 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
asg_arc_cut_complex_bub_links(sg, &bu, HARD_OL_DROP, HARD_OU_DROP, is_ou, b_mask_t);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL);
// prt_specfic_sge(sg, 10531, 10519, "--5--");
if(is_ou) {
if(ul_refine_alignment(uopt, sg)) update_sg_uo(sg, src);
@@ -2576,6 +2574,98 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
free(bu.a); free(ba.a);
}
int32_t gen_ext_tip(uint32_t v, asg_t *sg, ug_opt_t *uopt, uint8_t *ff, uint32_t min_ovlp, uint64_t *res, uint64_t *path_len)
{
uint32_t k, qn, tn, cc, ccu; int32_t r; asg_arc_t p, pmax;
ma_hit_t_alloc *s = &(uopt->sources[v>>1]); memset(&pmax, 0, sizeof(pmax)); pmax.ol = 0;
for (k = 0; k < s->length; k++) {
qn = Get_qn(s->buffer[k]); tn = Get_tn(s->buffer[k]);
r = ma_hit2arc(&(s->buffer[k]), sg->seq[qn].len, sg->seq[tn].len, uopt->max_hang, asm_opt.max_hang_rate, uopt->min_ovlp, &p);
if(r < 0) continue;
if((p.ul>>32) != v) continue;
if(p.ol < min_ovlp) continue;
if(p.ol > pmax.ol) pmax = p;
}
if(pmax.ol == 0) return 0;
tn = pmax.v>>1;
if((!(sg->seq[tn].del)) || (ff[tn])) return -1;///not a tip
get_R_to_U(uopt->ruIndex, tn, &cc, &ccu);
if(ccu == 1) return -1;
if((cc != (uint32_t)-1) &&
((!(sg->seq[cc].del)) || (ff[cc]))) {
return -1;///contained in an existing read
}
if((*path_len) >= sg->seq[v>>1].len) (*path_len) -= sg->seq[v>>1].len;
else (*path_len) = 0;
(*path_len) += (uint32_t)(pmax.ul) + sg->seq[pmax.v>>1].len;
uint32_t v0 = v; (*res) = pmax.v;
v = (pmax.v^1); s = &(uopt->sources[v>>1]); pmax.ol = 0;
for (k = 0; k < s->length; k++) {
qn = Get_qn(s->buffer[k]); tn = Get_tn(s->buffer[k]);
r = ma_hit2arc(&(s->buffer[k]), sg->seq[qn].len, sg->seq[tn].len, uopt->max_hang, asm_opt.max_hang_rate, uopt->min_ovlp, &p);
if(r < 0) continue;
if((p.ul>>32) != v) continue;
if(p.ol < min_ovlp) continue;
if(p.ol > pmax.ol) pmax = p;
}
if((pmax.v^1) == v0) {
ff[(*res)>>1] = 1; return 1;
}
(*res) = (uint64_t)-1; return -1;///not the longest overlap
}
void ug_ext_gfa(ug_opt_t *uopt, asg_t *sg, uint32_t max_len)
{
asg_arc_t *av; asg64_v res, idx; kv_init(res); kv_init(idx);
uint32_t v, nv, k, z, nvtx = sg->n_seq<<1, tip_n = uopt->tipsLen + 1, s, e;
uint8_t *ff; CALLOC(ff, sg->n_seq); int32_t rr; uint64_t v0, w, bn, plen;
asg_arc_t t, *p;
for (v = 0; v < nvtx; v++) {
if(sg->seq[v>>1].del) continue;
av = asg_arc_a(sg, v); nv = asg_arc_n(sg, v);
for (k = 0; k < nv && av[k].del; k++);
if(k < nv) continue;
for (z = 0, v0 = v, rr = 0, bn = res.n, plen = sg->seq[v0>>1].len; z < tip_n || plen < max_len; z++) {
rr = gen_ext_tip(v0, sg, uopt, ff, 2000, &w, &plen);
if(rr <= 0) break;
kv_push(uint64_t, res, (v0<<32)|w); v0 = w;
}
if(rr < 0 || (z >= tip_n && plen >= max_len)) {
for (z = bn; z < res.n; z++) {
ff[((uint32_t)res.a[z])>>1] = 0;
}
res.n = bn;
}
if(res.n > bn) {
bn <<= 32; bn |= (uint64_t)res.n;
kv_push(uint64_t, idx, bn);
}
}
if(idx.n > 0) {
for (k = 0; k < idx.n; k++) {
s = idx.a[k]>>32; e = (uint32_t)idx.a[k];
for (z = s; z < e; z++) {
v = res.a[z]>>32; w = (uint32_t)res.a[z];
assert((!sg->seq[v>>1].del) && (sg->seq[w>>1].del));
sg->seq[w>>1].del = 0;
rr = gen_spec_edge(sg, uopt, v, w, &t); assert(rr >= 0);
p = asg_arc_pushp(sg); *p = t;
rr = gen_spec_edge(sg, uopt, w^1, v^1, &t); assert(rr >= 0);
p = asg_arc_pushp(sg); *p = t;
}
}
free(sg->idx); sg->idx = 0; sg->is_srt = 0; asg_cleanup(sg);
fprintf(stderr, "[M::%s::] # tips::%u\n", __func__, (uint32_t)idx.n);
}
free(ff); kv_destroy(res); kv_destroy(idx);
}
bubble_type *gen_bubble_chain(asg_t *sg, ma_ug_t *ug, ug_opt_t *uopt, uint8_t **ir_het)
{
@@ -2586,17 +2676,17 @@ bubble_type *gen_bubble_chain(asg_t *sg, ma_ug_t *ug, ug_opt_t *uopt, uint8_t **
asg_t *copy_sg = copy_read_graph(sg);
ma_ug_t *copy_ug = copy_untig_graph(ug);
// fprintf(stderr, "0[M::%s]\n", __func__);
adjust_utg_by_primary(&copy_ug, copy_sg, TRIO_THRES, uopt->sources, uopt->reverse_sources, uopt->coverage_cut,
uopt->tipsLen, uopt->tip_drop_ratio, uopt->stops_threshold, uopt->ruIndex, uopt->chimeric_rate, uopt->drop_ratio,
uopt->max_hang, uopt->min_ovlp, &new_rtg_edges, &cov, uopt->b_mask_t, 0, 0);
ma_ug_destroy(copy_ug); copy_ug = NULL;
asg_destroy(copy_sg); copy_sg = NULL;
// fprintf(stderr, "1[M::%s]\n", __func__);
CALLOC(bub, 1); (*ir_het) = cov->t_ch->ir_het;
cov->t_ch->ir_het = NULL; cov->is_r_het = NULL;
identify_bubbles(ug, bub, (*ir_het), NULL);
// fprintf(stderr, "2[M::%s]\n", __func__);
kv_destroy(new_rtg_edges.a); destory_hap_cov_t(&cov);
if(asm_opt.purge_level_primary == 0) {///all nodes are het
uint32_t k;
@@ -2604,6 +2694,7 @@ bubble_type *gen_bubble_chain(asg_t *sg, ma_ug_t *ug, ug_opt_t *uopt, uint8_t **
if(IF_HOM(k, *bub)) bub->index[k] = bub->f_bub+1;
}
}
// fprintf(stderr, "3[M::%s]\n", __func__);
return bub;
}
@@ -9890,6 +9981,30 @@ static inline int usg_end(const usg_t *g, uint32_t v, uint64_t *lw)
return ASG_ET_MERGEABLE;
}
uint32_t usg_real_tip(usg_t *g, uint32_t v0, double rate)
{
usg_arc_t *av, *aw; uint32_t nv, nw, i, k, ov, ow;
av = usg_arc_a(g, v0);
nv = usg_arc_n(g, v0);
for (i = 0; i < nv; i++) {
if(av[i].del) continue;
aw = usg_arc_a(g, (av[i].v^1));
nw = usg_arc_n(g, (av[i].v^1));
for (k = ov = ow = 0; k < nw; k++) {
if(aw[k].del) continue;
if((aw[k].v>>1) == (v0>>1)) {
if(aw[k].ol > ov) ov = aw[k].ol;
} else {
if(aw[k].ol > ow) ow = aw[k].ol;
}
}
if(ow == 0 || ov >= ow) return 0;
if(ov > (ow*rate)) return 0;
}
return 1;
}
uint32_t usg_arc_cut_tips(usg_t *g, uint32_t max_ext, uint32_t ignore_ul, asg64_v *in)
{
asg64_v tx = {0,0,0}, *b = NULL;
@@ -9940,7 +10055,7 @@ uint32_t usg_arc_cut_tips(usg_t *g, uint32_t max_ext, uint32_t ignore_ul, asg64_
if(kv <= max_ext) {
ff = 0;
if(!ignore_ul) {
if(!ignore_ul) {///consider UL
for (i = pb; i + 1 < b->n; i++) {
p = get_usg_arc(g, ((uint32_t)b->a[i]), ((uint32_t)b->a[i+1])); assert(p);
if(p->ou > 1) {//ignore ou == 1
@@ -9961,10 +10076,30 @@ uint32_t usg_arc_cut_tips(usg_t *g, uint32_t max_ext, uint32_t ignore_ul, asg64_
}
}
// if((v>>1) == 139131) {
// fprintf(stderr, "+[M::%s::] v>>1::%u, v&1::%u, kv::%u, ff::%u, pb::%u, b->n::%u, w>>1::%u, w&1::%u\n", __func__, v>>1, v&1, kv,
// ff, pb, (uint32_t)b->n, ((uint32_t)b->a[b->n-1])>>1, ((uint32_t)b->a[b->n-1])&1);
// }
if((ff == 0) && (pb < b->n) && (!usg_real_tip(g, ((uint32_t)b->a[b->n-1]), 0.75))) {
ff = 1;
}
// if((v>>1) == 139131) {
// fprintf(stderr, "-[M::%s::] v>>1::%u, v&1::%u, kv::%u, ff::%u, pb::%u, b->n::%u, w>>1::%u, w&1::%u\n", __func__, v>>1, v&1, kv,
// ff, pb, (uint32_t)b->n, ((uint32_t)b->a[b->n-1])>>1, ((uint32_t)b->a[b->n-1])&1);
// }
if(ff == 0) {
for (i = pb; i < b->n; i++) usg_seq_del(g, ((uint32_t)b->a[i])>>1);
cnt++;
}
// if((v>>1) == 139131) {
// fprintf(stderr, ">[M::%s::] v>>1::%u, v&1::%u, kv::%u, ff::%u, pb::%u, b->n::%u, w>>1::%u, w&1::%u, del::%u\n", __func__, v>>1, v&1, kv,
// ff, pb, (uint32_t)b->n, ((uint32_t)b->a[b->n-1])>>1, ((uint32_t)b->a[b->n-1])&1,
// g->a[((uint32_t)b->a[b->n-1])>>1].del);
// }
}
b->n = pb;
}
@@ -13036,7 +13171,9 @@ void renew_usg_t_bub(ul_resolve_t *uidx, usg_t *ng, uint32_t *id_map, uint8_t *f
if(k + 1 == u->n) id_map[v>>1] |= p[(v)&1];
// dn++;
}
// fprintf(stderr, "+[M::%s] i::%u\n", __func__, i);
merge_hybrid_utg_content(u, uidx->l1_ug, uidx->sg, ng, &e);
// fprintf(stderr, "-[M::%s] i::%u\n", __func__, i);
ug->g->seq[i].len = u->len;
}
kv_destroy(e.a);
@@ -13917,6 +14054,11 @@ static void gen_scaffold_id(void *data, long i, int tid) // callback for kt_for(
uc_block_t *p, *n, *m0, *m1;
sv = (v&1?((ug->u.a[v>>1].a[0]>>32)^1):(ug->u.a[v>>1].a[ug->u.a[v>>1].n-1]>>32));
sw = (w&1?((ug->u.a[w>>1].a[ug->u.a[w>>1].n-1]>>32)^1):(ug->u.a[w>>1].a[0]>>32));
// if((((v>>1) == 15016) && ((w>>1) == 2778)) || (((v>>1) == 2778) && ((w>>1) == 15016))) {
// fprintf(stderr, "[M::%s::]\t%.*s(%c)(sv>>1::%u)(utg%.6ul(%c))\t%.*s(%c)(sw>>1::%u)(utg%.6ul(%c))\n", __func__,
// (int)Get_NAME_LENGTH(R_INF, (sv>>1)), Get_NAME(R_INF, (sv>>1)), "+-"[sv&1], sv>>1, (v>>1)+1, "+-"[v&1],
// (int)Get_NAME_LENGTH(R_INF, (sw>>1)), Get_NAME(R_INF, (sw>>1)), "+-"[sw&1], sw>>1, (w>>1)+1, "+-"[w&1]);
// }
sv = Get_READ_LENGTH(R_INF, (sv>>1)); sw = Get_READ_LENGTH(R_INF, (sw>>1));
max_ql = MIN(sv, sw); max_ql = -max_ql;
// if((((v>>1) == 235) && ((w>>1) == 153)) || (((v>>1) == 153) && ((w>>1) == 235))) {
@@ -13942,8 +14084,8 @@ static void gen_scaffold_id(void *data, long i, int tid) // callback for kt_for(
if(n->base || (!n->el) || (!n->pchain) || (n->pidx != (uint32_t)-1)) continue;
uw = (((uint32_t)(n->hid))<<1)|((uint32_t)(n->rev));
if(uw == w) {
// if((((v>>1) == 235) && ((w>>1) == 153)) || (((v>>1) == 153) && ((w>>1) == 235))) {
// fprintf(stderr, "+[M::%s] ul_id::%lu\n", __func__, a[k]>>32);
// if((((v>>1) == 15016) && ((w>>1) == 2778)) || (((v>>1) == 2778) && ((w>>1) == 15016))) {
// fprintf(stderr, "+[M::%s] ul_id::%lu, qs::%u, qe::%u\n", __func__, a[k]>>32, n->qs, n->qe);
// // print_ul_alignment(ug, &UL_INF, a[k]>>32, "+++");
// }
ql = ((int64_t)n->qs) - ((int64_t)p->qe);
@@ -13961,9 +14103,9 @@ static void gen_scaffold_id(void *data, long i, int tid) // callback for kt_for(
if(n->base || (!n->el) || (!n->pchain) || (n->aidx != (uint32_t)-1)) continue;
uw = (((uint32_t)(n->hid))<<1)|((uint32_t)(n->rev)); uw ^= 1;
if(uw == w) {
// if((((v>>1) == 235) && ((w>>1) == 153)) || (((v>>1) == 153) && ((w>>1) == 235))) {
// fprintf(stderr, "-[M::%s] ul_id::%lu\n", __func__, a[k]>>32);
// // print_ul_alignment(ug, &UL_INF, a[k]>>32, "---");
// if((((v>>1) == 15016) && ((w>>1) == 2778)) || (((v>>1) == 2778) && ((w>>1) == 15016))) {
// fprintf(stderr, "-[M::%s] ul_id::%lu, qs::%u, qe::%u\n", __func__, a[k]>>32, n->qs, n->qe);
// // print_ul_alignment(ug, &UL_INF, a[k]>>32, "+++");
// }
ql = ((int64_t)p->qs) - ((int64_t)n->qe);
if(ql > max_ql && n->qe < p->qe && n->qs < p->qs) {
@@ -13976,13 +14118,13 @@ static void gen_scaffold_id(void *data, long i, int tid) // callback for kt_for(
}
if(m0 && m1 && id != ((uint32_t)-1)) {
// if((((v>>1) == 235) && ((w>>1) == 153)) || (((v>>1) == 153) && ((w>>1) == 235))) {
// fprintf(stderr, "[M::%s] id::%u, v>>1::%u(%c), w>>1::%u(%c)\n",
// __func__, id, v>>1, "+-"[v&1], w>>1, "+-"[w&1]);
// }
z->ulid = (id<<1)|rev;
z->raw_sof = m0->qe;
z->raw_eof = m1->qs;
// if((((v>>1) == 15016) && ((w>>1) == 2778)) || (((v>>1) == 2778) && ((w>>1) == 15016))) {
// fprintf(stderr, "[M::%s] id::%u, v>>1::%u(%c), w>>1::%u(%c), raw_sof::%u, raw_eof::%u\n",
// __func__, id, v>>1, "+-"[v&1], w>>1, "+-"[w&1], z->raw_sof, z->raw_eof);
// }
}
}
@@ -14069,7 +14211,7 @@ ma_ug_t *convert_usg_t(usg_t *ng, ma_ug_t *ref)
return ug;
}
int32_t gen_spec_rc_edge(asg_t *rg, ug_opt_t *uopt, uint32_t v, uint32_t w, asg_arc_t *t)
int32_t gen_spec_rc_edge(asg_t *rg, ug_opt_t *uopt, uint32_t v, uint32_t w, asg_arc_t *t, ma_hit_t **te)
{
uint32_t k, qn, tn; int32_t r; ma_hit_t_alloc *s = &(uopt->reverse_sources[v>>1]); asg_arc_t p;
for (k = 0; k < s->length; k++) {
@@ -14078,7 +14220,7 @@ int32_t gen_spec_rc_edge(asg_t *rg, ug_opt_t *uopt, uint32_t v, uint32_t w, asg_
r = ma_hit2arc(&(s->buffer[k]), rg->seq[qn].len, rg->seq[tn].len, uopt->max_hang, asm_opt.max_hang_rate, uopt->min_ovlp, &p);
if(r < 0) continue;
if((p.ul>>32) != v || p.v != w) continue;
*t = p; t->ou = 0;
*t = p; t->ou = 0; if(te) (*te) = &(s->buffer[k]);
return 1;
}
return -1;
@@ -14086,10 +14228,10 @@ int32_t gen_spec_rc_edge(asg_t *rg, ug_opt_t *uopt, uint32_t v, uint32_t w, asg_
void push_direct_scaff_node(usc_t *psa, asg_t *ng, ug_opt_t *uopt, uint64_t vx, uint64_t wx)
{
uint64_t uol, dif; int32_t r; asg_arc_t t, *p;
uint64_t uol, dif; int32_t r; asg_arc_t t, *p; ma_hit_t *rc0, *rc1;
uol = psa->raw_sof - psa->raw_eof;
assert(uol >= (uint64_t)uopt->min_ovlp);
r = gen_spec_rc_edge(ng, uopt, vx, wx, &t);
r = gen_spec_rc_edge(ng, uopt, vx, wx, &t, &rc0);
if(r >= 0) {
if(uol > t.ol) dif = uol - t.ol;
else dif = t.ol - uol;
@@ -14097,11 +14239,17 @@ void push_direct_scaff_node(usc_t *psa, asg_t *ng, ug_opt_t *uopt, uint64_t vx,
if(dif > 16) r = -1;
}
}
// if((((vx>>1) == 3313901) && ((wx>>1) == 3555473)) || (((vx>>1) == 3555473) && ((wx>>1) == 3313901))) {
// fprintf(stderr, "[M::%s::]\t%.*s(%c)(sv>>1::%lu)\t%.*s(%c)(sw>>1::%lu)\tr::%d\n", __func__,
// (int)Get_NAME_LENGTH(R_INF, (vx>>1)), Get_NAME(R_INF, (vx>>1)), "+-"[vx&1], vx>>1,
// (int)Get_NAME_LENGTH(R_INF, (wx>>1)), Get_NAME(R_INF, (wx>>1)), "+-"[wx&1], wx>>1, r);
// }
if(r >= 0) {
p = asg_arc_pushp(ng); *p = t;
r = gen_spec_rc_edge(ng, uopt, wx^1, vx^1, &t);
r = gen_spec_rc_edge(ng, uopt, wx^1, vx^1, &t, &rc1);
assert(r >= 0); p = asg_arc_pushp(ng); *p = t;
add_ma_hit_t_alloc(&(R_INF.paf[Get_qn((*rc0))]), rc0);
add_ma_hit_t_alloc(&(R_INF.paf[Get_qn((*rc1))]), rc1);
} else {
push_scaff_node(R_INF.paf, vx, wx, psa->raw_sof-psa->raw_eof, ng);
@@ -14408,7 +14556,14 @@ void u2g_hybrid_detan_iter(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, uint
// ncut += ug_ext_strict(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64);
// prt_usg_t(uidx, ng, "ng.db");
if(ncut) {
usg_cleanup(ng); usg_arc_cut_tips(ng, max_ext, 1, ub);
usg_cleanup(ng);
// if((max_ext>>1) > 0) usg_arc_cut_tips(ng, (max_ext>>1), 1, ub);
usg_arc_cut_tips(ng, max_ext, 1, ub);
// if(ng->n > 139131) {
// fprintf(stderr, ">[M::%s::] ng->a[139131].del::%u\n", __func__,
// ng->a[139131].del);
// }
}
// prt_usg_t(uidx, ng, "ng2");
// u2g_hybrid_aln(uidx, ng, ob, ub);
@@ -16521,15 +16676,21 @@ ul_renew_t *ropt)
if(sg->seq[i].del) continue;
sg->seq[i].c = PRIMARY_LABLE;
}
// fprintf(stderr, "0[M::%s]\n", __func__);
hic_clean(sg);
// fprintf(stderr, "1[M::%s]\n", __func__);
ma_ug_t *init_ug = ul_realignment(uopt, sg, 0);
// fprintf(stderr, "2[M::%s]\n", __func__);
// exit(1);
filter_sg_by_ug(sg, init_ug, uopt);
// fprintf(stderr, "3[M::%s]\n", __func__);
// print_debug_gfa(sg, init_ug, uopt->coverage_cut, "UL.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0);
// print_ul_alignment(init_ug, &UL_INF, 47072, "after-0");
bub = gen_bubble_chain(sg, init_ug, uopt, &r_het);
// fprintf(stderr, "4[M::%s]\n", __func__);
// print_ul_alignment(init_ug, &UL_INF, 47072, "after-1");
ul_resolve_t *uidx = init_ul_resolve_t(sg, init_ug, bub, &UL_INF, uopt, r_het);
// fprintf(stderr, "5[M::%s]\n", __func__);
// print_ul_alignment(init_ug, &UL_INF, 47072, "after-2");
// exit(1);
// print_raw_uls_seq(uidx, asm_opt.output_file_name);
@@ -16556,7 +16717,7 @@ ul_renew_t *ropt)
asg_t *ng = renew_ng(uidx->uovl.h_usg, init_ug, sg, uopt, ropt->cov, ropt->ruIndex, 16);
destory_bubbles(uidx->bub); free(uidx->bub); ma_ug_destroy(init_ug);
destory_bubbles(uidx->bub); free(uidx->bub); ma_ug_destroy(init_ug); free(r_het);
destroy_ul_resolve_t(uidx); destory_all_ul_t(&UL_INF);
asg_destroy((*(ropt->sg))); (*(ropt->sg)) = ng;
gfa_ut.h
+2
View File
@@ -3,6 +3,7 @@
#include "Overlaps.h"
#include "hic.h"
typedef struct {
asg_t *g;
ma_hit_t_alloc *src;
@@ -35,5 +36,6 @@ void post_rescue(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc
// void print_raw_u2rgfa_seq(all_ul_t *aln, R_to_U* rI, uint32_t is_detail);
bubble_type *gen_bubble_chain(asg_t *sg, ma_ug_t *ug, ug_opt_t *uopt, uint8_t **ir_het);
void filter_sg_by_ug(asg_t *rg, ma_ug_t *ug, ug_opt_t *uopt);
void ug_ext_gfa(ug_opt_t *uopt, asg_t *sg, uint32_t max_len);
#endif
hic.cpp
+96 -10
View File
@@ -2142,6 +2142,7 @@ void print_hits_simp(ha_ug_index* idx, kvec_pe_hit* hits)
char dir[2] = {'+', '-'};
for (k = 0; k < hits->a.n; ++k)
{
if(((hits->a.a[k].s<<1)>>shif) == ((hits->a.a[k].e<<1)>>shif)) continue;
fprintf(stderr, "r-%lu-th\t%c\trs-utg%.6dl\t%lu\t%c\tre-utg%.6dl\t%lu\n",
hits->a.a[k].id,
dir[hits->a.a[k].s>>63], (int)((hits->a.a[k].s<<1)>>shif)+1, hits->a.a[k].s&idx->pos_mode,
@@ -7088,7 +7089,7 @@ void get_bub_graph(ma_ug_t* ug, bubble_type* bub)
adjecent = 0;
while (pre_id != v)
{
if(connect_bub_occ(bub, pre_id>>1, bub->check_het) > 0)
if(connect_bub_occ(bub, pre_id>>1, bub->check_het) > 0)///some bubbles between v and k
{
adjecent = 1;
break;
@@ -16480,7 +16481,7 @@ void optimize_u_trans(kv_u_trans_t *ovlp, kvec_pe_hit* hits, ha_ug_index* idx)
for (i = 0; i < ovlp->n; i++){
x = &(ovlp->a[i]);
if(x->qn > x->tn) continue;
if(x->f != RC_2 || x->del) continue;
if(x->f != RC_2 || x->del) continue;///no need to adjust RC_0/RC_1
occ = get_oe_occ(x->qn, x->tn, hits, idx) + get_oe_occ(x->tn, x->qn, hits, idx);///how many UL bridging qn and tn
kv_pushp(u_trans_t, k_trans, &p);
(*p) = (*x); p->nw = (x->nw*(1-(((double)(occ<<1))/((double)(hits->occ.a[x->qn]+hits->occ.a[x->tn])))));
@@ -16731,8 +16732,13 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o
write_hc_hits(&sl.hits, idx->ug, asm_opt.output_file_name);
}
sl.hits.uID_bits = idx->uID_bits; sl.hits.pos_mode = idx->pos_mode;
optimize_u_trans(&(idx->t_ch->k_trans), &sl.hits, idx);
filter_kv_u_trans_t(&(idx->t_ch->k_trans), idx->ug, 0.5);
/***debug***/
if(sl.hits.idx.n == 0) idx_hc_links(&(sl.hits), idx, NULL);
// optimize_u_trans(&(idx->t_ch->k_trans), &sl.hits, idx);
// filter_kv_u_trans_t(&(idx->t_ch->k_trans), idx->ug, 0.5);
/***debug***/
// print_kv_u_trans_t(&(idx->t_ch->k_trans), idx->ug);
// flter_by_cov(idx, &sl.hits, 2);
// update_hits(idx, &sl.hits, idx->t_ch->is_r_het);
@@ -16740,7 +16746,7 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o
////dedup_hits(&(sl.hits), sl.idx);
///write_hc_hits_v14(&sl.hits, asm_opt.output_file_name);
if(asm_opt.misjoin_len > 0)
if(asm_opt.misjoin_len > 0/** && (!(asm_opt.ar))**/)//disable it for the UL assembly
{
update_switch_unitig(idx->ug, idx->read_g, &(sl.hits), &(idx->t_ch->k_trans), 10, 20, asm_opt.misjoin_len, 0.15);
renew_idx_para(idx, idx->ug);
@@ -16782,13 +16788,15 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o
/*******************************for debug************************************/
mc_solve(NULL, NULL, &k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag,
(bub.round_id == 0? 1 : 0), s->s, 1, &bub,
&(idx->t_ch->k_trans), 0, 0/**(((bub.round_id+1) == bub.n_round)?1:0)**/);
&(idx->t_ch->k_trans), 0, /**(((bub.round_id+1) == bub.n_round)?1:0)**/0);
// mc_solve(NULL, NULL, &k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag,
// (bub.round_id == 0? 1 : 0), s->s, 1, NULL, &(idx->t_ch->k_trans), 0, 0);
// if((bub.round_id+1) == bub.n_round) {
// prt_debug_hic(asm_opt.output_file_name, idx->ug, idx, opt, &sl.hits,
// &(idx->t_ch->k_trans), &k_trans, &link, s->s, &(bub));
// }
/**
if((bub.round_id+1) == bub.n_round) {
prt_debug_hic(asm_opt.output_file_name, idx->ug, idx, opt, &sl.hits,
&(idx->t_ch->k_trans), &k_trans, &link, s->s, &(bub));
}
**/
/*******************************for debug************************************/
label_unitigs_sm(s->s, NULL, idx->ug);
@@ -16853,6 +16861,84 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o
return 1;
}
void bp_solve(ug_opt_t *opt, kv_u_trans_t *ref, ma_ug_t *ug, asg_t *sg, bubble_type *bub, double cis_rate)
{
uint64_t k, i, v, nvx, n[3], nv; u_trans_t *z;
kv_u_trans_t in; memset(&in, 0, sizeof(in));
ma_utg_t *u = NULL; asg_arc_t *av; double w;
ps_t *s = init_ps_t(11, ug->g->n_seq);
for (i = 0; i < ug->g->n_seq; i++) {
s->s[i] = 0;
if((ug->g->seq[i].del) || (IF_HOM(i, (*bub)))) {
continue;
}
u = &ug->u.a[i];
if(u->m == 0 || u->n == 0) continue;
for (k = n[0] = n[1] = n[2] = 0; k < u->n; k++) {
if(R_INF.trio_flag[u->a[k]>>33] == FATHER) n[1]++;
else if(R_INF.trio_flag[u->a[k]>>33] == MOTHER) n[2]++;
else n[0]++;
}
if(n[1] >= n[2] && n[1] >= n[0]) s->s[i] = 1;
else if(n[2] >= n[1] && n[2] >= n[0]) s->s[i] = -1;
}
kv_resize(u_trans_t, in, ref->n);
for (k = in.n = 0; k < ref->n; k++) {
if((IF_HOM(ref->a[k].qn, (*bub))) || (IF_HOM(ref->a[k].tn, (*bub))) || ref->a[k].del) {
continue;
}
kv_push(u_trans_t, in, ref->a[k]);
}
if(cis_rate > 0) {
nvx = ug->g->n_seq<<1;
for (v = 0; v < nvx; v++) {
if((ug->g->seq[v>>1].del) || (IF_HOM((v>>1), (*bub)))) continue;
av = asg_arc_a(ug->g, v); nv = asg_arc_n(ug->g, v);
for (i = 0; i < nv; i++) {
if((av[i].del) || (av[i].v>>1) == (v>>1) ||
(ug->g->seq[av[i].v>>1].del) || (IF_HOM((av[i].v>>1), (*bub)))) {
continue;
}
w = ((double)av[i].ol)*cis_rate;
if(w <= 0) continue;
kv_pushp(u_trans_t, in, &z);
memset(z, 0, sizeof((*z)));
z->f = RC_3; z->nw =-w; z->rev = (v^av[i].v)&1;
z->qn = v>>1; z->tn = av[i].v>>1;
z->qs = 0; z->qe = sg->seq[z->qn].len;
if(av[i].ol < sg->seq[z->qn].len) {
if(v&1) {
z->qe = z->qs + av[i].ol;
} else {
z->qs = z->qe - av[i].ol;
}
}
z->ts = 0; z->te = sg->seq[z->tn].len;
if(av[i].ol < sg->seq[z->tn].len) {
if(av[i].v&1) {
z->ts = z->te - av[i].ol;
} else {
z->te = z->ts + av[i].ol;
}
}
}
}
}
// fprintf(stderr, "[M::%s::] ==> nup::%lu\n", __func__, nup);
clean_u_trans_t_idx_adv(&in, ug, sg);
mc_solve(NULL, NULL, &in, ug, sg, 0.8, NULL, 0, s->s, 1, bub, ref, 0, 0);
label_unitigs_sm(s->s, NULL, ug);
destory_ps_t(&s); kv_destroy(in); kv_destroy(in.idx);
}
spg_t *hic_short_pre_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_opt_t *opt, kvec_pe_hit **rhits)
{
// double index_time = yak_realtime();
hic.h
+2
View File
@@ -8,6 +8,7 @@
#define RC_0 0
#define RC_1 1
#define RC_2 2
#define RC_3 3
hc_edge* get_hc_edge(hc_links* link, uint64_t src, uint64_t dest, uint64_t dir);
hc_edge* push_hc_edge(hc_linkeage* x, uint64_t uID, double weight, int dir, uint64_t* d);
@@ -111,5 +112,6 @@ void dedup_hits(kvec_pe_hit* hits, uint64_t is_dup);
void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, uint32_t is_poy, kvec_pe_hit **rhits);
spg_t *hic_pre_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, kvec_pe_hit **rhits);
void prt_bubble_gfa_adv(FILE *fp, bubble_type *bub, const char* utg_pre, const char* bub_pre, const char* chain_pre);
void bp_solve(ug_opt_t *opt, kv_u_trans_t *ref, ma_ug_t *ug, asg_t *sg, bubble_type *bub, double cis_rate);
#endif
horder.cpp
+12 -6
View File
@@ -960,7 +960,7 @@ int64_t cov_s_pos, int64_t cov_e_pos, uint64_t *s, uint64_t *e)
(*s) = (*e) = (uint64_t)-1;
res->n = 0;
get_hic_cov_interval(b, b_n, cutoff, &cov_s_pos, &cov_e_pos, res);
if(res->n == 0) return;
if(res->n == 0) return;///res keeps all intervals with >= cutoff coverage
int64_t max = -1, max_cur = 0;
int64_t max_s_idx, max_e_idx, cur_s_idx;
for (i = 0; i < res->n; i++)//all intervals have cov >= cutoff
@@ -1842,7 +1842,7 @@ void break_phasing_utg(ma_ug_t *ug, asg_t *rg, kvec_pe_hit *hits, kv_u_trans_t *
kv_pushp(h_cov_t, join, &p);
p->dp = (b_points->a[l].s<<32)|(ug->u.n-1);
p->s = rsi; p->e = ug->u.a[ug->u.n-1].len;
kv_push(uint64_t, dbg_N50, (uint64_t)(ug->u.a[ug->u.n-1].len)+((uint64_t)(1)<<63));
kv_push(uint64_t, dbg_N50, (uint64_t)(ug->u.a[ug->u.n-1].len)+((uint64_t)(1)<<63));///new nodes
}
}
pidx = idx;
@@ -1858,13 +1858,13 @@ void break_phasing_utg(ma_ug_t *ug, asg_t *rg, kvec_pe_hit *hits, kv_u_trans_t *
kv_pushp(h_cov_t, join, &p);
p->dp = (b_points->a[l].s<<32)|(ug->u.n-1);
p->s = rsi; p->e = ug->u.a[ug->u.n-1].len;
kv_push(uint64_t, dbg_N50, (uint64_t)(ug->u.a[ug->u.n-1].len)+((uint64_t)(1)<<63));
kv_push(uint64_t, dbg_N50, (uint64_t)(ug->u.a[ug->u.n-1].len)+((uint64_t)(1)<<63));///new nodes
}
}
if(de_u)
{
kv_push(uint64_t, dbg_N50, ug->u.a[b_points->a[l].s].len);
kv_push(uint64_t, dbg_N50, ug->u.a[b_points->a[l].s].len);///old nodes
update_nus(&(ug->u.a[b_points->a[l].s]), ug, join.a + pn, join.n - pn);
free(ug->u.a[b_points->a[l].s].a); free(ug->u.a[b_points->a[l].s].s);
memset(&(ug->u.a[b_points->a[l].s]), 0, sizeof(ug->u.a[b_points->a[l].s]));
@@ -1890,7 +1890,7 @@ void break_phasing_utg(ma_ug_t *ug, asg_t *rg, kvec_pe_hit *hits, kv_u_trans_t *
if(i < oug_n)
{
kv_pushp(h_cov_t, join, &p);
p->dp = (i<<32)|(m);
p->dp = (i<<32)|(m);///current id | updated id
p->s = 0; p->e = ug->u.a[i].len;
}
@@ -1963,7 +1963,11 @@ void break_phasing_utg(ma_ug_t *ug, asg_t *rg, kvec_pe_hit *hits, kv_u_trans_t *
free(k_trans->a); free(k_trans->idx.a); memset(k_trans, 0, sizeof(*k_trans));
k_trans->a = n_trans->a; k_trans->n = n_trans->n; k_trans->m = n_trans->m;
free(n_trans);
clean_u_trans_t_idx(k_trans, ug, rg);
clean_u_trans_t_idx_adv(k_trans, ug, rg);
// clean_u_trans_t_idx(k_trans, ug, rg);
uint64_t uID_bits, pos_mode;
for (uID_bits=1; (uint64_t)(1<<uID_bits)<(uint64_t)ug->u.n; uID_bits++);
pos_mode = ((uint64_t)-1)>>(uID_bits+1);
@@ -2069,6 +2073,8 @@ uint64_t min_ulen, double boundaryRate)
// fprintf(stderr, "consensus_break-s: %lu, e: %lu\n", cov_buf.a[i].s, cov_buf.a[i].e);
// }
/*******************************for debug************************************/
///res -> all low coverage intervals;
///cov_buf -> consensus coverage intervals;
get_read_breaks(&(ug->u.a[get_hit_suid(*hits, l)]), rg, &cov_buf,
&res, hits, l, k, ulen, &bs, &dp);
if(bs == (uint64_t)-1) fprintf(stderr, "ERROR-read\n");
inter.cpp
+1170 -126
View File
File diff suppressed because it is too large Load Diff
inter.h
+2
View File
@@ -2,6 +2,7 @@
#define __INTER__
#include "Overlaps.h"
#include "Process_Read.h"
#include "hic.h"
#define G_CHAIN_BW 16//128
#define FLANK_M (0x7fffU)
@@ -116,5 +117,6 @@ void gen_ul_vec_rid_t(all_ul_t *x, All_reads *rdb, ma_ug_t *ug);
void update_ug_arch_ul_mul(ma_ug_t *ug);
void print_ul_alignment(ma_ug_t *ug, all_ul_t *aln, uint32_t id, const char* cmd);
void clear_all_ul_t(all_ul_t *x);
void trans_base_infer(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, kv_u_trans_t *res, bubble_type *bub);
#endif
rcut.cpp
+66 -6
View File
@@ -3241,6 +3241,15 @@ mc_clus_t *init_mc_clus_t(const mc_opt_t *opt, mc_g_t *mg, bubble_type* bub, uin
return p;
}
void des_mc_clus_t(mc_clus_t *p)
{
if((!p)) return;
uint32_t k; free(p->lock);
for (k = 0; k < p->n_thread; k++) free(p->aux[k].vis.a);
free(p->aux); free(p->cc.ng.a); free(p->cc.nn.a);
free(p);
}
void mc_solve_core_adv(const mc_opt_t *opt, mc_g_t *mg, bubble_type* bub)
{
double index_time = yak_realtime();
@@ -3277,7 +3286,7 @@ void mc_solve_core_adv(const mc_opt_t *opt, mc_g_t *mg, bubble_type* bub)
// if(bp) mc_solve_bp(bp);
///mc_write_info(g, b);
mc_svaux_destroy(b);
mc_svaux_destroy(b); des_mc_clus_t(bc);
// if(bp) destroy_mc_bp_t(&bp);
fprintf(stderr, "[M::%s::%.3f] ==> Partition\n", __func__, yak_realtime()-index_time);
}
@@ -4464,6 +4473,28 @@ void dump_bubble_type(bubble_type* bub, FILE *fp)
dump_ma_ug_t(bub->ug, fp);
}
void dump_rid(All_reads *rdb, FILE *fp)
{
fwrite(&(rdb->total_reads), sizeof(rdb->total_reads), 1, fp);
fwrite(&(rdb->name_index_size), sizeof(rdb->name_index_size), 1, fp);
fwrite(rdb->name_index, sizeof((*(rdb->name_index))), rdb->name_index_size, fp);
fwrite(&(rdb->total_name_length), sizeof(rdb->total_name_length), 1, fp);
fwrite(rdb->name, sizeof((*(rdb->name))), rdb->total_name_length, fp);
}
void load_rid(All_reads *rdb, FILE *fp)
{
fread(&(rdb->total_reads), sizeof(rdb->total_reads), 1, fp);
fread(&(rdb->name_index_size), sizeof(rdb->name_index_size), 1, fp);
MALLOC(rdb->name_index, rdb->name_index_size);
fread(rdb->name_index, sizeof((*(rdb->name_index))), rdb->name_index_size, fp);
fread(&(rdb->total_name_length), sizeof(rdb->total_name_length), 1, fp);
MALLOC(rdb->name, rdb->total_name_length);
fread(rdb->name, sizeof((*(rdb->name))), rdb->total_name_length, fp);
}
void dump_debug_phasing(const char* fn, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate,
uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref)
{
@@ -4481,6 +4512,7 @@ uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *re
fwrite(&is_sys, sizeof(is_sys), 1, fp);///is_sys
dump_bubble_type(bub, fp);///bub
dump_kv_u_trans_t(ref, fp);///ta
dump_rid(&R_INF, fp);
fclose(fp);
free(buf);
@@ -4586,15 +4618,43 @@ uint32_t *renew_s, int8_t **s, uint32_t *is_sys, bubble_type **bub, kv_u_trans_t
CALLOC((*bub), 1); load_bubble_type(*bub, fp);///bub
CALLOC((*ref), 1); load_kv_u_trans_t(*ref, fp);///ta
load_rid(&R_INF, fp);///read id
fclose(fp);
free(buf);
}
void prt_rcut_res(const char* fn, int8_t *s, ma_ug_t *ug, asg_t *sg)
{
uint32_t i, k, flag = AMBIGU;
char *buf = (char*)calloc(strlen(fn) + 50, 1);
sprintf(buf, "%s.rcut.res.log", fn);
FILE *fp = fopen(buf, "w");
uint8_t *rs = NULL; MALLOC(rs, sg->n_seq);
memset(rs, AMBIGU, sg->n_seq*sizeof((*rs)));
for (i = 0; i < ug->g->n_seq; i++) {
if(ug->g->seq[i].del) continue;
flag = AMBIGU;
if(s[i] == 0) continue;
flag = (s[i] > 0? FATHER:MOTHER);
for (k = 0; k < ug->u.a[i].n; k++) rs[ug->u.a[i].a[k]>>33] = flag;
}
for (i = 0; i < sg->n_seq; i++) {
if(rs[i] == AMBIGU) continue;
fprintf(fp, "%.*s\t%u\n", (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i), rs[i]);
}
fclose(fp); free(rs);
free(buf);
}
void quick_debug_phasing(const char* fn)
{
kv_u_trans_t *ta; ma_ug_t *ug; asg_t *read_g; double f_rate;
uint32_t renew_s; int8_t *s; uint32_t is_sys, k; bubble_type *bub; kv_u_trans_t *ref;
uint32_t renew_s; int8_t *s; uint32_t is_sys; bubble_type *bub; kv_u_trans_t *ref;
load_debug_phasing(fn, &ta, &ug, &read_g, &f_rate, &renew_s, &s, &is_sys, &bub, &ref);
@@ -4602,10 +4662,10 @@ void quick_debug_phasing(const char* fn)
mc_solve(NULL, NULL, ta, ug, read_g, f_rate, NULL, renew_s, s, is_sys, bub, ref, 0, 0);
// mc_solve_core_adv(const mc_opt_t *opt, mc_g_t *mg, bubble_type* bub, kv_u_trans_t *ref)
for (k = 0; k < ug->g->n_seq; k++) {
fprintf(stderr, "utg%.6dl(len::%u), s[k]::%d\n", (int32_t)(k)+1, ug->g->seq[k].len, s[k]);
}
// for (k = 0; k < ug->g->n_seq; k++) {
// fprintf(stderr, "utg%.6dl(len::%u), s[k]::%d\n", (int32_t)(k)+1, ug->g->seq[k].len, s[k]);
// }
prt_rcut_res(fn, s, ug, read_g);
exit(1);
}
sketch.cpp
+1 -1
View File
@@ -461,7 +461,7 @@ void sf##_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is
uint32_t buf_p[256], min_s = (uint32_t)-1;\
tiny_queue_t tq;\
assert(len > 0 && (int64_t)(len) < (int64_t)((((uint64_t)1)<<PosBits)) && (int64_t)(rid) < (int64_t)((((uint64_t)1)<<RidBits)) && (w > 0 && w < 256) && (k > 0 && k <= 63));\
if (dbg_ct != NULL) dbg_ct->a.n = 0;\
if (dbg_ct != NULL) dbg_ct->a.n = 0;\
if (k_flag != NULL) {\
kv_resize_km(km, uint8_t, k_flag->a, (uint64_t)len);\
k_flag->a.n = len;\